www.biogeosciences.net/7/1237/2010/ © Author(s) 2010. This work is distributed under the Creative Commons Attribution 3.0 License.

Abstract. Terrestrial vegetation influences climate by modifying the radiative-, momentum-, and hydrologic-balance. This paper contributes to the ongoing debate on the question whether positive biogeophysical feedbacks between vegetation and climate may lead to multiple equilibria in vegetation and climate and consequent abrupt regime shifts. Several modelling studies argue that vegetation-climate feedbacks at local to regional scales could be strong enough to establish multiple states in the climate system. An Earth Model of Intermediate Complexity, PlaSim, is used to investigate the resilience of the climate system to vegetation disturbance at regional to global scales. We hypothesize that by starting with two extreme initialisations of biomass, positive vegetation-climate feedbacks will keep the vegetation-atmosphere system within different attraction domains. Indeed, model integrations starting from different initial biomass distributions diverged to clearly distinct climate-vegetation states in terms of abiotic (precipitation and temperature) and biotic (biomass) variables. Moreover, we found that between these states there are several other steady states which depend on the scale of perturbation. From here global susceptibility maps were made showing regions of low and high resilience. The model results suggest that mainly the boreal and monsoon regions have low resiliences, i.e. instable biomass equilibria, with positive vegetation-climate feedbacks in which the biomass induced by a perturbation is further enforced. The perturbation did not only influence single vegetation-climate cell interactions but also caused changes in spatial patterns of atmospheric circulation due to neighbouring cells constituting in spatial vegetation-climate feedbacks. Large perturbations could trigger an abrupt shift of the system towards another steady state. Although the model setup used in our simulation is rather simple, our results stress that the coupling of feedbacks at multiple scales in vegetation-climate models is essential and urgent to understand the system dynamics for improved projections of ecosystem responses to anthropogenic changes in climate forcing.


Introduction
Climate exerts a dominant control over distribution of terrestrial plants and surface properties (e.g.Woodward, 1987) that in turn affect the climate by changing the atmospheric water and energy budget.
It has been hypothesized that vegetation-climate feedbacks can explain regime shifts of vegetation on multidecadal time scales over the last century (e.g.Foley et al., 2005), as well as those that occurred in the Mid-Holocene (Kropelin et al., 2008;Brovkin and Claussen, 2008), as inferred from simulated multiple equilibria and desertification thresholds in the Sahel/Sahara region (Claussen and Esch, 1994;Brovkin et al., 1998;Liu et al., 2006;Dekker et al., 2007;Wang and Eltahir, 2000a, b;Zeng and Neelin, 2000).These multiple equilibria in climate and vegetation, in which an actual state depends on the history of the system, are caused by positive (amplifying) feedbacks between vegetation and climate.
There are numerous climate-relevant feedbacks between vegetation and atmosphere.
Biogeophysical feedbacks operate through modification of (i) surface albedo, (ii) evapotranspiration, and (iii) surface roughness by vegetation (e.g.Barry and Chorley, 2003).Although vegetation-climate feedbacks may originate at the local scale of vegetation, they might affect large-scale circulation patterns due to the impact on moisture recycling and temperature gradients.Vegetation-climate feedbacks can be both positive (reinforcing change) and negative (counteracting change) with different strengths, can interact at different scales and will depend on geographical position.
Whether positive biogeophysical vegetation-climate feedbacks can lead to regime shifts for vegetation and climate is debated in current literature.Bistable vegetation-climate states were found in models for the Amazon by Oyama and Nobre (2003) and for Central-Asia and North-Africa by Claussen (1998), while other modelling studies suggest an absence of multiple vegetation-climate equilibria in the climate-forest system (Brovkin et al., 2009;Levis et al., 1999).Also it has been claimed that the occurrence of multiple equilibria could be a model artefact, caused by using discrete vegetation classes (Kleidon et al., 2007).
A frequently applied method to test the resilience and bistability is to impose perturbations to the system.In this study we aim to assess global multiple stable equilibria in climate and vegetation which can be explained by vegetation-climate feedbacks.To simulate a continuous response of vegetation to climate, we used a simple vegetation model with one class of vegetation (Kleidon, 2006).
First, we checked whether the model has multiple vegetation-climate states.Second, we tested the hypothesis whether multiple stable equilibria are caused by biogeophysical feedbacks solely related to different initialisations of vegetation, excluding discrete vegetation classes.
We used the intermediate complexity model Planet Simulator (PlaSim) (Fraedrich et al., 2005a, c) and performed transient runs with minimum and maximum initial vegetation density, ran the model to steady state and made global maps of areas with multiple vegetation-climate equilibria.To test the stability of these equilibria, stepwise perturbations on vegetation states were performed.This enables us to construct geographically explicit susceptibility maps highlighting regions which are highly sensitive to small perturbations in the model.

