Introduction

The idea that the distribution of Earth’s biomes is governed by ‘climate envelopes’ has a long history (Von Humboldt and Bonpland 1807; Schimper 1903). According to this view, climate has a deterministic effect on vegetation types, implying that the locations of stable biome boundaries are predictable. Correspondingly, the distributions of the biomes on Earth, such as tropical forests and savannas, are expected to shift gradually with gradual climate change. In recent years, however, evidence has been growing that the effect of climate on tropical forests and savannas is non-linear (Reyer and others 2015). Notably, studies based on satellite data of tree cover showed that under a range of mean annual precipitation (MAP) high frequencies of both savanna (~20 %) and forest tree-cover levels (~80 %) can be found, whereas intermediate values are rare (Hirota and others 2011; Staver and others 2011b). This bimodality of tree cover in both satellite data and field data under given environmental conditions indicates that forest and savanna may represent separate basins of attraction, which is simply referred to as alternative stable states (Scheffer and Carpenter 2003; Hirota and others 2011). An important mechanism that could explain such bistability is a feedback between grass and fire (Staver and others 2011b; Murphy and Bowman 2012; Dantas and others 2016). In savannas, grasses act as fuel for fires (Bond 2008). Savanna trees are adapted to fires by allocating many resources to developing thick barks (Keeley and others 2011). However, this resource allocation implies investing less resources in leaf area than forest tree species do (Silva and others 2013). Because the open canopy in savannas allows fires to occur and the closed canopy in forests suppresses fire, both states are self-stabilizing (Murphy and Bowman 2012). Drought, however, may destabilize a forest by allowing it to burn. Burned forests can be invaded by grasses that facilitate recurrent fires (Brando and others 2014), which may trigger a transition to savanna. The reverse transition occurs if a fire-free interval exceeds a threshold, determined by available resources, above which forest-tree species can form a closed canopy (Hoffmann and others 2012). The theory of alternative stable states (Scheffer and others 2001) predicts that bistable systems can undergo a transition from one state to the other if environmental conditions (such as MAP) reach a so-called ‘tipping point.’ After a transition has occurred, the system may not easily recover to the original state even if the change in conditions is reversed, a phenomenon called hysteresis. Hysteresis implies that there is a range of conditions at which either state can be present and large stochastic disturbances could cause transitions between the states (the bistability range). Therefore, in a bistable forest-savanna system, boundaries are difficult to predict based on external conditions (Moncrieff and others 2016). Hence, bistability of forest and savanna on the one hand, and the climatic-deterministic explanation of forest-savanna distributions on the other, seem to represent two fundamentally incompatible views.

