A multi-model analysis of risk of ecosystem shifts under climate change

GVMs are developed to simulate changes in biomass, carbon turnover, water flows and ecophysiological functional strategies (woody or non-woody, evergreen or deciduous, needleleaf or broadleaf etc) on a coarse scale. While changes in these stocks and fluxes are interesting in themselves, they do not directly reveal the risk of shifts in ecosystems under pressure from climate change. Ecosystems are characterized by complex networks of interactions between species, communities and their local niches. However, here we suggest using changes in the biogeochemical state of the vegetated land surface, as simulated by GVMs, as a proxy for risk of such change.

Following [15], we argue that any of the following effects are indicative of risk that ecosystem structures will be affected: change in functional strategies of the vegetation present (ΔV); relative, changes in local biogeochemical carbon and water stocks and fluxes (c); absolute, local change in biogeochemical stocks and fluxes with respect to globally averaged changes in stocks and fluxes (g); changes in the relative (to one another) magnitude of key biogeochemical exchange fluxes (b); or change in interannual variability of biogeochemical stocks and fluxes (S(x,σ_x)). The variability term S(x,σ_x) is a normalized sigmoid function of the ratio of the respective components to their standard deviation in the reference period. We use the aggregate of these effects, based on the biogeochemical quantities listed in table S3 (where 'S' refers to the supplementary material, available at stacks.iop.org/ERL/8/044018/mmedia), as a proxy for the risk of ecosystem shift, Γ [15]. The multiple dimensions are first normalized using a sigmoid function, and then combined according to

$$\begin{eqnarray} &&\Gamma =[\Delta V+c S(c,{\sigma }_{c})+g S(g,{\sigma }_{g})+b S(b,{\sigma }_{b})]/4. \end{eqnarray} \tag{ 1 }$$

Each of the components c,b and g is constructed by first considering a set of biogeochemical fluxes and stocks, which are combined into a state vector for each grid cell, at each point in time, and then comparing the state vector for the same grid cell at a given time with that of a reference period. In the case that a GVM does not model changes in vegetation composition, the vegetation composition change component of the metric (ΔV) is not included in the calculation of Γ, resulting in an increased weighting of the other three components (c, b and g). A comprehensive description of the proxy can be found in [15].

Our correlation of Γ to risk of ecosystem change makes use of a space-for-time evaluation (i.e. using the difference between contemporaneous, but geographically separate states as a substitute for different points in time for the same location; see figure S1, available at stacks.iop.org/ERL/8/044018/mmedia). In the absence of both a well-founded theory of global-scale ecosystem changes under climate change, and a sufficient density and duration of instrumented test sites, space-for-time can give an indication of the degree of difference between ecosystem states implied by a given value of Γ [25]. We do not attempt to correlate Γ fully to risk of ecosystem change, but consider two thresholds for moderate and severe risk (Γ = 0.1 and Γ = 0.3 respectively), chosen in view of the space-for-time analysis.

Comparison of two identical vegetation states produces Γ = 0, whereas a replacement of a biome by a completely different biome produces a Gamma of nearly 1 (e.g. Γ = 0.98 for a change from rainforest to semi-desert) [14]. Analysis shows that when global vegetation modelled by a GVM (here, LPJmL) is mapped into 16 major biomes, the difference between any two of these biomes is never smaller than about Γ = 0.1 [14], which is the value characterizing the difference between closely related, but different biomes, e.g. the difference between a temperate coniferous and a mixed forest. Hence, we designate Γ = 0.1 as a substantial, but still moderate risk of shifts in ecosystem properties. Most between-biome differences are larger, with Γ = 0.3 characterizing, for example, the difference between a boreal evergreen forest and a temperate broadleaved deciduous forest (Γ = 0.32), or a warm wood- and shrubland and a tropical seasonal forest (Γ = 0.31). Since these biomes are substantially different, we designate changes of Γ > 0.3 as severe risk of ecosystem change. All other between-biome differences studied produce larger Γ [14]. For example, between a temperate grassland and Arctic tundra Γ = 0.57; between warm, woody savannah and a tropical rainforest Γ = 0.51; and between temperate and tropical vegetation Γ > 0.5. The table in figure S1 provides a comprehensive listing of Γ values between different biomes and is based on calculations performed using the LPJmL model [14]. It should be noted that a number of factors may combine to produce a particular value of Γ; the examples given should be used as an illustrative guide only.

