The impacts of modelling prescribed vs. dynamic land cover in a high-CO₂ future scenario – greening of the Arctic and Amazonian dieback

Terrestrial biosphere models are a key tool in investigating the role played by land surface in the global climate system. However, few models simulate the geographic distribution of biomes dynamically, opting instead to prescribe them using remote sensing products. While prescribing land cover still allows for the simulation of the impacts of climate change on vegetation growth and the impacts of land use change, it prevents the simulation of climate-change-driven biome shifts, with implications for the projection of future terrestrial carbon sink. Here, we isolate the impacts of prescribed vs. dynamic land cover implementations in a terrestrial biosphere model. We first introduce a new framework for evaluating dynamic land cover (i.e., the spatial distribution of plant functional types across the land surface), which can be applied across terrestrial biosphere models alongside standard benchmarking of energy, water, and carbon cycle variables in model intercomparison projects. After validating simulated land cover, we then show that the simulated terrestrial carbon sink differs significantly between simulations with dynamic vs. prescribed land cover for a high-CO₂ future scenario. This is because of important range shifts that are only simulated when dynamic land cover is implemented: tree expansion into the Arctic and Amazonian transition from forest to grassland. In particular, the projected change in net land–atmosphere CO₂ flux at the end of the 21st century is twice as large in simulations with dynamic land cover than in simulations with prescribed land cover. Our results illustrate the importance of climate-change-driven biome shifts for projecting future terrestrial carbon sink.

1 Introduction

Terrestrial biosphere models simulate the exchange of CO₂, water, and energy between the atmosphere and the land surface under both current and future climatic and socioeconomic conditions. The results from these models have led to the critical understanding that the terrestrial biosphere currently sequesters approximately a third of anthropogenic CO₂ emissions (Friedlingstein et al., 2022). Vegetation currently covers ∼ 75 % of global land area and plays a major role in regulating the land–atmosphere exchange of CO₂ (Fisher and Koven, 2020). Vegetation is a dynamic component of the global climate system that responds and feeds back to changes in environmental conditions at timescales ranging from seconds to centuries. Stomata, pores through which plants exchange CO₂ and water with the outside air, respond to changes in environmental conditions within seconds (Vialet-Chabrand et al., 2017). Over the course of a day, vegetation responds to diurnal cycles of energy. Seasonality determines annual cycles of plant phenology, and interannual climate variability regulates the land–atmosphere exchange of CO₂ over decades. Over decades to centuries, the geographical distribution of natural vegetation can also change as the geographical distribution of different biomes shifts with long-term climate variations (Pecl et al., 2017). These changes to vegetation have critical implications for the terrestrial carbon (C) sink and feedbacks to the climate system. Because changes to the geographical distribution of natural vegetation occur over such long timescales, terrestrial biosphere models generally opt to prescribe the geographical distribution of natural vegetation rather than simulate it dynamically. This avoids the challenges associated with accurately reproducing the geographical distribution of natural vegetation, which introduces another degree of uncertainty. However, when modelling the terrestrial biosphere over longer timescales, changes to the geographical distribution of natural vegetation are magnified and become critically important (Renwick and Rocca, 2015).

Terrestrial biosphere models generally represent vegetation using a set of plant functional types (PFTs) due to the intractability associated with modelling ecosystem processes for each plant species individually. This classification scheme allows terrestrial biosphere models to simplify the diversity in physiology and abiotic interactions across plant species by clustering plant species by their fundamental structure and function (Box, 1996; Reich et al., 2003). Models may choose to represent a set of PFTs that are distinguished based on leaf form (needleleaf or broadleaf), leaf phenology (evergreen or deciduous), stature (trees, grasses, or shrubs), photosynthetic pathway (C₃ or C₄), or geographical location (tropical, temperate, or boreal). When a terrestrial biosphere model prescribes land cover rather than representing it dynamically, the spatial distribution of PFTs across the land surface is specified using current observations of land cover based on remote sensing products (National Research Council, 2008). Remote sensing products are first reclassified into the PFTs that a given model represents (Hartley et al., 2017). Following this, historical information on land use change is incorporated: over the past century, crop and pasture areas have increased, while natural vegetation area has correspondingly decreased (Chini et al., 2021). The resulting land cover forcing allows for the simulation of how deforestation and the conversion of natural vegetation to crop and pasture, i.e., land use change, impact net terrestrial C sequestration (Houghton et al., 2012). Because the expansion of crop and pasture areas has been greater than climate-driven changes to the geographical distribution of natural vegetation over the historical period, prescribing land cover using a forcing that only accounts for land use change is a reasonable assumption for simulations over the historical period (although historical climate-driven changes to the geographical distribution of natural vegetation are not negligible, Pecl et al., 2017). However, prescribing land cover does not capture the shifting ranges of different biomes driven by climate change and its subsequent impacts on the net land–atmosphere CO₂ flux. Changes to the geographical distribution of natural vegetation should be especially important for high-CO₂ future scenarios with correspondingly strong climate warming because these scenarios exhibit the largest changes to vegetation productivity (Arora et al., 2020; Koven et al., 2022) and thus the largest potential changes to the ranges of different biomes.

In its simplest form, the dynamic behaviour of vegetation can be classified into two aspects: vertical and horizontal changes. Vertical changes in vegetation structure include changes to leaf area index, vegetation height, rooting depth, etc. given a certain spatial extent. Horizontal changes in vegetation structure describe changes in this spatial extent. Global change influences both aspects of the dynamic behaviour of vegetation: the response of plant growth to global change drivers such as CO₂ fertilisation and climate variability (Dusenge et al., 2019; Huang et al., 2018; Wu et al., 2011) and range shifts in different PFTs. The range of a given PFT is modulated by competitive interactions for space and resources with other PFTs and by the bioclimatic limits within which a PFT can exist. Both competitive interactions and bioclimatic limits are affected by climate change (Lenoir and Svenning, 2015; Thomas, 2010; Walther, 2010). While a terrestrial biosphere model with prescribed land cover can capture how global change modulates plant growth of a given PFT within its specified range (vertical changes), it cannot capture accompanying climate-change-driven range shifts in its PFTs (horizontal changes) aside from specified land use change. Conversely, a terrestrial biosphere model with dynamic land cover can capture both variation in plant growth (vertical changes) and range shifts (horizontal changes) driven by climate change alongside specified land use change.

The Global Carbon Project produces an annual quantification of the global C budget, which includes an estimate of the current terrestrial C sink (i.e., the positive global atmosphere–land CO₂ flux) from an ensemble of terrestrial biosphere models (Friedlingstein et al., 2022). In the most recent budget, only 3 out of the 16 contributing terrestrial biosphere models implemented dynamic land cover. Furthermore, only 3 out of the 11 Earth system models contributing to the Coupled Climate–Carbon Cycle Model Intercomparison Project (C⁴MIP) within the sixth phase of the Coupled Model Intercomparison Project (CMIP), which coordinates the analysis of C cycle interactions, included dynamic land cover (Arora et al., 2020). Here, we isolate the impacts of prescribed vs. dynamic land cover by comparing two versions of the Canadian Land Surface Scheme Including Biogeochemical Cycles (CLASSIC): a version with prescribed land cover and a version with dynamic land cover. All other aspects of the model are identical between the two versions, including land use change, thereby isolating the effects of land cover implementation (prescribed vs. dynamic). CLASSIC has been contributing to the Global Carbon Project since 2016, and CLASSIC is the land component of the Canadian Earth System Model (CanESM), which contributes to CMIP and the assessment reports of the Intergovernmental Panel on Climate Change (IPCC).