Experiments to test hysteresis are challenging, especially on landscape-scale systems such as forests and savannas (Bowman and others 2015). So far, experiments on alternative stable states are either conducted in laboratory microcosms (for example, Dai and others 2013) or on isolated ecosystems (Carpenter and others 2011). Moreover, theoretical studies of alternative stable states often use simple models that assume well-mixed systems without explicit consideration of spatial components (Scheffer 2009). However, the behavior of such models can change considerably when local spatial interactions are taken into account (Wilson and others 1996; Bel and others 2012; Van de Leemput and others 2015; Villa Martín and others 2015). Spatial interactions are usually modeled in reaction–diffusion equations where the relevant state variables (such as biomass) have local dynamics (‘reaction’) and exchange between neighboring patches (‘diffusion’). In the case of forest and savanna, this biomass exchange between neighboring cells would capture, in a simple manner, mechanisms by which forest patches increase tree cover in neighboring savanna patches (for example, by means of seed dispersal) and by which savanna patches decrease tree cover in neighboring forest patches (for example, by means of fire spread). If patches of forest and savanna spatially interact in such a way, then boundaries between alternative stable states may become unstable, meaning that the two states cannot coexist on a local scale. This impossibility of local coexistence of alternative stable states is a generic result in partial differential equations, regardless of the specific mechanisms that cause bistability (Murray 2002; Van de Leemput and others 2015). However, whether destabilization of local coexistence occurs depends on how strong these mechanisms are relative to spatial interactions. With an increase in spatial interaction strength, the range of environmental conditions at which local coexistence of alternative states is possible becomes smaller (Van de Leemput and others 2015). In cases with very strong spatial interactions, it becomes inevitable that at given conditions the most stable alternative state will eventually dominate. This is because even if the landscape is occupied by the least stable state, a local invasion by the most stable state is sufficient to trigger a domino effect of a spatial transition to that state (Pomeau 1986; Wilson and others 1996). This implies that in a gradient of environmental conditions, the boundary between the alternative stable states is highly predictable as it occurs at the conditions where both states are equally stable (Wilson and others 1996; Van de Leemput and others 2015). It also implies that there is no hysteresis in response to changing environmental conditions as transitions between alternative states always occur when the conditions of equal stability are crossed (Figure 1). Environmental conditions of equal stability are referred to as the Maxwell point (after Maxwell 1875). In cases where spatial interactions are so strong that hysteresis is eliminated, the system is said to follow the ‘Maxwell convention.’ With increasing spatial interaction strength, hysteresis becomes more narrow. If spatial interactions are negligible, the system is said to follow the ‘delay convention’ (Gilmore 1981; Fort 2013). In such cases the patches can be in either state independently of the state of neighboring cells and the patches can be in either state over a range of conditions. Moreover, hysteresis will make the conditions at which boundaries between alternative stable states are found unpredictable.

Fig. 1
figure 1

The theoretical effect of spatial interactions on critical transitions in a bistable system. a Without spatial interactions between patches in an ecosystem, the patches undergo critical transitions between stable states (indicated by solid lines; unstable equilibria are indicated with the dashed line) when conditions cross tipping points (vertical arrows). Thus, the system displays hysteresis and therefore follows the ‘delay convention.’ b An increase in spatial interactions narrows the range of conditions at which patches can be in alternative stable states. At other conditions the least stable state (indicated by dotted lines) is not resilient against perturbations. Although the patches undergo critical transitions at different environmental conditions than without spatial interactions between them, the system still displays hysteresis. c With strong spatial interactions between patches, the least stable state is not resilient against perturbations: the invasion of at least one patch of the most stable state results in dominance of that state in the landscape. Transitions in both directions occur at the Maxwell point (double arrow). Thus, the system does not display hysteresis and therefore follows the ‘Maxwell convention.’ For this example, we used the model from Noy-Meir (1975), but the qualitative result is independent of the bistable model used. Note that our simplified illustration ignores the effects of stochasticity, which allows for back-and-forth transitions to occur in between the tipping points

Full size image

Which of the two conventions applies to forest-savanna dynamics is unclear. Forest-savanna dynamics can be characterized as patch dynamics (Goetze and others 2006; Dantas and others 2013) acting on the scale of several hundreds of meters (Favier and others 2012). When patches do not interact, the delay convention would apply, which means that we may expect a roughly equal distribution of local coexistence of patches of forest and savanna along a range of environmental conditions, corresponding to the bistability range (see Appendix S1 for a modeling example). If, however, there are significant interactions between patches, the range where the two states can coexist locally would become smaller than the bistability range for isolated patches. The narrowing of the range of conditions with local coexistence depends on the spatial interaction strength, with the ultimate possibility that there will be only one point where a stable boundary between both states is possible. In this ultimate case the Maxwell convention would apply, which means that we would find a distinctly peaked distribution of local coexistence of patches of forest and savanna around the Maxwell point.