Planet Simulator
The numerical experiments were performed with a Model of Intermediate Complexity, the Planet Simulator (PlasSim) (Fraedrich et al., 2005a, c), which simulates the dynamics of the global atmosphere based on parameterizations of aerodynamics, vertical and horizontal diffusion, radiation, moist processes and the effects of clouds.
In this setup PlaSim is forced with a prescribed present-day climatology of sea surface temperatures (SST's), carbon dioxide concentration, solar radiation land and sea-ice (Rayner et al., 2003).As the dynamic vegetation model, SimBA (Simulator for Biospheric Aspects, Kleidon (2006), Table 1.
Numerical experiments: Names, description and simulation time of model runs.D eq,i and G eq,i are the equilibria biomasses of a cel (i) at the end of the D (Minimum) and G (Maximum) run respectively.see also www.mi-uni-hamburg.de/plasim) is used.Carbon uptake by gross primary production (GPP) depends on incoming photosynthetically active solar radiation, surface temperature and soil moisture status and is calculated as the minimum of a light limited or water limited GPP.The change in biomass per time step is calculated as net primary production minus respiration and litter production.

Name
The PlaSim modules of the atmosphere and vegetation (SimBa) interact via the radiative balance and hydrological cycle.The atmosphere provides solar radiation, temperature and precipitation and SimBa calculates soil water status and actual GPP affecting albedo, LAI and roughness length, feeding back to the atmosphere.

Numerical experiments
All simulations are run at T21 horizontal resolution, which represents a spatial resolution of about 5.6 • × 5.6 • latitude/longitude, 10 vertical layers and integrated at 40 min time steps.Model output is stored as daily averages.Biogeochemical effects of altered carbon cycling were not accounted for because all simulations were forced with a prescribed present-day ambient CO 2 level (360 ppm).For every single perturbation simulation, all ice free land cells were equally perturbed with constant amounts of biomass (Table 1).In total eight simulations were conducted: one simulation initialised as desert (0 kg C m 2 for all ice free land cells; D-run), one initialised with a global green cover (12 kg C m −2 for all ice fee land cells; G-run) and six runs starting from perturbations of the D and G equilibria, from now D eq and G eq , reached after 1000 years (Table 1 and Fig. 1).
The size of the six perturbation runs differed (Table 1).The D − and G + run have a perturbation of respectively −1 kg C m −2 and +1 kg C m −2 biomass.Within the range of the D eq and G eq , four perturbations were made.Because the D equilibrium is lower compared to the G different scale of axis.426 were run for 600 years, which was long enough to reach steady state conditions.We have explicitly chosen fixed perturbations and not proportional perturbations because we are interested in a map of potential multiple equilibria.Using proportional perturbations will have the big disadvantage that grid cells with small amounts of biomass will be perturbed with small perturbations.It is not likely that for these cells multiple equilibria will be found.
If the biomass at a single individual cell resulted in a negative biomass after perturbation, the carbon content of this cell was set to 0 kg m −2 .From the perturbation runs, the susceptibility index (S) for every cell (i) is calculated as the difference between the new (B i ) and old (Bo i ) equilibria biomass state averaged over the last 100 simulation years divided by its perturbation (P): The following cases can be distinguished: -S i = 0: The resilience of the system is so high that a perturbation imposed to the system does not change the original equilibrium biomass (B i equals Bo i ), indicating that the state is stable.
-0<S i ≤1: In this case, the perturbation imposed to the system will move back to the direction of the original equilibrium, but will not reach it.This means that there is another stable steady state in vicinity of Bo i , and that the perturbation overshoots the new state.Positive feedbacks exist but are not strong enough to keep the perturbation sustainable.
-S i >1: The resilience of the system is so low that a change in biomass induced by the perturbation is amplified indicating a net positive vegetation-climate feedback (Bi>Bo+P).
-S i <0: For this case the change in biomass between the new and old state is in the opposite direction of the perturbation, implying a negative biomass anomaly.
The rationale for introducing the susceptibility index S i is as follows.The classical Lyapunov's stability theory only quantifies whether the steady state is stable or not depending on the sign of eigenvalues of the linearized system.In our case, all biomass states before or after perturbation are stable in the Lyapunov's sense and therefore we need another indicator that quantifies the sensitivity of the system to shift from one state to the other.