We first demonstrate that CLASSIC with dynamic land cover successfully reproduces the current geographical distribution of natural vegetation when compared to remote sensing products following a statistical framework, which we adapt for the first time to include an evaluation of the geographical distribution of natural vegetation. This framework can be applied across terrestrial biosphere models alongside standard benchmarking of energy, water, and carbon cycle variables. We then analyse CLASSIC simulations of Shared Socioeconomic Pathway 585 (SSP5-8.5), which describes a “fossil-fuelled development scenario” from 2015–2100, and examine how climate change could influence the future geographical distribution of natural vegetation. We identify differences between projections with prescribed land cover vs. dynamic land cover and examine how different land cover implementations impact future terrestrial C sink.

2 Methods

2.1 CLASSIC overview

The Canadian Land Surface Scheme Including Biogeochemical Cycles (CLASSIC) (Melton et al., 2020; Seiler et al., 2021) is the successor to and is based on the coupled Canadian Land Surface Scheme (CLASS; Verseghy, 1991; Verseghy et al., 1993) and the Canadian Terrestrial Ecosystem Model (CTEM; Arora and Boer, 2005a; Melton and Arora, 2016). Older versions of CLASSIC (under the name CLASS-CTEM) have served as the land component in the family of Canadian Earth system models (CanESMs), including CanESM5, which contributes to CMIP (Swart et al., 2019).

The physical component of CLASSIC simulates fluxes of energy, momentum, and water (Verseghy, 1991; Verseghy et al., 1993). The structural attributes of vegetation are characterised by leaf area index (LAI), canopy mass, vegetation height, and rooting depth, all of which are dynamically simulated by the biogeochemical component of CLASSIC (described below). The soil profile is represented by 20 soil layers, starting with 10 soil layers of 0.1 m thickness followed by soil layers of increasing thickness up to a soil layer of 30 m thickness for a total depth of 61.4 m. The depth of permeable soil layers and thus the depth to bedrock soil layers vary geographically and are specified based on the SoilGrids250m dataset (Hengl et al., 2017). Soil temperature and soil moisture content (liquid and frozen) are simulated for each permeable soil layer. Where the climate permits snow to exist, the temperature, mass, density, and albedo of a single snowpack layer are simulated. The physical calculations yield net radiation, soil heat flux, latent and sensible heat fluxes, evapotranspiration, and runoff at the land–atmosphere boundary. Each grid cell is simulated independently, and there are no lateral transfers of energy or matter between grid cells.

The biogeochemical component of CLASSIC simulates the land–atmosphere exchange of CO₂ via photosynthesis, autotrophic respiration, heterotrophic respiration, land use change, and fire (Arora and Boer, 2005a). CLASSIC prognostically simulates the amount of C in vegetation, litter, and soil organic matter pools for each PFT and over the bare soil fraction in each grid cell. Vegetation C is represented by leaf, stem, and root components, each of which consists of structural and non-structural carbohydrate pools. Photosynthesis generates non-structural carbohydrates that are allocated between the non-structural leaf, stem, and root C pools. Autotrophic respiration occurs from the non-structural leaf, stem, and root C pools (Arora and Boer, 2005a). Non-structural C is converted to structural C within each vegetation component (Asaadi et al., 2018). Leaf, stem, and root turnover transfer C from the vegetation C pool to the litter C pool. In addition to normal leaf turnover, leaf turnover also occurs due to drought stress, cold stress, and shorter day lengths, affecting leaf phenology (Arora and Boer, 2005a). Land use change transfers C from the vegetation C pool to land use change product C pools (with turnover times corresponding to pulp and paper products and wood products), whereas fire emits C into the atmosphere and also transfers C from the vegetation C pool to the litter C pool (Arora and Boer, 2005b; Arora and Melton, 2018). Decomposition transfers C from the litter C pool to the soil C pool. Finally, heterotrophic respiration occurs from both the litter pool and the soil C pool (Melton et al., 2015). While CLASSIC includes a representation of nitrogen cycling (Asaadi and Arora, 2021; Kou-Giesbrecht and Arora, 2022), the interactions between C and nitrogen cycling are not considered in this study.

2.2 CLASSIC land cover implementation

Biogeochemical processes in CLASSIC are modelled for nine plant functional types (PFTs): needleleaf evergreen (NE) trees, needleleaf deciduous (ND) trees, broadleaf evergreen (BE) trees, deciduous broadleaf cold (DBC) trees, deciduous broadleaf dry (DBD) trees, C₃ crops, C₄ crops, C₃ grasses (C3Gs), and C₄ grasses (C4Gs). These nine PFTs map directly to four PFTs used for simulating physical processes (needleleaf trees, broadleaf trees, crops, and grasses). When prescribed land cover is implemented, the time-varying fractional coverage of each PFT in each grid cell is specified by a land use forcing (described below). Land use change is the only driver of variation in the fractional coverage of PFTs over time. When crop area increases, natural vegetation area proportionally decreases. When crop area decreases, natural vegetation area proportionally increases. When dynamic land cover is implemented, the time-varying fractional coverages of C₃ crops and C₄ crops in each grid cell are specified by a land use forcing, and the fractional coverages of natural PFTs in each grid cell evolve due to competition and mortality within the non-crop area of each grid cell. For both prescribed and dynamic land cover implementations, the fractional coverages of all PFTs and bare ground in a grid cell sum to 1:

$$\begin{array}{}\text{(1)}& {\sum }_{n=\mathrm{1}}^{N+\mathrm{1}}{f}_n=\mathrm{1},\end{array}$$

where f_n is the fractional coverage of a given natural PFT $n\in \left\{\mathrm{1},\mathrm{\dots },N\right\}$, and f_N+1 is the fractional coverage of bare ground in a given grid cell. Additionally, a land surface fraction for each grid cell is specified to exclude waterbodies (oceans and lakes), glaciers, and ice sheets.

When dynamic land cover is implemented, the fractional coverage of a natural PFT is the result of colonisation and mortality (described in detail in Melton and Arora, 2016, and in the supplementary information). Competition between natural PFTs is based on modified Lotka–Volterra equations which describe the interactions between two populations using first-order nonlinear differential equations. In this case, competition is for area in a grid cell rather than for population size. The colonisation rate of natural PFT n (c_n; d⁻¹) is determined by its net primary productivity and leaf area index (described in detail in Appendix A). The mortality rate of natural PFT n (m_n; d⁻¹) is the sum of intrinsic or age-related mortality, mortality due to reduced growth, mortality due to fire, and mortality when a PFT exists outside its bioclimatic limits. Intrinsic or age-related mortality is determined by its maximum age (described in detail in Appendix A). Mortality due to reduced growth is determined by its growth rate over the previous year. Mortality due to fire is described in detail in Arora and Melton (2018). Area burned, which generates bare ground, depends on a PFT-specific fire spread rate (where grasses have a higher fire spread rate than trees, and needleleaf trees have a higher spread rate than broadleaf trees), wind speed, and soil moisture as well as the probability of fire occurrence (which depends on the availability of vegetation biomass as a fuel source, the combustibility of the fuel source based on soil moisture, and the presence of an ignition source based on lightning and population density) and fire suppression (which depends on population density).

Mortality when a PFT exists outside its bioclimatic limits (m_bioclim,n; d⁻¹) ensures that a given natural PFT n does not venture outside of its bioclimatic envelope. Six bioclimatic indices are used to determine the spatial range of PFTs, representing the physiological limits to their survival: air temperature of the coldest month, air temperature of the warmest month, aridity index (ratio of potential evaporation to precipitation), growing degree days (cumulative number of days with air temperature above 5 °C), dry-season length (number of consecutive months with precipitation less than potential evaporation), and precipitation surplus (difference between precipitation and potential evaporation). Each bioclimatic index for each grid cell is updated annually on a 25-year timescale using exponential smoothing:

$$\begin{array}{}\text{(2)}& X\left(t+\mathrm{1}\right)=X\left(t\right)e^{-\mathrm{1}/\mathrm{25}}+x\left(t\right)\left(\mathrm{1}-e^{-\mathrm{1}/\mathrm{25}}\right),\end{array}$$

where x(t) represents a bioclimatic index at year t, and X(t) represents the smoothed bioclimatic index at year t. This accounts for time lags in the response of vegetation to climate change drivers (Wu et al., 2015). At the beginning of each time step, each grid cell is assigned with a small fractional coverage of each PFT (0.001). Whether the PFT persists is determined by $m_{\mathrm{bioclim},n}:m_{\mathrm{bioclim},n}=\mathrm{0.25}$ when any bioclimatic index is outside its corresponding bioclimatic limit for a given PFT, and $m_{\mathrm{bioclim},n}=\mathrm{0}$ when all bioclimatic indices are inside their corresponding bioclimatic limits for a given PFT. Bioclimatic limits for each natural PFT are given in Table B1. These were derived by first calculating the bioclimatic indices for each grid cell using the CRUJRA climate dataset averaged over 1900–1920 (Harris et al., 2020). Then, the 1 % and 99 % quantiles of these bioclimatic indices over all grid cells for each land cover type in the Collection 5 MODIS Global Land Cover Type International Geosphere-Biosphere Programme (IGBP) product (Friedl et al., 2010) (time-averaged between 2001–2010) were calculated and associated with a CLASSIC PFT. Note that this assumes that the current biome ranges are in equilibrium with the 1900–1920 climate due to the migration lag that occurs between climatic change and observed differences in established plant species ranges, especially in long-lived plant species such as trees (Corlett and Westcott, 2013). While this approach accounts for slower range shifts due to climatic change, it does not account for rapid range shifts due to fire disturbance, which have influenced range shifts over recent decades (Macander et al., 2022; Wang et al., 2020).

To simulate competition between PFTs, each natural PFT is given a dominance rank (i) where PFTs with higher dominance ranks invade PFTs with lower dominance ranks and bare ground. Tree PFTs are given higher dominance ranks than grass PFTs. Within tree PFTs and within grass PFTs, a PFT with a higher c_n is dominant over all PFTs with lower c_n. Each PFT is thus ranked from $\mathrm{1},\mathrm{2},\mathrm{\dots },i-\mathrm{1},i,i+\mathrm{1},\mathrm{\dots },N$, where the PFT with dominance rank 1 is the most dominant. Overall, the change in fractional coverage of a PFT with dominance rank i is due to (1) its colonisation into the areas of PFTs with lower dominance ranks ($i+\mathrm{1},\mathrm{\dots },N)$ and bare ground (N+1), (2) colonisation by PFTs with higher dominance ranks ($i,\mathrm{\dots },i-\mathrm{1})$ into its area, and (3) its mortality:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{f}_i}{\mathrm{d}t}}& =\left(c_i{f}_{i+\mathrm{1}}+\mathrm{\dots }+c_i{f}_N+c_i{f}_{N+\mathrm{1}}\right)\\ & -\left(c_{\mathrm{1}}{f}_i+\mathrm{\dots }+c_{i-\mathrm{1}}{f}_i\right)-m_i{f}_i.\end{array}\end{array}$$

The change in fractional coverage of bare ground is due to mortality of the PFTs and colonisation of bare ground by the PFTs:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{f}_{N+\mathrm{1}}}{\mathrm{d}t}}={\sum }_{n=\mathrm{1}}^N{m}_n{f}_n-{\sum }_{n=\mathrm{1}}^N{c}_n{f}_{N+\mathrm{1}}.\end{array}$$

The conceptual form of Eqs. (3) and (4) is different from the standard Lotka–Volterra equation and allows coexistence of PFTs (Arora and Boer, 2006).

2.3 Simulation descriptions

We use CLASSIC (offline) to simulate the principal aspects of the land surface energy, water, and C cycles at the global scale over the historical period (1851–2020) and over the future period (2015–2100) for Shared Socioeconomic Pathways 585 (SSP5-8.5) at a spatial resolution of 1°. Simulations are described in Table 1.

Table 1 Description of simulations. Simulations analysed (S1–S5) are indicated for the corresponding time period. Corresponding pre-industrial spin-ups and/or historical simulations (to initialise simulations of the future period) are also indicated.

n/a: not applicable.

To evaluate dynamic land cover over the historical period (1851–2020), we conducted three simulations: two simulations with prescribed land cover (S1 and S2) and a simulation with dynamic land cover (S3). These three historical simulations were initialised from corresponding pre-industrial spin-ups (described below). Simulations with prescribed land cover use land cover forcings specifying the fractional coverage of the nine CLASSIC PFTs that are derived from two remote sensing products: the Global Land Cover 2000 (GLC2000) product (S1) (Bartholomé and Belward, 2005) and the European Space Agency Climate Change Initiative (ESACCI) version 2 product (S2) (Defourny et al., 2023). The process of generating a land cover forcing specifying the fractional coverage of all PFTs has three steps. First, the present-day fractional coverages of the model's PFTs are obtained from a remote sensing product, which involves reclassifying the land cover classes in the remote sensing product into the model's PFTs. Wang et al. (2006) describe how the 22 land cover categories in the GLC2000 product are reclassified into the nine CLASSIC PFTs. L. Wang et al. (2023) show how the 37 land cover categories in the ESACCI product from 2018 are reclassified into the nine CLASSIC PFTs. Second, the fractional coverages of crop PFTs are replaced with values from a land use change forcing for a given year in the present day, and the fractional coverages of natural PFTs are adjusted accordingly such that the total fractional coverage of vegetation is unchanged in that year. Third, the time series of fractional coverages of crop PFTs from the land use change forcing is incorporated by adjusting the fractional coverages of natural PFTs over a given time period relative to the single-year land cover forcing generated in the previous step. These three steps are used to generate a land cover forcing specifying the fractional coverage of all model PFTs derived from a remote sensing product over a given time period. Simulations with dynamic land cover use land cover forcings that only specify the fractional coverages of crop PFTs.

S1–S3 follow the TRENDY protocol (used for contributions to the Global Carbon Project) (Friedlingstein et al., 2022). The TRENDY protocol uses the merged monthly Climate Research Unit (CRU) and 6-hourly Japanese 55-year Reanalysis (JRA-55) dataset from Harris et al. (2020), atmospheric CO₂ forcing from Lan et al. (2024), population density forcing from HYDE3.3 (Klein Goldewijk et al., 2017), and land use change forcing from the Land-Use Harmonization 2 (LUH2-)GCB2022 dataset (Chini et al., 2021; Hurtt et al., 2020; Klein Goldewijk et al., 2017). S1–S3 were initialised from corresponding pre-industrial spin-ups (Table 1). We conducted different pre-industrial spin-ups for each simulation that used the corresponding land cover implementation. Pre-industrial simulations used atmospheric CO₂, land cover, and population density forcings corresponding to the year 1851. Pre-industrial simulations repeatedly used meteorological forcings from 1901–1920. We ran the pre-industrial spin-up until the C pools came into equilibrium. A threshold of 0.05 Pg C yr⁻¹ for the global net atmosphere–land CO₂ flux was used to assess if equilibrium was achieved. Note that all simulations exclude Antarctica and Greenland.