In order to calculate Γ for each year, we calculate changes to the quantities listed in table S2 (available at stacks.iop.org/ERL/8/044018/mmedia) by comparing the running mean over a 30-year period centred on the year of interest and the average for the period 1980–2009, which is taken to be the baseline state of the ecosystem. The 30-year window ensures that year-to-year variability does not dominate the signal, favouring long-term shifts in the basic ecosystem properties. In addition, it allows for the required quantification of changes in the variability of the considered variables. Simulations were performed for the time frame 1980–2100. Where a contributing variable of Γ was not supplied by a model (see table S3), it was left out of the calculation. Where no dynamic vegetation composition was modelled, this component was not included in the calculation and the other components (carbon and water stocks and fluxes) were scaled up accordingly (i.e. in equation (1), ΔV = 0 and the factor in the denominator is 3). The risk of ecosystem changes is presented at different levels of global mean temperature change (ΔGMT, compared to 1980–2009 levels, which are in turn approximately 0.7 ° C above pre-industrial levels), thus contributing to the discussion about targets for limiting climate change.

This study reports results from seven GVMs, forced by the bias-corrected ISI-MIP climate data set [23] on a 0.5° × 0.5° grid (JeDi and JULES were simulated on a 1.25° × 1.85° grid) for four RCPs, with and without variable atmospheric CO₂ concentration. The bias correction method maintained the climate sensitivities of the GCMs (absolute change in monthly mean temperature and relative change in precipitation and other climate variables), whilst adjusting the absolute mean monthly climate variables to statistically match a historical data set. The method also corrects the daily variability of all climate variables to statistically match the observational data set. A detailed description of the bias correction method can be found in [26].

Each GVM performed a spin-up separately for each GCM, with the aim of bringing carbon and water pools into equilibrium for 1950 climate conditions. This spin-up was performed by recycling de-trended climate fields simulated by each GCM over the period 1950–1980 (JULES performed a single spin-up using HadGEM2-ES data). The length of the spin-up was determined individually according to the needs of each model. A CO₂ concentration of 280 ppm was adopted for the period until the year 1765. From 1765 CO₂ concentration was increased according to historical data until 2006 [27]. The atmospheric CO₂ concentration for 2006–2100 is prescribed by the four representative concentration pathways ([22] RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5) used to drive the GCMs. Table S3 gives details of the model runs performed by each GVM in combination with the GCMs, as well as ΔGMT for the period 2070–2099 compared to the 1980–2009 for each climate scenario. Further description of the climate data set can be found in [23].

The GVMs and their main characteristics are summarized in table 1. The Γ proxy was calculated using all available variables from each model. Four of these models (LPJmL, Hybrid, JeDi and JULES) simulate the changing vegetation composition due to changing climate and atmospheric conditions, albeit with unique classification of plant types into clusters (see the supplementary material for a full listing, available at stacks.iop.org/ERL/8/044018/mmedia). For JeDi, the vegetation changes, denoted ΔV, were not used to calculate Γ, since the large number of growth strategies considered were not easily categorized as plant functional types, as required by Γ. For the three models that provided dynamic vegetation composition (see table 1 for more details), Γ values were additionally calculated ignoring the ΔV component (results shown in the supplementary material), in order to limit the inconsistency in comparing these two classes of models.

In contrast to the coupled climate–land-surface modelling intercomparison as part of CMIP5, in this exercise the entire global surface was assumed to be covered by natural vegetation only (i.e. no anthropogenic land-cover), resulting in output of only climate driven rather than land-use-change-driven vegetation composition and biogeochemical changes. In this way we separate the pure climate impact from the extensive impacts of changing human agricultural practices and urbanization and focus on the vulnerability of natural ecosystems. For the two models where agricultural pastures was part of the default setup (ORCHIDEE and VISIT), results for the naturally vegetated portion of each cell were scaled up to 'fill' the cell. For display purposes, cells with <50% natural vegetation cover (based on the MIRCA 2000 data set of irrigated and rainfed crop areas [28]) are not shown on the world maps (coloured white). However, for the aggregation of the total land surface at risk of ecosystem change, all grid cells with >2.5% actual vegetation cover according to the GVMs are considered and weighted by the fractional area of natural vegetation. For grid cells where actual vegetation cover <2.5%, no Γ is calculated (coloured grey in the world maps) and the cells are not included in the aggregation calculations.