In this study, we aimed to identify climatic Maxwell points and test whether forests and savannas in South America and Africa tend to coexist at these conditions. If they tend to coexist around a Maxwell point, then this would indicate that shifts between forests and savannas in both directions can be expected to occur when climatic changes cross a Maxwell point, implying the absence of hysteresis. If forests and savannas tend to coexist roughly equally throughout the bistability range, then this would be consistent with the classical tipping point model (the delay convention as in, for example, Van Nes and others 2014; Staal and others 2015). In that case, shifts from forests to savannas and vice versa can be expected at different climatic conditions, implying hysteresis even in spatially interacting landscapes. Because our study differentiates between these two extremes, it improves our understanding of the effect of climate on forest-savanna boundaries and the catastrophic behavior of forest-savanna shifts in response to climate change.

Methods

The stability of a system state can be described as its ‘potential energy’ (Strogatz 1994), which is often visualized as the depth of the system’s stability landscape (Figure 2). The stability landscape of forest and savanna at given climatic conditions can be inferred from the frequency distribution of tree cover at these conditions (Livina and others 2010; Hirota and others 2011). Hence, the Maxwell point, at which forest and savanna have equal potential energy (Van de Leemput and others 2015), should be the conditions at which the savanna and forest peaks in the tree-cover frequency distribution are equally high (Figure 2). We thus attempted to empirically identify Maxwell points for climatic variables for the forest-savanna systems of tropical South America and Africa. We considered MAP and rainfall seasonality. Independently, we investigated the effect of these variables on the probability of local forest-savanna coexistence.

Fig. 2
figure 2

The probability density of tree cover at the Maxwell point. a Bistability of forest and savanna occurs in the range of a control parameter such as mean annual precipitation. Stable states are depicted by the solid line and the unstable state by the dashed line in the bottom plane. Stability is also visualized by constructing a stability landscape (after Scheffer and others 2001) based on the physical idea of potential energy (Strogatz 1994). This potential of states is represented by the depth of the valleys; the tendency of a system to seek a state of lower potential is illustrated by balls (the system state) rolling down valleys. b The depth of the two valleys is equal when the control parameter is at the Maxwell point. c Potential is inversely proportional to the logarithm of the system state’s probability density distribution in space or time (Gilmore 1979; Gardiner 1985). Thus, at the Maxwell point, the two modes in the frequency distribution of tree cover should have the same height. However, the fraction of points per attraction basin, which is resilience as defined by Hirota and others (2011), need not be the same for the two basins. For a description of the model, see Appendix S1

Full size image

We took tree-cover observations from South America and continental Africa between 15°N and 35°S (Hirota and others 2011) from the MODIS Vegetation Continuous Field (VCF) Collection 5 dataset from 2009 at 250 m resolution (DiMiceli and others 2011). We excluded human-used areas, open water, and bare ground using the 2009 ESA Globcover dataset (values 11–30 and ≥190) at 300 m resolution. However, this is an imperfect exclusion of all human influence such as deforestation. We also excluded mountain regions (≥1500 m elevation) using the Shuttle Radar Topography Mission (SRTM) elevation data at 1 km resolution assembled by Hijmans and others (2005). The resulting continental datasets consisted each of more than 100 million observations. Precipitation data were taken from the Climate Research Unit monthly dataset at 0.5° resolution (Mitchell and Jones 2005). Seasonality of rainfall was measured as Markham’s Seasonality Index (MSI; Markham 1970):

$$ MSI = 100\frac{{\sqrt {\left[ {\sum\nolimits_{m} {P_{m} } \sin \left( {m\frac{2\pi }{12}} \right)^{2} } \right] + \left[ {\sum\nolimits_{m} {P_{m} } \cos \left( {m\frac{2\pi }{12}} \right)^{2} } \right]} }}{{\sum\nolimits_{m} {P_{m} } }} $$

MSI considers each monthly precipitation P m value to be a vector, where its direction is determined by the time of year (in month number m) and its length by the amount of rainfall. For each year, monthly vectors are summed; the summed vector is standardized by dividing by total rainfall over the year. This generates an index that is independent of MAP and ranges from 0 to 100 %, where 0 implies an equal distribution of rainfall over the months and 100 % implies that all rainfall occurs in a single month. Climatic variables were calculated for the period 1961–2002. All variables were resampled to a consistent spatial resolution of 250 m.