Multiple equilibria of biomass and climate
The mean biomasses of all ice free land cells during the 1000 years are shown for the D and G numerical experiments (Fig. 1a).Clearly distinct global multiple equilibria are found, for G (Green run) with a stable global mean biomass at carrying capacity G eq of 6.6 kg C/m 2 (σ =0.031 kg C m −2 ) and for D (Desert run), a stable minimum global mean biomass at carrying capacity D eq of 5. D is 0.947 m yr −1 (σ =0.014 m yr −1 ) and for G is 1.01 m yr −1 (σ =0.013 m yr −1 ), and mean simulated global temperature for D 285.1 K (σ =0.14 K) and for G 285.8 K (σ =0.11 K).
To quantify the strength and spatial extent of the possible vegetation-climate feedbacks a global map of biomass differences between the G eq and D eq is shown (Fig. 3).Clear distinct states are located in the Northern Hemisphere especially at high latitudes, in the monsoon area of Africa and in south-eastern Asia.A positive anomaly may constitute a net positive vegetation-climate feedback.Regions that show negative anomalies were not expected because this would imply that higher initialized biomass results in lower simulated equilibrium biomass.Nevertheless clear negative anomalies were found in east North-America and South-East Asia.Overall, in most regions such as South America including the Amazon no large differences of biomass equilibria were found, suggesting that the positive vegetation-climate feedbacks are not strong enough to develop multiple equilibria.In total 38 cells (8% of all land cells) show no significant differences between the G and D run.
The biomass anomalies (G minus D) averaged over the last 100 simulation years per land cell give a clear correlation (r = 0.86) with simulated total precipitation anomalies above land (Fig. 4a), with a maximum increase of 0.5 m y −1 .Cells with negative anomalies in biomass also give negative anomalies in precipitation.Total precipitation is the sum of large scale and convective precipitation.On average, convective precipitation above land was higher for the G-simulation, while convective precipitation above the ocean was higher for the D-simulation.Above land, mean large scale precipitation was almost equal for the G and D run, respectively 0.248 m y −1 and 0.242 m y −1 .
For temperature (Fig. 4b), the northern and southern high latitudes give a clear positive correlation between biomass and temperature anomalies of land cells averaged over the last 100 simulation years (r = 0.84).As light limited gross primary production is reduced at lower temperatures, growing seasons are smaller resulting in limited biomass development for the D-simulation.At low and mid latitudes (45 N-45 S) a negative correlation is found between biomass and temperature anomalies (r = −0.5).The largest positive biomass anomalies are found in Central Africa caused by higher precipitation which in turn is negatively correlated to temperature.

Perturbations
After 1000 years, the G and D multiple equilibria were disturbed with six different perturbations (Table 1) and were run till steady state (Fig. 1b). Figure 5 shows the global mean multiple equilibria in steady state for all land cells, averaged over the last 100 simulation years for all eight simulations.As expected the D − and G + runs did not show any significant difference with the D and G equilibria for both biomass and climate.This implies that all global mean biomass states lower as D eq are attracted by D eq , and that all global mean biomass states higher as G eq are attracted by the G eq .
The D + , D ++ and G − and G −− runs, however, show significant differences with the D eq and G eq values for both biomass and precipitation.For temperature, significant differences were only found for the D ++ and G −− runs.These multiple equilibria are within the range of the minimum and maximum stable states for global mean biomass and climate.The new steady states in biomass, precipitation and temperature are also plotted in the phase space figures of Fig. 2.