We quantify the risk of ecosystem change using the ecosystem shift proxy Γ described in [15]. In order to understand the drivers of the risk of ecosystem change, we first consider the relative contribution of changes to carbon stocks, carbon fluxes, water fluxes and vegetation composition changes (where available). This is done by calculating Γ as given in equation (1) for each GVM only for the chosen biogeochemical or vegetation properties. An example of the results are shown in figure 1, where the colour denotes the dominant component of Γ in 2084 for the median of the ecosystem-model ensemble for all RCP 8.5 GCM runs, and the intensity of the colour denotes the magnitude of that component (see also figures S13–S15, available at stacks.iop.org/ERL/8/044018/mmedia). It should be noted, that these results do not represent the output of any one GVM, and that results across the GVMs vary both in magnitude and spatial distribution. In most cases, changes to carbon fluxes between ecosystem and atmosphere (red colours in figure 1) mainly drive the risk of ecosystem change, which arises from the climate sensitivity of photosynthesis, respiration and plant water-use efficiency to atmospheric CO₂ concentrations increases. Steady increases with ΔGMT in net primary production across all models except Hybrid, and total vegetation carbon across all models, is driven largely by CO₂ fertilization of photosynthesis [46]. This is confirmed by the reduced risk of ecosystem shift when atmospheric CO₂ levels are held constant at present-day levels (see figures S2–S4, available at stacks.iop.org/ERL/8/044018/mmedia). This difference is more pronounced at low temperatures (e.g. LPJmL projects 5% natural vegetation at risk of severe ecosystem change with fixed, present-day CO₂ compared to 16% with changing CO₂ concentration at ΔGMT = 2 °C) compared to high temperatures (27–32% in both cases for LPJmL at ΔGMT = 4 °C). The Amazon is an exception, where changes to water fluxes and carbon stores also play a prominent role in most ecosystem models (see figures S13–S15). For ORCHIDEE changes in water fluxes dominate in the far northern latitudes and for both ORCHIDEE and VISIT, changes in water fluxes also dominate the region in and around the Democratic Republic of the Congo.

Zoom In Zoom Out Reset image size

No cells are dominated by vegetation change (ΔV, green) since four of the seven GVMs considered do not include dynamic vegetation composition changes. However, vegetation changes dominate Γ for some regions of the northern Boreal forests according to the three models with dynamic GVMs. Even when only these three models with dynamic vegetation are considered, only 4% of the land surface is dominated by vegetation changes (see figure S14), most notably due to shifts in the Boreal treeline northwards in Canada and Russia, greening in the Sahel, and transitions between grass and shrub biomes in western Russia and eastern India. The relatively small contribution of ΔV to the overall value of Γ is further confirmed by the 15% reduction in land surface at risk of severe change when the vegetation component of Γ is accounted for. However, comparison of the solid and dashed (with and without ΔV respectively) curves in figure 3 shows this effect does not result in a systematic offset between those models with dynamic vegetation and those without. Furthermore, figure S7 (available at stacks.iop.org/ERL/8/044018/mmedia) shows maps of Γ for one scenario run (HadGEM2-ES RCP 8.5), where Γ is calculated without (top) and with (bottom) the ΔV component. There does not appear to be a qualitative difference in the pattern of risk of ecosystem change, whilst there is a slight tendency for the Γ values to be higher when ΔV is not included, possibly resulting from a slower response of vegetation composition than biogeochemical ecosystem properties.

In order to assess the global risk of ecosystem change under different warming scenarios, we calculate the fraction of natural vegetation globally (by surface area) where the risk of ecosystem shift exceeds the moderate or severe threshold (Γ ≥ 0.1 and Γ ≥ 0.3 respectively) as a function of ΔGMT. We emphasize that the use of a single proxy may mask compensatory changes in its different components and information on processes that cause the change in Γ, but the aggregation facilitates evaluation and intercomparison. We find that ΔGMT is a reasonable proxy for the regional climate changes that drive shifts in the biogeochemical processes represented in the GVMs. In figure 2 the fraction of global natural vegetation at risk of change for the different RCP pathways is shown for single GCM–GVM combinations. Each point represents a comparison between the biogeochemistry of the 30-year period centred on the year shown and the baseline conditions (1980–2009). The Γ pathway is clearly dependent on the emissions scenario in the plot against time. However, when plotted against ΔGMT, the results are relatively independent of the RCP scenario (see also figure S8, available at stacks.iop.org/ERL/8/044018/mmedia), despite the strong effects of elevated CO₂ concentrations and potential inertia in the system.

Zoom In Zoom Out Reset image size

