**Alternative Vegetation States in Tropical Forests and Savannas: The Search for Consistent Signals in Diverse Remote Sensing Data**

**Sanath Sathyachandran Kumar 1,\*, Niall P. Hanan 1, Lara Prihodko 2, Julius Anchang 1, C. Wade Ross 1, Wenjie Ji <sup>1</sup> and Brianna M Lind <sup>1</sup>**


Received: 25 January 2019; Accepted: 1 April 2019; Published: 4 April 2019

**Abstract:** Globally, the spatial distribution of vegetation is governed primarily by climatological factors (rainfall and temperature, seasonality, and inter-annual variability). The local distribution of vegetation, however, depends on local edaphic conditions (soils and topography) and disturbances (fire, herbivory, and anthropogenic activities). Abrupt spatial or temporal changes in vegetation distribution can occur if there are positive (i.e., amplifying) feedbacks favoring certain vegetation states under otherwise similar climatic and edaphic conditions. Previous studies in the tropical savannas of Africa and other continents using the MODerate Resolution Imaging Spectroradiometer (MODIS) vegetation continuous fields (VCF) satellite data product have focused on discontinuities in the distribution of tree cover at different rainfall levels, with bimodal distributions (e.g., concentrations of high and low tree cover) interpreted as alternative vegetation states. Such observed bimodalities over large spatial extents may not be evidence for alternate states, as they may include regions that have different edaphic conditions and disturbance histories. In this study, we conduct a systematic multi-scale analysis of diverse MODIS data streams to quantify the presence and spatial consistency of alternative vegetation states in Sub-Saharan Africa. The analysis is based on the premise that major discontinuities in vegetation structure should also manifest as consistent spatial patterns in a range of remote sensing data streams, including, for example, albedo and land surface temperature (LST). Our results confirm previous observations of bimodal and multimodal distributions of estimated tree cover in the MODIS VCF. However, strong disagreements in the location of multimodality between VCF and other data streams were observed at 1 km scale. Results suggest that the observed distribution of VCF over vast spatial extents are multimodal, not because of local-scale feedbacks and emergent bifurcations (the definition of alternative states), but likely because of other factors including regional scale differences in woody dynamics associated with edaphic, disturbance, and/or anthropogenic processes. These results suggest the need for more in-depth consideration of bifurcation mechanisms and thus the likely spatial and temporal scales at which alternative states driven by different positive feedback processes should manifest.

**Keywords:** Savanna; alternative stable states; MODIS VCF; land surface temperature; albedo

#### **1. Introduction**

Global distributions of vegetation structure, species composition, and functional composition are governed primarily by climate and phylogenetic processes [1–5]. In the seasonal tropics, precipitation patterns and ecohydrological interactions are the principal drivers of large scale vegetation patterns [6–8], with increasing precipitation generally corresponding to increasing tree density, canopy cover, and tree height [8–10]. While continental-scale average tree density, cover, and height may increase with precipitation, at fine spatial scales, local factors including soil type and hydrological interactions, fire and herbivory, and agriculture and wood harvest lead to considerable spatial heterogeneity in vegetation structure [3,7,11–15]). Thus, landscapes with similar climate are not necessarily similar in vegetation structure (e.g., tree cover, density, and height), and patch mosaics emerge at different spatial scales, reflecting different edaphic conditions and disturbance histories.

More interesting in the context of this paper, however, is the degree to which spatial variability in a landscape is (a) proportional to the frequency, intensity, and time since disturbance (e.g., a harvest event that removes a fraction of tree cover in a particular location) and subsequent regrowth ("linear patch dynamics"), or is (b) amplified by positive feedbacks that lead to alternative state dynamics in response to disturbance events ("bifurcating patch dynamics"). The classic example of an amplifying feedback producing bifurcating patch dynamics (between forest and savannas in mesic systems, and savanna and grasslands in drier savannas) is the so-called fire-trap, where an initial loss (or gain) of tree cover promotes (or suppresses) grass growth, providing more (or less) fuel for fires and increased (or decreased) tree mortality in a continuing cycle of tree loss (or gain) [16–23].

Numerous studies have examined satellite observations in search of empirical support for the prevalence of alternative vegetation states [18,19,24–27] in spatial and temporal data [23,28]. The majority of studies that presented evidence for alternate states were based on interpretation of frequency distributions of tree cover over large spatial regions with similar mean annual rainfall (implicitly assuming similar local edaphic and ecological conditions). Regions following linear patch dynamics would be expected to result in unimodal histograms while, in sharp contrast, regions following bifurcating patch dynamics would result in bi-modal or multimodal histograms, with each mode defining an alternative state caused by the presence of a positive feedback. The earlier studies [18,19,24–27] used tree cover estimates from the MODIS vegetation continuous field product (VCF) [29,30] and found that bimodal and multimodal tree cover distributions were common, particularly in the semi-arid and mesic tropical savanna regions, suggesting the presence of alternative states. In these studies, potential spatial variability in ecological, edaphic, and anthropogenic factors within rainfall zones were mostly ignored. Further, the MODIS VCF retrieval was based on classification and regression tree (CART) algorithms and trained using pre-averaged binning techniques that have been linked to possible artifacts (artificial clumping of predictions) that could be wrongly interpreted as alternative states not present in reality [31–33]. Thus, the true extent and prevalence of bifurcating patch dynamics in the tropical savannas remains unknown.

Remote sensing data are invaluable for our understanding of savanna ecology and the importance of positive feedbacks at landscape, continental, and global scales. However, to advance our understanding of alternative vegetation states, it is critical that we identify not only the specific processes potentially leading to alternative states but also the spatial scale at which the emergent patterns should be detectable and thus the appropriate spatial resolution for data used in its detection. We suggest that the controversy and uncertainty in the extent to which VCF-based analyses have proven the existence of alternative states in tropical savannas can be resolved by systematic comparisons among diverse remotely sensed data streams which, in theory, should also respond to changing tree cover. We posit that bifurcations (i.e., discontinuities) in the frequency distributions of tree cover should also logically be expressed in other remotely-sensed data streams, including albedo and surface temperature. More particularly, not only should the bimodal distributions in the VCF tree cover be detected in other physically based remote sensing variables, but those bifurcations should exhibit spatial consistency (i.e., occur in the same geographic locations). Consistency among diverse data-streams would provide strong evidence supporting earlier conclusions from analysis of the VCF product and remove potential bias associated with any one product. Inconsistencies, on the other hand, would require reevaluation of the earlier results and more in-depth consideration of the

mechanisms of bifurcation and thus the likely spatial and temporal scales at which alternative states driven by different positive feedback processes should manifest. In this work, we examine MODIS derived satellite VCF, albedo, and land surface temperature to determine the degree of similarity in the occurrence of multimodality and inference of alternative stable states across Sub-Saharan Africa (Figure 1). For this analysis, the Sub-Saharan study area includes the African continent from 22◦ North to 35◦ South, including Madagascar but excluding the Arabian Peninsula (Figure 1). We ask the following questions: (1) are the apparent forest-savanna and savanna-grassland bifurcations detected in earlier VCF analysis also present in albedo and surface temperature data?; and (2) to what extent are emergent bifurcations in VCF, albedo, and land surface temperature dependent on the spatial scale of analysis?

We report results of three distinct multi-scale analyses (at the continental scale, based on zones of similar rainfall, and at the landscape scale) designed to quantify the presence of bimodality as a diagnostic of forest-savanna and savanna-grassland bifurcations. Further, we undertake a comprehensive analysis to identify evidence of multiple vegetation states across different data streams [VCF, near-infrared (NIR) albedo, and LST] at the same geographic locations, noting that the independent (i.e., non-VCF) products are derived using fundamental physical principles with relatively little empirical model fitting, while VCF estimates rely on an empirical CART approach. For this analysis, the Sub-Saharan study area includes the African continent south of the Sahara Desert, including Madagascar (Figure 1). The dominant (~66%) ecosystems in this region are tropical and subtropical grassland, savanna, and shrublands.

**Figure 1.** Terrestrial ecosystems [34] over the Sub-Saharan study area.

#### **2. Data**

#### *2.1. Precipitation Data*

We use the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) [35]. CHIRPS is a 30+ year (1981 onwards to current) quasi-global (50◦S–50◦N) rainfall dataset defined at monthly time-steps for 0.05◦ × 0.05◦ spatial grids. CHIRPS combines satellite based precipitation estimates with in-situ station data to create high quality reliable rainfall time series with reduced systematic bias relative to other available global precipitation datasets [35–37], particularly over Sub-Saharan Africa. The mean annual precipitation (MAP) used in this study is computed using 30 years of CHIRPS monthly precipitation data (1981 to 2011; Figure 2).

**Figure 2.** Long-term mean annual precipitation [MAP (mm yr<sup>−</sup>1)] in Sub-Saharan Africa using 30 years (1981–2011) of data from Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS).

#### *2.2. MODIS Vegetation Continuous Fields*

The most recent MOD44B collection 6 VCF data available at a spatial scale of 250 m are used in this study [29,30]. VCF estimates are generated using a CART based algorithm and MODIS optical bands spanning blue to shortwave infrared wavelengths. A suite of metrics derived using annual multi-spectral reflectance observations and vegetation indices are used as predictor variables. The CART is trained using estimates of tree cover, non-tree vegetation, and non-vegetated area observations from field data that are up-scaled using Landsat data [29,30]. The VCF product also includes standard deviation estimates for retrieved tree cover and non-vegetated percent cover along with quality and cloud or water status information. In this study, a spatially explicit VCF map for the study region is created using only the highest quality non-cloudy 2005 data over Sub-Saharan Africa (Figure 3A). The 2005 MOD44B 250 m data are nearest neighbor resampled to a spatial resolution of

1 km to match the resolution of other MODIS data streams (NIR albedo and LST) used in this study using the nearest neighbor method. We use the nearest neighbor approach since averaging pixels tends to reduce spatial heterogeneity and thus suppress identification of alternative vegetation modes.

#### *2.3. MODIS NIR Albedo*

The albedo defines the ratio of radiant energy scattered away from a surface in all directions to the down-welling irradiance incident upon that surface. The albedo over vegetated regions is a function of the local illumination conditions that includes solar geometry and atmospheric conditions, the reflectance properties of the vegetation (leaf optical properties, leaf area index, and distribution of leaf angles), and the underlying soil [38,39]. The MODIS albedo retrieval algorithm first characterizes the surface anisotropy by combining multi-date multi-band atmospherically corrected surface reflectance data to derive the albedo (MCD43B3 V5, 1 km spatial and 16 day temporal resolutions) defined in three spectral regions: visible (VIS: 0.3–0.7 μm), near-infrared (NIR: 0.7–5.0 μm), and shortwave (SWIR: 0.3–5.0 μm) albedo. Changes in the VIS and NIR albedo have also been related to season and vegetation growth [40], while the SWIR albedo is more influenced by non-vegetative components of the land cover. The VIS and NIR albedo are highly correlated; we use the NIR albedo for this study [39], as it is less affected by the atmosphere. The MCD43B3 albedo data include a solar angle specific black sky albedo, solar angle independent diffuse white sky albedo, and associated quality flags. The median NIR albedo map is computed using 45 eight day composites over Sub Saharan Africa for the year 2005. The median (instead of the mean) is used because it is less sensitive to outliers in the data sets (Figure 3B).

#### *2.4. MODIS Land Surface Temperature*

Land surface temperature is a function of both net radiation and the fluxes of latent and sensible heat that are mediated by vegetation cover [41–47]. Observations [46,47] suggest that land surface temperatures depend on the abundance and type of vegetation, with regions with higher vegetation density—such as forests—showing lower surface temperatures than sparsely vegetated surfaces. Hence, surface temperature is expected to reflect underlying tree cover distributions. The observations in the thermal infrared (TIR: 10–13 μm) bands in conjunction with observations in the middle infrared (MIR: 3 μm) are used here to deduce both the land surface temperatures and emissivity using physically based models. The MYD11A2 [48] V6 daytime LST (1 km spatial resolution and eight day composite) includes LST observations and emissivity estimates. LST observations for Sub Saharan Africa in the year 2005 are used to compute the median daytime LST for analysis in this study. The median (instead of the mean) is preferred to reduce sensitivity to outliers (Figure 3C).

**Figure 3.** Illustration of all data streams used in this study. Spatially explicit maps of: (**A**) tree cover estimates from MODerate Resolution Imaging Spectroradiometer (MODIS) vegetation continuous fields (VCF), (**B**) median 2005 near-infrared (NIR) white sky albedo, and (**C**) median 2005 land surface temperature (LST, K).

#### **3. Methods**