To examine prescribed vs. dynamic land cover over the future period (2015–2100) for SSP5-8.5 (Riahi et al., 2017), we conducted two simulations: one simulation with prescribed land cover from the ESACCI-derived land cover product (S4) and one simulation with dynamic land cover (S5). These future simulations for SSP5-8.5 were initialised from the end of their corresponding historical simulations, which were in turn initialised from their corresponding pre-industrial spin-ups (described below). S4 and S5 use bias-corrected meteorological forcings from CanESM5 (Buchner and Reyer, 2021). They use atmospheric CO2, land use change, and population density forcings from CMIP and the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) (Warszawski et al., 2014). Our study uses a single meteorological forcing (from the corresponding model CanESM5). The effects of different meteorological forcings have been explored elsewhere (Arora et al., 2023), as have the effects of bias correction (Seiler et al., 2024). S4 and S5 were initialised from the end of their corresponding historical simulations, which were in turn initialised from their corresponding pre-industrial spin-ups (Table 1). We conducted different pre-industrial spin-ups for each simulation. Pre-industrial simulations used atmospheric CO₂, land use change, and population density forcings corresponding to the year 1851. Pre-industrial simulations used pre-industrial meteorological forcings provided by ISIMIP.

2.4 Model evaluation

We evaluated historical simulations using the Automated Model Benchmarking R (AMBER) package developed by Seiler et al. (2021), which quantifies model performance in reproducing observation-based datasets using a skill score system that is based on the International Land Model Benchmarking (ILAMB) framework (Collier et al., 2018). Five scores compare model output to a reference dataset for a given variable, assessing the simulated time-mean bias (S_bias), monthly centralised root mean square error (S_rmse), seasonality (S_phase), interannual variability (S_iav), and spatial distribution (S_dist) in comparison to the observation-based dataset. Scores are dimensionless and range from 0 to 1, where higher values indicate better model performance. The overall score for each variable (S_overall) is calculated as

$$\begin{array}{}\text{(5)}& S_{\mathrm{overall}}=\mathrm{mean}\left(S_{\mathrm{bias}},S_{\mathrm{rmse}},S_{\mathrm{phase}},S_{\mathrm{iav}},S_{\mathrm{dist}}\right).\end{array}$$

For variables with more than one observation-based dataset, we calculated a benchmark score for the observation-based datasets, which quantifies how well the independently derived observation-based datasets agree with each other, thereby providing an estimate of the uncertainty in the observation-based datasets themselves. The benchmark score was calculated by iteratively comparing pairs of independently derived observation-based datasets, where one observation-based dataset is treated as model output, and the other observation dataset is treated as the reference dataset. If the two observation-based datasets were the same, the benchmark score would be 1. The calculation of each score and benchmark scores is described in detail in Seiler et al. (2022).

We adapted this statistical framework, which operates at a grid cell scale, to evaluate land cover fraction for each natural PFT within a grid cell. For simulations with dynamic land cover, we evaluated the land cover fraction for each natural PFT against both the GLC2000-derived land cover product and the ESACCI-derived land cover product. We compared the GLC2000-derived land cover product and the ESACCI-derived land cover product to obtain a benchmark score that quantifies the uncertainty in these land cover products. Note that a qualitative assessment of dynamic land cover at a spatial resolution of 2.81° in CLASSIC was presented in Melton and Arora (2016).

Additionally, we evaluated the standard set of energy, water, and C cycle variables following Seiler et al. (2022). Specifically, we evaluated aboveground biomass C (AGB), surface albedo (ALBS), area burned (BURNTAREA), soil C (CSOIL), vegetation C (CVEG), fire emissions (FIRE), gross primary productivity (GPP), soil heat flux (HFG), latent heat flux (HFLS), sensible heat flux (HFSS), leaf area index (LAI), runoff (MRRO), soil moisture (MRSLL), net biome productivity (NBP), net ecosystem exchange (NEE), ecosystem respiration (RECO), net surface longwave radiation (RLS), net surface radiation (RNS), net surface shortwave radiation (RSS), and snow water equivalent (SNW). Observation-based datasets used for evaluation are summarised in Table B2.

3 Results

3.1 Historical simulations and evaluation

CLASSIC with dynamic land cover successfully reproduces the global area of trees, grasses, and total natural vegetation over the historical period, falling within the range of ESACCI- and GLC2000-derived land cover products as shown in Fig. 1. Simulated global tree, grass, and natural vegetation areas are closer to the ESACCI-derived land cover product than the GLC2000-derived land cover product for the present day (2000–2020). Over the historical period, simulated global natural vegetation area decreases due to increasing crop area, i.e., land use change (Fig. B1). Land use change is the only driver of variation in natural vegetation area in the ESACCI- and GLC2000-derived land cover products. CLASSIC with dynamic land cover simulates a smaller decrease in natural vegetation area than that attributable to land use change. This is because decreasing natural vegetation area due to land use change is partially offset by increasing natural vegetation area driven by stimulated plant growth from CO₂ fertilisation and climate change over the historical period.

Figure 1 Natural vegetation area over the historical period. (a) Global area of natural vegetation, trees, and grasses simulated by CLASSIC with dynamic land cover in comparison to ESACCI- and GLC2000-derived land cover products. Natural vegetation fraction per grid cell in the (b) ESACCI-derived land cover product and (c) GLC2000-derived land cover product and (d) simulated by CLASSIC with dynamic land cover. Simulations and land cover products are averaged over 2000–2020 in (b)–(d). Natural vegetation includes seven natural PFTs (described in Sect. 2.2), excluding crop area and bare ground.

CLASSIC with dynamic land cover also broadly reproduces the latitudinal distribution of the natural PFTs in comparison to ESACCI- and GLC2000-derived land cover products as shown in Fig. 2. However, there remain some differences between simulations and observations. CLASSIC does not simulate needleleaf evergreen trees between 20 and 40° N, whereas observations suggest that this PFT exists at these latitudes (e.g., high-elevation tropical and subtropical coniferous forests). Rather, CLASSIC simulates broadleaf deciduous cold trees at these latitudes, whereas observations suggest that this PFT declines between 40 and 20° N. Furthermore, CLASSIC slightly overestimates C₄ grass area and slightly underestimates C₃ grass area in the Northern Hemisphere. Reasons for these disagreements are discussed below.

Figure 2 Latitudinal distributions of each natural plant functional type (PFT) simulated by CLASSIC with dynamic land cover in comparison to ESACCI- and GLC2000-derived land cover products. (a) Evergreen needleleaf, (b) deciduous needleleaf, (c) evergreen broadleaf, (d) deciduous broadleaf cold, (e) deciduous broadleaf dry, (f) C₃ grass, and (g) C₄ grass. Simulations and land cover products are averaged over 2000–2020.

Both CLASSIC with prescribed land cover (simulated with ESACCI- and GLC2000-derived land cover products) and CLASSIC with dynamic land cover capture the primary C, water, and energy cycle variables (Fig. B2) reasonably well. The agreement between CLASSIC simulations and observations is quantified by scores of model performance in reproducing observation-based datasets for each biogeochemical and biophysical variable. These scores indicate that CLASSIC with dynamic land cover performed similarly to CLASSIC with prescribed land cover (Fig. B3). Figure 3a and b show model scores (averaged across all CLASSIC simulations, with both prescribed and dynamic land cover) plotted against benchmark scores for biogeochemical and biophysical variables, respectively. Figure 3 shows the correspondence between model performance in reproducing observation-based datasets (i.e., model scores) and the uncertainty in the observation-based datasets themselves (i.e., benchmark scores). Model scores for both versions of CLASSIC (prescribed and dynamic land cover) and benchmark scores are similar, falling close to the 1:1 line, which indicates that model performance is comparable to the uncertainty in the observation-based datasets themselves. Both versions of CLASSIC perform slightly better for biophysical variables (Fig. 3b) than for biogeochemical variables (Fig. 3a). Finally, CLASSIC with both prescribed and dynamic land cover simulates the terrestrial C sink over the historical period reasonably well in comparison to other models in the Global Carbon Project (Friedlingstein et al., 2022) (Fig. B4).

