A model investigation of vegetation-atmosphere interactions on a millennial timescale

Introduction Conclusions References

Recently, Davin and Noblet-Ducoudre (2010) found that the biogeophysical (i.e.albedo change and evapotranspiration) impact of global deforestation leads to global cooling of 1 K because the albedo effect is dominant over evapotranspiration.Lawrence et al. (2010) and Zhongfeng et al. (2009) found that for land cover change, surface hydrology and its interaction with vegetation and soils also has large impacts on the climate.Changes in land cover can also have remote effects: Chen et al. ( 2001); Snyder et al. (2004) and Avissar et al. (2004) have found that tropical deforestation is likely to induce changes in atmospheric circulation, and that these changes may have consequences on precipitation and temperature patterns on a global scale.
In view of the aforementioned climate effects of land cover change there is a need to understand the impact of large perturbations to land cover.Are there interactions between forests and the atmosphere (i.e.vegetationprecipitation feedback, vegetation-temperature feedback) that gets severely altered for large scale land cover change?N. Devaraju et al.: A model investigation of vegetation-atmosphere interactions Specifically, if deforested areas are abandoned and forests are allowed to re-grow, will the climate system go back to its natural state?How long will it take for the climate system to reach a new equilibrium in such a case?What determines the time scale?We will address these questions in this paper using a dynamic vegetation model coupled to an atmospheric general circulation model.Our model has representation for dynamic vegetation, interaction between atmosphere and land surface via temperature precipitation, soil moisture and plant water uptake but it lacks detailed representation for nitrogen cycle, evolution of atmospheric CO 2 , disturbance such as fires and pests, anthropogenic land cover change, and management of forests.Therefore, our investigation mainly aims to find out if interaction among atmosphere, land biophysics and dynamic natural potential vegetation will result in multiple equilibriums.
Multiple steady states in the atmosphere biosphere system can arise as a consequence of positive feedbacks between atmosphere and land surface.Claussen (1994Claussen ( , 1998) ) find that multiple vegetation-climate states can exist in the presentday semi-desert regions of northern Africa and western Asia.A modeling study by Levis et al. (1999) and Brovkin et al. (2003) find only one possible stable steady state in the climate system in the high latitudes while bi-stable vegetationclimate states are found in models for the Amazon by Oyama and Nobre (2003), and for Central-Asia and North-Africa by Claussen (1998).Wang and Eltahir (2000) study the stability of similar states in response to disturbances.Brovkin et al. (1998) uses a conceptual model to compare desert versus "green" equilibrium under parameter estimates typical of current climate and of mid-Holocene climate, respectively.Kleidon et al. (2007) find that multiple steady states occur in their earth system model of intermediate complexity only if vegetation is represented by a few vegetation classes.With an increased number of classes, the difference between the numbers of multiple steady states diminishes and disappears completely when vegetation is represented by 8 classes or more in their model.The two major positive vegetation feedbacks that result in the emergence of multiple steady states are related to the temperature limitation on productivity in high latitudes and the water limitation in the tropics and subtropics.
The possibility of multiple forest-climate states in the climate system has been investigated recently by Brovkin et al. (2009) in the Earth System Model of the Max Planck Institute for Meteorology (MPI-ESM).In their study, when the geographic distributions of vegetation in forest-only and grassland-only simulations are allowed to interactively respond to climate, both forest and grassland simulations converge to essentially the same climate state as in the control simulation after subsequent 500 years of interactive climatevegetation dynamics.This convergence suggests an absence of multiple climate-forest states in MPI-ESM.Dekker et al. (2010), using an Earth Model of Intermediate Complexity, PlaSim, find that model integrations starting from dif-ferent initial biomass distributions diverge to clearly distinct climate-vegetation states in terms of climatic (precipitation and temperature) and biotic (biomass) variables.Their simulations suggest that the boreal and monsoon regions have low resilience, i.e. unstable biomass equilibrium with positive vegetation-climate feedbacks in which the biomass change induced by a perturbation is further enhanced.Large perturbations trigger an abrupt shift of the system towards another steady state, and hence, Dekker et al. (2010) stress the importance of coupling at multiple scales in vegetation-climate models and indicate the urgent need to understand the system dynamics for improved projections of ecosystem responses to anthropogenic changes in climate forcing.
In the present study, we investigate the possibility of multiple states in the climate system by carrying out transient simulations using the NCAR atmospheric general circulation model CAM2 (Community Atmosphere model) coupled to a dynamic vegetation model and a mixed layer ocean model.Our objective is to use another model to simulate the vegetation dynamics to understand whether, when and how positive climate-vegetation feedbacks will lead to multiple steady states.We perform a millennial time scale simulation using a comprehensive coupled atmospheric general circulation model and a terrestrial biosphere model to investigate the possibility of multiple states.We show that the climatevegetation system has no multiple steady states in this modeling study but it takes nearly a millennium for the system to return to the initial equilibrium state due to long time scales of vegetation dynamics in the northern high latitudes.As a caveat, we note that our simulations are idealized and intended to understand the earth system behavior for possible multiple climate-forest states.This study is not intended to realistically represent current or future land cover change.