Past studies that suggested empirical evidence of alternate vegetation states analyzed tree cover distributions in rainfall zones defined using MAP intervals over large geographic areas (see, for example, the rainfall intervals in Figure 2). Here, we report results of three distinct multi-scale analyses designed to quantify the presence of bimodality as a diagnostic of savanna bifurcations: (1) a non-spatially explicit continental-scale analysis, (2) a non-spatially explicit rainfall zone analysis, and (3) a spatially explicit landscape-scale analysis. Further, we undertake a comprehensive analysis to identify evidence of multiple states across different data streams (VCF, NIR albedo, and LST) at the same geographic locations, noting that the independent (i.e., non-VCF) products are derived using fundamental physical principles with relatively little empirical model fitting, while VCF estimates rely on an empirical CART approach. We also note that while MODIS VCF is known to be unreliable for low tree cover values, distinctions of high and low tree cover may be more reliable [31–33]. Quantitative evidence for the presence or absence of unimodality (i.e., multimodality) can be inferred statistically using Hartigan's dip test [48] on a population of observations [49,50]. This work uses the Hartigan's dip test for populations of observations at multiple scales as evidence for or against unimodality. Hartigan's dip test is sensitive to skew [24,51], thus the VCF data are arcsine transformed to reduce skew. Hartigan's dip test *p* values of less than 0.05 and 0.01 are considered to be moderately and highly significant indicators of multimodality, respectively. The following sections detail the multi-scale analysis approach followed in this work.

#### *3.1. Non-Spatially Explicit Continental Scale Analysis*

For this initial analysis, observations over all of Sub-Saharan Africa are examined both qualitatively and quantitatively. Histograms of each data set (VCF, albedo, and LST) are used to quantitatively infer the presence or absence of multiple vegetation states using the Hartigan's dip test. Scatter plots of each data set with long-term MAP as the dependent variable are also visually inspected for evidence of discontinuities indicative of multimodality.

#### *3.2. Non-Spatially Explicit Rainfall Zone Analysis*

Rainfall zones defined by ranges of MAP at continental scales are often spatially extended and, in some cases, spatially disconnected across the African continent (Figure 1). We examine the histograms of MODIS VCF, albedo, and LST in 200 mm rainfall zones (e.g., 0–200 mm/yr, 200–400 mm/yr, etc.) for all of Sub-Saharan Africa, replicating earlier analyses that used VCF tree cover stratified by rainfall to infer bifurcation dynamics [18,19,24–27]. We anticipate that bimodality detected in VCF data should also be present in the albedo and LST measurements. Hartigan's dip test p values are tabulated for each precipitation zone, and similarities among data sets are used to identify consistent/inconsistent signals of alternative vegetation states.

#### *3.3. Spatially Explicit Landscape-Scale Analysis*

We finally conduct a finer-scale, spatially explicit analysis to detect evidence for bimodality at the landscape level (local spatial extents of 15 km × 15 km). The premise is that bimodality in savanna structure (i.e., tree cover) detected at these local scales will provide more compelling evidence for positive feedbacks and bifurcations occurring within the same local climate and broad edaphic, biotic, and anthropogenic environments. This is in contrast to the continental scale rainfall zone analyses (above) where observed bimodality might reflect much coarser scale differences in soils, topography, and anthropogenic activity not reflected in coarse scale MAP stratifications. Strong evidence for alternative states is inferred if landscapes with multimodality are spatially consistent among independent data sets (i.e., occur in the same geographic locations). Conversely, the lack of spatial consistency might be cause for caution in interpreting apparent bimodality as evidence for the existence of alternative vegetation states in those locations. To test for spatial consistency at finer pixel level

scales, spatially explicit signals of multimodality for each data product for each pixel are inferred using the dip test over its surrounding n × n pixel window (n = 15 × 15 = 225, 1 km2 pixels), ignoring pixels labeled as cloudy and pixels labeled as water in the VCF data. The number of locations that shows evidence for multimodality for each product on its own and in conjunction with other products are tabulated. Spatial consistency between two variables is inferred by generating a confusion matrix and interpreting the Cohen's kappa [52] agreement/disagreement metric.

#### **4. Results**

#### *4.1. Non-Spatially Explicit Continental Scale Analysis*

Figure 3 shows maps of MODIS VCF, NIR albedo, and LST analyzed in this study. All datasets exhibit broadly similar spatial patterns at continental scales, with the drylands distinct from the mesic savannas and forest zones. Spearman's correlation between the MODIS VCF-NIR albedo is moderately high at −0.43, while it is higher for MODIS VCF-LST at −0.71 and LST-NIR albedo at 0.64. However, the histograms (Figure 4) are markedly different for each of the data streams. The MODIS VCF and LST show at least two distinct modes, while NIR albedo appears unimodal visually. The MODIS VCF shows distinct peaks in the frequency of tree cover at <20% and ~80%, while the MODIS LST shows three distinct peaks between 300 K and 315 K. Hartigan's dip test indicates that all histograms show highly significant (*p* value << 0.001) statistical evidence for multimodality.

**Figure 4.** Histograms showing data distributions for (**A**) tree cover, (**B**) NIR albedo, and (**C**) LST. A total of 21,081,051 1 km<sup>2</sup> locations with MAP <sup>≥</sup> 100 mm are used in the histogram analysis, with 20 equal bins spanning the range of each dataset. The red line traces the density (>512 equal bins) for visual identification of multiple modes not easily seen in the histogram. Multimodal distributions are visually evident in the VCF and LST, but less so in the NIR albedo. Hartigan's dip test results indicate highly significant departure from unimodal distribution for all data streams.

Figure 5 shows scatter plots illustrating the relationships between the MODIS data streams and MAP. MODIS VCF shows a distinct bimodality (disconnected islands of red, orange, and yellow at a similar MAP in Figure 5A) between low (<20%) and high (~80%) tree cover, with bimodality for MAP regions >1200 mm/yr. However, such bimodality is not visually evident in the NIR albedo, which shows a broad exponential decline with MAP associated with increasing vegetation cover but without visually distinct discontinuities (Figure 5B). LST shows a decreasing trend with increased MAP (Figure 5C), as more trees tend to reduce the daytime land surface temperature, but the bimodality observed in the tree cover data is not visually evident in the LST data. These observations corroborate past researcher observations of bimodality in VCF, but they also suggest that the discontinuities seen in the VCF tree cover are not readily apparent in other related MODIS data at the 1 km scale.

**Figure 5.** Scatter plots of each data stream with respect to long term MAP. (**A**) VCF tree cover, (**B**) NIR albedo, and (**C**) LST. The density of points is color-coded to show increasing point-density in warmer (red) colors.4.2. Non-Spatially Explicit Rainfall Zone Analysis

Similar to past work [26], we analyze tree cover distributions in 200 mm MAP precipitation intervals, providing 22 MAP zones between 0 and 4400 mm MAP. Hartigan's dip test results applied to all MODIS data over each zone are shown in Table 1. Individual histograms for each precipitation interval are shown in supplemental Figures S1–S3 for reference. It can be inferred from Table 1 that the VCF tree cover shows evidence for multimodality over the entire MAP range. Visual examination of the tree cover data (Figure 5) shows distinct bifurcation between 1200 and 2000. However, the dip test (Table 1) and Figure S1 suggest multimodality over almost the entire range of MAP. The dip test results for LST also indicate the presence of multimodality over most of the range of MAP, while multimodality in the albedo data is indicated in relatively fewer MAP ranges at low and intermediate rainfall.

**Table 1.** Hartigan's dip test for multimodality in three independent remote sensing datasets assessed within 200 mm MAP zones of Sub-Saharan Africa. Statistical significance is color coded, with red highly significant (*p* < 0.01), yellow moderately significant (*p* < 0.05), and gray not significant. These data organized by rainfall zones suggest multimodality is common, however, see text and Figure 4 for analysis of the spatial organization of bifurcation in these datasets.


To further understand the results in Table 1, we show the tree cover, albedo, LST histograms, and associated spatial maps for the regions of the continent with mean annual rainfall between 1600–1800 mm MAP, where all three datasets show evidence for bifurcations (Figure 6). The histogram for the VCF clearly shows a marked bimodality between the highest (~80%) and lowest (<20%) VCF values. However, the potential problems in using such MAP intervals and histograms to analyze for savanna bifurcations is evident from the spatially explicit map in Figure 6. In particular, when we examine these data spatially, it becomes evident that the apparent bifurcations in the VCF histogram occur in distinctly different regions. This suggests that histograms based on rainfall zones are multimodal not because of local-scale feedbacks and emergent bifurcations (the definition of alternative states), but because of regional scale differences in woody dynamics associated with edaphic, disturbance, and anthropogenic regimes.

**Figure 6.** Histograms and maps of (**A**) tree cover (%), (**B**) NIR albedo, and (**C**) LST (K) in the 1600–1800 mm/yr mean annual precipitation zones of Africa, showing that bimodality seen in the tree cover datasets occurs in spatially disconnected regions. While the histograms in this MAP interval are visually and statistically bimodal, the apparent bifurcations occur in spatially distant locations on the fringes of the moist tropical forests of West Africa (e.g., the low tree cover mode in red) and the fringes of the Congo Basin (e.g., the high tree cover mode in green). Similar separation of low and high tree cover modes is apparent between eastern and western Madagascar. The LST and NIR albedo multimodality, while much less pronounced than in the tree cover data, is also geographically separated between West Africa and Congo Basin locations.

#### *4.2. Spatially Explicit Landscape-Scale Analysis*

Figure 7 shows results of our analysis of multimodality within 15 × 15 km (i.e., 225 pixels) moving windows across all of Africa, with red and yellow colors showing locations with statistical evidence for multimodality assessed using the Dip test. The VCF data have the highest number of locations across most of Sub-Saharan Africa, with evidence for bimodal or multimodal distributions in most areas at this scale. At this scale, the occurrence of multimodality in the albedo and LST data streams is much less common, with a general lack of spatial consistency between VCF, albedo, and LST datasets (Table 2). Quantitatively, the percent agreement of multimodality among all three data

streams (VCF, LST, and albedo) conveys that only 0.52% of the locations show multimodality in the VCF data.

Analysis of Cohen's kappa (Table 2) quantifies this apparent lack of spatial consistency in the locations where multimodality is observed. The exclusion of VCF to examine the degree of spatial coherence between LST and albedo marginally improves the spatial consistency (Table 2), but there is still considerable disagreement between LST and albedo in the locations with multimodality.