Figure 3 Model performance in reproducing observation-based datasets (model scores) relative to the uncertainty in the observation-based datasets themselves (benchmark scores) for (a) biogeochemical variables, (b) biophysical variables, and (c) land cover fractions of each natural plant functional type (PFT). Model scores are averaged across simulations of CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing and a GLC2000-derived land cover forcing) and dynamic land cover for biogeochemical and biophysical variables. Whiskers indicate the maximum and minimum model scores. Abbreviations for C, water, and energy variables and observation-based datasets are described in the Methods section. Abbreviations for PFTs are evergreen needleleaf (EN), deciduous needleleaf (DN), evergreen broadleaf (EB), deciduous broadleaf cold (DBC), deciduous broadleaf dry (DBD), C₃ grass (C3G), and C₄ grass (C4G).

When dynamic land cover is implemented in CLASSIC, the ability of CLASSIC to reproduce land cover fractions of each natural PFT can be assessed alongside C, water, and energy variables as shown in Fig. 3c. Additionally, benchmark scores can be calculated for the land cover fractions of each natural PFT by comparing the ESACCI-derived land cover product to the GLC2000-derived land cover product. Model scores tend to be lower than benchmark scores, especially for tropical PFTs (broadleaf deciduous dry trees and C₄ grasses). For these PFTs uncertainty is introduced because remote sensing products do not distinguish between cold and dry broadleaf deciduous trees and between C₃ and C₄ grasses. This distinction is made during the reclassification of the land cover categories into the nine CLASSIC PFTs and introduces some subjectivity (L. Wang et al., 2023). Specifically, deciduous trees above 30° in both hemispheres are assumed to be cold deciduous, deciduous trees below 20° in both hemispheres are assumed to be dry deciduous, and the fraction of trees that are cold vs. dry deciduous varies linearly between 0 and 1 with latitude for deciduous trees between 20 and 30°. The separation of grasses into C₃ and C₄ is based on the time-invariant global product from Still et al. (2003).

3.2 Future simulations

For SSP5-8.5, CLASSIC with dynamic land cover simulates the increasing global area of trees, grasses, and total natural vegetation (Fig. 4a). Global natural vegetation area decreases in the ESACCI-derived land cover product due to land use change (increasing crop area), which is again the only driver of variation in natural vegetation area in the ESACCI-derived land cover product (Figs. 4a, B1). Decreasing natural vegetation area due to land use change occurs primarily in South America, Africa, and Asia (Fig. 4b). CLASSIC with dynamic land cover simulates increasing natural vegetation area because the loss of natural vegetation area due to land use change is offset by increasing natural vegetation area, especially at high latitudes (Fig. 4c).

Figure 4 Change in natural vegetation area for SSP5-8.5 (2015–2100). (a) Change in global area of natural vegetation, trees, and grasses simulated by CLASSIC with dynamic land cover in comparison to the ESACCI-derived land cover product. Change in natural vegetation fraction per grid cell (b) in the ESACCI-derived land cover product and (c) simulated by CLASSIC with dynamic land cover. Differences reflect the change between the averages calculated over 1995–2015 and 2080–2100. Only grid cells with a significant temporal trend are shown (P<0.05; assessed with a Mann–Kendall trend test). Figure B5 shows the change in fraction per grid cell for each PFT simulated by CLASSIC with dynamic land cover. Figure B13 shows the absolute value of the natural vegetation fraction per grid cell averaged over 1995–2015 and 2080–2100 simulated by CLASSIC with prescribed land cover and CLASSIC with dynamic land cover.

At high latitudes, CLASSIC with dynamic land cover simulates increasing areas of needleleaf evergreen trees, broadleaf deciduous cold trees, and C₃ and C₄ grasses, expanding into what was previously bare ground (Fig. B5). This corresponds to increasing temperatures at high latitudes (Fig. B6). In particular, C₄ grasses, which are currently primarily found at low latitudes, expand into high latitudes due to their dependence on temperature (Luo et al., 2024). At low latitudes, CLASSIC with dynamic land cover simulates decreasing natural vegetation area in South America but variably decreasing and increasing natural vegetation area in Africa and Asia. In the Amazon, the area of broadleaf evergreen trees decreases, while the areas of C₃ and C₄ grasses increase (Fig. B5). This corresponds to decreasing precipitation in this region (Fig. B6). Conversely, in Africa and Asia, the area of broadleaf evergreen trees increases, while the areas of C₃ and C₄ grasses decrease (Fig. B5). This corresponds to increasing precipitation in this region (Fig. B6). As an exception, in the Sahel, the areas of C₃ and C₄ grasses increase, expanding into what was previously bare ground (Fig. B5). This also corresponds to increasing precipitation in this region (Fig. B6).

CLASSIC with prescribed land cover (simulated with an ESACCI-derived land cover product) and CLASSIC with dynamic land cover both simulate increasing net primary productivity (NPP) at higher latitudes for SSP5-8.5 (Fig. 5a–b). NPP increases to a greater extent at higher latitudes in the simulation with dynamic land cover because of increasing natural vegetation area, especially that of broadleaf trees at the expense of needleleaf trees (Figs. 5c–f and 6a–c). NPP changes are similar at lower latitudes in both versions of CLASSIC (Fig. 5a–b). Both versions of CLASSIC project decreasing NPP in the Amazon (Fig. 5a–b). NPP decreases in simulations with prescribed land cover due to decreasing precipitation in the Amazon (Fig. B6). NPP decreases slightly less in the simulation with dynamic land cover because decreasing tree area and NPP are offset by increasing grass area and NPP in the Amazon (Figs. 5c–f and 6d–f). In the Sahel, NPP increases slightly in the simulation with dynamic land cover but decreases slightly in the simulation with prescribed land cover (Fig. 5a–b). NPP decreases in simulations with prescribed land cover due to land use change in the Sahel (Fig. 6h). NPP increases slightly in the simulation with dynamic land cover because of increasing grass area, outweighing the effects of NPP in the Sahel (Figs. 5c–f and 6g–i). Both versions of CLASSIC simulate similar increases in vegetation C (Figs. B7 and B8). CLASSIC with dynamic land cover simulates a smaller decrease in soil C due to increased natural vegetation area at high latitudes (Figs. B7 and B9).

CLASSIC dynamically simulates fire (Arora and Melton, 2018). When prescribed land cover is implemented, fire reduces vegetation biomass density, whereas when dynamic land cover is implemented, fire both reduces vegetation biomass density and creates bare ground that can then be colonised by a different plant functional type. CLASSIC with dynamic land cover simulates slightly higher area burned and fire CO₂ emissions than CLASSIC with prescribed land cover at the global scale due to higher natural vegetation area (Fig. B10), but spatial patterns are relatively similar between both versions of CLASSIC (Figs. B11 and B12). In particular, neither version of CLASSIC simulates substantial changes in high-latitude fires, which CLASSIC underestimates (Arora and Melton, 2018). CLASSIC with dynamic land cover simulates lower fire in the Amazon and higher fire in the Sahel, corresponding to lower and higher natural vegetation area in the Amazon and Sahel, respectively. This suggests that there are minor differences in fire occurrence between different PFTs, which should be improved in future model development.

Figure 5 Change in net primary productivity (NPP) simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). (a, b) Total NPP, (c, d) tree NPP, and (e, f) grass NPP. Differences reflect the change between the average over 1995–2015 and 2080–2100. Only grid cells with a significant temporal trend are shown (P<0.05; assessed with a Mann–Kendall trend test).