For each continent, we investigated whether the 250 m tree-cover data indicate bistability of forest and savanna within ranges of MAP (500–2500 mm) and MSI (10–70 %). We inferred the climatic ranges of bistability by analyzing the tree-cover probability densities using potential analysis (Livina and others 2010), which can be used to construct stability landscapes (Figure 2). We applied this potential analysis to 1 ‰ of the tree-cover data points from each continent (Hirota and others 2011). We smoothed the densities of tree cover with the MATLAB kernel smoothing function ksdensity (with a bandwidth according to Silverman’s rule of thumb; Silverman 1986). We inferred stable and unstable states of tree cover (minima and maxima in the potentials) for moving windows of the climatic variables (20 mm increments for MAP and 1 % increments for MSI), where we applied Gaussian weights to the climatic variables (standard deviation was 5 % of the entire range of the climatic variable) (Hirota and others 2011). As the MODIS VCF product used here is not without bias (Hanan and others 2014, 2015), it precludes detailed interpretations for specific tree-cover values. However, it has recently been validated that observed forest-savanna bimodalities are no artifact of the dataset (Staver and Hansen 2015), implying that the bias is smaller than the signal. The observed bimodalities have also been shown to be consistent with independent Landsat data (Xu and others 2015), providing sufficient support for the use of the dataset for our purpose.

To estimate Maxwell points of forest and savanna, we randomly took from each continent a sample of 1000 tree-cover data points for classes of MAP (500–2500 mm; bin size 100 mm with 20 mm increments) and MSI (10–70 %; bin size 6 % with 1 % increments). We determined for each sample whether the distribution was bimodal by fitting both a uni- and bimodal distribution using latent class analysis with the MATLAB function gmdistribution. We adopted a conservative approach in determining whether the distribution was bimodal by considering three criteria: the bimodal distribution should have both (1) a lower Akaike information criterion and (2) a lower Bayesian information criterion than the unimodal distribution, and (3) there should be a significant (α = 0.05) deviation from a unimodal distribution according to Hartigan’s dip test (Hartigan and Hartigan 1985). Furthermore, because we were interested in forest-savanna bimodality, we only considered bimodal distributions of which the modes were separated by a tree-cover value of 0.50, which is the approximate boundary between both basins of attraction. We calculated for each bimodal distribution the ratio of the height of the forest mode over that of the savanna mode (the bimodal ratio; Zhang and others 2003). Heights were determined after smoothing each distribution, again using the MATLAB function ksdensity. The climatic value at which the two peaks are on average equally high (that is, the mean bimodal ratio crosses 1) was used as estimator of the Maxwell point (Figure 2). To generate confidence intervals for the Maxwell point, we did the sampling and analysis 1000 times; we only report bimodal ratios when bimodality was found in at least 10 % of the samples.

At a Maxwell point, there may be stable local coexistence of both states independent of the strength of spatial interactions (Figure S1; Van de Leemput and others 2015). We determined for each continent and for each climatic variable if there was a relation with the probability that savanna and forest coexist locally. We determined local forest-savanna coexistence by testing if a grid cell at 0.5° resolution (about 50 km) contained a bimodal distribution of tree-cover data points at a scale of 250 m. From each 0.5° cell, we took a random sample of 1000 out of approximately 25,000 tree-cover data points. For each sample, we determined whether the distribution was bimodal, following the same procedure as for the analyses on tree cover sampled on a continental scale. However, we excluded 0.5° grid cells that were either at least 50 % occupied by human-used areas, water and bare ground, or located above 1500 m elevation. The proportion of grid cells at certain climatic conditions that was significantly bimodal we call the probability of local coexistence at these conditions. To test whether the probability of local coexistence varied within the bistability range (as inferred with potential analysis), we binned the grid cells (bin sizes 100 mm for MAP and 3 % for MSI) and performed χ 2 tests for independence. If a significant effect (α = 0.05) of the climatic variable on local forest-savanna coexistence was found, we performed post hoc pairwise χ 2 tests between the bin with highest probability of local forest-savanna bimodality and the other bins. We thus estimated the 95 % confidence interval for the conditions at which the probability in local bimodality peaks (also see Appendix S2). To map the distribution of forest, savanna and bimodal areas, we used the dip test to determine whether the mode in each unimodal cell was significantly (α = 0.05) in the savanna range (mode <50 % tree cover), significantly in the forest range (mode ≥50 % tree cover) or significantly in neither.