**Figure 7.** Spatially explicit signals of local bifurcations (non-unimodality) calculated within window sizes of 15 × 15 km for data shown in Figure 3, (**A**) tree cover, (**B**) albedo, (**C**) LST. Red and yellow colors show 15 × 15 km landscapes with high and moderate statistical evidence for local bifurcations (i.e., Hartigan's unimodality dip test is rejected in these locations).

**Table 2.** Tests for spatial consistency in detections of multimodality among three remote sensing datasets, showing the number of 15 km × 15 km locations across Sub-Saharan Africa with evidence for multimodality individually (columns 2–4) and in conjunction with other products (columns 5–8, where, for example, VCF and LST quantify locations scored as multimodal in both VCF and LST). Only locations that have a MAP ≥ 100 mm are considered. Statistical significance is based on Hartigan's dip test *p* value < 0.05. The Cohen's kappa (κ), a measurement of agreement between data sets, is given in parenthesis. Kappa values close to unity indicate strong agreement, while values close to zero imply low or poor agreement, and negative values imply strong disagreement. The percent agreement of multimodality among all data streams over the entire continent is also tabulated.


#### **5. Discussion**

In this study, we examine the distributions of MODIS VCF tree cover estimates in Sub-Saharan Africa alongside data on albedo and LST to determine the degree of similarity among datasets in the occurrence of multimodality and inference of alternative stable states. While we argue that albedo and LST should be closely correlated with tree cover, we also recognize that the darkening and cooling effect of trees in a landscape may also occur with dense herbaceous vegetation, potentially reducing the effectiveness of both NIR-albedo and LST as indices of tree cover.

The above caveat notwithstanding, our results (Figures 4–7) show clearly that bimodality is common in MODIS VCF when visualized at a continental scale, when sorted into MAP bands (consistent with past studies, [18,19,24–26]; Table 2), and when analyzed at landscape scales (15 × 15 km blocks; Table 2). However, analysis of NIR albedo and LST datasets suggest multimodality is much less common in these datasets. Examination of the spatial distribution of high and low tree cover within an example MAP interval (Figure 6) highlights the potential dangers in the approach. In particular, while Figures 4 and 5 suggest strong bimodality in tree cover, examination of the spatial distribution of the data reveals that the

high and low modes are geographically separated, indicating regional differences in tree cover related to possible edaphic or anthropogenic impacts on mean tree cover rather than positive feedbacks giving rise to bifurcations at local scales under similar edaphic conditions.

Past studies presenting empirical evidence for alternative states at different levels of MAP did not examine independent data sets, nor did they consider the spatial organization of tree cover modes within each MAP band. Based on the conclusion that detection of multimodality in spatially extensive MAP intervals may not be appropriate for the diagnosis of bifurcations at regional scales, we carry out spatially explicit analyses designed to detect bifurcations at landscape scales (15 × 15 km) in the three independent data streams. Multimodality is suggested across most of Sub-Saharan Africa in the VCF dataset but exists in much more constrained regions in the NIR albedo and LST datasets. Metrics of spatial consistency in multimodality among data-streams based on Cohen's kappa (Table 2) indicate that VCF is distinct from LST and albedo (i.e., very little spatial consistency), with more agreement (but still considerable differences) in the locations of bimodality in the albedo and LST data streams. These results support concerns raised elsewhere [31–33] that the VCF product may not be appropriate for the diagnosis of alternative states (i.e., bifurcations) because of the calibration approach utilized and the discontinuities inherent in CART-based model-fitting approaches.

#### **6. Conclusions**

We posed the following two questions: (1) are the apparent forest-savanna and savanna-grassland bifurcations detected in earlier VCF analysis also present in albedo and surface temperature data?; and (2) to what extent are emergent bifurcations in VCF, albedo, and land surface temperature dependent on the spatial scale of analysis? Our results show that apparent bifurcation in the VCF data, diagnosed via statistical tests for multi-modality, far exceeds bifurcation detections in albedo and land surface temperature data. Further, our multi-scale analysis highlights that the "rainfall zone" analysis approach adopted by earlier authors (i.e., aggregating continental data using rainfall bins; Section 4.2) tends to identify multimodality associated with regional differences in tree cover rather than detecting landscape-scale bifurcations. Our explicit landscape scale analysis finds that, while apparent bifurcations are ubiquitous in the VCF data, they are far less common in both the albedo and the land surface temperature data.

The lack of consistent signatures of bimodality in satellite estimates of tree cover, land surface temperature, and albedo highlights the need for caution in analyses suggesting that alternative states are widespread in tropical savannas. A more considered review of the associated ecological theory points to several explanations for this difference of opinion: (1) while positive feedbacks relating to fire, tree-grass competition, and herbivory can—in theory—lead to the emergence of areas of high and low tree cover under similar climate [12,18,53], in biologically and edaphically complex landscapes, the processes potentially leading to bifurcations may be buffered by other interactions [54]; (2) in connected landscapes, spatial interactions may reduce the climate space where bifurcations emerge, making bifurcations less common than expected based on simple non-spatial models [18,55]; and (3) the positive feedbacks potentially driving bifurcations in tree cover are processes that operate at distinct spatial scales that may not be captured by the relatively coarse grain of many satellite tree cover datasets. For example, while facilitation and competition among plants in arid and semi-arid systems is a well-known cause of the patterned vegetation states known as tiger-bush (i.e., stripes of dense vegetation separated by stripes of bare soil [56]), the processes involved (i.e., redistribution of water and root competition) occur at relatively fine spatial scales (<~100 m). We should therefore not anticipate being able to detect this type of alternative vegetation state in data with a spatial grain >~50 m. While remote sensing data are invaluable for our understanding of savanna ecology at landscape, continental, and global scales, to advance our understanding of alternative vegetation states, it is critical that we identify not only the specific processes potentially leading to alternative states but also the spatial scale at which the emergent patterns should be detectable.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/11/7/815/s1, Figure S1: Histograms of MODIS tree cover (VCF) over regions stratified by mean annual precipitation (MAP in mm/y). Each histogram was generated using 20 equally spaced bins spanning the range of the data in each region Red lines trace the kernel density estimates to help visualize the shape of the histograms., Figure S2: Histograms of MODIS NIR albedo over regions stratified by mean annual precipitation (MAP in mm/y). Each histogram was generated using 20 equally spaced bins spanning the range of the data in each region Red lines trace the kernel density estimates to help visualize the shape of the histograms, Figure S3: Histograms of MODIS Land Surface Temperature over regions stratified by mean annual precipitation (MAP in mm/y). Each histogram was generated using 20 equally spaced bins spanning the range of the data in each region Red lines trace the kernel density estimates to help visualize the shape of the histograms.

**Author Contributions:** S.S.K., N.P.H. and L.P. conceived research, S.S.K. conducted the analyses with inputs from J.A., C.W.R., W.J. and B.M.L. All authors contributed towards writing the manuscript.

**Funding:** Research was supported by the US National Aeronautic and Space Administration (NASA) as part of the NASA Carbon Cycle Science program (Grant # NNX17AI49G).

**Acknowledgments:** We thank the editor and the two anonymous reviewers for their insightful suggestions that greatly improved our manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Healthy Park Needs Healthy Vegetation: The Story of Gorongosa National Park in the 21st Century**

#### **Hannah Herrero 1,\*, Peter Waylen 2, Jane Southworth 2, Reza Khatami 2, Di Yang <sup>3</sup> and Brian Child <sup>2</sup>**


Received: 6 December 2019; Accepted: 29 January 2020; Published: 3 February 2020

**Abstract:** Understanding trends or changes in biomass and biodiversity around conservation areas in Africa is important and has economic and societal impacts on the surrounding communities. Gorongosa National Park, Mozambique was established under unique conditions due to its complex history. In this study, we used a time-series of Normalized Difference Vegetation Index (NDVI) to explore seasonal trends in biomass between 2000 and 2016. In addition, vegetation directional persistence was created. This product is derived from the seasonal NDVI time series-based analysis and represents the accumulation of directional change in NDVI relative to a fixed benchmark (2000–2004). Trends in precipitation from Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) was explored from 2000–2016. Different vegetation covers are also considered across various landscapes, including a comparison between the Lower Gorongosa (savanna), Mount Gorongosa (rainforest), and surrounding buffer zones. Important findings include a decline in precipitation over the time of study, which most likely drives the observed decrease in NDVI. In terms of vegetation persistence, Lower Gorongosa had stronger positive trends than the buffer zone, and Mount Gorongosa had higher negative persistence overall. Directional persistence also varied by vegetation type. These are valuable findings for park managers and conservationists across the world.

**Keywords:** remote sensing; vegetation dynamics; vegetation persistence; conservation; savannas; Africa

#### **1. Introduction**

Savannas are commonly defined as grassland with scattered trees, but in practice are a mix of grass, shrub, and trees [1–3] which comprise over 55% of the southern African landscape, and approximately 20% of Earth's landcover [1]. Savannas are controlled by key factors: precipitation, herbivory, fire, and humans [4–7]. Savannas are water-limited systems, therefore precipitation is the main driving factor for the ecosystem [5,8]. Up to a precipitation threshold of 750 mm/year, grasses dominate, between 750-950 mm/year, more mixed systems occur, and above 950 mm/year, woodlands dominate [5,7–9]. Decreases in precipitation alone can lead to declines in vegetation, and when coupled with rising temperatures the resulting lower availability of soil moisture causes even further declines in vegetation [10–12]. Savannas are an important conservation landscape in southern Africa. Understanding and conserving these systems is vital as they support high faunal and floral biodiversity. Savannas are also an important biome affecting the carbon cycle, supporting high human populations, and contributing 14% of global net primary productivity [8,13–17]. Conservation areas are a key

resource in southern Africa both in terms of their ecological importance (biodiversity, biomass, carbon, etc.) and as socioeconomic drivers for the surrounding communities [18].

Remote sensing of these savanna systems is a challenge due to the heterogeneity of savanna landscapes, but given their ecological and socioeconomic importance, they are a target for development of innovative image analysis techniques. One such approach to evaluating savannas with imagery includes continuous time series analysis [19,20]. The normalized difference vegetation index (NDVI) can be used as a proxy for vegetation biomass, health, and abundance [20–26]. The value of this vegetation index may be enhanced with new time series measures to provide better understanding of trends in vegetation health over time and space [19]. One such measure, derived from NDVI, is that of vegetation persistence, which has demonstrated its usefulness in drylands globally [20]. This is an NDVI-based time series approach to compare vegetation greenness in a given season and year to a baseline period. This methodology allows the user to detect inter-annual changes of NDVI in a spatially explicit manner and identify key trends on the landscape, with a statistical significance attached to each pixel [19,20].

Gorongosa National Park (GNP) in central Mozambique contains two important ecosystem types—savanna and montane rainforest. This paper examines the complex savanna system of the lower part of GNP (herein referred to as "Lower Gorongosa") and compares it to the montane rainforest of Mount Gorongosa [27]. Like many mountainous areas [28], the latter provides a haven for biodiversity much of which remains to be discovered, as evidenced by the extraordinary discoveries from the rapid biodiversity collection event held by E.O. Wilson in 2011 [29]. GNP is chosen for analysis because of the interesting conditions under which it was established and also due to the fact that it has been restored in the last decade. GNP was originally established in 1920 under Portuguese colonial rule, and the park served as headquarters for both sides at some point during both the War for Independence and the Civil War, which lasted from 1964 to 1994. As a result, between 90% and 99% of all mammals in the park were extirpated by crossfire, lack of sustenance, and poaching, and this lack of mammals over time may affect vegetation [30]. Following peace, a variety of illegal hunting activities continued for years in the park, pressuring already drastically decreased populations. In 2004, a relationship developed between the Mozambican government and an American NGO initiating the Gorongosa Restoration Project, the first conservation collaboration of its kind in the world. This project has three main goals: science and conservation, sustainable tourism, and human development. These goals include the restoration of the ecosystem to a pre-war baseline of wildlife [31]. A significant management presence began in the park in 2008 with the aim to restock animals, which might then affect vegetation. The park is a natural experiment for the impact of herbivores on vegetation. Several studies note that the absence of large herbivores enables trees to attain maturity here in GNP, unlike other places in southern Africa where bush densification has been the trend [32,33]. The primary focus of restoration has been wildlife, however, a healthy park also needs healthy vegetation to sustain its wildlife populations. Therefore, the purpose of this paper is to evaluate trends in vegetation biomass during pre- (2000–2008) and post- (2009–2016) management periods in GNP.

The objective of this study is to understand the trends in vegetation in Lower Gorongosa and to compare them to the buffer zone and Mount Gorongosa. The questions we addressed are: 1. Is there a difference in the precipitation between the period since restoration management has been imposed and the prior period? 2. Has there been a trend in the seasonal precipitation time period between 2000 and 2016? 3. Can trends in vegetation density and biomass be detected in the period 2000 to 2016? 4. Does the vegetation directional persistence methodology detect significant changes in NDVI as an indicator of vegetation change, in the period 2000 to 2016?

#### **2. Materials and Methods**

#### *2.1. Study Area*

The GNP occupies approximately 4000 km<sup>2</sup> in central Mozambique (18.2◦S and 34.0◦E) (Figure 1). Annual minimum and maximum temperatures on average range between 15 ◦C and 30 ◦C, in the dry and wet seasons respectively. The ecosystem of the original area of the park, Lower Gorongosa (LG), is markedly different to that of the more newly incorporated Mount Gorongosa (MG), which lies to the northwest at 700 m elevation. Based on the park management report, LG contains various savannas ranging from "open" savanna (floodplain grassland, 21% of LG) in the central part of the park, to "mixed" savanna (grass-shrub-open tree, 44% of LG), and "closed" savanna (thick woodlands dominated by several different tree species, 35% of LG). Montane rain forest vegetation dominates the flanks of MG, while short grassland/rock covers the summit. The protected area on MG is critical to the health of both LG and people living in the surrounding area, as it constitutes the dominant source of stream flow. The areas surrounding parks, known as buffer zones, provide a comparison to vegetation trends occurring inside the park because they are not as protected as parks, but have similar environmental conditions affecting them. They therefore help better understand how the park is functioning in the larger context of the landscape. Based on literature indicating that a buffer zone analysis should encompass approximately 3–4 times the size of the park [34], a 20 km buffer zone around the park was considered in this work. There is some limited human activity in this zone where communities live. Six study areas are addressed: (1) LG (first as a whole), and then further divided into land cover types- (2) open savanna, (3) mixed savanna, and (4) closed savanna, (5) a 20 km buffer zone around LG and (6) MG.

**Figure 1.** Study area map of low Lower Gorongosa, Mount Gorongosa and buffer area shown on a true color composite Landsat image. In the top inset map, open savanna, mixed savanna, and closed savanna are shown.

#### *2.2. Vegetation Change Characteristics*

This section addresses factors used to measure vegetation change and related drivers.

#### 2.2.1. Mean NDVI

The Normalized Difference Vegetation Index (NDVI) is used to establish trends in vegetation health and density over time. The metric is a band ratio between red (R) (highly absorptive in healthy vegetation) and near infrared (NIR) (highly reflective in healthy vegetation) energy [35]:

$$NDVI = \frac{NIR - R}{NIR + R} \tag{1}$$

NDVI ranges from −1 to +1. Higher values indicate the presence of higher vegetation density, and values below 0 indicate no vegetation [35]. NDVI composite products from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, i.e., MOD13Q1.006 Terra Vegetation Index, were used in this research due to their fine temporal resolution. The data span March 2000 to November 2016, have a spatial resolution of 250 m, and are a 16-day composite of the highest NDVI value per pixel. The data are atmospherically corrected and masked for clouds, cloud shadows, water, and heavy aerosols. Furthermore, as water produces negative NDVI values, the Lake Urema system in LG is also masked out and removed from further analysis. Data are aggregated into seasonal composites based on per-pixel seasonal mean values for the following seasons; December, January, February (wet), March, April, May (wet), June, July, August (dry), and September, October, November (dry). Data are further aggregated by mean pixel value for the seasons into different periods corresponding to different on-the-ground management practices in GNP: 2000–2004 (no management), 2005–2008 (management contracts negotiated), 2009–2012 (beginning of on the ground management), and 2013–2016 (fully implemented on the ground management). These can be grouped broadly into 2000–2008, pre-management and 2009–2016, post-management [29]. The analysis is separated into each of the six study areas. Spatial extents of the three savanna types in LG (open, mixed, and closed) were provided by the park management team.

#### 2.2.2. Seasonal Precipitation Totals

As savannas are water limited systems, monitoring trends in precipitation is critical in understanding changes in vegetation. Estimates of monthly precipitation totals from 2000 to 2016 were extracted from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) gridded dataset (0.05◦ resolution). Given the coarse spatial resolution of the data and the relatively small study area, a mean value for each time period was extracted over the whole study area [36]. Total seasonal precipitation for the park is created for each year from the monthly estimates observations. The annual migration of the Intertropical Convergence Zone (ITCZ) is considered to be the principal driver of intra-annual precipitation [37,38]. Given the delay between the arrival of precipitation and green up of vegetation, seasonal estimates of vegetation are lagged by one month [39]. Therefore, the seasonal precipitation-vegetation pairs are respectively as follows: (NDJ-DJF; FMA-MAM; MJJ-JJA; ASO-SON). Mean annual precipitation (2000–2016) is 1075 mm, and individual annual totals range from 665 mm in 2015 to 1415 mm in 2001.