Figure 6 Change in natural vegetation area and net primary productivity (NPP) at high latitudes, the Amazon, and the Sahel by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). (a) Map of high latitudes. (b) Change in global area of natural vegetation, broadleaf trees, needleleaf trees, and grasses in high latitudes. (c) Change in NPP in high latitudes. (d) Map of the Amazon. (e) Change in global area of natural vegetation, broadleaf trees, needleleaf trees, and grasses in the Amazon. (f) Change in NPP in the Amazon. (g) Map of the Sahel. (h) Change in global area of natural vegetation, broadleaf trees, needleleaf trees, and grasses in the Sahel. (i) Change in NPP in the Sahel.

Net biome productivity (NBP), i.e., the net land–atmosphere CO₂ flux (which includes photosynthesis, autotrophic and heterotrophic respiration, fire emissions, and land use change emissions), increases in CLASSIC simulations of SSP5-8.5 with both prescribed and dynamic land cover (Fig. 7a). The increase in NBP is almost 2-fold greater in the simulation with dynamic land cover (5.28 Pg C yr⁻¹ averaged over 2080–2100) than in the simulation with prescribed land cover (2.29 Pg C yr⁻¹ averaged over 2080–2100). At the global scale, this increase in NBP is driven by increasing NPP (Fig. 7b–c), which is primarily due to CO₂ fertilisation in both versions of CLASSIC. In CLASSIC with dynamic land cover, there is a stronger NPP increase at high northern latitudes (Fig. 7c). NPP increases at low latitudes in both versions of CLASSIC because increasing NPP in Africa and Asia outweighs decreasing NPP in South America. The change in NPP at low latitudes is similar between both versions of CLASSIC despite changing tree vs. grass areas as described above due to offsetting effects. Surface albedo decreases due to snow loss driven by increasing temperature in both versions of CLASSIC (Fig. 7d), but this decrease is slightly stronger at high northern latitudes in CLASSIC with dynamic land cover due to the expansion of natural vegetation into higher latitudes (Fig. 7e). Evapotranspiration increases in both versions of CLASSIC (Fig. 7f), but this increase is slightly stronger at high northern latitudes in CLASSIC with dynamic land cover due to the expansion of natural vegetation into higher latitudes (Fig. 7g).

Figure 7 Change in (a) net biome productivity (NBP), (b, c) net primary productivity (NPP), (d, e) surface albedo, and (f, g) evapotranspiration simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). For the latitudinal distributions, CLASSIC simulations were averaged over 2080–2100. Thin lines indicate the annual value, and thick lines indicate the moving average over 10 years.

4 Discussion

Dynamic land cover is rarely implemented in terrestrial biosphere models (e.g., in only 3 out of 11 models contributing to C⁴MIP within the sixth phase of CMIP), and the simulated geographical distribution of natural vegetation has not yet been evaluated robustly as it lacked the methodology to do so. The Global Carbon Project uses a statistical framework to evaluate the ability of its models to reproduce observations of the C, water, and energy cycles – the International Land Model Benchmarking (ILAMB) framework (Collier et al., 2018) – which has since been expanded to include a method to quantify the uncertainty in the observation-based datasets themselves (Seiler et al., 2022). We build on this to introduce a framework with which the dynamically simulated geographical distribution of natural vegetation can be quantitatively evaluated against observations. We use this framework to show that CLASSIC with dynamic land cover successfully simulates the geographical distribution of its natural PFTs for the present day, reproducing observations (Figs. 1–3). This framework can be applied across terrestrial biosphere models to evaluate dynamic land cover alongside evaluations of the principal C, water, and energy cycle processes.

Future projections for Shared Socioeconomic Pathway 5-8.5 (fossil-fuelled development scenario, SSP5-8.5) differ significantly between simulations with prescribed and dynamic land cover because in the latter the geographical distribution of natural vegetation is able to respond to changing environmental conditions. The net land–atmosphere CO₂ flux (or net biome productivity, NBP) is nearly twice as large in the simulation with dynamic land cover (5.28 Pg C yr⁻¹ averaged over 2080–2100) than in the simulation with prescribed land cover (2.29 Pg C yr⁻¹ averaged over 2080–2100) (Fig. 7). This difference is attributed to three important range shifts that are only simulated when dynamic land cover is implemented: (1) the expansion of trees into high latitudes, which increases NBP; (2) the recession of trees accompanied by the expansion of grasses in the Amazon, which decreases NBP; and (3) the expansion of grasses into the Sahel, which increases NBP (Figs. 5, 6, and 7). These responses occurred over the historical period as well but to a much smaller extent (Fig. 1). We now discuss each of these range shifts in greater detail.

“Greening of the Arctic” describes the expansion of vegetation into high-latitude regions, which occurs as previously uninhabitable regions become habitable due to increasing temperatures (Chen et al., 2011; Elmendorf et al., 2012; Myers-Smith et al., 2020; Pearson et al., 2013; Piao et al., 2020; Tape et al., 2006). While terrestrial biosphere models with prescribed land cover capture increasing plant growth at high latitudes within their present-day range due to increasing temperatures (J. Wang et al., 2023), these models are not able to simulate range shifts. Our results suggest these range shifts could make a sizeable contribution to the terrestrial C sink via C sequestration in vegetation and soils at high latitudes due to the northward expansion of trees and grasses. Fire alongside rising temperature is known to drive a transition from needleleaf trees to broadleaf trees (Hisano et al., 2021). Because broadleaf trees are less flammable than needleleaf trees, this transition could reduce fire activity (Baltzer et al., 2021; Johnstone and Chapin, 2006; Mack et al., 2021). These processes are not represented explicitly here but warrant further study, especially given their importance and significant biases across land surface models at high latitudes (Braghiere et al., 2023; Wang et al., 2021). Although the greening of the Arctic drives increased terrestrial C sequestration at high latitudes, it also drives other important feedbacks to the climate system, such as reduced albedo (Chapin et al., 2005) and enhanced evapotranspiration (Swann et al., 2010).

“Amazonian dieback” describes the combination of deforestation and the increased mortality of trees in the Amazon rainforest due to changing climate in the region (Malhi et al., 2008; Parry et al., 2022). Tree mortality occurs as previously habitable regions become inhabitable due to the sensitivity of trees to drought and heat stress (Aleixo et al., 2019; Phillips et al., 2009). While the direct effects of reduced precipitation and elevated temperature on plant growth as well as deforestation are captured by terrestrial biosphere models with prescribed land cover, these models are incapable of simulating transitions from forest to grassland. Fire modulates these transitions, as grasses rapidly regenerate post-fire and are more flammable than trees, thereby intensifying fire activity (Davidson et al., 2012). Finally, Amazonian dieback self-amplifies via feedbacks to the climate system, as decreased forest cover decreases evapotranspiration, which could amplify regional drought (Zemp et al., 2017). The effects of these feedbacks and the feedbacks associated with the greening of the Arctic are not evaluated in our offline simulations and warrant further study.

Finally, the “greening of the Sahel” describes the expansion of vegetation into the Sahara, which occurs as previously uninhabitable regions become habitable due to increasing precipitation (Brandt et al., 2015; Herrmann et al., 2005; Olsson et al., 2005). This range shift is not captured in terrestrial biosphere models with prescribed land cover. However, there is large uncertainty in future projections of precipitation in this region (Monerie et al., 2020). While less explored than the previous two effects, greening of the Sahel is regionally important and could make a non-trivial contribution to the global C, water, and energy cycles, warranting further study.

