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

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.


Introduction
Global distributions of vegetation structure, species composition, and functional composition are governed primarily by climate and phylogenetic processes [1][2][3][4][5].In the seasonal tropics, precipitation patterns and ecohydrological interactions are the principal drivers of large scale vegetation patterns [6][7][8], with increasing precipitation generally corresponding to increasing tree density, canopy cover, and tree height [8][9][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][12][13][14][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][17][18][19][20][21][22][23].
Numerous studies have examined satellite observations in search of empirical support for the prevalence of alternative vegetation states [18,19,[24][25][26][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][25][26][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][32][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.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.

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][36][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).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][36][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).

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 250m 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

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.

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).

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][42][43][44][45][46][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).
tends to reduce spatial heterogeneity and thus suppress identification of alternative vegetation modes.

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 nonvegetative 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).

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][42][43][44][45][46][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).

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][32][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.

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.

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][25][26][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.

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 km 2 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.

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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 14 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.

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 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 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.

.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 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.
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.

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 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 tropical forests 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.

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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 10 of 14 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.

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

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][25][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][32][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.

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.

Figure 3 .
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).

Figure 3 .
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).

Figure 4 .
Figure 4. Histograms showing data distributions for (A) tree cover, (B) NIR albedo, and (C) LST.A total of 21,081,051 1km 2 locations with MAP >= 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 4 .
Figure 4. Histograms showing data distributions for (A) tree cover, (B) NIR albedo, and (C) LST.A total of 21,081,051 1 km 2 locations with MAP ≥ 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 .
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

Figure 5 .
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

Figure 7 .
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 >= 100mm 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.

Figure 7 .
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.