Analysis of seasonal precipitation totals proceeds in several ways. Totals are expressed as being above or below the long-term (2000–2016) mean and the subsequent counts of the number of above/below mean years in each time period recorded. To investigate if total seasonal precipitation has increased/decreased between the two management periods, a hypergeometric test was used. The hypergeometric distribution determines the probability of experiencing a given number of years with values above (or below) mean characteristic at random [40]. It therefore highlights whether the observed number of years (above/below mean) is significantly greater than/less than, would be expected at random. Years during which a seasonal total fall above the mean are considered "successes" and during which it falls below the mean are considered "failures." The probability distribution of the number of successes (failures) is given as:

$$p(\mathbf{x}) = \frac{\binom{k}{\mathbf{x}} \binom{N-k}{n-\mathbf{x}}}{\binom{N}{n}} \tag{2}$$

where *N* is the number of years in the record (i.e., *N* = 17 during 2000–2016), *n* is the sample size (the number of years in a particular period during which a set of management practices is followed), and *k* is the total number of observations in the entire record with the desired property (success or failure) [40]. The sample time periods for the hypergeometric distribution was chosen based on the main management time frame of 2000–2008, where there was not a management presence on the ground, and 2009–2016, where there was a heavy management presence on the ground.

Precipitation values were also graphed as a time series, normalized for seasonality, anomalies were determined, and linear trends were added. Precipitation is also quantitatively compared to mean NDVI, as savanna vegetation-precipitation association is well established [5,7–9,12,20–26].

#### 2.2.3. Vegetation Directional Persistence

Newly developed directional persistence (DP), a continuous time series metric, is calculated on a per-pixel basis and is the cumulative direction of change over time in NDVI calculated in comparison to a benchmark value [19,20]. This approach illuminates trends in vegetation over time and space [19,20]. We implemented the DP analysis using the time-series data from the MOD13Q1.006 Terra Vegetation Index products outlined above in Section 2.2.1. To calculate DP, per-pixel seasonal mean NDVI over the period 2000–2004 was calculated as the benchmark. Then, all subsequent seasonal values 2005–2016 (*N* = 12) were compared to it at a pixel level, generating the equivalent of a random walk in each pixel. Under the null hypothesis of persistent vegetation, i.e., no change in vegetation, there is an equal chance of an observation above the benchmark (success) or below it (failure) [19]. For a "success" in a particular year the pixel receives a value of +1, if a "failure", it receives a −1. These values are accumulated through time, resulting in a net value of DP in each pixel. In the case of a N of 12, a pixel value of +12 would indicate that every year had a NDVI value higher than the benchmark. Likewise, a pixel value of −12 would indicate that every year had a NDVI value lower than the benchmark. The DP values can be assigned to statistical significance levels based on the hypergeometric distribution. Given the preliminary nature of this exploration and the comparatively short (2005–2016) period of evaluation, a 0.10 level of significance, or critical DP of +/−6, is chosen for testing. The respective mean directional persistence, and the percentage of pixels returning significant (positive or negative) persistence are calculated for each of the six study areas.

#### **3. Results**

#### *3.1. Precipitation Anomolies and Trends*

Figure 2 highlights trends in seasonal precipitation, with the blue and red bars indicating the years with precipitation above and below the overall mean, respectively. Each season except ASO, returns approximately equal numbers of years above and below the long-term mean. However, it appears that there are more below-mean years later in the time series. Additionally, the magnitude of the negative anomalies are high compared to those of the positive anomalies. There were some periods of consecutive seasonal rainfall below the mean, such as in ASO from 2004 to 2010, however application of the hypergeometric test revealed no statistically significant differences at the 0.1 level in the number of positive and negative anomalies in any season or time period. This means that the number of years with rainfall above or below the mean were not statistically different from what we might expect at random when we compare the pre and post management periods.

**Figure 2.** Seasonal precipitation anomalies 2000–2016 for November, December, January (NDJ), mean 550 mm; February, March, April (FMA), mean 384 mm; May, June, July (MJJ), mean 80 mm; August, September, October (ASO), mean 61 mm. The numbers in the graphs are count data for above and below mean precipitation years in the time period.

The time series in Figure 3 suggests a decreasing trend in precipitation over time, particularly in the latter part of the time series. In Figure 3a, when seasonal anomalies are represented by their standard deviates, there is a statistically significant decline over this time period at the 0.1 level of significance. In Figure 3b, where seasonal totals were separated out, there was a statistically significant decline over time in February, March, April (FMA) and May, June, July (MJJ) at the 0.05 significance level.

**Figure 3.** *Cont.*

**Figure 3.** Precipitation time series 2000-2016 for (**a**) continuous raw seasonal totals, seasonal anomalies (mm), and seasonal anomalies expressed as standard deviates (z) and (**b**) separated seasonal anomalies for November, December, January (NDJ), February, March, April (FMA), May, June, July (MJJ), and August, September, October (ASO).

#### *3.2. Changes of NDVI*

Mean NDVI values computed by season (MAM, JJA, SON, DJF) and time period (2000–2004, 2005–2008, 2009–2012, and 2013–2016) are mapped across the LG and MG (Figure 4). This clearly shows the role of seasonality, with higher NDVI values in the wet season than in the dry season. Even within the wet/dry seasons the impact of soil moisture accumulation (MAM is greener than DJF), and depletion (JJA is greener than SON) is apparent. Each season also appears to have some variability of NDVI during the time periods. In general, the first period (2000–2004) had the highest NDVI across the study areas, the most recent period (2013–2016) had the lowest NDVI across the study areas, and the two periods between (2005–2008 and 2009–2012) were variable.

Daily NDVI values, extracted by pixel from MODIS data from 2000-2016, are examined to give further insight into the vegetation dynamics of LG, revealing several important findings, including notable shifts in both the magnitude and timing of anomalies in each period (Figure 5a). There was an overall decline in dry season NDVI after the first time period (2000–2004), which may have even greater negative effects on flora and fauna as this is already the most difficult part of the year for them to survive. In Figure 5a, during the initial period (2000–2004), the data supports that the landscape was greener than the long-term mean during the dry season. During the second period (2005–2008), which was a time of less precipitation, the data shows that it was much less green during the dry season than the long-term mean. During the third time period (2009–2012), NDVI was greener than the long-term mean during the wet season. In the fourth time period (2013–2016), it was less green in the first half of the wet season, and greener during the second half of the wet season than the long-term mean. Figure 5b highlights the magnitude of the positive and negative anomalies (seen in Figure 5a) in both absolute number and percentage.

**Figure 4.** Mean NDVI in Lower Gorongosa and Mount Gorongosa by season: March, April, May (MAM), June, July, August (JJA), September, October, November (SON), December, January, February (DJF) and time periods: 2000–2004, 2005–2008, 2009–2012, 2013–2016.

In order to identify any shifts in the timing of the highest and lowest values of NDVI, the dates during which the seasonally averaged NDVI exceeds (or falls below) the value of the upper (or lower) 10% of long-term NDVI values are recorded (Figure 5c) as an objective indicator of green-up and brown-down. Both the upper and lower limits have moved later in the season, indicating that shifts in climate (precipitation coming later) are already noticeable across this region. Several differences emerge in comparison to long-term characteristics: (a) high levels of greening are attained earlier in the first and second periods (by as much as 30 days), and later in the fourth period (by about 10 days); (b) high levels of NDVI persist until the end of May, except in the second period when declines are noted about half a month earlier; (c) the timing of the lowest values of NDVI show considerably less inter-period variability than those of high levels; (d) the lowest values of NDVI are generally experienced in mid-September, starting slightly earlier in the first and second periods; and (e) the landscape starts to return NDVI values higher at the end of October in all periods. Post management (2009–2016) there seems to have been less variability in these dates (Figure 5c). Shifts in season length are important for future vegetation response and health, with related implications on wildlife and people on this landscape. Specifically, they can have very real consequences for people's livelihood and wildlife health.

**Figure 5.** *Cont.*

**Figure 5.** Daily NDVI (2000-04, 2005-08, 2009-12, 2013-16) of Lower Gorongosa (**a**) positive and negative anomalies compared to the long term mean (**b**). NDVI anomalies displaying the frequency of anomaly both as an absolute count and as a percentage of days exhibiting negative and positive anomalies from the long term mean in each of the four periods (left to right and top to bottom 2000–2004, 2005–2008, 2009–2012, 2013–2016), and histograms of the distribution of the magnitudes of those anomalies. (**c**) The dates on which the uppermost (0.732) and lowermost (0.343) 10% of NDVI values occurred during the entire period of study and period of interest. The day count on the *x*-axis starts on 1 March, given our seasons for the year start with March–April–May, then June–July–August, then September–October–November, then December into January–February of the following year.

Time series of mean NDVI across the six study areas highlight several important findings. All data reveal decreasing trends, though not all are statistically significant (Figure 6). There are also some similarities between the LG and the 20 km buffer zone. However, the rate of decline within the buffer zone is steeper than that of LG, suggesting that the conservation and management practices within the park have been effective in mitigating some potential decline [34].

The overall higher NDVI exhibited by MG than LG is attributable to the differing land cover. LG contains varying types of drier savanna, and MG supports montane rainforest. There is a more pronounced and significant decline in NDVI on MG (Figure 6a). This may be due in part to the rainforest of MG supporting higher NDVI and the impacts of human-driven processes such as deforestation being particularly pronounced. Given the land cover types, there is also much less interannual variability in MG when compared to LG. Regardless, the same trends of declining NDVI and precipitation are seen across this landscape.

NDVI was also graphed by land cover type within LG and MG. All land cover types were clearly separable by mean NDVI values. The lowest NDVI values occur in the open savanna vegetation class (0.36–0.69). These increase to 0.41-0.76 in mixed savanna and 0.51–0.81 in closed savanna. The highest NDVI occurs in the heavily forested landscape on MG (0.62–0.82) even though there is a significant decline over time (Figure 6a). Given that these three savanna land cover classes within LG are respectively dominated by grassland, shrubland, and woodland, this separation of vegetation abundance (NDVI) makes sense biologically. There is a significant decline in the NDVI of open grassland savanna in LG (Figure 6b) and greater interannual variability than in closed woodland because grasslands are more dependent, in the short term, on precipitation than woodlands are. The smallest and non-significant decline is found in the closed savanna land cover class in LG (Figure 6c) as woody plants tend to be hardier/less vulnerable, which is also supported by the literature [5,8,20,41].

**Figure 6.** *Cont.*

**Figure 6.** Time series of mean NDVI seasonally (December–January–February (DJF); March–April–May (MAM); June–July–August (JJA); September–October–November (SON)) from 2000–2016 for (**a**) Mount Gorongosa (significant trend), (**b**) open savanna (significant trend), (**c**) and closed savanna (no significant trend) other non-significant trends (in mixed savanna, LG, and the buffer) are not shown.