These responses are consistent with those of other terrestrial biosphere models that implement dynamic land cover. Both Scholze et al. (2006) and Alo and Wang (2008) simulated increasing tree area at high latitudes and the replacement of trees with grasses in the Amazon. Sitch et al. (2008), who compared four models with dynamic land cover, found similar results, although the magnitudes of the responses differed between models. Port et al. (2012) conducted and examined a fully coupled simulation with a land model with dynamic land cover and similarly observed both greening of the Arctic and Amazonian dieback. Because there is large uncertainty in climate projections in the Sahel region, greening of the Sahel is often but not always simulated (Sitch et al., 2008). Our results follow the consistent pattern simulated by these other models, but our unique experimental design allows us to isolate the impact of dynamic land cover implementation (vs. prescribed land cover implementation) and compare it against the effects of other global change drivers.

5 Conclusions

Our results illustrate the potential effects of range shifts alongside variation in plant growth that is driven by climate change. While all terrestrial biosphere models account for the influence of global change drivers such as CO₂ fertilisation, increasing temperatures, and variable precipitation regimes on plant growth, terrestrial biosphere models with prescribed land cover cannot capture range shifts that are driven by climate change. However, few terrestrial biosphere models incorporate dynamic land cover. Importantly, the shifting ranges of biomes have critical feedbacks to climate change and critical consequences for biodiversity and human well-being (Pecl et al., 2017) alongside their consequences for the global C cycle.

Appendix A: Detailed description of dynamic land cover in CLASSIC

When dynamic land cover is implemented, the fractional coverage of a natural PFT is the result of colonisation and mortality (described in detail in Melton and Arora, 2016). The colonisation rate of PFT n (c_n; d⁻¹) is determined by its net primary productivity (NPP):

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{c}_n=& min\left[max\left[{\mathit{\lambda }}_{\mathrm{1},n},{\mathit{\lambda }}_{\mathrm{2},n}\right],\mathrm{0.1}\right]\\ & {\mathrm{NPP}}_n{\displaystyle \frac{\mathrm{1}}{s_{\mathrm{sap},n}max\left[\mathrm{0.25},min\left[C_{\mathrm{veg},n},\mathrm{5.0}\right]\right]}},\end{array}\end{array}$$

where NPP_n is the NPP of PFT n, and C_veg,n is the vegetation C biomass of PFT n. s_sap,n is a factor for converting vegetation C biomass to seedling C biomass of PFT n (unitless) (Table B3). λ_1,n and λ_2,n determine the fraction of NPP that is used for spatial expansion within a given grid cell of PFT n (unitless):

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1},n}=\left\{\begin{array}{l}\mathrm{0},{\mathrm{LAI}}_n\le {\mathrm{LAI}}_{\min ,n}\\ \mathrm{0.1}\frac{{\mathrm{LAI}}_n-{\mathrm{LAI}}_{\min ,n}}{{\mathrm{LAI}}_{\max ,n}-{\mathrm{LAI}}_{\min ,n}},{\mathrm{LAI}}_{\min ,n}\\ <{\mathrm{LAI}}_n<{\mathrm{LAI}}_{\max ,n}\\ \mathrm{0.1},{\mathrm{LAI}}_{\max ,n}\le {\mathrm{LAI}}_n,\end{array}\right.\text{(A3)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{2},n}=\left\{\begin{array}{l}\mathrm{0},{\mathrm{LAI}}_n\le \mathrm{0.25}{\mathrm{LAI}}_{\min ,n}\\ \cosh \left(\mathrm{0.115}\left({\mathrm{LAI}}_n-\mathrm{0.25}{\mathrm{LAI}}_{\min ,n}\right)\right)\\ -\mathrm{1},\mathrm{0.25}{\mathrm{LAI}}_{\min ,n}<{\mathrm{LAI}}_n.\end{array}\right.\end{array}$$

LAI_min,n and LAI_max,n are the minimum and maximum leaf area index thresholds of PFT n (m² m⁻²), respectively (Table B3).

The mortality rate of PFT n (m_n; d⁻¹) is the sum of intrinsic or age-related mortality (m_intr,n; d⁻¹), mortality due to reduced growth (m_ge,n; d⁻¹), mortality due to fire (m_dist,n; d⁻¹), and mortality when a PFT exists outside its bioclimatic limits (m_bioclim,n; d⁻¹):

$$\begin{array}{}\text{(A4)}& {\displaystyle }{m}_n=m_{\mathrm{intr},n}+m_{\mathrm{ge},n}+m_{\mathrm{dist},n}+m_{\mathrm{bioclim},n},\text{(A5)}& {\displaystyle }{m}_{\mathrm{intr},n}=\mathrm{1}-\exp \left({\displaystyle \frac{-\mathrm{4.605}}{A_{\max ,n}}}\right),\text{(A6)}& {\displaystyle }{m}_{\mathrm{ge},n}={\displaystyle \frac{m_{\mathrm{ge},\max ,n}}{\mathrm{1}+\left(\mathrm{300}{\mathrm{m}}^{\mathrm{2}}\mathrm{kg}{\mathrm{C}}^{-\mathrm{1}}\right)\left(max\left[\mathrm{0},{\mathrm{ge}}_n\right]\right)}},\text{(A7)}& {\displaystyle }{m}_{\mathrm{dist},n}={\mathit{\varsigma }}_{\mathrm{r},n}{\displaystyle \frac{f_{\mathrm{burned},n}}{f_n}},\end{array}$$

where A_max,n is the maximum age of PFT n (year) (Table B3), and $m_{\mathrm{ge},\max ,n}$ is the maximum growth-related mortality rate occurring when no growth occurs in PFT n (d⁻¹) (Table B3). ge_n is the growth efficiency over the previous year of PFT n (kg C m⁻²):

$$\begin{array}{}\text{(A8)}& {\mathrm{ge}}_n={\displaystyle \frac{\mathrm{\Delta }{C}_{\mathrm{s},n}+\mathrm{\Delta }{C}_{\mathrm{r},n}}{L_{\max ,n}}},\end{array}$$

where ΔC_s,n is the stem C increment over the previous year of PFT n (kg C m⁻²), ΔC_r,n is the root C increment over the previous year of PFT n (kg C m⁻²), and L_max,n is the maximum LAI over the previous year of PFT n (m² m⁻²). ς_r,n is the susceptibility to stand-replacing fire of PFT n (fraction) (Table B3). f_burned,n is the fractional area burned within a day of PFT n and depends on wind speed, soil moisture, and the probability of fire occurrence (which depends on the availability of vegetation biomass as a fuel source, the combustibility of the fuel source based on soil moisture, and the presence of an ignition source based on lightning and population density) and suppression (which depends on population density). This is described in detail in Arora and Melton (2018).

Six bioclimatic indices are used to determine the spatial range of PFTs, representing the physiological limits to their survival: air temperature of the coldest month, air temperature of the warmest month, aridity index (ratio of potential evaporation to precipitation), growing degree days (cumulative number of days with air temperature above 5 °C), dry-season length (number of consecutive months with precipitation less than potential evaporation), and precipitation surplus (difference between precipitation and potential evaporation). Bioclimatic indices for each grid cell are updated annually on a 25-year timescale using exponential smoothing:

$$\begin{array}{}\text{(A9)}& X\left(t+\mathrm{1}\right)=X\left(t\right)e^{-\mathrm{1}/\mathrm{25}}+x\left(t\right)\left(\mathrm{1}-e^{-\mathrm{1}/\mathrm{25}}\right),\end{array}$$

where x(t) represents a bioclimatic index at year t, and X(t) represents the smoothed bioclimatic index at year t. This accounts for time lags in the response of vegetation to climate change drivers (Wu et al., 2015). At the beginning of each time step, each grid cell is assigned with a small fractional coverage of each PFT (0.001). Whether the PFT persists is determined by $m_{\mathrm{bioclim},n}:m_{\mathrm{bioclim},n}=\mathrm{0.25}$ when any bioclimatic index is outside its corresponding bioclimatic limit for a given PFT, and $m_{\mathrm{bioclim},n}=\mathrm{0}$ when all bioclimatic indices are inside their corresponding bioclimatic limits for a given PFT. Bioclimatic limits for each PFT are given in Table B1. These were derived by first calculating the bioclimatic indices for each grid cell using the CRUJRA climate dataset averaged over 1900–1920 (Harris et al., 2020). Then, the 1 % and 99 % quantiles of these bioclimatic indices over all grid cells for each land cover type in the Collection 5 MODIS Global Land Cover Type International Geosphere-Biosphere Programme (IGBP) product (Friedl et al., 2010) were calculated and associated with a CLASSIC PFT. Note that this assumes that the current biome ranges are in equilibrium with the 1900–1920 climate.

Each natural PFT is given a dominance rank (i), where PFTs with higher dominance ranks invade PFTs with lower dominance ranks and bare ground. Natural PFTs cannot invade crop area. Tree PFTs are given higher dominance ranks than grass PFTs. Within tree PFTs and within grass PFTs, a PFT with a higher c_n is dominant over all PFTs with lower c_n. Each PFT is thus ranked from $\mathrm{1},\mathrm{2},\mathrm{\dots },i-\mathrm{1},i,i+\mathrm{1},\mathrm{\dots },N$, where the PFT with dominance rank 1 is the most dominant. Overall, the change in fractional coverage of a PFT with dominance rank i is due to (1) its colonisation into the areas of PFTs with lower dominance ranks ($i+\mathrm{1},\mathrm{\dots },N)$ and bare ground (N+1), (2) colonisation by PFTs with higher dominance ranks ($i,\mathrm{\dots },i-\mathrm{1}$) into its area, and (3) its mortality:

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{f}_i}{\mathrm{d}t}}& =\left(c_i{f}_{i+\mathrm{1}}+\mathrm{\dots }+c_i{f}_N+c_i{f}_{N+\mathrm{1}}\right)\\ & -\left(c_{\mathrm{1}}{f}_i+\mathrm{\dots }+c_{i-\mathrm{1}}{f}_i\right)-m_i{f}_i.\end{array}\end{array}$$

