Trends and Drivers of Terrestrial Sources and Sinks of Carbon Dioxide: An Overview of the TRENDY Project

The terrestrial biosphere plays a major role in the global carbon cycle, and there is a recognized need for regularly updated estimates of land‐atmosphere exchange at regional and global scales. An international ensemble of Dynamic Global Vegetation Models (


Introduction
The land biosphere plays an important role in the global carbon cycle, as a natural carbon sink for anthropogenic CO 2 and for emissions from land-use and land cover changes (LULCC).The magnitude of the natural land carbon sink (S LAND ) has increased along with atmospheric CO 2 and is equivalent to almost one third of total anthropogenic CO 2 emissions, thus providing an important ecosystem service by mitigating climate change (Friedlingstein et al., 2022a).Currently, LULCC emissions (E LUC ) contribute to about 11% of total anthropogenic CO 2 emissions (Friedlingstein et al., 2022a) at global scale with important regional variability.Realistic estimates of future land carbon sinks and sources are vital to quantify the remaining carbon budget, that is, the CO 2 that society can still emit to stay within the targets of 1.5°C or well below 2°C considered by the UN Paris Agreement.This information is critical to understand the efficiency of land ecosystems to sequester carbon in the future (Arora et al., 2020).While the global magnitude of the net land-atmosphere carbon exchange is well known and observable (Ruehr et al., 2023), key challenges remain in identifying its geographical location, ecosystems and processes (e.g., split between S LAND and E LUC ) responsible for variability and trends, and the timescales at which they act (Pongratz et al., 2021).Hence the need for process-based models to inform on attribution.
Different families of land models are used in carbon cycle assessments, including here in the TRENDY project.Contemporary Dynamic Global Vegetation Models (DGVMs) have developed over the past 30 years, and have their roots in the disciplines of biogeography, biogeochemistry, forest gap dynamics, and land surface modeling (LSMs) (Fisher et al., 2014;Prentice et al., 2007).We use the broad term DGVM as it encapsulates dynamic modeling of vegetation and associated biogeochemical cycling at global scales.Several of the TRENDY DGVMs are also LSMs.LSMs are important tools to study the physical and biogeochemical land surface processes at subdaily timescales and can be coupled into Earth System Models.During the 1990s, LSMs were compared against each other for the first time in the Project for Intercomparison of Land Surface Schemes (PILPS) (Henderson-Sellers et al., 1993) in terms of their ability to simulate physical processes.The Coupled Model Intercomparison Project (CMIP) also started at around this time, focusing on coupled ocean-atmosphere general circulation models, and later the Coupled Climate-Carbon Cycle MIP (C4MIP) (Friedlingstein et al., 2006).The first biogeochemistry-focused MIPs were designed to synthesize results on land carbon dynamics, provide consensus, identify uncertainties, and quantify the role of individual processes to the land carbon cycle.These included the Net Primary Production MIP (Cramer et al., 2001) and two concurrent projects: the Vegetation-Ecosystem Modeling and Analysis Project (VEMAP) (Melillo et al., 1995), and the Carbon Cycle Model Linkage Project (CCMLP) (Heimann et al., 1998;McGuire et al., 2001).These early MIPs laid the foundation for subsequent DGVM intercomparisons, based on the idea of provision of a common simulation protocol, timelines, and driving data sets.McGuire et al. (2001) pioneered the factorial simulation approach whereby models were driven with a different combination of time-varying and fixed environmental forcing drivers.This approach aims to attribute change in any model simulated quantity to the underlying forcing drivers and their interactions.In addition, identification of sources of model uncertainty allows to quantify knowledge and data gaps, and thus guide future model improvements as well as inform empirical, monitoring and observational needs.
The TRENDY Project (Sitch et al., 2015), uses the same factorial approach.TRENDY is an international MIP which contributes regional and global Net Biome Production (NBP) estimates to the annual Global Carbon Budget (GCB) assessments (e.g., Friedlingstein et al., 2022a), with data available to the wider scientific community for numerous additional analyses on land biogeochemistry.TRENDY models contributed to the first GCB assessment (Le Quéré et al., 2009) and the Regional Carbon Cycle Assessment and Processes project, phase 1 (RECCAP-1) (Sitch et al., 2015).These early studies focused on diagnosing the natural land carbon sink, and in more recent years the protocol has expanded to include LULCC to estimate associated E LUC .As a community TRENDY-v1 started in 2009, and has been repeated every year since 2013 (TRENDY-v2), through to TRENDY-v11 in 2022 using updated and extended forcing data sets.Throughout this period the protocol has been developed and refined.
The aim of this paper is to document the recent protocol as used for TRENDY v11, including an updated description of the standard set of factorial simulations, underlying environmental forcing data sets up to the year 2021, and to present results for the contemporary period beyond those presented in GCB2022 (Friedlingstein et al., 2022a), and over 10 geographical regions used in RECCAP2.

Data Sets
A description of the forcing data sets is given below and summarized in Table S1 of Supporting Information S1.This includes a general description valid for TRENDY versions since Sitch et al. (2015) and where necessary specific additional information for TRENDY-v11/GCB2022.
Our general modeling approach is to use the most up-to-date historical data sets available and pre-processed by May/June each year to be able to run the DGVMs in July/August.Our choice of climate forcing is based on our need to run DGVMs with high-temporal fidelity (diel cycle) observation-based climate over long, century timescales with low latency.We need gridded forcing data covering the period 1900 until the end of the previous year, that is, for TRENDY-v11 until the end of 2021.In each new annual cycle of TRENDY these data are updated by the addition of 1 year, and retrospectively, as new historical meteorological station data become available and/or blending and interpolation methodologies change.Given the multi-decade to century timescales of carbon legacy fluxes associated with land use change (e.g., soil carbon stock change with deforestation), our LULCC forcing also needs to be global, spatially explicit (transformed onto a common grid with the climate forcing) and cover multiple centuries, that is, prior to the satellite era, for example, from 1700 onwards.Some underlying data sets are either not available every year, for example, data from the global Forest Resources Assessments (FRA), and/or are not updated in time for the GCB annual assessments, for example, the Food and Agricultural Organization of the United Nations (FAO) land-use statistics.Thus, we need to make extrapolations in time, and in some GCB years, for simplicity, FAO data were not updated at all.

Global Atmospheric CO 2
A 1700-2021 monthly time-series, derived by combining ice core CO 2 data with monthly observations of CO 2 concentration from NOAA stations since 1958, is prepared for the annual GCB assessments.Data from March 1958 onwards are based on monthly averages from Mauna Loa (MLO) and South Pole (SPO) atmospheric monitoring stations, provided by NOAA's Earth System Research Laboratory http://www.esrl.noaa.gov/gmd/ccgg/trends/.When no SPO data are available (including for the entire period 1959-1975 and some months after 1975), the SPO time series is gap-filled based on the trend in divergence of SPO from MLO across years and the climatological mean monthly difference between SPO from MLO during the period 1976-2014.Data prior to March 1958 are estimated with a cubic spline fit to CO 2 concentrations from air trapped in ice cores (Joos & Spahni, 2008).Annual mean fields are generated from these monthly data to force the DGVMs.

Climate Forcing
The Climatic Research Unit (CRU) monthly historical climate data at 0.5°spatial resolution are based on a synthesis of meteorological station data ("CRU TS") and are available for the period 1901-2021, or to the previous year of the simulations for more recent releases.These data are selectively merged with sub-daily climate data from a new generation reanalysis product, JRA-55 (Japanese 55-year Reanalysis), which are available for the period 1958-2021.The blended product yields a composite 0.5 degree CRUJRA 6-hourly historical climate forcing data set from 1901 to 2021 for the following climate variables: minimum, maximum, and mean air temperature (t min , t max , tmp); precipitation (pre), specific humidity (spfh), pressure (pres), wind speed (u and v-components, ugrd, vgrd) and downward longwave radiation (dlwrf) (Figure 1).The model resolution of JRA is T319 Gaussian grid (approximately 0.56°), which is close to 0.5°spatial resolution of the CRU data set.
The procedure for blending JRA with CRU data is described as follows.All JRA-55 data are regridded to the CRU 0.5°grid using NCL routines (based on ESMF routines) and masked to give a land-only (excluding Antarctica) data set.The time series for each variable is obtained by scaling the 6-hourly data simulated in the Japanese Reanalysis (JRA; Kobayashi et al., 2015) so that their monthly means or totals match those in the CRU TS v4.03 (and later versions) data set (Harris et al., 2020).This procedure, termed "alignment" hereafter, is applied to T mean , T max , T min , Spfh and Pre: T min and T max are aligned with CRU TS Mean Temperature (TMP) and Diurnal Temperature Range; T mean is aligned with TMP; Spfh is aligned with CRU TS Vapor Pressure (VAP, converted to Specific Humidity); Pre is aligned with CRU TS Precipitation (PRE).The other five variables (pres, ugrd, vgrd, dlwrf, dswrf) in the JRA-55 data are not aligned with CRU TS.For the time period 1958-2021 for which JRA-55 data are available they are blended with the CRU TS data using the procedure described above.For the time period 1901-1957, for which the CRU TS data are available but the JRA-55 data are not, randomly chosen (but now fixed -to allow repeatability) years from JRA-55 are chosen from the period 1958-1967 and assigned to the years 1901-1957, thus providing high-frequency variability in the earlier period.Each month in that period is then aligned with the appropriate CRU TS monthly value from 1901 to 1957, as above.

Revised Radiation Fields
Vegetation productivity is sensitive to both the quantity and quality of light.The latter relates to the scattering of incoming shortwave radiation through the atmosphere during cloudy skies and/or through aerosols in the atmosphere.Scattered light, that is, diffuse light conditions, provides a more uniform distribution of light through the vegetation canopy and can thus enhance photosynthesis in shaded leaves in the sub-canopy layers which are typically light-limited (Roderick et al., 2001).The impact of diffuse radiation on the simulated global land carbon sink was first evaluated in the framework of the JULES land model (Mercado et al., 2009).Mercado et al. (2009) found that the use of diffuse radiation for photosynthesis enhances the global land carbon sink by approximately one-quarter between 1960 and 1999.As more DGVMs include the ability to use both direct and diffuse radiation within their photosynthesis modules, a time series of diffuse radiation has been provided as part of the TRENDY protocol since TRENDY-v10/GCB2021 (Friedlingstein et al., 2022b).In TRENDY-v11, CABLE-POP, CLM5.0,JULES-ES, and YIBs used the diffuse fraction information.A separate study using two DGVMs and this diffuse radiation time series was also published (O'Sullivan et al., 2021).Here we summarize the diffuse forcing data set with more detailed information in O' Sullivan et al. (2021).
Six hourly fields of the diffuse fraction of surface shortwave fluxes were generated over the period 1901-2021 (Figure S1 in Supporting Information S1).A detailed evaluation of the data set against in situ radiation data is performed in O' Sullivan et al. (2021).Radiative transfer calculations are based on monthly averaged distributions of tropospheric and stratospheric aerosol optical depth, and 6-hourly distributions of cloud fraction.The procedure followed is described in the Methods section of Mercado et al. (2009), but with updated input data sets.Surface radiative fluxes account for aerosol-radiation interactions from both tropospheric and stratospheric aerosols, and aerosol-cloud interactions from tropospheric aerosols, except mineral dust.Tropospheric aerosols are also assumed to exert interactions with clouds.The radiative effects of those aerosol-cloud interactions are assumed to scale with the radiative effects of aerosol-radiation interactions of tropospheric aerosols, using regional scaling factors derived from HadGEM2-ES.
The time series of speciated tropospheric aerosol optical depth is taken from the historical and RCP8.5 simulations by the HadGEM2-ES climate model (Bellouin et al., 2011).To correct for biases in HadGEM2-ES, tropospheric aerosol optical depths are scaled over the whole period to match the global and monthly averages obtained over the period 2003-2020 by the CAMS Reanalysis of atmospheric composition (Inness et al., 2019), which assimilates satellite retrievals of aerosol optical depth.
The time series of stratospheric aerosol optical depth is taken from the climatology by Sato et al. (1993), which has been updated to 2012.The stratospheric aerosol optical depth for years 2013-2021 was assumed to be the same as for a normal year (2010), which saw background stratospheric aerosol concentration, with no large volcanic eruption or lingering effects of prior volcanic eruptions.This assumption is supported until at least 2016 included by the Global Space-based Stratospheric Aerosol Climatology time series (1979-2016;Thomason et al., 2018).Solar radiation is assumed to be 100% diffuse in cloudy skies.Atmospheric constituents other than aerosols and clouds are set to a constant standard mid-latitude summer atmosphere, as their variations do not affect the diffuse fraction of surface shortwave fluxes.Each grid is partitioned into cloudy and cloud free fractions at 6-hourly intervals using JRA-55 data, scaled to match CRU monthly observations.

Land Use Change
The earliest approaches to simulate the effect of land use change in DGVMs (McGuire et al., 2001) used gridded annual data sets of cropland coverage and a rule-based approach, adopted from book-keeping models (Houghton et al., 1983), to simulate the fate of deforested carbon over land where natural vegetation is converted to cropland.Newer data sets now include, in addition to changes in cropland area, changes in pasture and grazing area and detailed sub-grid transitional information between vegetation cover classes which allows to represent practices such as shifting cultivation.Models have also improved their representation of land-use change processes, including incorporating wood harvest, shifting cultivation, grazing, crop harvesting, and cropland management, although the individual capabilities vary between models (Arneth et al., 2017;Pongratz et al., 2017).The use of agricultural area change (including cropland and pasture) data sets is an indirect way to infer change in area of natural ecosystems, for example, forest loss and gain due to deforestation and secondary forest regrowth.The advantage of using information on changes in agricultural area and wood harvest as compared to using directly information on area and state of forest and other natural vegetation types is that anthropogenic and natural drivers are cleanly separated.Time-varying, spatially explicit data sets that are used to represent land use change in models are based on annual land use and population data from History Database of the Global Environment v3.3 (HYDE 3.3), land use states and transitions produced by updates to the Land Use Harmonization v2 model (LUH2-GCB; Chini et al., 2021), and Food and Agriculture Organization (FAO) wood harvest data.LUH2 has been developed for the purpose of driving land models as part of the CMIP endeavor and has been extended in scope to match progress in LULCC description in the DGVMs.

HYDE3.3
HYDE combines historical population estimates and model algorithms to generate spatially explicit, temporally varying maps of land use categories (Klein Goldewijk, Beusen, et al., 2017), for the period 10,000 Before the Common Era (BCE) to 2017 Common Era (CE) period.The land use maps are based on an allocation algorithm that uses country totals from the Food and Agricultural Organization of the United Nations statistical data (FAOSTAT) of "Arable land and permanent crops" and "Permanent meadows and pastures" (or "grazing land" in HYDE) for the 1961 to present period (FAO, 2020a), and other historical census data supplemented with percapita land use estimates for the pre-1960 period.HYDE uses the European Space Agency Climate Change Initiative (ESA-CCI) (Plummer et al., 2017) medium-resolution land cover data set (ESA, 2016;Harper et al., 2023), to spatially allocate the FAO land use areas for the present-time , which are gradually replaced by a combination of normalized proxies (climate, soil fertility, distance to water, slope, population density) to allocate land use area for the past.
To further inform land cover change in DGVMs, the grazing land use category is subdivided based on the intensity of use; grazing in rangelands (extensive grazing on the drier natural ecosystems), managed pastures (intensive grazing or mowing, on any natural vegetation type) and converted rangeland (forests converted to rangelands in areas with low human population density).ESA CCI land cover data are used to constrain the spatial allocation of cropland and grazing land on an annual basis for the period 1992-2017 (Klein Goldewijk, Dekker, et al., 2017).FAO agricultural statistics are used until 2017, and for TRENDY-v11 we extrapolated these data using the trend from 2012 to 2017 forward in time (see LUH2-GCB section).
However, using FAO statistics and its extrapolation may not be appropriate for capturing trends and variability in countries that have changed reporting methodologies without retrospective data revision or experienced recent rapid land-use change.For example, China (Yu et al., 2022) and Brazil (Rosan et al., 2021), have witnessed a recent increase in reforestation/afforestation and deforestation, respectively.This is a major source of uncertainty in GCB assessments (Bastos et al., 2020;Rosan et al., 2021).For Brazil, in TRENDY-v11 we have replaced the FAO state-level data for Brazil's cropland and grazing land with those derived from the in-country land cover data set MapBiomas (collection 6) for the 1985-2020 period (Souza et al., 2020).Note, ESA-CCI spatially explicit land cover maps are still used to spatially allocate the total country area following HYDE methodology.For Brazil, an estimate of land use change for the year 2021 is based on the trend 2015-2020 from the total areas calculated from MapBiomas.The pre-1985 period is scaled with the per capita numbers from 1985, so ensuring a smooth transition.The land use areas in the HYDE framework are computed at a 5-arc min spatial resolution and aggregated to a 0.25-degree resolution for use in LUH2.

LUH2-GCB
Gridded HYDE data on agricultural area (cropland and grazing area), population data until 2017 (2020 for Democratic Republic of the Congo, and 2022 for Brazil) and national wood harvest production data from FAO (2020b) are then used as input to the process that generates the updated Land-Use Harmonization 2 (LUH2-GCB) data set (Chini et al., 2021;Hurtt et al., 2019aHurtt et al., , 2019bHurtt et al., , 2020)).LUH2 harmonizes multiple datastreams into a single 0.25-degree spatial resolution data set of land-use transitions for 12 land-use states, the age and area of recovering secondary lands, wood harvest area, and agricultural management practices such as irrigated area and synthetic nitrogen fertilizer application (Hurtt et al., 2020).The 12 LUH2 land-use states are cropland (C3 annuals, C4 annuals, C3 perennials, C4 perennials, and C3 N-fixers), grazing land (managed pasture and rangeland), urban land, primary land (both forested and non-forested), and secondary land (forested and non-forested).
In TRENDY-v11, for all countries except Brazil and the Democratic Republic of the Congo, after the year 2017, LUH2-GCB linearly extrapolates on a per gridcell basis the cropland, pasture, and urban data, based on the trend over the previous 5 years (2012)(2013)(2014)(2015)(2016)(2017), to generate data for years 2018-2022.For gridcells within the Democratic Republic of the Congo this is needed only from 2020 until the year 2022.For gridcells within Brazil, LUH2-GCB uses the annual HYDE data as it is provided for the full length of the simulation.Land-use states are available for the period 850-2022 (Figure 2a), and land-use transitions for the period 850-2021 (Chini et al., 2021).Extrapolation using the trend over the previous 5 years is a pragmatic approach to enable necessary extension of the LULCC data sets.After successful adoption of mapbiomas for Brazil, we envisage the adoption of low latency, time-varying Earth Observation based maps of LULCC moving forward.
For the wood harvest inputs to LUH2-GCB, we use annual FAO wood harvest data from 1961 to 2019, and linearly extrapolate to the year 2022 because FAO updates were not available (Figure 2b).Wood harvest patterns are constrained with Landsat-based tree cover loss data (Hansen et al., 2013).HYDE population data are used to extend this time series back in time prior to 1961, along with data from Zon and Sparhawk (1923) which is used to specify national wood harvest demands in the year 1920.Some DGVMs use in addition to land-use transition data, directly information on wood harvest.LUH2-GCB data were provided in three separate files (for the states, transitions, and management data layers).A full set of LUH2-GCB diagnostics are given in Table S1 of Supporting Information S1. Figure 3 shows improved diagnostics in the LUH2-GCB2022 data set when compared with reference data for two important regions/ variables: cropland in the Democratic Republic of the Congo, and for recent primary forest loss over the satellite era in Amazonian Brazil.

Application of LUH2 in DGVMs: Conversion to Pasture/Rangeland
The LUH2 methodology uses the cropland, managed pasture, rangeland, and urban layers from HYDE.Natural vegetation is cleared for croplands.DGVM groups in the past have requested more information on whether natural vegetation is lost in conversion to managed pasture or rangeland, since only the former is typically associated with a change in land cover from natural vegetation to managed grassland.The LUH2 simple guidelines provided at https://luh.umd.edu/data.shtml(Ma et al., 2020) state: "all natural vegetation should be cleared for managed pasture, and only cleared for rangeland if it is forested."Using this rule/guideline, together with the forest layer from LUH2, gives maps of forest area, carbon density, and carbon emissions.Models can have static or dynamic natural vegetation but all use prescribed geographical distribution of cropland and grazing land (defined as the sum of managed pasture and rangeland).These LULCC data sets are common to all the DGVMs, but their implementation differs depending on the land-use processes and functionalities in individual models (Arneth et al., 2017; Table 1).LUH2-GCB provides fractional gridcell information on land-use states and their transitions that enables DGVM simulation of sub-grid transitions and gross LUC carbon fluxes, for example, associated with deforestation and forest regrowth from shifting cultivation within the same gridcell, and permanent deforestation and forest (re-)growth due to afforestation/reforestation.However, some DGVMs are only able to simulate net E LUC , and thus rely on the subset of net agricultural area changes based on HYDE, whereas others can simulate gross E LUC , that is, including sub-grid cell transitions/shifting cultivation, which require in addition the land-use transition information from LUH2 (Table 1).

Nitrogen Deposition and Fertilizer
A major effort has been made to also represent nutrient cycles in DGVMs, first with the inclusion of a Nitrogen (N) cycle (Thornton et al., 2007;Zaehle & Friend, 2010) and more recently a phosphorous cycle (Goll et al., 2017;Y. P. Wang et al., 2007, 2010, X. Yang et al., 2014).Broadly, high latitude ecosystems with lower nitrogen mineralization rates in cold environments may be nitrogen limited, whereas highly weathered tropical soils are thought to be primarily P limited.In a modeling study Zaehle et al. (2010), simulated a reduction in the global CO 2 fertilization effect in land ecosystems counteracted by a reduced positive climate-carbon cycle feedback with inclusion of a nitrogen cycle.N deposition fields can also alleviate potential N limitation.Carbon-nitrogen cycle DGVMs (Table 1) therefore need historical gridded forcing data on nitrogen deposition and fertilizer application.
Nitrogen fertilizer input data sets are available via the NMIP and NMIP2 projects (Tian et al., 2018(Tian et al., , 2019(Tian et al., , 2022)).Manure is an organic fertilizer (animal waste put on fields), and is important from the N-cycle perspective, because it is one of the natural fertilizer sources.However, since manure is organic N, it must be taken from other organic N sources in a model which causes problems with mass balance without representation of other relevant processes.Doing this incorrectly will influence the simulated C and N cycle.For TRENDY, DGVMs are thus recommended not to include manure application over croplands.Artificial inorganic N fertilizer application started after 1910 (the Haber-Bosch process was developed in the early 20th century), and thus DGVMs assume zero N fertilizer for years 1700-1910.Recently gridded N inorganic fertilizer application rates have been updated to 2020 (Tian et al., 2022).For TRENDY-v11 DGVMs assume N fertilizer application rates remain unchanged in years 2020-2021 (Tian et al., 2022).
For nitrogen deposition (Ndep), DGVMs use historical N-deposition rates for the period 1850-2014 and then transition to the future RCP8.5 N-deposition scenario for the period 2015-2021 (Hegglin et al., 2016).N deposition is available only from 1850, and DGVMs assume N deposition at the 1850 value for years 1700-1850.

Additional Data Sets
Each DGVM is free to use its own data sources for vegetation and soil types, soil textural properties, soil permeable depth and topography data (e.g., to simulate lateral water flows and soil erosion) as well as for the land mask.Additional data are required for fire-enabled DGVMs including, for example, spatio-temporal maps of lightning ignition and population density, the latter used to estimate human impacts on wildfire regimes.Given the uncertainty around lightning data sets (Kaplan & Lau, 2021), scaling factors, for example, proportion of cloud-cloud versus cloud-ground strikes, and potential need for model recalibration, and the fact that in TRENDY we want models to supply their best C-cycle representation, modeling groups are free to choose their lightning data set.Gridded population density is based on HYDE3.3.

Models
Seventeen DGVMs participated in TRENDY-v11 (Table 1, Table S2 in Supporting Information S1).While DGVMs are process-based, bottom-up models, the processes included, and their parameterizations vary across a Available but not provided to TRENDY.b Very rudimentary, represents some aspects of N Limitation on photosynthesis driven by soil processes.c Basic representation of demography (yearly cohorts, self thinning).
models.The key biological and land-use processes included in the individual DGVMs are listed in Table 1 (adapted from Table A1 in Friedlingstein et al. (2022a)).DGVMs run at a spatial resolution: 0.5 × 0.5 (longitude × latitude) degrees (or at a coarser resolution if necessary; ideally at 0.5°).

Factorial Simulations and Protocol
Initially all DGVMs conduct a spin-up with CO 2 and land cover at their year 1700 levels.While models are free to choose their land cover (typically consisting of the plant functional types they represent) the crop and pasture area distribution correspond to the year 1700.Some DGVMs further use 1700 harvest from LUH2 (Table 1).The CO 2 concentration is set to 276.59 ppm.In the absence of pre-industrial climate data, the climate data from the years 1901-1920 are used repeatedly until model pools reach equilibrium.Whether the climate data is cycled or randomly drawn is thereby model dependent.Subsequently to the spin-up, DGVMs perform four factorial simulations, S0-S3, from 1700 to 2021, as follows, and summarized in Table 2: S0: A control pre-industrial simulation in which all forcings stay as in the spin-up.S0 allows to diagnose any "cold start" issues and model drift.S1: In this simulation observed global CO 2 evolves over the historical period, the climate and land cover forcings stay at their pre-industrial levels like in simulation S0.Note: some models also use transient N deposition rates and/or population data in this simulation (see below).S2: In this simulation, CO 2 and climate evolve over the historical period, while the land cover stays at its preindustrial level.For the period 1700-1920, the climate data from 1901 to 1920 is used repeatedly and from 1921 to 2021 the climate data for the corresponding years is used.S3: All forcings (CO 2 , climate and land use) evolve over the historical period in this simulation.The meteorological data is used in the same way as in the S2 simulation.
These runs are used for attribution of changes in model's state variables and fluxes to the three global change drivers (CO 2 , climate change, and land use change) (Table 2).The S1 simulation mainly quantifies the effect of the anthropogenically increased CO 2 (Walker et al., 2020) but for some models also secondary effects following increases in N deposition and changes in fire dynamics following population increases (see below).The difference between S2 and S1 simulations quantifies the impact of climate variability and change, and the difference between S3 and S2 simulations yields the impact of land-use and land cover change on the evolution of the land carbon sink.Since the S2 simulation uses land-cover corresponding to the pre-industrial state (year 1700), the environmental change effects (CO 2 , climate) act on a larger area of natural land than in reality, that is, those areas affected by land-use change since the pre-industrial period.Thus, the simulated land sink (S LAND ) in the S2 simulation is likely overestimated globally given the net conversion of natural vegetation (e.g., forests) to agriculture, where the latter has faster carbon turnover times and a smaller sink potential.This process was first termed the land-use amplifier (Gitz & Ciais, 2003;Sitch et al., 2005), and more recently the "Loss of Additional Sink Capacity" (LASC) (Gasser & Ciais, 2013;Pongratz et al., 2014).Thus, DGVM land-use and land cover change emissions (E LUC ) derived from differencing S3 and S2 simulations, implicitly include this LASC term, which is substantial and accounts for ∼40% of the present-day DGVM E LUC (Obermeier et al., 2021).This estimate was based on additional DGVM factorial simulations in TRENDY-v8 (GCB2019) to estimate E LUC under pre-industrial and present-day environmental conditions, and thus used to diagnose the LASC (Obermeier et al., 2021).The inclusion of the LASC complicates comparison of E LUC estimates from DGVMs with other methodologies (e.g., bookkeeping and NGHGIs) which do not take LASC into account.Global Biogeochemical Cycles 10.1029/2024GB008102 In addition to the forcings described above, models that include a representation of the terrestrial nitrogen cycle (see Table 1) require pre-industrial or time varying N deposition and N fertilizer application rates as follows: S0: This simulation uses N deposition and fertilizer application rates corresponding to their pre-industrial level (year 1700).S1: Time varying N deposition rates but pre-industrial N fertilizer application rates are used.S2: Time varying N deposition rates but pre-industrial N fertilizer application rates are used S3: Time varying N deposition and N fertilizer application rates are used.DGVMs that represent fire use time-evolving population density, itself linked to fire ignition and suppression, in simulations S1-S3, and population density for the year 1700 in simulation S0.Our simplified logic is that LULCC and its direct consequences (Nfert) go together in E LUC , and all other indirect environmental changes (Ndep, population, climate, CO 2 ) in S LAND .

Criteria for Budget Inclusion
DGVMs are included in GCB assessments if they fulfill a set of three criteria for minimum model realism: 1. Steady state after spin-up.Steady-State defined as an offset of global NBP <0.10 PgC/yr, drift <0.05 PgC/yr per century (i.e., first is the average over 1700-2021, second is the NBP ∼ year slope × 100).Diagnosed from S0 run.2. The global net annual land flux (S LAND E LUC ) is a carbon sink over the 1990s and 2000s as constrained by global atmospheric and oceanic observations (Keeling & Manning, 2014).Diagnosed from S3 run.3. Global net annual land use flux (E LUC ) is a carbon source over the 1990s, based on results of book-keeping models (IPCC, 2007).Diagnosed from S3-S2 runs.

Benchmarking
DGVMs are evaluated and compared in the ILAMB benchmarking system and summary statistics given for each model (in summary table/ figures in the supplementary material of the annual GCB assessments).This allows to document model improvement each year, and to identify possible issues/model deficiencies to aid model development.We currently do not use the benchmarking results as criteria for budget inclusion or model weighting, but rather to encourage modeling groups to continue to improve their models.In future years and with consultation among participating groups, benchmarking could play a role in model inclusion or weighting.The selection of tests, metrics and performance criteria are subjective, and observation-based data uncertainty would further need to be included (Seiler et al., 2022).A medium-term objective will be to use benchmarking to weight model performance for the Multi-Model Ensemble Mean (MMEM), with aim to help reduce the Budget Imbalance, which is the remaining unexplained uncertainty in the global carbon budget.

Outputs
For GCB assessments four aggregated NBP values (PgC/yr) are requested from all simulations (S0 to S3): annual global NBP, annual northern extra tropics NBP (>30°N), annual tropical NBP (30°S < latitude < 30°N) annual southern extra-tropics NBP (<30°S); for the years 1700-2021.An extended list of gridded output variables for carbon and water cycling, and energy balance is requested (see Tables S3-S6 of Supporting Information S1), which are then available to the research community for spin-off studies.Variables are listed as either, Level 1: essential, which all DGVM groups are strongly encouraged to deliver, and Level 2 variables: as desirable for additional in-depth analysis/studies.To enable alignment with national inventories we need DGVMs to output NBP on a Plant Functional Type basis where possible (Grassi et al., 2018).This is to separate fluxes on managed versus unmanaged land.

Benchmarking Carbon and Hydrological Cycles
In GCB we provide a summary benchmarking table for the key carbon and hydrological cycle variables.An exhaustive benchmarking analysis has been performed by Seiler et al., 2022 on the earlier TRENDY-v9 with focus on uncertainty in models and observation-based data sets.We do not expect major differences in results and therefore give a broad overview of model performance for a small set of key carbon stocks and fluxes: biomass, soil carbon, and GPP; and water cycle variables: river runoff and evapotranspiration.In TRENDY, DGVMs have been evaluated within the ILAMB benchmarking tool (Collier et al., 2018).
Here we summarize the model-data evaluation for the S3 simulation, and the full suite of TRENDY ILAMB benchmarking results can be accessed below, which updated each GCB cycle provides an essential service to the modeling community: https://gws-access.jasmin.ac.uk/public/jules/TRENDY2022/GCP_220523_vn2.6_global/.Using four regional and global biomass data sets (Santoro & Cartus, 2021;Saatchi et al., 2011;Thurner et al., 2014; GlobalCarbon unpublished data) the MMEM performance is in reasonable agreement, in terms of spatial variation in biomass across biomes.The MMEM tends to under-estimate biomass against the ESA-CCI benchmark year 2010 (Figure 4) but overestimate the GlobalCarbon benchmark (see ILAMB link), despite the ESA-CCI benchmark only estimating aboveground biomass, while the other benchmarks estimate total biomass.This highlights how uncertainty across observationbased products needs to be considered in benchmarking analyses (Seiler et al., 2022).Against the regional biomass data sets, the MMEM overestimates biomass in the tropics, and in the boreal and tundra regions, and underestimates biomass in mid-latitude northern and temperate regions.H. Yang et al. (2023) suggests DGVMs may overestimate biomass in tropical forests in the absence of forest degradation processes (e.g., through selective logging and forest fires).
The ensemble performs well for soil carbon against the one global soil carbon data set (HWSD, Wieder, Boehnert, & Bonan, 2014;Wieder, Boehnert, Bonan, & Langseth, 2014; Figure 5), albeit overestimating in high latitudes, in Alaska and Siberia, against both the global and regional (see ILAMB link) benchmark data sets.Two DGVMs simulate very high soil stocks in the high northern latitudes.However, a note of caution, a second global soil database (https://soilgrids.org)also has higher soil C stocks in high latitudes.
DGVMs generally perform well against fluxnet-derived and FLUXCOM GPP for spatial variance over 1980-2014 (Figure 6), although the benchmark data are biased toward mid-latitude sites.The MMEM slightly underestimates the data-derived peak in amplitude in the seasonal cycle of GPP (not shown -see ILAMB link).
DGVMs perform very well in simulating spatial patterns of river runoff (Figure 7), particularly its interannual variability (see ILAMB link), however, tend to underestimate runoff in major river basins, for example, the Amazon, and in S.E.Asia, possibly due to the lack of lateral flow processes and river routing in most models.DGVMs also perform well against GLEAM evapotranspiration product, in particular the seasonal cycle (see ILAMB link).ILAMB also allows the evaluation of relationships between variables, for example, GPP versus ET (Figure 8), which is also available for each of the 10 main RECCAP2 regions (see 10 regional and global options under Global-Land): https://gws-access.jasmin.ac.uk/public/jules/TRENDY2022/GCP_220523_vn2.6/.

Global NBP and Environmental Drivers
The long-term evolution of the net land-atmosphere exchange and the attribution to individual drivers is shown in Figure 9. Since the 1950s, rising atmospheric CO 2 concentrations have increased carbon uptake to an average of 3.8 ± 0.8 PgC/yr over 2012-2021.The natural land sink is reduced by the effects of climate variability and change ( 0.6 PgC/yr ± 0.5 PgC/yr).Over the same period, E LUC has remained roughly constant ( 1.6 ± 0.5 PgC/yr).Hence, the land is a net sink of carbon globally (1.7 ± 0.6 PgC/yr) over the 2012-2021 period.On inter-annual timescales, climate impacts global land carbon fluxes primarily via ENSO variability (Keeling et al., 1995;Peylin et al., 2005;Piao et al., 2019;Tian et al., 1998).

NBP and Process Attribution
In recent decades, the simulated net land-atmosphere exchange is mainly driven by the CO 2 response of photosynthesis and biomass production (Walker et al., 2020;Figures 9, 10i, 10j, and 11).Further, increased nitrogen deposition may have also contributed to enhanced plant growth (Townsend et al., 1996).The models simulate a small reduction in global burnt area over 2012-2021 (Figure S2 in Supporting Information S1), in agreement with satellite-based data sets, which attribute changes primarily to cropland expansion and its intensification (Andela et al., 2017).However, the ability of recent high-resolution remote-sensing to detect widespread small fires cautions the validity of trends based on earlier coarse-resolution products (Chuvieco  (Jung et al., 2017;Tramontana et al., 2016).Compares time-means over the period 1980-2014. et al., 2022;Ramo et al., 2021).LULCC has reduced global fire emissions by ∼0.3 PgC/yr compared to the start of the 20 th century.The specific contribution of recent changes in climate-alone act to increase fire risk and increase global burnt area since around 1980 (Jones et al., 2022), but have a small impact on simulated carbon emissions (Figure S2 in Supporting Information S1).Indeed, while carbon emissions from fires have reduced in savannahs and increased in forests and approximately balanced globally (Zheng et al., 2021).There is no clear consensus as to how rising CO 2 concentrations influence global burnt area, but CO 2 -induced gains in biomass stocks lead to  Global Biogeochemical Cycles 10.1029/2024GB008102 increased fire emissions when vegetation is eventually burned.A more comprehensive assessment of fire-enabled DGVMs has been conducted by FireMIP (Hantson et al., 2016(Hantson et al., , 2020;;Li et al., 2019;Rabin et al., 2017).

Ecosystems Contribution to NBP
Biome contribution to mean, trend, and interannual variability in the natural sink is shown in Figure 12 and how their contributions vary per decade in Figures S3 and S4 of Supporting Information S1.Three biomes (tropical forest, semi-arid and grasslands/crops) are the most important on all three timescales, with the greatest contribution of the semi-arid biome to the interannual variability, as found in Ahlström et al. (2015), Piao et al. (2019), andPoulter et al. (2014).IAV in the carbon cycle is driven by ENSO climate variability, and it is notable that the last decade has been characterized by La Niñas (2011,2020,2021), punctuated by one large 2015/16 El Niño, in contrast with an El Niño dominated decade of the 1990s.The IAV contribution does not appear to vary significantly by decade (Figure S3 in Supporting Information S1).The relative biome contribution to the long-term trend also does not change (Figure S4 in Supporting Information S1), except for a decrease in contribution of the tropical forest between 1982 and 2001, and 2002-2021.This may indicate a decline in the tropical intact forest sink strength, which has been observed in forest plot networks (Brienen et al., 2015;Hubau et al., 2020), and broadly attributed to drought-impacts on forest function.DGVMs are likely to underestimate drought impacts, as they lack explicit drought mortality, however, will simulate reductions in GPP with water deficit, albeit only few DGVMs explicitly represent plant hydraulics, a field of active recent development (Eller et al., 2020;S. Wang et al., 2020;Y. Wang et al., 2020).

Regional Contribution to NBP
Regional attribution (RECCAP2-regions) is shown in Figure 13, and in Table 3 for the decades 2000-2009 and 2012-2021.A positive effect of increasing atmospheric CO 2 on the simulated land carbon sink is broadly consistent with empirical understanding and measurements at smaller scales albeit measurement of the CO 2 effect at regional to global scales is not currently possible (Walker et al., 2020).Boreal and temperate regions are simulated to have increased carbon stocks in response to moderate warming, for example, with longer growing seasons and enhanced nutrient cycling and N deposition (Piao et al., 2007), whereas low latitude regions are negatively impacted by climate, suggesting that warming combined with drought reduces plant productivity.The largest LULCC emissions are in the tropical regions over the past decade.Land abandonment and preservation, reforestation efforts and improved land management practices particularly in the second half of the 20th century (Nabuurs et al., 2003) led to carbon sinks in regrowing temperate forests, although there is evidence of potential sink saturation and subsequent declines during the 21st century, for example, in Europe (Nabuurs et al., 2013;Winkler et al., 2023).In contrast, DGVMs simulate a LULCC source in mid-latitude regions, which may reflect that forest-age, demography and wood harvest is not properly represented in most DGVMs, and emphasizes the importance of simulating gross transitions.This may lead to an underestimate in the DGVM simulated net land carbon sink in these regions.Russia and North America are simulated to have the largest net land carbon sink with positive effects of elevated CO 2 and climate on carbon sequestration.Typically, the tropical regions experience the largest gross impacts of CO 2 , LULCC, and climate, which counterbalance to give a more moderate net land-atmosphere exchange estimate.The budgets estimated by DGVMs reported here can be compared with the corresponding best estimates currently being developed in the RECCAP2 activity.
Although climate variability and change have the smallest impact on the mean land sink in most regions over the 2012-2021 decade it is the most important driver of regional trends in the 21st Century (Figure 14).Model results suggest a small reduction in land sink capacity in tropical South America and South-East Asia, the former driven by negative impacts of climate change offsetting decadal reductions in deforestation.In contrast climate variability and change has the largest positive contribution to the large land sink trend in Africa and North America.

Perspectives
TRENDY aims to further develop and apply model benchmarking in the future, aiming to weight modeled NBP by model performance across a range of observational-based benchmarking tests; that will require consideration of uncertainties in the observation data products.Following the findings of O'Sullivan et al. ( 2022), further internal-cycling variables will be requested and analyzed to better isolate uncertainty and sensitivity in process representations among DGVMs.These variables include plant allocation and turnover to and from individual tissue pools, and input and loss from individual litter and soils pools, and separate C fluxes associated with individual disturbance agents (e.g., fire, drought and background mortality).Further efforts are needed to reduce the large uncertainty in the CO 2 fertilization effect across DGVMs (Chen et al., 2022;Haverd et al., 2020;Huntingford & Oliver, 2021;Medlyn et al., 2015); empirical efforts with Free-Air-Carbon Dioxide Enrichment experiments in mature temperate forest (BiFOR) and soon the Amazon forest (AmazonFACE) promise large knowledge gains and data sets to constrain DGVMs.DGVM benchmarking is of crucial importance.Modelers need to actively engage with data set providers to understand the most appropriate data sets to include in the benchmarking exercise and discuss how to interpret uncertainty in these products.In TRENDY, DGVMs are forced with one climate data set (CRUJRA) and one land use change reconstruction.Studies with single DGVMs or other types of carbon cycle/book-keeping models have quantified the uncertainty in land carbon cycle variables related to alternative climate forcing data sets (Hardouin et al., 2022;Wu et al., 2017) and land-use data sets (Bayer et al., 2017(Bayer et al., , 2021;;Ganzenmüller et al., 2022), and found that these can be large (Hartung et al., 2021).A systematic assessment using multiple climate forcing data sets is needed across the range of TRENDY models to quantify the impact of climate forcing uncertainty.Many DGVMs (including those coupled within host Earth System Models) contribute to other MIP activities, such as ISIMIP/CMIP-land simulations (e.g., Hardouin et al., 2022), CMIP ESM simulations, the N2O Model Intercomparison (NMIP, Tian et al., 2018Tian et al., , 2019) ) and FireMIP (Rabin et al., 2017).So, improved benchmarking under TRENDY (a) helps to constrain these other assessments including future climate change work, (b) helps to direct model development that feeds into CMIP7, and (c) supports protocol and forcing data sharing among multiple assessments to improve comparability and reduce modeling teams workload when contributing to more than one assessment, for example, GCP-C (Friedlingstein et al., 2022a), GCP-CH 4 (Saunois et al., 2020) GCP-N 2 O (Tian et al., 2020).
Few DGVMs represent ecosystem demography and climate-stress induced mortality, and they vary widely in which land management practices (like wood harvesting, tillage, irrigation) are considered.There is currently a large uncertainty in simulated biomass and its change between DGVMs (H.Yang et al., 2020).Remote-sensing offers great potential to better constrain and calibrate models for biomass, disturbance effects, regrowth and recovery in different forest types (Heinrich et al., 2021(Heinrich et al., , 2023).These will be key processes for model development to robustly represent the impact of climate change and extremes on the land sink (e.g., response to and recovery from heat extremes, drought, and fire), and the impact of land management, on the near-term evolution of the land sink and on land-use fluxes.This information is key to the periodic UNFCCC Global Stocktakes and to inform the remaining human emissions to avoid dangerous climate change.Additionally, TRENDY model outputs have been extensively used to assess retrospectively the impacts of extreme events or other events of interest, for example, El Niño events (Bastos et al., 2018).These typically lag the events of interest by several months if not years.A new initiative is to conduct fast-track assessments of the impact of regional climate extremes and megafires on the carbon cycle in near real time.A new protocol based on Bastos et al. (2020) is currently being developed to complement annual TRENDY, whereby DGVMs are run at set points during a given year (or when required), with forcing updated from ERA5 reanalysis (Hersbach et al., 2020), with an assessment latency of ∼3 months.DGVMs will be applied for an "S2"-style (CO 2 and Climate only, with fixed present-day LUC) simulation, and can be restarted from the end of the previous simulation (to avoid computationally expensive simulations).The additional S2 simulation with ERA5 can also be contrasted with the DGVM simulation with CRUJRA to quantify the impact of climate-driven uncertainty, although the protocols differ considerably.Rosan et al. (2021) demonstrated the benefit of adopting country-specific LULCC data sets and expertise, particularly those based on low latency, spatio-temporal Earth Observation data.We envisage extending this approach, first for those countries with the largest contribution to ELUC and eventually all countries.There are Global Biogeochemical Cycles 10.1029/2024GB008102 recognized issues relating to the LULCC data used for China (X.Wang et al., 2024;Yu et al., 2022) and this represents a major uncertainty in GCB.This in part relates to changing methodologies in reporting agricultural areas to FAO through time (Yu et al., 2022).We aim to adopt these alternative data sets in the future.
A final medium-term objective is to combine the land GHG budgets for CO 2 , CH 4 , and N 2 O in near-real-time assessments, for example, to contribute to initiatives like the WMO for a GHG monitoring system (WMO, 2023).

Summary
The TRENDY project has completed its 11th annual iteration and delivered natural land sink and land-use emission estimates to the GCP-Global Carbon Budget 2022.The project has served the global land modeling community for over a decade, providing a forum to evaluate, intercompare and guide future model development, and its data have been used in over 100 peer-reviewed publications (https://blogs.exeter.ac.uk/trendy/publications/) to enhance our understanding of land biogeochemical cycles.In this study we have described the latest TRENDY project protocols and input data sets and provided summary information on the terrestrial carbon cycle  across large regions and their process attribution, and individual model performance across a range of benchmarking tests.The largest gross fluxes are simulated in the tropical regions, with the negative effects of land-use change, through deforestation, plus a small negative contribution of climate variability and change, that counteract simulated positive effects of CO 2 fertilization on the net land-atmosphere exchange.In contrast, in extratropical regions CO 2 fertilization dominates and is supplemented by a small positive effect of climate change, extending high-latitude growing seasons, with a smaller contribution of land-use change on the net landatmosphere exchange.

Figure 1 .
Figure 1.Temporal variation of all meteorological driving variables (global means over land presented).

Figure 2 .
Figure 2. (a) Historical land use and land cover change from LUH2 based on HYDE3.3 from 1900 to 2021.(b) Temporal variation in global wood harvest over the same time period.

Figure 8 .
Figure 8. Relationships between variables from the benchmark data sets (a) GPP Fluxcom versus ET GLEAM (b) GPP versus ET Multi-Model Ensemble Mean (MMEM) (c) data-derived (black) versus MMEM (brown) relationship between GPP and ET.Colors on (a) and (b) indicate the fraction of gridboxes in each bin.Error bars on (c) indicate ±1 standard deviation of GPP at gridboxes with ET in a given bin.

Figure 9 .
Figure 9. Timeseries of global NBP attributed to processes CO 2 , Climate and Land Use and Land Cover Changes and the sum of the three (NET).Shown is the multi-model ensemble mean.

Figure 10 .
Figure 10.Dynamic Global Vegetation Model Multi-Model Ensemble Mean and 1 standard deviation carbon fluxes for the present-day (2012-2021).

Figure
Figure Dynamic Global Vegetation Model Multi-Model Ensemble Mean and 1 standard deviation carbon stocks in vegetation and soils for the present-day (2012-2021).

Figure 12 .
Figure 12.Decomposition of natural land sink into mean, trend, and interannual variability for different biomes, following Ahlström et al. (2015) over the period 1982-2021.Box-plots show Multi-Model Ensemble Mean, interquartile range, and range.

Figure 13 .
Figure 13.Regional mean net land-atmosphere exchange (symbols) and process attribution (bars) for the 10 main RECCAP2 land regions over 2012-2021.

Table 2
Factorial Dynamic Global Vegetation Model Simulations With Varying (Observed Historical Evolution), Pre-Industrial (Constant) Forcing for Environmental Conditions (Climate, Atmospheric CO 2 Concentrations, NitrogenDeposition and  Fertilizer Application, Land Use and Land Cover Changes)