Results

In a range of MAP, the tree-cover frequency distribution was bimodal, corresponding with forest and savanna states (Figure 3a–b). This bimodality of tree-cover values of 80 and 10 % was present in the MAP ranges approximately 1100–2000 mm (South America) and approximately 1300–2000 mm (Africa). Along the bimodal ranges, the ratios of the height of the forest modes over the savanna modes (bimodal ratio) of the samples increased, indicating increasing stability of forest relative to savanna with increasing precipitation. Forest and savanna modes were equally high at 1760 mm MAP for South America (Figure 3c) and 1580 mm MAP for Africa (Figure 3d), which we consider an estimation of the Maxwell points of equal stability of forest and savanna. We could not explain the apparent decrease in forest stability in Africa around 2000 mm, but it could be an artifact of deforestation in the wet area of coastal West Africa (Hansen and others 2013), which may have been incompletely excluded from the data. We therefore assume in our analysis of the probabilities of coexistence of forest and savanna that their bistability range in Africa extends up to 2000 mm MAP.

Fig. 3
figure 3

Bimodality of tree cover depending on the mean annual precipitation (MAP), related to Maxwell points and the probability of local forest-savanna coexistence. Potential analysis for South America (a) and Africa (b) indicates stable (solid dots, local minima in the potential) and unstable (open dots, local maxima in the potential) states. Maxwell points (red dashed lines) for South America (c) and Africa (d) are estimated as the MAP where the height of the forest mode divided by the savanna mode crosses 1 (shown with 95 % bootstrap confidence interval; note the log scale). Probabilities of local forest-savanna coexistence, estimated as the proportion of bimodal grid cells at 0.5° resolution for South America (e) and Africa (f) are given as running means. Shaded areas indicate the 95 % confidence interval for the location of the peak in the probability of local coexistence. Bin size was 100 mm with 20 mm increments for c–f (Color figure online)

Full size image

We tested for bimodality within 0.5° grid cells to find the probability of local coexistence of forest and savanna. For South America (Figure 3e), there was a significant peak in the probability of local coexistence between 1300 and 2000 mm MAP (Appendix S2). This range contained the estimated Maxwell point at 1760 mm MAP. For Africa (Figure 3f), the probability of local coexistence also peaked significantly between 1500 and 2000 mm MAP (Appendix S2). This range contained the estimated Maxwell point at 1580 mm MAP.

Also in a range of Markham’s Seasonality Indices (MSI, seasonality in precipitation), the tree-cover frequency distribution sampled per continent was bimodal (Figure 4a–b). This bimodality was present in the MSI ranges 12–55 % (South America) and 12–35 % (Africa). Along the bimodal ranges, the height of the forest mode over the savanna mode of the samples decreased, indicating decreasing stability of forest relative to savanna with increasing rainfall seasonality. These bimodal ratios crossed 1 at MSI = 50 % for South America (Figure 4c) and at MSI = 24 % for Africa (Figure 4d), which we considered estimations of Maxwell points of equal stability of forest and savanna.

Fig. 4
figure 4