The change in fractional coverage of bare ground is due to mortality of the PFTs and colonisation of bare ground by the PFTs:

$$\begin{array}{}\text{(A11)}& {\displaystyle \frac{\mathrm{d}{f}_{N+\mathrm{1}}}{\mathrm{d}t}}={\sum }_{n=\mathrm{1}}^N{m}_n{f}_n-{\sum }_{n=\mathrm{1}}^N{c}_n{f}_{N+\mathrm{1}}.\end{array}$$

Appendix B: Supplemental figures and tables

Figure B1 Crop area over the historical period (1850–2015) and SSP5-8.5 (2015–2100) from the Land-Use Harmonization 2 (LUH2-)GCB2022 dataset (Chini et al., 2021; Hurtt et al., 2020; Klein Goldewijk et al., 2017).

Figure B2 Latitudinal distributions of (a) gross primary productivity (GPP), (b) aboveground biomass C, (c) soil C, (d) area burned, (e) fire CO₂ emissions, (f) soil moisture, (g) surface net latent heat flux, (h) surface net sensible heat flux, (i) net surface radiation, and (j) surface albedo simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing and a GLC2000-derived land cover forcing) and CLASSIC with dynamic land cover in comparison to present-day observations. CLASSIC simulations were averaged over 2000–2020. Observation-based datasets are described in the Methods section. Figure B2 shows the scores of model performance in reproducing observation-based datasets for each variable.

Figure B3 Scores of model performance in reproducing observation-based datasets of principal C, water, and energy variables as well as land cover fractions of each natural plant functional type (PFT) for simulations of CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing and a GLC2000-derived land cover forcing) and dynamic land cover. Abbreviations for C, water, and energy variables and observation-based datasets are described in the Methods section. Abbreviations for PFTs are evergreen needleleaf (EN), deciduous needleleaf (DN), evergreen broadleaf (EB), deciduous broadleaf cold (DBC), deciduous broadleaf dry (DBD), C₃ grass (C3G), and C₄ grass (C4G).

Figure B4 Net biome productivity (NBP) simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing and a GLC2000-derived land cover forcing) and CLASSIC with dynamic land cover. Boxes indicate estimates from the suite of models in the Global Carbon Project, where dark horizontal lines indicate the mean, and the extent indicates 1 standard deviation (Friedlingstein et al., 2022). Thin lines indicate the annual value, and thick lines indicate the moving average over 10 years.

Figure B5 Fraction per grid cell in the historical period (averaged over 1995–2015), in SSP5-8.5 (averaged over 2080–2100), and change between the historical period and SSP5-8.5 for each plant functional type (PFT) simulated by CLASSIC with dynamic land cover. (a) Evergreen needleleaf, (b) deciduous needleleaf, (c) evergreen broadleaf, (d) deciduous broadleaf cold, (e) deciduous broadleaf dry, (f) C₃ grass, and (g) C₄ grass. Differences reflect the change between the average over 1995–2015 and 2080–2100.

Figure B6 Temperature and precipitation in (a, c) the historical period (averaged over 1995–2015) and (b, d) SSP5-8.5 (averaged over 2080–2100) and (e, f) temperature and precipitation difference between the historical period and SSP5-8.5.

Figure B7 (a, b) Change in vegetation C and (c, d) soil C simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). For the latitudinal distributions, CLASSIC simulations were averaged over 2080–2100. Thin lines indicate the annual value, and thick lines indicate the moving average over 10 years.

Figure B8 Change in vegetation C simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). (a, b) Total vegetation C, (c, d) tree vegetation C, and (e, f) grass vegetation C. Differences reflect the change between the average over 1995–2015 and 2080–2100.

Figure B9 Change in soil C simulated by (a) CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and (b) CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). Differences reflect the change between the average over 1995–2015 and 2080–2100.

Figure B10 Change in (a, b) area burned and (c, d) fire CO₂ emissions by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). For the latitudinal distributions, CLASSIC simulations were averaged over 2080–2100. Thin lines indicate the annual value, and thick lines indicate the moving average over 10 years.

Figure B11 Change in area burned simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). (a, b) Total area burned, (c, d) tree area burned, and (e, f) grass area burned. Differences reflect the change between the average over 1995–2015 and 2080–2100.

Figure B12 Change in fire CO₂ emissions simulated by CLASSIC with prescribed land cover (with an ESACCI-derived land cover forcing) and CLASSIC with dynamic land cover for SSP5-8.5 (2015–2100). (a, b) Total fire CO₂ emissions, (c, d) tree fire CO₂ emissions, and (e, f) grass fire CO₂ emissions. Differences reflect the change between the average over 1995–2015 and 2080–2100.

Figure B13 Natural vegetation fraction per grid cell (a, b) in the ESACCI-derived land cover product and (c, d) simulated by CLASSIC with dynamic land cover over (a, c) the historical period (averaged over 1995–2015) and (b, d) for SSP5-8.5 (averaged over 2080–2100).

Table B1 Parameters associated with bioclimatic indices.

Table B2 Observation-based datasets used in AMBER. Described in detail in Seiler et al. (2022). Datasets are globally gridded or in situ (specified).

Table B3 Parameters associated with the calculation of colonisation and mortality rates for each PFT.