RCP 2.6, for which the climate and atmospheric CO₂ concentrations flatten out and slightly decline after 2050, is an exception to the RCP-independent ecosystem-model response, where the natural vegetation land area at risk of change continues to increase despite no further temperature rise (see also figure 1 in [23] and figure S8). This result suggests that adjustment of the ecosystem lags behind the changing climatic conditions, reflecting the long residence time of carbon stocks [2, 40]. This effect, manifested as a sharp up-turn in the fraction of natural vegetation at risk of change once ΔGMT ceases to rise, is more pronounced for the Γ > 0.1 curves, where slowly adjusting carbon stocks are sufficient to push Γ over the small changes threshold. Furthermore, the fall-off in atmospheric CO₂ concentration in RCP 2.6 results in a net CO₂ flux from vegetation to the atmosphere, despite only moderate climate change [41], which also registers as an increase in Γ through increased net primary production and ecosystem to atmosphere carbon flux. One could also suspect that inherently slow adjustment of vegetation composition plays a role, however the results in figure 3 exhibit no systematic offset between those models with and without (solid and dashed curves respectively) dynamic changes to vegetation composition, suggesting that vegetation changes are not the main driver of these changes.

Zoom In Zoom Out Reset image size

According to the multi-model ensemble depicted in figure 3, where results from all scenario combinations of RCP and GCM are included in the calculation, the median fraction of natural vegetation at risk of severe ecosystem change approximately doubles between ΔGMT = 2 and 3 ° C (from 13% (8%–20%) to 28% (20%–38%); bracketed values are inter-quartile ranges). However, agreement on the extent of natural vegetation at risk of severe change across the models is limited to a monotonically increasing trend as a function of ΔGMT. JULES projects the lowest fraction of naturally vegetated land area at risk of severe change, partially resulting from only small changes in vegetation composition compared to Hybrid and LPJmL. The highest projection comes from Hybrid, which probably results from lower mortality rates, especially in the northern latitudes, and significant reduction of vegetation carbon in the Amazon due to increased vapour pressure deficits arising from heat stress [46]. At lower temperatures significant uncertainty arises from both the GCM forcing and the GVMs, whereas at higher temperatures (where there are fewer available climate scenarios) the ecosystem models dominate the uncertainty budget (see also figure 5).

The spatial distribution of sensitivity of the biosphere to global warming is investigated in figure 4 using the median ΔGMT at which Γ first exceeds the severe change threshold (0.3) across all GVMs for RCP 8.5. Together with the far northern boreal forests of Canada and Russia, some regions of the tundra and shrublands of the Tibetan Plateau, where warmer winters result in longer growing seasons, are projected to be at risk of severe change as early as ΔGMT ≈ 1.5 °C. By ΔGMT ≈ 3 °C, the savannah regions in the Horn of Africa, along with large boreal forest regions across northern Canada and northern Russia, the Amazon forest and the grasslands of south-eastern India are projected to be at risk of severe change.

Zoom In Zoom Out Reset image size

However, individual maps of Γ for each GCM–GVM combination reveal large differences in the spatial distribution and intensity of ecosystem shifts (see figure S9–S12 for a full set of world maps, available at stacks.iop.org/ERL/8/044018/mmedia). In general, the uncertainty in the spatial distribution of risk of severe change coming from the GVMs is greater than from the GCMs (see section 2.4). For example, the projections of changes in the savannah region in the Horn of Africa range from little to small change for JeDi with the HadGEM2-ES RCP 8.5, to values above 0.3 for LPJmL for the same climate forcing. Projections of the spatial extent and magnitude of Γ in the Tibetan Plateau also vary significantly across the GVMs, despite a similar projection of global land fraction at risk, ranging from widespread, severe change for SGVM for IPSL-CM5A-LR RCP 8.5, through to only small risk of ecosystem change according to ORCHIDEE for the same forcing. These regional differences highlight the importance of such analysis to identify where the differences in process implementation across the GVMs lead to the greatest discrepancy in projections of ecosystem change. Global aggregations such as reported in figure 3 should therefore be treated cautiously, as they can obscure the fact that these global values arise from significantly different spatial distributions of change.