Model description
The coupled climate-vegetation model used is Integrated Biosphere Simulator (IBIS2) (Foley et al., 1996;Kucharik et al., 2000) coupled to Community Atmosphere Model 2.0 of NCAR (CAM2) (Collins et al., 2004).We used the finite volume configuration of the model with a horizontal resolution of 2 • latitude and 2.5 • longitude.There are 26 vertical levels.For this study, we coupled CAM2 to a mixed layer ocean and thermodynamic sea ice model, which allows for interactive surface for the ocean and sea ice components of the climate system.
Land surface biophysics, terrestrial carbon flux and global vegetation dynamics are represented in a single, physically consistent modelling framework within IBIS2.IBIS2 simulates a variety of ecosystem processes including energy, water, and carbon dioxide exchanges between soil, plants and the atmosphere, physiological processes of plants (photosynthesis and respiration), soil biogeochemistry, vegetation phenology including budburst, senescence, and dormancy of vegetation, plant growth and competition, nutrient cycling and soil physics.The coupled simulation of surface water, energy and carbon fluxes are performed on hourly time steps and integrated over the year to estimate annual water and carbon balance.The annual carbon balance of vegetation is used to predict changes in the leaf area index and biomass for each of 12 plant functional types, which compete for light and water using different ecological strategies.IBIS2 also simulates carbon cycling through litter and soil organic matter.When driven by observed climatological datasets, the model's nearequilibrium runoff, net primary productivity (NPP), and vegetation categories show a fair degree of agreement with observations (Foley et al., 1996;Kucharik et al., 2000).

Experiments
We first performed a long simulation with prescribed climatological sea surface temperatures (SST) and a round number of present day CO 2 concentration of 400 ppm to calculate the implied ocean heat transport which is needed to use CAM2 in slab ocean configuration.For this climatological-SST simulation, we used a soil carbon spin up factor of 40 and ran the model until biomass and soil carbon reached quasi equilibrium.A soil carbon drift of less than 0.1 Gt-C per year is used to define the quasi equilibrium state.The implied ocean heat transport is calculated after soil carbon reached equilibrium.Then, we used this spun-up state of the biosphere model in a 200 year mixed layer simulation until the soil carbon again reached quasi equilibrium.From this state, we performed two 1000-year mixed layer simulations as described below.
(a) a control simulation corresponding to present day conditions, (b) a land cover change simulation, denoted by "LCC", where we do not allow tree plant functional types (PFTs) to exist globally for 100 years; only grasses and shrubs are allowed.The biomass in tree leaves and fine roots is immediately (year 1) transferred to the litter pool, and the stem biomass becomes litter on a time scale of 10-50 years depending on the tree plant functional type.After 100 years, all PFTs are allowed to exist for the subsequent 900 years.By 100 years, majority of the biomass would be transferred to the litter pool and hence we choose this time scale for the deforested period.We choose 900 years for the re-growth period because that is the time it took for the coupled vegetation-climate system in this model to reach a quasi equilibrium state.The climate statistics presented below are the averaged values over the last 100 years of model simulations.The statistical significance here is tested using the student t-test with correction for lag-1 autocorrelation (Zwiers and von Storch, 1995).
Since our model is not a comprehensive earth system model with deep-ocean and ocean carbon cycle, we do not have the formulation for tracking the carbon and accounting for the carbon budget in this study.We have prescribed atmospheric CO 2 at a constant level (400 ppm) throughout the simulation period and hence it serves as an infinite reservoir for carbon for the terrestrial biosphere.In effect, we account for feedbacks due to biophysical effect of land cover change (e.g.albedo, evapotranspiration, roughness length changes) but omit the biogeochemical change since atmospheric CO 2 does not vary in response to land cover change.