Bimodality of tree cover depending on rainfall seasonality, related to Maxwell points and the probability of local forest-savanna coexistence. Potential analysis for South America (a) and Africa (b) indicates stable (solid dots, local minima in the potential) and unstable (open dots, local maxima in the potential) points. Maxwell points (red dashed lines) for South America (c) and Africa (d) are measured as the Markham’s Seasonality Index (MSI) where the height of the forest mode divided by the savanna mode crosses 1 (shown with 95 % bootstrap confidence interval; note the log scale). Probabilities of local forest-savanna coexistence, estimated as the proportion of bimodal grid cells at 0.5° resolution for South America (e) and Africa (f) are given as running means. The shaded area in e indicates the 95 % confidence interval for the location of the peak (there was no significant effect in f). Bin size was 6 % with 1 % increments for c–f (Color figure online)

Full size image

For South America (Figure 4e), there was a significant peak in the probability of local coexistence between MSI = 46–55 % (Appendix S2). This range contained the estimated Maxwell point at MSI = 50 %. For Africa (Figure 4f), the probability of local coexistence did not peak significantly within the inferred bistability range (Appendix S2).

Local coexistence of forest and savanna was mostly found at boundaries of the Amazon and Congo rainforests (Figure 5). We show the MAP and MSI at these locations in Figure S3 and provide Google Earth files of our results as Appendix S6.

Fig. 5
figure 5

The observed distribution of unimodal forest, unimodal savanna, local forest-savanna coexistence, and unimodal savanna or forest in South America and Africa (between 15°N and 35°S) within 0.5° cells (n = 1000 for each cell). Excluded cells are either at least 50 % occupied by human-used areas, water and bare ground or are located above 1500-m elevation. In South America, 14 % of the 3612 non-excluded cells are classified as containing local forest-savanna coexistence, which are mostly found around the Amazon rainforest. In Africa, 3 % of the 4720 non-excluded cells are classified as containing local forest-savanna coexistence, which are mostly found around the Congo rainforest. ‘Unimodal savanna or forest’ are unimodal cells with its 95 % confidence interval in both the savanna (tree cover below 0.50) and forest domains (tree cover at least 0.50). Digital files of these maps, suitable for display in Google Earth, are provided as Appendix S6 (Color figure online)

Full size image

Discussion

Our results confirm that tropical forests and savannas are bistable over a broad range of climatic conditions (Hirota and others 2011; Staver and others 2011b; Murphy and Bowman 2012). We also found local coexistence of patches of forest and savanna over a range of climatic conditions. However, the probability of local coexistence was not always evenly distributed along the range of climatic conditions at which we inferred bistability. This effect can be interpreted in the light of dynamical systems theory (Appendix S1; Van de Leemput and others 2015), stating that spatial interactions between patches narrow down the range of conditions with stable coexistence of the two states. Eventually, the most stable state would dominate if we assume that local invasions of forests and savannas are always possible and the system is not dominated by stochasticity. Spatial interactions between patches could thus eliminate (part of) the hysteresis. If this is the case, we would expect a significant decline of the probability of forest-savanna coexistence further away from a Maxwell point, even within the range of conditions at which both states are stable when they are in isolation. This effect becomes more pronounced with increasing spatial interaction strength (Figure S2).

Within the range of MAP at which we inferred bistability of forest and savanna, we observed a significant effect of MAP on the probability of local forest-savanna coexistence for both South America and Africa. In the case of rainfall seasonality (MSI), we observed a significant effect on local coexistence only for South America. The climatic conditions at which the probability of local coexistence peaked coincided with the conditions where the forest and savanna modes in the tree-cover frequency distribution were equally high. Based on dynamical systems theory (Gilmore 1979; Gardiner 1985; Livina and others 2010), we reason that such conditions correspond to the Maxwell point, where two alternative states are equally stable (Van de Leemput and others 2015; Villa Martín and others 2015). We used a simple model to illustrate the well-known principle (Murray 2002) that if spatial interactions between patches of two alternative stable states are very strong relative to the states’ local self-stabilizing mechanisms (that is, the Maxwell convention applies to the system), then local coexistence of the states under constant conditions is only stable at the Maxwell point. However, in our empirical results, the confidence intervals for the peaks in the probability of local coexistence were consistently wider than those for the Maxwell points. Thus, the probability of local forest-savanna coexistence peaked at a wider range of conditions than merely those at which the two states were estimated to be equally stable.