Susceptibility map of vegetation-climate feedbacks
As a final step we calculated the susceptibility index (S i ) of the vegetation-climate susceptibility to perturbations with Eq. 1.In Fig. 6  ) per cell between G and D equilibrium (B i, G eq ) and (B i, D eq ) against precipitation (m y −1 ) and temperature anomalies for every land cell, averaged over the last 100 years.comparable distributions (see insets Fig. 6a and b).Cells with no significant (ns) differences between the old and new biomass states after perturbations, i.e. S i =0, were found for the positive perturbations in 12% of the land cells and for the negative perturbations in 18% of the land cells.These cells are located as well in deserts as in regions with high biomasses.
A S i larger than 1, which means that the perturbation is enforcing the vegetation-climate feedback caused by a net positive climate-vegetation feedback, i.e. a low resilience, is found for the positive perturbation in 21% of the land cells and for the negative perturbation in only 10%.For the positive perturbations, regions with S i >1 were found at high latitudes in the Northern Hemisphere, across Siberia, in the north eastern North America and Central Europe and at low latitudes in Central Africa and Asia.For the Southern Hemisphere only regions with S i >1 were found in South Africa.Note that not only the percentage land cover of S i >1 but also the regions differ substantially between the positive and negative perturbations.Because D eq has more regions with S i >1 (21%) then the G eq (10%), we conclude that the D eq has a larger sensitivity and lower resilience against perturbations as the G eq .
In both S i maps, S i <0 was simulated for 29% of the land cells.Regions with high negative S i are neighbouring regions with positive S i >1.Negative S i indexes cannot be caused by single vegetation-climate cell interactions alone, changes in spatial patterns of atmospheric circulation due to neighbouring cells are required to influence neighbouring cells.This implies that biogeophysical feedbacks at small scales (cell) affect large scale circulation patterns and in turn influence cell interaction further away.Values (Mean ± standard deviations) are derived over last 100 years of the simulations.Different letters denote significant differences between runs at p <0.05.Arrows in Fig. 5a show the change from perturbed state towards the new state.

Discussion
The coupled dynamic atmosphere-vegetation earth system model PlaSim possesses multiple vegetation-climate equilibria which depend on the initial vegetation state.Despite interannual variations in precipitation and temperature caused by internal variability of the atmospheric processes, the steady states for climate and biomass for the desert and green initialisation are clearly distinct.We conclude that multiple equilibria can be found, which are caused by the interplay of different positive and negative vegetation-climate feedbacks affecting temperature and precipitation.These equilibria differ not only in mean climatic state, but also in response to the perturbations.The global average surface temperature of ice free land cells is in the G steady state 0.7 K higher than in the D steady state.This annual temperature increase is comparable to the study of Brovkin et al. (2009), who found an increase of 0.7 K and a decline of 0.6 K in a forest and grassland simulation and comparable to the study of Gibbard et al. (2005) who found an increase of 1.3 K due to global replacement of current vegetation by trees.Note that our simulations did not account for interactive feedbacks with the ocean or the carbon balance.Besides, in our case the difference is between two steady states and not between climate enforced by forest and grassland changes.
The different perturbation simulations showed interesting behaviour on a global mean scale.From an ecological point of view, equilibrium points can be stable or unstable (e.g.DeAngelis and Waterhouse, 1987).The D − and G + perturbations show that indeed a minimum and maximum stable state of vegetation and climate are found, indicating a stable equilibrium in the D − and G + direction.The other 4 perturbations, however, have shown that different stable equilibria in global mean biomass and climate were simulated between the minimum and maximum stable states indicating also stable states.So, the global G and D equilibria are two stable equilibria separated by a number of other stable steady states emerged at the global scale from smaller scale interactions.We conclude that the global perturbations trigger the vegetation-climate system to new steady states between the G and D equilibria (Fig. 2).
The S i -maps show regions with high and low resilience.Regions with low resilience are found for S i >1, in which the change of biomass induced by the perturbation is enforced.For S i <0, if the change in biomass between the new and old biomass state is in the opposite direction of the perturbation.
From the S i -maps we conclude that more cells have a S i >1 if the D eq is perturbed by positive perturbations as compared to the negative perturbations from the G eq meaning that the D eq have more cells with low resilience.At the grid scale, the S i -map shows regions with high S indexes, i.e. low resilience, caused by a net positive vegetation-climate feedback, mainly for the Northern Hemisphere at high latitudes above 40 • N, for the monsoon areas of Africa and Asia and for the Southern Hemisphere in South Africa.For high latitudes, the light limited GPP is lower for the D run as for the G run, because simulated surface temperature is lower and simulated growing season is shorter for the D run.Because simulated GPP is reduced at low temperatures in the vegetation model, almost no build up of biomass is simulated in the D run.Functional relationships between temperature and GPP for boreal regions are often used to take into account reduced photosynthetic production due to lower temperatures (e.g.Kellomaki and Vaisanen, 1997;Bergh et al., 1998).Simulated temperature anomalies at high latitudes between G eq and D eq is up to 7 K comparable to temperature differences between bare ground and forest earth found by Gibbard et al. (2005), although we did use prescribed SST's.The S i -map shows for boreal conditions clear multiple equilibria due to strong net positive feedback strengths.
In the monsoon areas of Africa and Asia multiple equilibria were also found.Here clear positive correlations are simulated between anomalies in biomass and precipitation and negative correlations between anomalies in biomass and temperature.Although the albedo effect would lead to warming, an overall cooling of up to 3 K is found.This cooling effect by the vegetation is caused by increased moisture recycling resulting in increased evapotranspiration (Barry and Chorley, 2003) and in turn increased precipitation (up to 0.5 m/y) for some grid cells in Central Africa.A general cooling in the tropics (up to 4 K) due to vegetation was earlier found by for instance Fraedrich et al. (1999;2005b).In our experimental set-up, with prescribed SST's and sea ice, there is no feedback to the ocean, and therefore lower land temperatures directly decrease the gradient between ocean and land and subsequently decrease the strength and length of the monsoon.This can be the cause that multiple equilibria for the Sahel Sahara were not found in our study while it was simulated by e.g.Wang and Eltahir (2000a, b) and Zeng and Neelin (2000).The Amazon region mainly simulated stable conditions.In comparison to the African region, average simulated precipitation for the Amazon is far more (about 2.0 m/y) suggesting that precipitation is large enough to regenerate biomass with the dynamic SimBa model.In contrast, two stable equilibria for the Amazon were found by for instance Oyama and Nobre (2003), which were caused by large decrease of local precipitation in the deforested situation (Malhi et al., 2008) while other model studies did not show multiple equilibria for the Amazon (e.g.Cowling et al., 2008).
Interesting to observe are the negative S i indexes.These cannot be caused by single vegetation-climate cell interactions alone, but need to be caused by changes in spatial patterns of atmospheric circulation over neighbouring cells.Circulation patterns changes of surface pressure and near surface wind speeds are observed (results not shown) in North-East America, North-Europe, Siberia and Central Asia, while no differences were found for South-America and Asia.
For mid-latitudes, the S i -map, especially with negative perturbations, show zonely oriented synoptic patterns of positive and negative S i at a 1500-2000 km scale for Europe and Asia that could be associated with recycling patterns.For the tropics no clear S i <0 pattern was observed.Although impacts on cell interaction on large scale circulation patterns can be explained in a physical way due to changes in pressure, wind and temperature, the model-resolution is to coarse to go in more detail.Moreover, a global perturbation of all cells can change circulation patterns more easily as with regional perturbed experiments.
We conclude that small-scale vegetation-climate feedbacks can affect large scale circulation patterns constituting an extra feedback loop between small-scale and large scale vegetation-climate interactions.Coupling of feedbacks at multiple scales is essential to understand the climate system, especially in studies with anthropogenic land-use change (afforestation, deforestation, and woody biofuel plantation).This study also implies that it is reasonable to assume that the current vegetation-climate state is one of several possible equilibria meaning that several ensemble runs with different initial biomasses should be used to analyse further climate change.
Fig. 1.(a) Simulated mean global biomass with the G (maximum) and D (minimum) initialization of vegetation.G eq and D eq are mean global equilibrium densities.(b) Simulated mean global biomass with G-perturbations and (c) D-perturbations.Note different scale of axis.