Results
The temporal evolution of annual land mean key climate and terrestrial carbon cycle variables is shown in Fig. 1.In response to the instantaneous deforestation in LCC simulation, the land surface temperature decreases by 3 K (1.7 K over globe) averaged over the first 100 years and land mean precipitation decreases by 14.4 % (3.7 % over globe).The sign of the changes are suggestive of the dominance of albedo effect (deforestation increases the surface reflectivity) over changes in evapotranspiration and the associated cloud effects in this model.Globla mean temperature change (cooling) simulated in our study is comparable to another global land cover change study that included only biogeophysical effect (1.3 K in Gibbard et al., 2005) but higher than obtained in Bala et al. (2007) because this later study also included the warming effect of increased atmospheric CO 2 from deforestation (the carbon cycle effect).
When tree-PFTs are allowed from the deforested state, temperature and precipitation and other variables show a fast recovery for about 100 years (Fig. 1) followed by a slow recovery for the subsequent centuries .After 700 years of model integration with interactive climate-vegetation dynamics from the deforested state, the land mean annual surface temperature and precipitation in the LCC simulation slowly converges to essentially the same climate state (Fig. 1a, b) as that of control simulation.Global mean surface temperature in LCC, averaged over the last century, is higher than the control state by only 0.08 • C. Since the standard deviation in the control is of the same magnitude (0.07 • C in control and 0.09 • C in LCC), we conclude that the climate of the last hundred years of LCC are indistinguishable from the control.
However, the time series of carbon cycle variables specifically the carbon stocks (Fig. 1c, d), suggest that the terrestrial biosphere is still converging towards the new equilibrium state.Immediately following deforestation, Net Primary Productivity (NPP) declines by 22 % from 79 to 62 Gt-C per year since grasslands have, on an average, less NPP than forests in this model (Fig. 1e).The "step-like" behavior of NPP at year 100 is due to "step-like" increase in gross primary productivity since tree PFTs have higher LAI than grasses and shrubs.We used an exponential fit (Eq. 1, given below) and find that time scale of this rapid evolution is only about 2 years.There is also a slow recovery component which has a time scale of about 600 years.There is no succession dynamics in the model and NPP is mainly dependent on the biomass in leaves which has a time scale of 2 years in our model.We find that the time scale of slow component is dictated by the recovery of NPP in high latitudes which is similar to the vegetation fraction dynamics there as discussed below (Fig. 2).
Net Ecosystem Exchange (NEE) represents the net land carbon uptake and is defined as NPP minus soil respiration which includes both microbial decomposition and root respiration.It shows declines and increases of similar magnitude to NPP after instantaneous deforestation and immediately following re-growth of forests (Fig. 1f).The sharp decline of NEE by 18 Gt-C immediately after deforestation is mainly due to the sharp decline in NPP.The subsequent increase in NEE is mainly due to decrease in soil respiration (Fig. 1g).Re-growth of forests at year 100 leads to a jump in NEE because of the step-function-like increase in NPP.Subsequently, increase in soil respiration leads to the decline in NEE to control simulation levels.The time scale of NEE recovery is about 100 years which is primarily dictated by the decomposition time scale of dead stem biomass in the litter pool (Fig. 1h).
After deforestation, there is a large decline of about 800 Gt-C in biomass (Fig. 1c).Soil carbon increases (Fig. 1d) by only about 120 Gt-C during this period because most of the dead biomass is transferred to the litter pool which transfers only a fraction to soil carbon pool and the rest decomposes from the litter pool.When forests are allowed to re-grow, there is a rapid recovery in biomass; it increases from less than 100 Gt-C to more than 800 Gt-C within 100 years.After this period, it asymptotically reaches the control simulation after several centuries.Soil carbon stock also shows a rapid recovery but it has about 20 Gt-C more than control even after 900 years of re-growth.It should be noted that there is no conservation of total carbon stock (biomass plus soil carbon) in this study because atmospheric CO 2 is prescribed, providing an infinite reservoir for carbon to the terrestrial biosphere.
The differences in annual precipitation values over land between the two runs at year 100 after trees allow growing are low (0.3 mm d −1 , Fig. 1b).However the differences in biomass are large (Fig. 1c).This implies that the moisture recycling due to vegetation changes is very low in this model Biogeosciences, 8, 3677-3686, 2011 www.biogeosciences.net/8/3677/2011/setup.We also find similar evolution for tropical region which implies relatively small differences for the tropical regions in moisture recycling between control and LCC.Therefore, our model simulates low moisture recycling to vegetation changes and hence has low sensitivity to precipitationvegetation feedbacks.This may be one of the causes for not simulating multiple climate-forest states in the model.The time evolution of global and regional tree cover area fraction is shown in Fig. 2 We have fitted the tree cover fraction of LCC to an exponential form: where f is the fractional area of tree cover; A 0 ,A 1 ,A 2 are constants, t is time in years; and t 1 and t 2 are time scales in years.This exponential fit is applied to tree cover evolution simulated for global, tropical, mid latitude and high latitude regions.Table 1 lists the time scales thus obtained.When tree cover is allowed to recover, we see that the dynamics of tree cover converges faster in tropics (∼4.5 years).In mid latitudes, the dynamics is dominated by a 12 year fast time scale.In the high latitudes, a long time scale of 553 years dictates the evolution of tree cover suggesting that the long time scale in the global tree cover comes from the dynamics of boreal tree cover evolution.
To investigate the dynamics in the boreal region further, we plot the key climate and carbon cycle variables for the latitudinal belt 50 • N to 90 • N in Fig. 3.The evolution of all the variables are similar to the global-mean (Fig. 1) but the time scales are much longer.As found in other studies (Kleidon et al., 2007;Bala et al., 2006;Bonan, 2008), it is likely that the feedback between temperature (snow cover related albedo) and vegetation (tree cover extent) is primarily the cause for the longer time scales that is simulated for the high latitudes (Fig. 3).
The spatial pattern of temperature change in LCC relative to control averaged over the last 100 years shows that the significant changes are seen in high latitudes (Fig. 4).Over the globe, temperature changes over 11.8 % (12.0 % over land) during the winter (DJF) and 13.8 % (18.8 % over land) in the summer (JJA) are statistically significant at 5 % level.The LCC simulation produces in general more precipitation over land than the control simulation during both the DJF and JJA seasons (Fig. 4c, d), and there is a reduction in the tropics.Precipitation changes over globe are significant at the 5 % level over 5.6 % (8.3 % over land) and 7.59 % (9.4 % over land) in DJF and JJA seasons, respectively.In summary, Fig. 4 shows that spatial pattern of the physical (abiotic) climate system in LCC is almost indistinguishable from the control.
Figure 5 shows NPP, biomass and soil carbon in the LCC and control simulations and their difference averaged over the last 100 years.In the control and LCC simulations, NPP and biomass are larger in thickly vegetated regions of the world such as Amazon, central Africa, South and Southeast Asia, Europe and eastern North America.There is a large bias in simulated NPP and biomass in Australia when www.biogeosciences.net/8/3677/2011/Biogeosciences, 8, 3677-3686, 2011 Table 1.Values of the constants in the exponential fit (Eq. 1) for the dynamics of fractional tree cover (f ) plotted in Fig. 2  compared to a standard suite of global products characterizing the NPP based on satellite observations (Fig. 5 in Running et al., 2004).Soil carbon is also higher in places where NPP and biomass have larger values.In addition, soil carbon has maxima in colder places such as Siberia, Alaska and the Himalayas.The spatial pattern of NPP, biomass and soil carbon is similar in both control and LCC simulations.Differences in these variables are small and mostly not significant except soil carbon differences in some high latitude locations in Siberia and Alaska: the differences in NPP, biomass and soil carbon are significant over 3.9, 21.3, and 25.2 % of the land regions.Since the magnitudes of the differences are small, we conclude that the spatial distribution of the terrestrial carbon fluxes and carbon stocks (biotic component) in LCC in the last 100 years are also almost indistinguishable from control.
The areal extent (in percentage) of simulated dominant vegetation types in the last 100 years in LCC and control cases are given in Table 2.We found that the model simulates the locations of major biomes reasonably well; tropical evergreen and deciduous forests in the tropics, temperate forests in the mid latitudes and boreal forests in the high latitudes.However, there were major biases such as the simulation of tropical and temperate forests in desert regions such as Saudi Arabia, Australia and northern Africa and, consequently, an underestimation of deserts, Tundra and polar desert regions.We identified that a warm and wet bias in the control simulation was the main cause for these biases: the global and annual mean surface air temperature and precipitation in the control simulation are higher by 1.3 • C and 10 % when compared to NCEP reanalysis (Kistler et al., 2001) and GPCP (Adler et al., 2003), respectively.We use kappa statistics (Monserud and Leemans, 1992) to compare vegetation distributions in LCC and control.Kappa takes on a value of 1 with perfect agreement.It has a value close to zero when the agreement is approximately the same as would be expected by chance.Global comparison of vegetation distributions of LCC and control gives a kappa value of 0.93 (excellent agreement).Except for Tundra, Kappa for major biomes between LCC and control suggests either excellent or very good agreement (Table 2).The low value of Kappa for Tundra is partly related to the smallest area occupied by Tundra in the control experiment (changes are relatively bigger for Tundra) and partly because the vegetation dynamics has not yet reached equilibrium in the high latitudes (Fig. 2d).This Kappa statistics suggests that the spatial distribution of the vegetation state in LCC in the last 100 years is also almost indistinguishable from control.