At ΔGMT = 4 °C (reached only using RCP 8.5; see table S3), the uncertainty in Γ arising from the GVMs is on average approximately twice as large as the uncertainty arising from the GCMs (shown as the standard deviation in Γ across GCMs and GVMs in the top and bottom panels of figure 5 respectively). The standard deviation across the ecosystem models is particular high in south-eastern North America, Turkey, south-east Australia, the Tibetan Plateau, south-east India, Canada and southern South America. In many of these regions, LPJmL and Hybrid project relatively large changes in vegetation composition, which cannot be mimicked in those models without dynamic vegetation composition. Additionally, differences across impact models are considerably more pronounced at ΔGMT = 4 °C compared to 2 ° C (see figure S16, available at stacks.iop.org/ERL/8/044018/mmedia), in line with the conclusions drawn from figure 2. Major differences in the behaviour of stomatal conductance in response to vapour pressure deficit may contribute to spread in the ecosystem models, in particular in tropical forests. Mortality is also handled differently across the models, which contributes to the factor of two difference in carbon residence times across the models [46], and contributes particularly strongly to the uncertainty in Γ projections in the northern latitude boreal forests. In addition, the shifts in this region are driven strongly by the response, among other processes, of water-use efficiency to atmospheric CO₂ concentration, which is shown in figure S3 (available at stacks.iop.org/ERL/8/044018/mmedia) to vary greatly across the ecosystem models. Despite these differences in the magnitude of Γ, at ΔGMT = 4 °C approximately 60% of all combinations of ecosystem model and GCM project Γ > 0.3 across the northern Amazon forest, southern India and the Tibetan Plateau (see figure S17 for the percentage of model runs agreeing on Γ > 0.3 at different levels of global warming, available at stacks.iop.org/ERL/8/044018/mmedia).

Zoom In Zoom Out Reset image size

Uncertainty in the projections of risk of ecosystem change arising from the GCMs dominates the Sahel region, where the median temperature at which the median Γ first exceeds 0.3 is 2.5 °C < ΔGMT < 3.5 °C and water fluxes dominate Γ (see figure 1). Γ in the monsoon region of India is also dominated by uncertainty arising from the GCMs, where the relative change in discharge compared to present day is also projected to increase by over 30% based on multi-model projections conducted as part of ISI-MIP [49]. It is interesting to note that very few of the regions of projected risk of severe ecosystem change correspond to regions projected to get drier under climate change. This most likely arises from the increased water-use efficiency of plants under elevated atmospheric CO₂, which help to counter this effect.

In many cases, the regions projected to be threatened by severe ecosystem change at ΔGMT > 2 °C coincide with regions that harbour exceptional biodiversity according to 'The Global 200' (compiled by [42] and comprising 142 terrestrial regions, see figure S18, available at stacks.iop.org/ERL/8/044018/mmedia). This set of regions was selected by analysing patterns of biodiversity to identify distinctive and irreplaceable biomes and biogeographic realms, whilst relying on the discretion of the authors to make a final selection. The assessment accounts for species richness and endemic species, as well as unusual ecological or evolutionary phenomena and rarity of global habitats. Regions of overlap between projected risk of severe ecosystem changes and regions of exceptional biodiversity are myriad, including the montane grass and shrublands of the Tibetan Plateau Steppe and the Kamchatka Taiga and Grasslands in north-eastern Russia. Several moist forest regions also overlap with regions of projected severe change, including the Guayanan Highland Moist forests and Amazon forests, and the moist and dry forests of north-eastern India. The Cerrado woodlands south-east of the Amazon basin, constituting one of the largest savanna forest complexes in the world, coincide with an area where the median projects risk of severe ecosystem changes at 3 °C < ΔGMT < 4 °C. In a similar temperature range, the Acacia Savannas in the Horn of Africa and the Northern Australian and Trans-Fly Savannas also overlap with model predictions of risk of severe ecosystem shifts.

On the regional scale, we can relate Γ to local changes in precipitation and temperature patterns. Taking the much-studied moist Amazon forest as an example [43, 44], in figure 6 we plot the percentage of ecosystem models that agree on a risk of moderate (left) and severe (right) ecosystem change in the Amazon region for each year of each combination of RCP and GCM. Unlike in [45], where a strong trend in the plant functional type composition change with both regional precipitation and temperature change was projected, no trend is visible here. Note that relative precipitation changes here are limited to ≤20%. The strong trend in agreement with increasing regional absolute temperature change is most likely driven by increased vegetation carbon stocks due to CO₂ fertilization effects (see figure 1 and [46]). It is also important to note that Γ is not only a measure of dieback or greening (although both Hybrid and LPJmL project a decline in Amazon rainforest trees [46]), but rather of the aggregated biogeochemical biomes changes. A similar approach to regional analysis of ecosystem shifts could be extended to develop regional ecosystem impact functions, facilitating inclusion of risk of ecosystem change in globally aggregated and sectorally integrated models.

Zoom In Zoom Out Reset image size