In Figure 7a, seasonal precipitation and lagged NDVI are shown. Dry season NDVI is much more sensitive (steeper slope) to small increments in precipitation than the rainy season NDVI. In Figure 7b, there is a significant relationship between precipitation and vegetation in terms of interannual variability (e.g., if there was a wetter than normal wet season versus an unusually high NDVI). This relationship was significant across all seasons of comparison (precipitation to NDVI respectively): FMA to MAM (r-squared = 0.466, slope = 0.000416, and is significant at the 0.01 level); MJJ to JJA (r-squared = 0.146, slope = 0.002981, and is significant at the 0.10 level); ASO to SON (r-squared = 0.570, slope = 0.005185, and is significant at the 0.01 level); and NDJ to DJF (r-squared = 0.247, slope = 0.000441, and is significant at the 0.05 level), thus supporting the notion of savanna as a water limited system and precipitation as the main driver of vegetation growth [5].

**Figure 7.** Precipitation and lagged NDVI (respectively) (**a**) plotted by season together with different colors for: February, March, April (FMA) to March, April, May (MAM); May, June, July (MJJ) to June, July, August (JJA); August, September, October (ASO) to September, October, November (SON); and November, December, January (NDJ) to December, January, February (DJF) (**b**) by season individually with fitted regression lines.

#### *3.3. Tendencies in Vegetation Directional Persistence*

Maps of DP for these study areas (Figure 8) show declines in vegetation, particularly during the dry season and in central LG (which is the open savanna cover, Figure 1) relative to the initial baseline period. This area corresponds with open savanna (grassland/floodplain) land cover. There is a dominating negative trend in vegetation persistence during the dry season, particularly in the central part of LG and on MG. During the wet season fewer significant values (positive or negative) of DP prevail. Again similar tendencies are found both inside LG and in the buffer area, though the buffer area has greater negative values of DP.

**Figure 8.** Directional persistence (2005–2016 compared to initial baseline period 2000–2004 at a pixel level) in Lower Gorongosa (including the three savanna types, open, mixed, and closed), Mount Gorongosa, and the buffer zone highlighted by season (MAM = March, April, May, JJA = June, July, August, SON = September, October, November, DJF = December, January, February).

The seasonal properties of mean directional persistence of all six study areas are presented in Figure 9. Across all study areas, the absolute value of mean persistence is larger during the dry season (more negatively deviated from zero), than the wet season (less positively deviated from zero). This suggests the tendency for NDVI in the dry season to be more negative through time than the baseline value to a greater magnitude than in the wet season. When comparing LG, the buffer zone, and MG, LG had the most positive/least negative mean DP values across seasons. This indicates that, of the six study areas, LG had higher percentages of positive vegetation persistence over time. MG with negative mean DP values across all seasons, shows the most marked decline in vegetation persistence over time. Within the three land cover types in LG, closed savanna had the most positive/least negative

mean DP values, whereas open savanna had the most negative DP values (Figure 9). OS appears to suffer most in terms of vegetation declines, as found elsewhere in the literature.

**Figure 9.** Mean directional persistence values for Lower Gorongosa (LG), the buffer zone (B), Mount Gorongosa (MG), and the three Lower Gorongosa land cover types: open savanna (O), mixed savanna (M), and closed savanna (C) by season (December-January-February (DJF); March-April-May (MAM); June-July-August (JJA); September-October-November (SON)).

Percentage significant positive and negative DP (+/−6) is also calculated for each study area (Figure 10). LG has higher percentages of significant positive DP and lower percentages of significant negative DP than the buffer zone. MG returns higher percentages of significant negative DP and lower percentages of significant positive DP than any other study area. In the dry season, as much as 70% of the MG landscape has a significant negative DP value. Of the three LG land cover types, open savanna exhibits the highest percentages of significant negative DP. In the dry season, as much as 75% of the OS landscape has a significant negative DP value. OS also has the lowest percentages of significant positive DP across all seasons.

**Figure 10.** Percentages of significant directional persistence (positive and negative) across the landscapes for: Lower Gorongosa (LG), the buffer zone (B), Mount Gorongosa (MG), and the three Lower Gorongosa land cover types: open savanna (O), mixed savanna (M), and closed savanna (C) by season (December–January–February (DJF); March–April–May (MAM); June–July–August (JJA); September–October–November (SON)).

#### **4. Discussion**

This study evaluates the relationship between savanna vegetation and precipitation, the role of protected areas in vegetation conservation, and the importance of distinguishing vegetation type for evaluation of its health/persistence. We found that there was a decline in precipitation in Gorongosa National Park over the 2000 to 2016 period, which could be linked to a decline in NDVI. However, we found that there were differences in decline of NDVI amongst savanna land cover classes within LG, with OS having the steepest decline. Compared to the buffer zone, inside the protected area of

LG had less vegetation decline. MG had the greatest decline of all six study areas in terms of NDVI, as well as the highest negative DP. This may be due to more localized drivers outlined below, but with further conservation efforts, we hope to see more positive trends in the future.

Degradation has been of concern for several decades in savannas [42,43]. Degradation can include several phenomena, including a decrease in biomass of vegetation (seen in this study in and around GNP) and bush encroachment. Bush encroachment is a type of degradation that results in an increase in lower quality short woody plants [12,44,45]. Important drivers of savanna vegetation and therefore potential degradation include, decline in precipitation (already seen here), changes in herbivory (increases to a certain population level may mitigate degradation, but beyond that level, may actually exacerbate it), decreased fire frequency, and humans [4–7,42–47]. Under continued climate change and rewilding, these systems will need to be carefully monitored for further degradation [15,48]. One type of bush encroachment reported in the park, particularly in the vulnerable floodplain grasslands, is that of *Mimosa pigra*. This is due to the previous absence of herbivory (almost no herbivores survived the civil war, allowing encroachers to increase). Now with the heavy reintroduction of herbivores, the park is seeing a decrease in these woody encroachers [84]. However, given the relatively recent rewilding, our longer-term study is seeing the trend of degradation in OS. Though a constant fire regime (the fire regime remains consistent pre- and post-war) has enabled large trees to be able to survive and thrive here [29,32,33]. A common way to mitigate the effects of some drivers on the savanna protected area landscape is to implement buffer zones. Buffer zones are areas surrounding parks that act as a more neutral space for activity than the broader landscape. In buffer zones there may be some human activity, but generally it is limited [49]. Creation of these areas increases the distance between protected zones and potential disturbances, limits edge effects, and thereby increases the conservation value of the protected area [50,51].

Jung et al. (2010) [11] point to a global decline in precipitation leading to reduced soil moisture and evapotranspiration. This study found some decrease in precipitation in this region of southeastern Africa during the 21st century [10,44], linked to a decline in vegetation health [10–12], [52,53]. The observed decline in vegetation health (using NDVI as a proxy) further strengthens the link between vegetation biomass in savannas and its main driver, precipitation [5,7–9,12,20–26]. Declines in precipitation and vegetation health suggest that management has not yet had a significant impact on vegetation [10,53] and is still subordinate to the effects of precipitation.

Global greening has been detected over recent decades by remote sensing [53–60]. The most important factors influencing greening appears to be an increase in carbon dioxide and nitrogen levels, land use, temperature, and weather [52,61–63]. Some limited studies have found an opposite trend, specifically that of forest browning [53,64]. This research matches more with the browning literature, with decreases in vegetation amount over time. Negative directional persistence in the forests of MG is consistent with literature on browning hotspots globally, especially in subequatorial Africa [10]. de Jong et al., (2013) stress the importance of the scale of the study, and, at a localized scale, direct human interventions, such as deforestation in this landscape, may have more effect than climate. Since rainforests are less climatically driven than savannas [9,54], the significant decline in NDVI on MG may reflect its recent incorporation into the national park and on-going struggles with deforestation for conversion to agriculture and timber harvesting [29]. Therefore, this study provides a local exception to the literature on global greening [53–60]. de Jong et al., (2013) report much less greening in more open savanna cover types, and Southworth et al., (2015) [20] note significant negative trends in dryland vegetation of Mozambique from 1982 to 2010.

These results also highlight the importance of understanding specific vegetation types that make up the heterogeneous savannas [44], whose impacts may otherwise have been masked by contrasting trends. The various rates of decline of NDVI across differing vegetation types underline the importance of understanding the vegetation structural groups within LG and MG [31,65–68]. While all vegetation types examined exhibit declining NDVI, particularly in the dry season, the decline is not equal. In LG, open grassland shows the greatest and most significant decline, whereas the declines of woody

vegetation cover types (mixed savanna and closed savanna) are not statistically significant. This finding further supports the proposal that open grasslands are the most vulnerable type of savanna because of their strong water-dependence and vulnerability to encroachment [5,8,10,20,41].

A principal purpose of protected areas is to conserve vegetation. Placement of the park on the landscape and proximity to humans can be as important in achieving this goal as biophysical variables, such as climate and soil [34,69–73]. The NDVI time series and directional persistence demonstrate that despite evidence of a strong link between precipitation and vegetation [11], protected areas can be effective in conserving vegetation. Vegetation inside protected areas tends to fare better than in surrounding areas [34]. Management practices of the LG protected area seem to be conserving vegetation better than in the buffer zones, as suggested by the smaller declines in NDVI compared to the surrounding buffer areas. MG is unique in the sense that management practices are not yet fully implemented and some deforestation persists, with concomitant detrimental impacts on NDVI as reported globally [20,24,25,74–84].

More locally, this study will help the park management team identify areas for conservation and those most vulnerable to vegetation decline like the grasslands and MG, both of which are key to tourism—the socioeconomic driver of this region. While precipitation data correspond very well with NDVI, a limitation of this study is that the data are resolved at a coarse spatial scale (0.05◦) relative to the study area [36], and that gridded datasets are dependent on the availability of local observations and satellite data, which may not be robust in this system. Another important note concerning DP values is that this is an NDVI comparison to the baseline in number of years above or below the baseline, but this measure does not account for magnitude of difference.

In the context of the broader landscape, this park seems to be maintaining vegetation health, and proving its effectiveness as a protected area in this respect. There were similar, but more negative trends detected in the buffer zone. Going forward, the park management team's work with community projects should also have a positive effect in the buffer zone. Though, socioeconomically speaking, this new model of conservation is already mutually beneficial for conservation [29]. One effort to combat deforestation on MG and stimulate the local economy is the production of shade-grown coffee [29]. This year, the shade-grown coffee being grown on MG is being sold globally. The park is currently expanding this coffee by 100 hectares per year, and by the year 2025, they will have developed 1000 hectares of coffee plantation. This translates to over 5000 hectares of restored forest cover on and around MG [29].

Significant mammal restocking, both herbivore (i.e., antelope) and carnivore (i.e., painted wolves), is also underway in the park, and more is planned in the future. Given the importance of herbivory as a driver in savanna landscapes and possible bush encroachment, the monitoring and measurement of vegetation cover change in this region is essential. The park can serve as an interesting comparison study to many other parks in southern Africa that are already suffering from bush encroachment and degradation due to overstocking [42–46]. The role of GNP, both as a crucial conservation area in this landscape, and as an important area for monitoring and research given its unique historical trajectory, cannot be overstated. Future monitoring and research will be of key importance to the management of African parks.

#### **5. Conclusions**

This study finds a declining trend in precipitation in GNP from 2000–2016, though the pre- (2000–2008) and post- (2009–2016) management precipitation regimes are not found to be statistically different. There is a corresponding decrease in NDVI over time across the entire study area. Different vegetation types display varying rates of decline in NDVI with grasslands and rainforest exhibiting the steepest (and significant) decline. Woody savanna vegetation experienced the least decline in NDVI over time.

The directional persistence metric evidences significant declines in vegetation biomass 2005–2016 compared to the 2000–2004 baseline. Negative persistence is strongest on MG, possibly reflecting continuing deforestation issues and the area's recent incorporation under protection. Open grassland savanna in LG also experienced significant negative persistence. LG experiences smaller declines in vegetation biomass and smaller values of directional persistence compared to the buffer zone, suggesting that the park is functioning effectively in conserving vegetation. Overall, there have been declines in vegetation in and around GNP during the 21st century. These declines may be prompted by a decline in precipitation, a shift in the seasonal timing of this precipitation, and direct human action. However, LG fares better than the surrounding area, further strengthening the justification for protected areas. The "laboratory" provided by this particular park and its unique history will be useful not only to scientists, but also to managers of protected areas, particularly savannas, as they can separate out the differing impacts of natural changes (such as precipitation which are outside manager control) and manmade changes (such as linked to deforestation, management rules etc.).