There can be several explanations for the finding that coexistence occurs at a wider range of conditions than the Maxwell point (Appendix S5). Firstly, it is possible that spatial interactions between patches are often too weak to establish the most stable state; secondly, climatic conditions of equal stability may differ between locations due to other environmental factors; thirdly, fine-scale spatial heterogeneity increases the occurrence of coexistence; and fourthly, tree cover is often not in equilibrium, due to either natural or human disturbances.

The first explanation for the wider range of conditions at which coexistence occurs implies that we cannot fully confirm that the Maxwell convention applies to tropical forests and savannas. The results therefore do not refute that transitions between the two states may occur at tipping points.

The second explanation relates to the fact that the Maxwell points that we determined represent the average value for many locations. In reality, stability of forest and savanna depends on other environmental factors than only precipitation. Therefore, variation in other relevant factors (such as soil type) would also cause variation in the Maxwell point for a certain climatic variable (such as MAP). This variation in the environment would also inflate our estimates of hysteresis and may potentially give a false impression of alternative stable states (and its prediction of sudden ecological change) (Van Nes and others 2014).

The third explanation relates to local spatial heterogeneity, which also expands the range of climatic conditions at which local coexistence could be found (Appendix S1; Favier and others 2004, 2012; Van Nes and Scheffer 2005). Examples of such heterogeneity include local differences in water availability, differential responses of species to water stress and differences in soil properties (Levine and others 2016). Examples include gallery forests, which occur in savanna landscapes where local water availability is high. Environmental variability on a scale in the order of the tree-cover data (250 m) could thus deterministically cause local bimodality of tree cover. However, there is also evidence that local heterogeneity in soils may be enhanced by the vegetation itself and therefore be not only cause, but also effect of vegetation patchiness. It has, for instance, been shown that expansion of forest trees over savannas and grasslands increases soil fertility (Silva and others 2008; Silva and Anand 2011; Paiva and others 2015). More fertile soils in turn favor the development of a closed canopy, which suppresses fire (Hoffmann and others 2012). Thus, a positive feedback between forest trees and soil fertility may reinforce forest-savanna bistability (Hoffmann and others 2012; Silva and others 2013; Staal and Flores 2015) and therefore enhance coexistence in the absence of predetermined environmental heterogeneity.

The fourth explanation is that the distribution of forests and savannas is often not in equilibrium, because these systems are dynamic and always subject to disturbances. These disturbances can not only be anthropogenic (such as deforestation, which has been imperfectly excluded from the data) but also natural. The latter is indicated, for instance, by the large variability in tree cover observed both in the field and from satellites (Figure S4; Sankaran and others 2005; Hirota and others 2011; Ratajczak and Nippert 2012; Van Nes and others 2012). Stochasticity may occasionally push the system away from the most stable state and can thus be expected to expand the range of climatic conditions at which forest-savanna coexistence is found. The closer conditions are to the Maxwell point, the longer it takes for a system to move to its equilibrium (Scheffer and others 2015b; Van de Leemput and others 2015). Thus, even when forest and savanna can be alternative stable states, part of the forest-savanna coexistence we observe may be transient ecosystem configurations. Note, however, that these effects of temporal as well as spatial variability in conditions do not affect whether or not the probability of coexistence differs within inferred bistability ranges.