Fig. 2 .
Fig. 2. Trajectories of G and D runs over the first 300 years of biomass (kg C m −2 ) and precipitation (m y −1 ) (a) and biomass and temperature (b).Dots represent yearly global averaged biomass, precipitation and temperature values.Clear global different D and G equilibria in biomass, precipitation and temperature are found.Different global mean steady states were found between the two equilibria D eq and G eq .

Figure 4 :
Figure 4: Biomass anomalies (kg C m -2 ) per cell between G and D equilibrium (B i,Geq ) and 441 (B i,Deq ) against precipitation (m y -1 )and temperature anomalies for every land cell, averaged 442 over the last 100 years.443 444

Fig. 5 .
Fig. 5. Mean global equilibria of biomass (kg C m −2 ) (a), precipitation (m y −1 ) (b) and temperature (K) (c) of all 8 runs.Values (Mean ± standard deviations) are derived over last 100 years of the simulations.Different letters denote significant differences between runs at p <0.05.Arrows in Fig.5ashow the change from perturbed state towards the new state.

Fig. 6 .
Fig. 6.Susceptibility map (S indexes) from D to G with positive perturbations (a) and from G to D with negative perturbations (b); ns means no significant differences between the old and perturbed biomass state, meaning S i =0.Insets show histograms of land cover distributions of S i .