**Author Contributions:** Conceptualization, H.H., J.S., and B.C.; methodology, J.S. and P.W.; software, R.K.; formal analysis, H.H., R.K., P.W., and D.Y.; data curation, R.K., H.H., and D.Y.; writing—original draft preparation, H.H. and J.S.; writing—review and editing, J.S., P.W., R.K., H.H., D.Y. and B.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** These authors would like to thank the Gorongosa National Park Management/Science Team for their input on vegetation covers in the park and background information.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Functional Phenology of a Texas Post Oak Savanna from a CHRIS PROBA Time Series**

#### **Michael J. Hill 1,2,\*, Andrew Millington 2, Rebecca Lemons <sup>1</sup> and Cherie New <sup>1</sup>**


Received: 19 September 2019; Accepted: 13 October 2019; Published: 15 October 2019

**Abstract:** Remnant midwestern oak savannas in the USA have been altered by fire suppression and the encroachment of woody evergreen trees and shrubs. The Gus Engeling Wildlife Management Area (GEWMA) near Palestine, Texas represents a relatively intact southern example of thickening and evergreen encroachment in oak savannas. In this study, 18 images from the CHRIS/PROBA (Compact High-Resolution Imaging Spectrometer/Project for On-Board Autonomy) sensor were acquired between June 2009 and October 2010 and used to explore variation in canopy dynamics among deciduous and evergreen trees and shrubs, and savanna grassland in seasonal leaf-on and leaf-off conditions. Nadir CHRIS images from the 11 useable dates were processed to surface reflectance and a selection of vegetation indices (VIs) sensitive to pigments, photosynthetic efficiency, and canopy water content were calculated. An analysis of temporal VI phenology was undertaken using a fishnet polygon at 90 m resolution incorporating tree densities from a classified aerial photo and soil type polygons. The results showed that the major differences in spectral phenology were associated with deciduous tree density, the density of evergreen trees and shrubs—especially during deciduous leaf-off periods—broad vegetation types, and soil type interactions with elevation. The VIs were sensitive to high densities of evergreens during the leaf-off period and indicative of a photosynthetic advantage over deciduous trees. The largest differences in VI profiles were associated with high and low tree density, and soil types with the lowest and highest available soil water. The study showed how time series of hyperspectral data could be used to monitor the relative abundance and vigor of desirable and less desirable species in conservation lands.

**Keywords:** savanna; post oak; vegetation index; ecosystem function; phenology; encroachment; evergreen; deciduous

#### **1. Introduction**

Midwestern oak savannas originally occupied the transition zone between the true prairie grasslands and the eastern deciduous forests in North America. They stretched from eastern Texas in the south to northern Minnesota and occupied as much as 20 M ha in a continuous band [1]. Despite their association with rolling landforms on light sandy soils, most of these savannas were cleared for grazing and arable agriculture in the 19th century, leaving only about 12,000 ha of highly fragmented remnants [2]. Much of the remaining southern remnants are post oak woodlands and savanna grasslands in varying states of intactness, conservation status, or extirpation. Post oak savanna and savanna grasslands occupy approximately 2.7 M ha of eastern Texas [3]. Other studies place the total area in Texas at around 4 M ha [4,5]. In eastern Texas, Post Oak Motte and Woodland occupies about 38% and Savanna Grassland occupies about 58% of this vegetation type, with other minor types making up the remainder [3].

The Gus Engeling Wildlife Management Area (GEWMA) represents a relatively intact example of post oak savanna and savanna grassland and includes hardwood bottomland forest and riparian areas [5]. As a result of fire suppression, the upland oak savanna dominated by post oak (*Quercus stellata*), red oak (*Q. falcata*), blackjack oak (*Q. marilandica*), and bluejack oak (*Q. incana*) has thickened to exclude understory grasses such as little bluestem (*Shizachyrium scoparium*) in many areas, and open areas have suffered from woody encroachment by native evergreen species such as eastern red cedar (*Juniperus virginiana*; [6]) and yaupon (*Ilex vomitoria*). The GEWMA is under an active program of land cover management aimed at the regeneration of open savannas by clearing sections of oak and replanting little bluestem grassland.

For the past 15 years, the spaceborne hyperspectral remote sensing of vegetation has depended upon two experimental sensors, EO-1 Hyperion (Earth Observer 1 Hyperion) and CHRIS/PROBA (Compact High-Resolution Imaging Spectrometer/Project for On-board Autonomy) for data at a spatial resolution (30–34 m pixels) for field and landscape-scale studies [7–11]. Although the CHRIS sensor only covers the 400–1000 nm wavelength range, it also provides nadir images and two forward and backward-looking off-nadir angles. It has been applied to a wide range of studies on classification [12], vegetation properties [13,14], vegetation mapping [15,16], water [17–19], canopy function [20], canopy radiative transfer [21,22], and vegetation indices [23].