In South America, the effects of MAP and MSI on the probability of local forest-savanna coexistence within the inferred bistability ranges were both strongly significant, whereas in Africa only MAP had an effect and it was weakly significant. The weak correspondence between climate and forest-savanna coexistence in Africa could relate to any of the four explanations listed above. However, we conjecture that in Africa there are weak spatial interactions between forest and savanna patches relative to the local self-stabilizing mechanisms of these patches (Appendix S1). Specifically, the growth of grasses is faster and responds more strongly to MAP in Africa than in South America (Lehmann and others 2014). Indeed, fire frequencies are higher in Africa than in South America (Lehmann and others 2011). We therefore hypothesize that the grass-fire feedback is faster in African savannas than in South American ones, which could explain why we found that forest-savanna boundaries more strongly follow the Maxwell point in South America. There, the estimated Maxwell point for MSI coincided very well with the independently observed peak in the probability of local forest-savanna coexistence. This peak was also more narrow and pronounced for MSI than for MAP, which suggests that rainfall seasonality more strongly affects the South American forest-savanna distribution than average rainfall does. By contrast, our results suggest that average rainfall more strongly affects the African forest-savanna distribution than rainfall seasonality does. This contrast resonates with the yet unanswered question why savannas in South America can persist under wetter conditions than in Africa (Lehmann and others 2014). Part of the explanation can be that African forests are better adapted to dry conditions than South American ones. Indeed, analyses of satellite radar data, which can indicate canopy mortality, showed that African forests are resilient against occasional drought (Asefi-Najafabady and Saatchi 2013), whereas Amazon forests have higher drought sensitivity (Saatchi and others 2013).

Most cases of local forest-savanna coexistence that we observed occurred at the boundaries of the Amazon and Congo rainforests, whereas it was relatively rare at other locations. This is consistent with the modeling result that local coexistence of forests and savannas is not always stable, even in the inferred bistability range. However, it is possible that our conservative approach in detecting statistically significant bimodality on the basis of 0.5° grid cells has contributed to its scarcity. Furthermore, many 0.5° cells in the Brazilian savanna (‘cerrado’) were excluded from the analysis due to human impacts, which may have caused an underestimation of the occurrence of local forest-savanna coexistence in this region. On the other hand, we may have included human impacted areas at the boundaries of the rainforests and thus overestimated the coexistence at these boundaries. Local forest-savanna coexistence was about twice as probable in South America as in Africa. This may have been caused by a difference in the size of the continents’ major rainforest, as the circumference of the Amazon relative to the area of South America is about twice as large as the circumference of the Congo forest relative to the area of Africa.

We have presented a novel way of identifying the relative stability of forests and savannas, building upon modeling results (Appendix S1; Van de Leemput and others 2015; Villa Martín and others 2015). The climatic conditions of average equal stability, at which a boundary between alternative stable states in theory stabilizes, can be detected from peak heights in bimodal frequency distributions of tree cover sampled on a continental scale. This differs from the definition of resilience by Hirota and others (2011). According to that definition, resilience of a state is the probability that an observation of tree cover falls in the state’s basin of attraction, where the cut-off between savanna and forest was set at 60 % tree cover. Our approach does not require the definition of a cut-off between basins of attraction for inferring stability.

In summary, our results support previous studies that have suggested that tropical forests and savannas can be alternative stable states, thereby opposing climate determinism in the distribution of these biomes (for example, Hirota and others 2011; Staver and others 2011a; Staver and others 2011b; Favier and others 2012; Murphy and Bowman 2012; Staal and Flores 2015). Thus, our results indicate that changes in precipitation regimes could cause spatial shifts of the boundaries between forests and savannas when tipping points are crossed. However, both theoretical and empirical results of our study also indicate that forests and savannas may have less hysteresis than suggested by previous studies, which did not account for spatial interactions between patches of the two alternative stable states. An implication of our results is then that invasion of a low tree-cover state into a forest, facilitated for instance by logging, could cause the spread of a savanna-like landscape into the forest. Thus, understanding and preventing forest-savanna transitions requires the consideration of different scales, as changes in climate and in tree cover could together cross critical points (Staal and others 2015) after which undesirable biome shifts become inevitable (Scheffer and others 2015a).