Discussion
In this study, we examined the possibility of multiple climate-forest states using a terrestrial biosphere model IBIS2 coupled to NCAR CAM2.For this purpose, we performed a global deforestation experiment for 100 years followed by re-growth of forests.We find that there are no multiple climate-forest states in our model; the simulation with deforestation followed by re-growth converges to the control climate.Our conclusion is similar to a recent modeling study (Brovkin et al., 2009) which also found an absence of multiple climate-forest states in the Earth System Model of the Max Planck Institute for Meteorology, but different from another study by Dekker et al. (2010) which found multiple equilibrium states.In our study, we find that the climate system takes about 700 years to come back to the original natural state.This is despite the fact that our simulation did not have representation for deep ocean feedbacks; we have used only a mixed layer ocean model for representing the interaction between the atmosphere and oceans.The millennial time scale for recovery in our model is dictated by the vegetation dynamics in the high latitudes.Brovkin et al. (2009) uses a tiled structure of the land surface with 8 PFT's and two types of bare surface but no biogeochemical effects of vegetation cover changes are included Fig. 5. Comparison between LCC and control simulations for terrestrial biosphere variables.The first and second column panels represent last 100-year mean from the control and LCC simulations, respectively.The third column panels represent the difference between LCC and the control simulations.Hatched areas are regions where changes are statistically significant at the 95 % confidence level.Significance level is estimated using a Student-t test with a sample of 100 means and standard error corrected for serial correlation (Zwiers and von Storch, 1995).
in their model (JSBACH).Trees and grasses compete for free available area.IBIS2 uses 12 PFTs with grass layer underneath the tree canopy.The biogeochemical effects of vegetation cover changes are included.IBIS2 also has representation for soil biogeochemistry.Both models have differing representation for competition among PFTs.Whether vegetation dynamics has multiple states or not is apparently independent on these minor differences in the formulations between IBIS2 and JSBACH.However, multiple vegetation states be may be inherent to vegetation models with different formulations that are probably not be represented in IBIS2 and JSBACH (Dekker et al., 2010).
The dynamic vegetation model used in this study has some limitations: IBIS2 does not have representation for nitrogen and other nutrient cycles (Cramer et al., 2001;McGuire et al., 2001;Bala et al., 2007).IBIS2 model, in its current form, does not include a dynamic fire module (Foley et al., 1996).It does not account for changes in pest attack or grazing by animal in a changed climate.Suitable climatic conditions are sufficient for the existence of plant functional types in IBIS2: seed dispersal mechanisms which are crucial for reestablishment of forests are not represented.The real world, therefore, could have multiple climate-forest states, and this present modeling study is unlikely to have represented all the essential ecological processes (such as altered fire regimes, seed sources and seedling establishment dynamics) for the re-establishment of major biomes.
The atmospheric conditions affect vegetation productivity in terms of available light, water and heat, and different levels of vegetation productivity can result in differing energy and water partitioning at the land surface, thereby leading to different atmospheric conditions.Kleidon et al. (2007) find that if there are less number of discrete PFTs the climate conditions do not overlap and multiple vegetation equilibria can result.With an increased number of classes the difference between the numbers of multiple steady states diminishes and disappears completely when vegetation is represented by 8 classes or more in their model.Therefore, it is possible that our model with 12 PFTs has no multiple steady states.An analysis of the sensitivity of multiple equilibriums to the number of PFTs and other parameters in the model is beyond the scope of this paper and our future studies will investigate this.
Major limitations of this modelling study are the lack of deep-ocean dynamics, dynamic sea ice, and representation of biogeochemical effects of vegetation cover changes and ocean carbon cycle.As discussed before, the effect of changes in atmospheric CO 2 from deforestation and interactive vegetation is not modeled in this study.For instance, we have overestimated the amount of cooling immediately after deforestation on shorter time scales because we have not included the effect of CO 2 emission from deforestation.However, on centennial to millennial time scale the ocean biogeochemistry could buffer most of the atmospheric CO 2 changes induced by an altered land cover.Future modelling on the investigation of multiple climate-forest states should use coupled climate and carbon cycle models that will have realistic representation of these long term feedbacks in the climate system.
The results discussed in this paper are from a single modelling study.Climate model differ in their representation of physical processes and hence models show differing sensitivities to climate perturbations.Many of the local and remote effects depend heavily on the model structure and the simulated effects are therefore subject to have a wide range.Therefore, an intercomparision of multiple models in the future will be required to investigate robustness in the behaviour climate system.

Fig. 2 .
Fig. 2. Evolution of global, tropical (20 • S to 20 • N), mid latitude (20 • N to 50 • N), and high latitude (50 • N to 90 • N) annual mean fractional area of tree cover (Tropical evergreen, Tropical deciduous, Temperate and Boreal forests) in the control and LCC simulations.

Fig. 4 .
Fig. 4. Mean differences between LCC and control simulations for surface Temperature ( • C) in (a) winter (DJF), (b) summer (JJA), and precipitation (mm day −1 ) in (c) winter (DJF), and (d) summer (JJA).Mean differences are obtained by averaging over the last 100year.Hatched areas are regions where changes are statistically significant at the 95 % confidence level.Significance level is estimated using a Student-t test with a sample of 100 annual means and standard error corrected for temporal serial correlation(Zwiers and von  Storch,1995).

Table 2 .
Areal extent of dominant vegetation types and Kappa statistics.