Time series of hyperspectral data have proven valuable in vegetation classification and the detection of changes in functional processes. However, dense time series over a single site are relatively rare. The GEWMA site, with its thickening of post oak woodland and woody encroachment by evergreens, provides an opportunity to use CHRIS data to explore the extent to which suites of functionally based vegetation indices (VI) can distinguish differences in vegetation canopy properties and behaviors over a full seasonal growing cycle. A very large number of VIs have been constructed to target various aspects of photosynthetic function, canopy structure, fractional cover of photosynthetic and non-photosynthetic components, and canopy water content [24]. In general, narrow-band VIs acquired from visible and near infrared wavelengths can provide quantitative information about photosynthetic pigments such as chlorophyll, carotenoids and anthocyanins ([25–27], photosynthetic radiation use efficiency (RUE; [28,29]), canopy senescence using the chlorophyll to carotenoid ratio [30], and subtle differences in photosynthetic capacity based on many variants and more complex surrogates of the Normalized Difference Vegetation Index (NDVI; see [24]). Imagery from sensors such as Hyperion with short-wave infrared bands can be used to derive VIs with sensitivity to canopy water content and dynamics [31] and the fractional cover of photosynthetic and non-photosynthetic canopies and bare soil [32].

Previous work has explored the application of selected VIs derived from EO-1 Hyperion imagery to the identification of vegetation states in tropical savannas [33] and North American savannas and grasslands [34]. However, the image acquisition frequency in these studies was insufficient for the characterization of phenology. In this study, a sequence of 17 CHRIS images was acquired between June 2009 and October 2010, providing an opportunity to explore variation in canopy dynamics among deciduous and evergreen trees and shrubs and savanna grassland in seasonal leaf-on and leaf-off conditions for a valuable remnant of the threatened midwestern oak savannas. The objectives of the study were as follows:


#### **2. Materials and Methods**

#### *2.1. Study Area*

The Gus Engeling Wildlife Management Area (GEWMA) is located about 29 km northwest of Palestine, near Bethel in Anderson County, Texas (Figure 1a). The GEWMA was acquired by the State of Texas between 1950 and 1960, and comprises 4434 ha of prairie, post oak savanna, hardwood bottomlands, and riparian areas including pitcher plant bogs ([3,5,35]. The area is regarded as being less disturbed than much of Texas, since soils are too poor for arable agriculture. However, it was grazed as open range by cattle and hogs from the mid-1800s, becoming severely overgrazed prior to acquisition by the State [5]. The landscape has a rolling terrain with elevation ranging from 75 to 150 m and vegetation cover dominated by a mix of deciduous and evergreen trees and shrubs (Figure 1a,b). The climate is characterized by hot summers and mild winters, with average monthly maximum temperatures approaching 35 ◦C in August and average monthly minima approaching 0 ◦C in January (Figure 2). Average monthly rainfall reaches a peak in May and October; however, rainfall can be highly variable. Very heavy rainfall in May 2010 resulted in the flooding of the lower parts of GEWMA and the surrounding villages and towns (Figure 2).

The soils range from loamy sands and sandy loams on the terraces and lowlands to fine sands on the upland post, black and bluejack oak woodlands (Table 1). The Ecological Site Description (ESD) provided with soil maps by the Natural Resources Conservation Service (NRCS) of the United States Department of Agriculture (USDA) provides descriptions of original reference vegetation and state and transition frameworks for vegetation succession in response to disturbance, fire exclusion and management of herbivory (Table 1; Soil Survey Staff, NRCS, 3 October 2017). However, the detailed vegetation mapping in [5] and in Phase 2 of the Texas Vegetation Classification Project [3] (Figure 3b; Table 2) provides a better indication of the vegetation on GEWMA at the time of this study.

**Figure 1.** (**a**) Location of Gus Engeling Wildlife Management Area (GEWMA) in Anderson County, near Bethel, Texas showing the elevation of the county landscape. (https://tpwd.texas.gov/huntwild/ hunt/wma). (**b**) National Agricultural Imagery Program (NAIP) images from 2009 and 2012 show the canopy cover in leaf-off and leaf-on states, the classification window (red box) used later in the study, and the structured clearing of thickened woodland to re-establish savanna in the north-west of the area (2012 image).

**Figure 2.** Average monthly and actual daily maximum and minimum temperatures and precipitation for Palestine, Texas for 2009 and 2019. Note the flooding rains in late May 2010. Letters indicate the timing of Compact High-Resolution Imaging Spectrometer (CHRIS) acquisitions with Y/N indicating whether the image was used in the analysis or not.

The vegetation canopy structure on GEWMA is shown in Figure 1b contrasting leaf-off and leaf-on images from 2009 and 2012. Note that site managers started to clear part of the thickened woodland in 2010 in order to restore some of the original oak savanna parkland environment. The result of this is visible in the 2012 image. A simplified version of the Phase 2 land cover map identifies four major and three minor vegetation associations on GEWMA: savanna grassland, post oak motte and woodland, seasonally flooded hardwood forest, floodplain hardwood evergreen forest; and juniper woodland, mesquite woodland, and post oak/red cedar motte and woodland (Figure 3b). The species composition of these main types is shown in Table 2, indicating the diversity of the genera and species involved. The proportions of these vegetation sites associated with each major ESD and NRCS soil types are shown in Table 1 (the minor land cover types with evergreen shrubs are all aggregated into one for this purpose). The savanna grassland, post oak motte and woodland are strongly concentrated on Darco, Kirvin, Teneha and Tonkawa fine sands, whilst floodplain hardwoods and seasonally flooded hardwoods are primarily concentrated on Thenas fine sandy loam, Naconiche loamy fine sand, and Nahatche and Pluck soils (Table 1).

Observations from field work carried out in November 2009 indicated that the upland oak savanna was thickened and dense with little or no understory (Figure 4a), the savanna grassland (Figure 4b,c) had areas of bare soil, thin areas of little bluestem and scattered evergreen shrubs, and the lowlands and terraces contained dense stands of hardwoods with evergreen shrub understories (Figure 4c background, Figure 4d roadside). The map developed in [5] identifies natural vegetation types with more spatially explicit detail but does not identify invasive evergreens aside from long-established eastern red cedar and mesquite stands.


**Table 1.** Dominant ecological site descriptions (provisional status), soil types and major Texas Land Cover Project classes for GEWMA.

<sup>1</sup> certain soil polygons do not have an Ecological Site Description (ESD) association: LuA (Lufkin fine sandy loam), W (Water); <sup>2</sup> POW – post oak woodland; SG – savanna grassland; HF – hardwood forest; SFH – seasonally flooded hardwood forest; JM – juniper or mesquite woodland.

**Figure 3.** (**a**) Major soil types by soil polygon from the United States Department of Agriculture (USDA) Natural Resources Conservation Service (NRCS) for GEWMA; see Table 1 for descriptions (minor soil types in the southern area are shown as blank for simplicity). (**b**) Major land cover classes from the Texas Vegetation Classification Project ([2] for GEWMA (minor classes covering <200 pixels in the dissected southern area have been excluded for simplicity).

**Table 2.** Major species associated with the most important land cover types from the TNRIS classification and major encroaching evergreen species mapped in this study [3,5]. Species are listed from major to minor; however, woodlands and forests are very heterogeneous.


**Figure 4.** (**a**) Thickened post oak woodland on the sandy upland ArD or DaD soil type with virtually no understory. (**b**) Savanna grassland with sparse grasses, yaupon and eastern red cedar. (**c**) Degraded savanna grassland with bottomland hardwood forest in the background. (**d**) Dense mixed oak woodland with shrubby understory by the road. (Photos: M.J.H and A.M.).

#### *2.2. Remote Sensing Data and Processing*

#### 2.2.1. Aerial Photography

To aid with the geometric correction of CHRIS data, characterization of the vegetation cover, and identification of evergreen vegetation, a National Agricultural Imagery Program (NAIP) 1 m resolution airborne image with red, green, blue and near infrared bands was acquired for January 2009 coinciding with the winter leaf-off period for deciduous trees. The four-band image was subjected to an unsupervised isodata classification to create 25 classes. These classes were then aggregated to five vegetation classes (water or shadow, deciduous trees, leaf/grass understory; bare or sparse herbaceous vegetation, and evergreen trees and shrubs) based on a visual interpretation of the NAIP image (Figure 5a). Potential for error was greatest in distinguishing deciduous tree cover from understory since the separation of trees from understory under leaf-off conditions was dependent upon trunk/shadow mixtures. Varying the aggregated class allocation for the most marginal class resulted in unacceptably high "omission" differences of -86% for deciduous tree cover and a "commission" difference of 140% for understory, which did not match our knowledge of the site based on measurements of canopy trees along two transects in central GEWMA, as well as qualitative surveys throughout the area. Using the most conservative aggregation for evergreens changed the evergreen percentage from around 20% to around 10%. However, this excluded classes with clear red edge signatures indicative of photosynthetic vegetation. Since this classification was undertaken at the NAIP pixel resolution of 1 m, but subsequent analysis was undertaken at a resolution of multiple CHRIS pixels, minor accuracy errors in this classification had very little impact on the results of the study.

**Figure 5.** (**a**) Classification of the winter 2009 NAIP image into five classes specifically identifying evergreen trees; and (**b**) Fishnet polygons corresponding to a 5 × 5 grid of CHRIS Project for On-Board Autonomy (PROBA) pixels showing the proportion of evergreen trees.

#### 2.2.2. CHRIS PROBA Data

The ESA sensor CHRIS/PROBA was launched on 22 October 2001 [8]. The CHRIS sensor on the PROBA-1 platform was designed to capture high-resolution images at multiple view angles in 62 narrow reflectance channels for wavelengths between 400 and 1000 nm [8] (Table 3). The sensor could be tasked in several modes which varied the swath width, pixel resolution and number of spectral channels. For this study, CHRIS data were acquired in Mode 1 which provided 62 spectral channels and five view angles at a nominal pixel resolution of 34 m (Table 3). The CHRIS data were acquired on 17 dates between June 2009 and October 2010 over the GEWMA study area (Table 4).

The CHRIS images were processed using the CHRIS Toolbox [36] in version 4.7 of the BEAM software package to process data from a variety of European Space Agency sources [37] and in ENVI®(ENvironment for Visualizing Images). Following the procedure in [36], the standard noise reduction, cloud identification and masking, and atmospheric correction with spectral calibration, aerosol optical thickness retrieval and column water vapor retrieval were carried out to obtain surface reflectance images. The atmospheric correction module of the CHRIS Toolbox [38] is based on MODTRAN4 (MODerate resolution atmospheric TRANsmission) [39]. Although the full set of multi-angular images was acquired and processed for almost all dates, only nadir images were used in this study, as the focus was on the temporal domain rather than the angular domain.

Nadir surface reflectance images were converted to tiff format, exported and imported into ENVI®. Only those images without cloud over the GEWMA site were subjected to further processing, and the 20 December image was excluded due to a targeting error. This left 11 useable images for further processing (Table 4). Images for each date were individually geometrically corrected using ground control points from the NAIP image. Due to the nature of the GEWMA environment with limited manmade features, CHRIS images were warped to an 18 m pixel resolution with a registration accuracy of approximately 1 pixel. Although this leads to the oversampling and smearing of the detected signal, this was not important in the context of this study, in which patch-based temporal profiles were the major target. Finally, the individually registered images were co-registered using Texas Road centerlines with an accuracy of approximately 2 pixels (36 m).





<sup>1</sup> Sen - senescing.

#### 2.2.3. Vegetation Indices

A total of 51 narrow band VIs was calculated from the 62 channel surface reflectance images based on the formulae documented in [24] with wavelengths adjusted for CHRIS band centers. Given the nature of the study site, with very dense tree cover except for certain soils where savanna grassland occurred, the identification of the key VIs focused first on known pigment-sensitive wavelength combinations found to be effective in previous studies with Hyperion imagery [33,34], and then on whether the many different narrow-band indices creating variants on the NDVI [38] and targeting the red edge provided any additional sensitivity to augment NDVI. The VIs that utilized blue wavelengths below 480 nm often failed due to low-quality blue bands. Based on [33,34], visual assessment of images and pixel histograms, eight VIs were selected as indicators of vegetation canopy properties (Table 5). The NDVI was ultimately selected as the indicator of greenness/photosynthetic capacity since no advantage was gained from the other greenness VI variants. Although band centers were well aligned with narrow band reference values, CHRIS band widths at 6–12 nm can only broadly approximate the narrow bands used in the original derivation of some of these indices. For example, [28] used an average band width of 2.8 nm in deriving the PRI2. Nevertheless, the image value ranges for highly precise VIs such as PRI and PRI2 matched the ranges observed in the original studies [28,29].


**Table 5.** Vegetation indices selected as indicators. Formulae depict band reflectance (R) wavelengths in nanometers.

<sup>1</sup> From [24]; <sup>2</sup> band centers based on [41].

#### *2.3. Analysis*

#### 2.3.1. Fishnet Polygon Extraction of Tree Cover Density and VI Profiles

In order to capture the broad patterns of vegetation canopy dynamics with such heterogeneous species mixtures (Table 2), the analysis was based on averaging across pixels at a patch scale. A fishnet polygon layer was constructed with a polygon cell size of 5 × 5 CHRIS pixels making a 90 × 90 m cell. Each fishnet polygon was attributed a unique number. Pixel counts of tree and other cover classes were then extracted from the classified NAIP image (with each fishnet polygon cell containing 8100 NAIP pixels), converted to cover class proportions and assigned to each fishnet polygon cell (Figure 5b; percentage of evergreen class shown). The fishnet was then merged with the soil polygon layer for GEWMA to create a combined layer with vegetation cover fractions and soil type for each fishnet polygon. The fishnet was then used to extract VI values for each fishnet polygon from each VI time

series creating a database with cover fraction, soil type, and VI value for each 90 × 90 m polygon across GEWMA. To minimize the effects of uncertainty in the attribution of aggregated classes in the NAIP image classification, the time series profiles of each VI were then examined for average deciduous and evergreen tree cover at 10% intervals. Associations with the Texas Land Cover classes were also examined.

#### 2.3.2. Combined Indicators of Vegetation Function

Based on image coverage and quality, two dates (September 2009—leaf-on—and 25 January 2010—leaf-off) were selected to explore information provided by combining VIs. The data were windowed to focus on the northern half of GEWMA, where a consistent gradient from the upland sandy savannas to the lowland hardwood floodplain and riparian areas could be sampled. The six best VIs from the phenological analysis were stacked (Anthocyanin Reflectance Index 2 (ARI2), Carotenoid Reflectance Index (CRI), NDVI, Photochemical Reflectance Index 2 (PRI2), Plant Senescence Reflectance Index (PSRI) and Water Index (WI)) for each date. An isodata unsupervised classification was carried out on the VI stacks for each date to create a 20-class image. The classifications were then aggregated to eight classes based on the dendrograms and class distances. The mean and standard deviation of VI index values, NAIP classification cover percentages and dominant soil type percentage were extracted for each classified image.

#### **3. Results**

The extracted time series profiles showed major differences in relation to the major land cover types and deciduous tree density levels and in relation to the proportion of evergreen trees and shrubs present.

#### *3.1. Overall VI Profiles*

The comparison of profiles for the selected VIs showed a clear separation between areas with trees and those without during leaf-on. For most leaf-on dates, savanna grassland VI profiles were separated from the 40% and 80% tree cover profiles by more than one standard deviation (Figure 6). The differences among the deciduous tree canopies largely disappeared during leaf-off.

Although increments in VI response to tree density were evident, there was a substantial overlap in standard deviations, with CRI and ARI1 and ARI2 exhibiting a higher variation in values within a tree density level and between dates. In general, magnitudes of pigment-sensitive VIs during leaf on could be ranked 80% > 40% hardwood > 40% post oak >> savanna/grassland (for PRI2, low values represent a high response). The largest range of response was exhibited by CRI during leaf-on, but the rankings among tree densities were the same. The PSRI exhibited little sensitivity to hardwood tree density, but a consistently higher response for post oak at 40% tree cover and savanna/grassland (Figure 6g). The compressed and low values for CRI, ARI and ARI2 for day 576 (July 30) are likely anomalous as they are not supported by parallel behaviors in NDVI, WI, PRI or PRI2. The WI showed large seasonal variation in canopy water content, but the main difference in profiles was between the savanna grassland and the areas with trees (Figure 6h).

#### *3.2. Response of VIs to Evergreen Density*

The effect of evergreens on the VI profiles was examined for the DaD soil type dominated by post oak woodland with total tree cover of >80% (see Table 1) and evergreen tree cover in increments of 10% up to the maximum observed of 40% (Figure 7). The pigment-sensitive VIs exhibited small increases in values in the leaf-off period between 350 and 430 days (mid-December 2009–mid-March) in response to increasing evergreen cover, except for the ARI (Figure 7). Although the increments due to evergreen density were evident, standard deviations within date and density levels were relatively large and exhibited a significant overlap. The PRI2 and NDVI exhibited lower values at 40% tree cover than at 80% tree cover in the leaf-on period of 2010 (Figure 7a,b). The PSRI and WI provided

the clearest response to increasing evergreen percentage, with PSRI declining and WI increasing incrementally with each 10% increase in evergreens (Figure 7c,d). The responses to increasing the evergreen percentage were similar for the other soil types, but the magnitude of the difference between zero and 40% evergreens was generally smaller. Evergreen percentages were very low on the ArD soil type except when adjacent to bottomlands.

**Figure 6.** Overview of seasonal phenology of vegetation indices that are indicators of photosynthetic activity and pigments averaged across soil polygons with predominant vegetation types of post oak, bottomland hardwoods and grassland at tree cover percentages of 80, 40 and <10 and no evergreens. The bars denote standard deviations. The gap in the ARI2 profile is due to an unrecoverable error in the calculation for 19 October (see Table 1 for land cover association with soil types). (**a**) PRI; (**b**) PRI2; (**c**) NDVI; (**d**) CRI; (**e**) ARI; (**f**) ARI2; (**g**) PSRI; and (**h**) WI.

#### *3.3. Responses within Tree Density Levels among Soil Types*

Available soil water (ASW) is a major variable among dominant soil types on GEWMA, with upland sandy soils having 6–10% ASW and loamy bottomland soils having approximately 14–16% ASW (Table 1). In contrast with the other profile analyses, differences in PRI2 and WI among the different combinations of soil type and tree density mostly exceeded one standard deviation indicating the importance of ASW in canopy processes. The comparisons between selected soils at the extremes of ASW content showed major differences during the leaf-on period among the three soil types predominantly associated with grassland savanna or woodland at tree densities of <10% for PRI2, and smaller differences during the leaf-off period (Figure 8a). There were also differences between the ArD and Na at 80% tree cover in the leaf-on period, but these disappeared in the leaf-off period. There were small differences in WI between ArD and both DkF and DaD soils in the leaf-on and leaf-off periods, and between ArD and Na in the leaf-on period (Figure 8b), but the relative differences were smaller than for PRI2.

**Figure 7.** Variation in responses of the Vegetation Indices (VIs) most sensitive to increasing levels of evergreen tree cover within the DaD soil type dominated by post oak woodland. The bars denote standard deviations. (**a**) PRI2; (**b**) NDVI; (**c**) PSRI; and (**d**) WI.

**Figure 8.** Variation in VI indicators across extremes of deciduous tree cover (<10% and >80%) without evergreens, and across soil types in terms of available soil water (see Table 1). The bars denote standard deviations. (**a**) Photosynthetic activity, PRI2; and (**b**) vegetation water content, WI.

#### *3.4. Sensitivity of VIs to Spatial Variability in Leaf-on and Leaf-o*ff *Periods*

Images from 6 September 2009 and 25 January 2010—peak leaf-on and leaf-off conditions, respectively—were chosen to examine the sensitivity to vegetation and canopy variation and underlying soil differences that could be obtained by combining the VIs by classification. The 2009 date was chosen because site managers commenced clearing parts of the upland woodland in mid-2010 (see Figure 1) and consistent vegetation cover between leaf-on and leaf-off was required. The images from these dates clearly reflect the differences in leaf-on and leaf-off canopies (Figure 9) and highlight the presence of evergreen vegetation in January (Figure 9b,d,f,h). Whilst the NDVI shows limited within-canopy sensitivity, PRI2 shows a marked difference between evergreen-dominated and deciduous oak-dominated areas in September (Figure 9a,c). The upland post oak woodland and savanna environment associated with drier soil types (Table 1) is sharply delineated in January NDVI and ARI2 (Figure 9b,f). The PSRI and WI show the same patterns as the pigment-specific indices with a similar sharp delineation of the upland post oak woodland and savanna in January in both VIs (Figure 10) and the presence of evergreen vegetation in January (Figure 10b,d).

The image classification results for September and January produced distinctly different class patterns, with the leaf-on September image highlighting the grassland/savanna areas, variation in the post oak woodland canopy, and no clear distinction of hardwoods on bottomland (Figure 11a). By contrast, the results for January clearly distinguished the post oak woodland and savanna grassland on ArD and DaD soil types, areas of evergreen vegetation and the flooded stream bottomland associated with the Na soil type (Figure 11b).

The extracted profiles for each classification illustrate the difference in VI sensitivity, class association with tree and evergreen density, and class association with soil type (Figure 12). Available soil water is a major variable among dominant soil types on GEWMA, with upland sandy soils having 6–10% ASW and loamy bottomland soils having approximately 14–16% ASW (Table 1). Within the image, the dominant soil types were DaD > ArD > DkF > Na, thus representing a contrast between drier soils supporting post oak woodland and savanna grassland, and bottomland loams supporting hardwood forest (Figure 12). In September, most of the area was occupied by classes 1, 2, 6, 7 and 8. Classes 1 (78%) and 2 (61%) were strongly associated with elevations between 100 and 129 m, while classes 7 (55%) and 8 (76%) were strongly associated with elevations between 78 and 100 m. Neither the evergreen vegetation nor the understory fraction varied much between classes, but bare soil was most associated with class 1 and least associated with class 8. In all cases, VI class values showed either a regular positive progression from class 1 to class 8 for ARI2, CRI, NDVI, PRI2 (negative PRI2 represents a higher value) and WI, and a negative progression for PSRI (Figure 12). These progressions represented a regular increase in pigments, photosynthetic activity and canopy water content, and a regular decrease in brownness/yellowness with class number. Class 1 occupied 40% of soil type ArD, while class 2 occupied 20% of each of ArD, DkF and FuB.

In January, most of the area was occupied by classes 1, 3, 4, 5, 6, and 7 (Figure 12b). Bare areas were associated with classes 1 and 3, and evergreen vegetation was associated with classes 7 and 8. The class association with elevation was more marked than in September, with between 88–100% of classes 5–8 being associated with elevations between 78–115 m, and 63% of class 2 being associated with elevations between 78 and 85 m. Classes 7 and 8 had the highest NDVI, CRI and ARI2 and the lowest PRI2 values as well as the highest WI and lowest PSRI. Among classes 1, 3 and 4, associated with the deciduous post oak woodland and savanna grassland, there were lower values for WI for classes 3 and 4 with greater understory cover than class 1 with more bare cover. Class 1 occupied >40% of the Na soil type, class 3 occupied between 30 and 40% of soil types ArD, DaD and FuB, and class 6 occupied around 30% of the minor soil types KnE, DkF and PeC (Figure 12b).

**Figure 9.** Photosynthesis and pigment-sensitive VIs at leaf-on and leaf-off for deciduous trees. (**a**) NDVI 6 September 2009; (**b**) NDVI 25 January 2010; (**c**) PRI2 6 September 2009; (**d**) PRI1 25 January 2010; (**e**) ARI2 6 September 2009; (**f**) ARI2 25 January 2010; (**g**) CRI 6 September 2009; (**h**) CRI 25 January 2010.

**Figure 10.** Dryness and moisture-sensitive VIs at leaf-on and leaf-off for deciduous trees. (**a**) PSRI 6 September 2009; (**b**) PSRI 25 January 2010; (**c**) WI 6 September 2009; (**d**) WI 25 January 2010.

**Figure 11.** Classification of the northern section of GEWMA based on six VIs aggregated to eight classes. (**a**) September 6, 2009: deciduous leaf-on; (**b**) January 25, 2010: deciduous leaf-off.

**Figure 12.** Vegetation Index profiles of cover types from classified NAIP imagery and proportion of NRCS Web Soil Survey soil types associated with (**a**) September 2009 and (**b**) January 2010 classifications shown in Figure 11. Numbers at the top indicate the percentage of the area occupied by each class. Numbers on the soil type legend indicate the percentage of the area occupied by each major soil type (minor soil types were excluded).

#### **4. Discussion**

This study on a remnant post oak woodland area has shown that the major differences in spectral phenology as defined by VIs are associated with deciduous tree density, the density of evergreen trees and shrubs—especially during deciduous leaf-off periods—broad vegetation types, and soil types. Although major differences in VIs at peak leaf-on and leaf-off were clearly associated with high and low evergreen density and extremes of ASW associated with soils, combining VIs in a simple classification also revealed variation in canopy water content and photosynthetic capacity within ostensibly uniform post oak savanna on the sandy soils with low ASW. The photosynthetic capacity of evergreens during deciduous leaf-off was clearly indicated by NDVI, PRI1 and PRI2 profiles with obvious growth advantages for further encroachment. Although tree species vary between post oak savanna and bottomland hardwoods, VI phenology was remarkably similar at the same tree densities. Heterogeneous and shared species composition post oak woodland and hardwood forests limited patch-scale discrimination with VI profiles.

The temporal pattern of spectral VIs was dominated by the deciduous tree leaf-on/leaf-off cycle of post oak woodland and hardwoods, and spatial variation in amplitude was primarily determined by tree density from >80% cover down to <10%. There was a greater range in the amplitude of NDVI and CRI than PSRI—an indicator of the chlorophyll/carotenoid ratio [30]—or ARI and ARI2 during leaf-on in response to the tree densities and major vegetation types. The precipitation records from September 2009 to October 2010 indicate that the site did not undergo any unusual moisture stress during the study. In addition, due to a flooding event in May of 2010 with substantial cloud in the preceding couple of months, 2010 spring green-up was not captured. There are significant differences in pigment magnitudes [42], anthocyanin synthesis [43] and rates of decay during senescence associated with stress or different mixes of deciduous species [44,45]. Such differences were either minimal or undetectable in this study. Since deciduous species are highly intermingled in both post oak savanna and bottomland hardwood forest (see Table 2 [3,5]; the pixel resolution of CHRIS and the 90 m resolution of the fishnet used in the analysis meant that patch-scale species differences were likely to be distinguishable).

At the landscape scale, the upland post oak savanna was clearly delineated in the leaf-off VI images, and there was a strong association with soil type and elevation. However, in leaf-on images, PRI2 and WI showed most sensitivity to the differences between upland and bottomland in soil properties, elevation and vegetation type, and to variation within the upland post oak savanna woodland and savanna grassland. The PRI and PRI2 are sensitive to photosynthetic RUE [28,29] and evergreens exhibit reduced midday PRI levels compared to deciduous species [29]. In this study, lower values for PRI and higher values for PRI2 were observed for 40% evergreen tree cover compared to 0–30% evergreen cover at 80% total tree cover in the peak leaf-on period of 2010. This suggests that evergreens suffer some disadvantage in photosynthesis when competing with deciduous trees during the leaf-on period. However, this is a minor factor when compared with the opportunity for photosynthesis and growth in the 3–4 month deciduous leaf-off period. An examination of the weather during the study period (Figure 2) and the long-term climate indicates that the temperature and precipitation are favorable for the growth of evergreens, especially coniferous *Juniperus* spp., while deciduous trees are in a leaf-off state.

The largest differences in pigment/senescence and photosynthetic-sensitive VIs occurred when comparing tree canopies among the sandy soil types with the lowest ASW at tree densities of <10%, and between the soil types with the lowest and highest ASW for tree densities >80%. This suggests that water relations are not only the most important driver of the vegetation types between upland woodland and bottomland hardwood forest at the landscape scale, but also potentially influence the vigor and location of savanna tree species at the patch scale. While differences in leaf-off classes derived from combined VIs were obviously associated with the presence or absence of an evergreen canopy, leaf-on classes associated with deciduous cover of >50% and ≤20% bare (classes 5–8) showed a significant variation in WI that is indicative of different canopy water contents. Studies have shown significant differences in water use attributes such as access to soil water, stomatal conductance, water use efficiency, and leaf-specific hydraulic conductance among *Quercus* spp. [46,47]. The post oak woodland contains a mix of *Quercus* spp., but tree-by-tree identification for any local site was not available or collected for this study. However, spatial variation in class distribution in the water-limited upland savanna environment in September 2010 may be indicative of both local variation in *Quercus* species composition and the vigor of trees associated with differences in water relations attributes.

The observation and analysis of woody encroachment across North America suggest that the potential for "release" of woody encroachers favored by fire suppression is greater at sites with higher precipitation [48], and by inference might be better on bottomland soils with higher ASW rather than on upland soils with lower ASW. Certain *Quercus* species have exhibited higher pre-dawn and midday water potentials than eastern red cedar [49,50] which is known to use water year-round, and almost 100% of precipitation reaches the soil surface [51]. In this study, higher densities of evergreens were clearly negatively associated with the soils with the lowest ASW, higher elevation (compare Figures 3 and 5) and those most dominated by thickened post oak woodland. Evergreens were mainly distributed among a range of soils in the bottomlands as an understory to hardwood canopies and along small tributaries of the main stream, Catfish Creek, at relatively lower elevations than the adjacent terrain. Among the native evergreens, eastern red cedar may have the most potential to encroach on low ASW soils since it has a record of highly successful encroachment in grasslands where it has an advantage in terms of rooting depth [52] and in upland oak forests in Oklahoma [53]. At the time of this study, the upland post oak savanna area was densely thickened with no grasses and a dense thatch of leaf litter (Figure 4a). The combination of high oak density, low ASW, rapid soil drainage at higher elevations, and the advantage in water relations attributes over eastern red cedar may explain its almost complete absence from the upland post oak woodlands.

Close to the end of the study period, a section of the northern upland post oak savanna was partially cleared to enable the planting of grasses and regeneration of a savanna landscape as part of an on-going management plan by the Texas Parks and Wildlife Department (Figure 1b). Small areas of GEWMA had been replanted with *Shizachyrium scoparium* prior to the commencement of this study. Airborne and satellite remote sensing could play a key role in the monitoring and management of GEWMA by tracking the development of replanted areas, identifying evergreen woody encroachment in newly planted areas, selecting areas for spraying, burning and/or mechanical control of eastern red cedar, and monitoring the long-term health of the important remnants of original savanna grassland, post oak woodland and hardwood forest (e.g., [54]). Such monitoring could be carried out with relatively frequent imagery at a 10–20 m resolution from Sentinel 2a and 2b. A prior study simulated Sentinel 2 with Hyperion imagery across a range of sites in North America and showed that suites of VIs could provide clear insights into differences among vegetation states and processes [34]. Various baseline conditions and approaches for the restoration of oak savannas [55,56] could be informed with suitable passive and active remote sensing data (e.g., [34,57]). This study used a now almost 20-year-old sensor with relatively low signal to noise and 400–1000 nm spectral coverage to detect a variety of differences seasonal phenology due to plant functional type, vegetation type, soils, and tree density. Advances in the hyperspectral detection of plant species diversity and function [58–62] provide exciting potential for monitoring important repositories of biodiversity with heterogeneous species composition from space or airborne platforms.

#### **5. Conclusions**

This study showed that functional VIs can identify differences in vegetation types that are strongly associated with soil type and landform and indicative of differences in physiological function and plant functional type. The study applied a relatively rare time series of CHRIS imagery at an acquisition frequency that captured variation in canopy dynamics both between and within major vegetation types and tree densities. Although the discrimination of known species-level differences in pigment dynamics and soil water access could not be detected with precision using the CHRIS data due to heterogeneous species mixtures and the patch-scale analysis applied, the potential for this with airborne

sensors or new spaceborne hyperspectral sensors was clearly indicated. Although no explicit difference in photosynthetic function among deciduous canopies was detected in this study, nevertheless, the VIs sensitive to carotenoids, anthocyanins and photosynthetic efficiency can discriminate among species and landscape components when retrieved with better-quality sensors. In the meantime, Sentinel 2 with its additional bands in the pigment-sensitive wavelengths at a 20 m pixel resolution, may be used effectively for current monitoring. Although evergreen encroachment on rare oak savanna conservation lands may be easily monitored with leaf-off aerial photography, hyperspectral imagery can provide additional insights into the phenological behaviors that relate to physiological function and competitive advantage. This information could be helpful in the management and control of evergreen shrubs on bottomland areas [63] that present a substantial challenge due to their competitive advantages over deciduous oaks and hardwoods in the leaf-off period.

**Author Contributions:** Conceptualization, M.J.H. and A.M.; methodology, M.J.H.; field observations and photos A.M. and M.J.H.; formal analysis, M.J.H., R.L. and C.N.; data curation, M.J.H.; writing—original draft preparation, M.J.H.; writing—review and editing, A.M.; project administration, A.M. and M.J.H.

**Funding:** The involvement of the senior author was partially supported by grants from NASA (Contract #NNX09AQ81G and #NNX10AH20G).

**Acknowledgments:** We thank the managers of GEWMA for providing access to the site for field observations. Tony Fillipi assisted with field work. Field work costs were covered by funds from Texas Agriculture Experiment Station, Texas A&M University to A.M.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
