Process-oriented analysis of dominant sources of uncertainty in the land carbon sink.

The observed global net land carbon sink is captured by current land models. All models agree that atmospheric CO2 and nitrogen deposition driven gains in carbon stocks are partially offset by climate and land-use and land-cover change (LULCC) losses. However, there is a lack of consensus in the partitioning of the sink between vegetation and soil, where models do not even agree on the direction of change in carbon stocks over the past 60 years. This uncertainty is driven by plant productivity, allocation, and turnover response to atmospheric CO2 (and to a smaller extent to LULCC), and the response of soil to LULCC (and to a lesser extent climate). Overall, differences in turnover explain ~70% of model spread in both vegetation and soil carbon changes. Further analysis of internal plant and soil (individual pools) cycling is needed to reduce uncertainty in the controlling processes behind the global land carbon sink.

Introduction

Over the last 60 years, there has been a continuous rise in anthropogenic CO2 emissions. Around the equivalent of a quarter of these emissions have been taken up by the land biosphere (known as the natural land sink), acting as a strong negative feedback to mitigate climate change[1]. To be able to project the global carbon cycle (and hence the climate response) in the future, we need to understand the underlying processes and their timescales that drive the contemporary land sink. There are many distinct but interdependent mechanisms that regulate the flow of carbon into, through, and out of the land. The timescales at which these processes (e.g. photosynthesis, allocation, plant growth, litterfall, plant mortality and soil turnover) act on the carbon cycle range from days to centuries, and the interplay between these determines the changes in land carbon storage[2].

Recent work suggests the global net land sink is located in northern latitudes[3]. The attribution of drivers is highly uncertain, but a combination of increasing atmospheric CO2 concentrations[4], reactive nitrogen (N) deposition[5] and atmospheric warming[6] are likely responsible. Tropical lands are probably nearer net zero carbon sinks due to large land-use and land cover changes (LULCC) carbon losses counteracting the ‘natural’ sink[7,8]. Atmospheric observations also show that the Amazon forest is currently carbon neutral (‘natural’ sink = LULCC-driven source)[9], and this provides further evidence of a limited net tropical sink. Yet, the largest gross fluxes between the land and atmosphere are often found in tropical regions[10], and so changes in tropical ecosystem functioning can have significant global impacts.

Process-based dynamic global vegetation models (DGVMs) that simulate processes of carbon uptake and release can help to elucidate the roles of individual drivers (rising atmospheric CO2, changes in climate, nutrient deposition, and LULCC), attribute the change to processes, and quantify regional sinks over timescales needed given the multitude of timescales over which terrestrial processes act (i.e. beyond the period of reliable empirical and remote-sensed data).

DGVMs are used to estimate the natural land sink (hereafter simply land sink) as part of the global carbon budget (GCB[8]). Globally, the DGVM multi-model mean estimate of the global land sink is consistent with the global carbon budget residual land sink (the difference between fossil and land-use emissions and atmospheric and ocean sinks, see Table 5 in ref. 8), however, there is a significant spread across models, and here we show the spread widens at regional scales, when quantifying changes in vegetation and soil carbon, or attributing changes in internal processes to external drivers. Therefore, while model ensembles are helpful in analysing global-scale processes, they must be interpreted in the context of their process representation and unique and collective biases[2,4].

These deficiencies are of high significance when trying to understand past changes but also when DGVMs are used to make predictions of future carbon cycling[11], as limited trust can be placed in projections when certain fundamental processes and carbon-climate feedbacks are not fully captured[2]. Therefore, it is of high priority to identify the leading sources of uncertainty in DGVMs and understand the relationship between drivers and processes on varying spatial scales. In this study, we use a triple (3D-matrix) approach to identify where models agree and, just as importantly, disagree, and thus guide future modelling efforts:

What are the (1) external drivers (concurrent rises in atmospheric CO2 and N deposition, climate and LULCC), (2) main regions (tropics and extra tropics) and (3) processes (production vs turnover) primarily responsible for the changes in the global net land carbon sink?

We use the suite of 18 DGVMs from the GCB2021 (TRENDYv10; ref. 8) to quantify changes in net carbon exchange and carbon stocks over the period 1959–2020. TRENDYv10 provides a set of simulations to attribute these changes to drivers, and we use a process attribution framework to decompose changes in carbon stocks into those driven by productivity and turnover separately. This framework enables us to express productivity and turnover-induced changes in carbon stocks in units of PgC, which allows for a direct comparison between both processes, which have units of carbon per unit time and time, respectively (see Methods).

Results

Drivers of global and regional land sinks

The global net land sink is derived from the difference between the fossil fuel emissions (EFOS) and the CO2 accumulation in the atmosphere (GATM) and uptake by the oceans (SOCEAN). We refer to this estimate of the net land sink as the ‘observed’ GCB budget constraint (= EFOS-GATM-SOCEAN; see methods and ref. 8) which grew from 0.2 ± 0.4 (mean ± std. dev) PgC yr−1 in the 1960s to 1.7 ± 0.6 PgC yr−1 in the decade 2011–2020, with DGVMs (S3 simulation; see Methods) capturing the increase (−0.1 ± 0.6 to 1.6 ± 0.5 PgC yr−1; Fig. 1). DGVMs suggest the enhanced net sink over the past 60 years has mainly been driven by rising atmospheric CO2 concentrations and nitrogen deposition (1.2 ± 0.2 PgC yr-1 over 1960–1969 rising to 3.5 ± 0.8 PgC yr-1 over 2011–2020) (Fig. 1b), with the sink partially offset by relatively constant net LULCC emissions of 1.3 ± 0.5 PgC yr−1 (predominantly arising from tropical deforestation and shifting cultivation[12]; Fig. 2 and Supplementary Fig. 1). DGVMs indicate that climatic variability drives the large year-to-year changes (±2 PgC yr−1) in the land sink (Fig. 1), with long-term (multi-decadal) climate trends reducing the land sink since the 1980s (−0.4 ± 0.5 PgC yr−1 over 1980–2020) (Fig. 1b).

Fig. 1

Global net land sink and attribution to drivers.

a Net annual land carbon sink (PgC yr−1) as estimated by dynamic global vegetation models (DGVMs) with all drivers varying (red) and the top-down Global Carbon Budget (GCB) constraint (black) and b the decomposition of the DGVM net sink into contributions from rising atmospheric CO2 concentrations and N deposition (blue), changes in climate (yellow), and land-use and land cover change (grey). Thick lines show a locally fitted regression. Shading around DGVM estimates corresponds to 1σ, and the uncertainty on the GCB constraint is taken from ref. 8.

Fig. 2

Driver attribution to spatial changes in total ecosystem carbon.

Maps show the multi-model mean a net change in ecosystem (vegetation and soil) carbon (kgC m−2) from 1959–2020 and b the contribution of each driver to overall change; CO2 and N deposition (green), climate (red), and land-use and land cover change (blue). Stippling in panel a indicates <80% of models agree on the direction of change. The colours in panel b are calculated by assigning a red-green-blue (RGB) value to each grid depending on the relative magnitude of change due to each driver. Transparency is determined by the magnitude of the net change in panel a.

The CO2 and N deposition-driven land carbon sink is located in northern and tropical forests (Fig. 2). LULCC-induced carbon losses occur across the globe but are most apparent in tropical latitudes and certain hotspot regions across the northern hemisphere (China, USA and West Eurasia), which recently have been (within our study period) densely forested (Fig. 2 and Supplementary Fig. 1). There are signs of European carbon sink in part driven by a/reforestation or natural land (re)establishment. The impact of changes in climate can also be detected across the globe, with vegetation and soil carbon losses in the Amazon, the Sahel and South Africa, and climate-induced carbon gains in east Brazil, Australia, and across the high northern latitudes (Fig. 2 and Supplementary Fig. 1). In summary, the DGVM ensemble agrees on the large sink located in the world’s forests. However, the models are not in full agreement in the direction of change in ecosystem carbon across much of the globe, in particular in regions with competing CO2 and LULCC effects (Fig. 2a).

Key processes and uncertainties behind the sink

The top-down global GCB budget constraint combined with the bottom-up DGVM estimates give high confidence the land has been a net sink of carbon over the last 60 years, however, the relative contribution of each driver (e.g. CO2, N deposition, climate, LULCC), and ecological attribution (vegetation, soil) of the sink remains elusive.

The DGVM multi-model mean suggests a net global increase in both vegetation (ΔC = 28 ± 26 (mean ± std. dev) PgC) and soil (ΔC = 21 ± 32 PgC) carbon stocks, although the uncertainty is large (Fig. 3a, b), with two models (LPX-Bern and YIBs) simulating a net loss of vegetation carbon and three models (CLASSIC-N, ISAM and LPJ-GUESS) simulating a net loss of soil carbon (Fig. 3c).

Fig. 3

Temporal changes in global vegetation and soil carbon stocks.

Time-series show the change in global a vegetation () and b soil () carbon stocks from 1959–2020 due to each of the three external drivers (CO2 and N deposition, climate, land-use, and land cover change). Lines represent the mean of the dynamic global vegetation models (DGVMs) and shading the ±1σ of the DGVMs. The DGVM output is first smoothed using a fourth-order spline. The cumulative net (sum of three drivers) change in global carbon stocks by 2020 is shown (red crosses show each model and red circle shows the model mean). c Shows the change in vegetation and soil stocks for each of the 18 models and the grey region is the Global Carbon Budget net land sink constraint (see Methods).

This relatively high uncertainty is in part due to large opposing fluxes (Fig. 3a) driven by increasing atmospheric CO2 and N deposition (carbon increase) vs LULCC (carbon decrease at global scale), with soil carbon gains (from rising CO2) also counteracted by negative fluxes driven by climate change (Fig. 3b). In general, there are no observational constraints on long-term (pre-satellite era) changes in global vegetation or soil carbon, and so reducing model uncertainty is challenging. Therefore, although the DGVMs generally capture the global net carbon sink, the attribution to vegetation or soil stocks is completely unconstrained (Fig. 3c) and is driven by alternate model structures and parameterisations between DGVMs.

DGVMs have inherently different baseline productivity, allocation, and turnover rates as well as varying degrees of sensitivity of processes to environmental change. To gain a deeper understanding of the processes driving carbon sink changes, we use our attribution framework (see Methods) to associate the modelled changes in vegetation and soil carbon to changes in inputs (NPP for vegetation and litterfall for soil) and outputs (turnover rates for vegetation and soil). For each of these processes, we identify their changes due to CO2, climate change or LULCC, according to the models.

The ensemble-simulated increase in global vegetation carbon over the last 60 years (ΔC = 28 ± 26 PgC; Fig. 4a) was driven by enhanced plant productivity in response to rising CO2 concentrations and N deposition, with small (but highly uncertain) losses due to changes in vegetation turnover ( = 87 ± 26 PgC, = −3 ± 21 PgC; Fig. 4a). Although the central estimate of the net CO2 and nitrogen deposition biomass response appears relatively well constrained at global scale ( = 84 ± 31 PgC), model estimates range from a gain of 30 PgC (YIBs) to a gain of 150 PgC (LPJ) in biomass due to rising atmospheric CO2 (Supplementary Fig. 2).

Fig. 4

Process and driver attribution of changes in global vegetation and soil carbon stocks.

Change in global a vegetation () and b soil () carbon stocks over 1959–2020 (PgC). The contribution to net changes in carbon stocks (green bars) from changes in inputs (net primary productivity for vegetation () and vegetation to soil flux for soil (), red bars), outputs/turnover ( for vegetation and for soil, orange bars), and the interaction term ( for vegetation and for soil, blue bars) are shown. The bars depict the multi-model mean with the range as ±1σ of the models. The arrows show the direction of change in carbon stocks due to each process. The panels from left to right show the changes due to all drivers varying (ALL), changes in atmospheric CO2 and N deposition, climate (CLIM), and land-use and land cover change (LULCC).

Model spread is driven by a combination of uncertainty in the response of both NPP (e.g. driven by uncertainty in leaf-level photosynthetic response, scaling to canopy and landscape scales, carbon allocation and nutrient limitation) and vegetation turnover (τ) (e.g. driven by the treatment of nutrient limitation coupled with differences in allocation and in some cases forest-stand packing constraints[13]) to rising atmospheric CO2 and N deposition (Fig. 5 and Supplementary Fig. 3). There is some confidence in the direction of change in global and regional NPP, albeit not in the magnitude, whereas models disagree on the direction of Δτ due to rising CO2 (Fig. 4). CO2 and N deposition-driven Δτ depends on changes in stand dynamics (in some models) and plant allocation. For example, an increase in short-lived root production, at the expense of wood growth, to alleviate nutrient limitations on plant growth can reduce τ. This increase in root allocation also changes vegetation nitrogen demand (because of the substantially higher C:N in wood compared to roots), which feeds back onto the whole-plant productivity response to rising CO2. Conversely, some models increase their wood allocation fraction when NPP (or production) increases, which will increase τ[14].

Fig. 5

Attribution of uncertainty to processes in modelled changes in carbon stocks.

The relative uncertainty (defined as the standard deviation among model estimates) in the change in global vegetation () and soil () carbon stocks resulting from each of the driving terms (Eqs. 9 and 10 in Methods). The four terms are baseline input (baseline productivity ) for vegetation and baseline litterfall/mortality () for soil), change in inputs (change in productivity ) for vegetation and change in litterfall/mortality () for soil), baseline turnover ( for vegetation and for soil), and change in turnover ( for vegetation and for soil). For each of these terms in Eqs. 9 and 10, we calculate the standard deviation in and using the multi-model mean values of all other terms in the equations and the individual model values for that term.

This increase in vegetation stocks has been partially offset by emissions associated with LULCC from vegetation ( = −58 ± 21 PgC; Fig. 4a), predominantly in the tropics (Supplementary Fig. 1), which increases vegetation turnover (loss of biomass to the atmosphere or wood products). About half ( = −23 ± 22 PgC, ~40%; Fig. 4a) of the net LULCC-driven global vegetation carbon losses are from a reduction in inputs (NPP) to the land, as crops or pastures often have lower productivity than the forests that they replace. The other half ( = −29 ± 18 PgC; Fig. 4a) is attributable to changes in turnover related to forests containing woody biomass, which has a slower turnover than leaves and fine roots. The turnover effect is further aggravated through the phenomenon that due to their longer-lived biomass, forests take up more carbon under rising CO2 levels, but this sink is lost by clearing for agricultural use. This is known as the “loss of additional sink capacity” (LASC)[15,16] and accounts for the impact of environmental (CO2 and climate) changes on carbon uptake of deforested land (found in S3 simulation) compared to potential vegetation (found in S2 simulation). Ref. 16 attributes ~40% of DGVM estimated LULCC vegetation and soil carbon losses to the LASC, albeit with a different methodology and study period.

All models simulate a net loss of global biomass following LULCC, although the magnitude is not well constrained. Model spread is driven by differences in LULCC process representation in models[7]. We find that models that include wood harvest (routine harvest of established managed forests), grazing harvest, or shifting cultivation simulate larger C losses than models without these land management processes (Supplementary Fig. 4), in line with earlier studies showing the high importance of land management on vegetation carbon stocks[17]. Further, in northern ecosystems, models do not agree on the direction of or Δτ following LULCC (Supplementary Fig.  5a). In addition to the land management processes mentioned above, uncertainty in large-scale regional changes in C is in part driven by smaller-scale regional compensations, where losses in Russia and USA are somewhat countered by a forest regrowth sink in Eurasia (Supplementary Fig. 1). Therefore, the balance between these two opposing fluxes determines net northern changes in C from LULCC, with uncertainties in both carbon loss fluxes and regrowth uptake[18,19]. In addition, spatial patterns of northern LULCC losses and gains are not entirely consistent between models (Supplementary Fig. 6), which adds to uncertainty in local, regional, and hemispheric net LULCC-driven changes in vegetation carbon.

Changes in climate have a much smaller but similarly uncertain influence on vegetation processes at a global scale, with no agreement among models on the direction of change in biomass driven by shifts in NPP ( = −1 ± 14 PgC over 1959–2020) or vegetation turnover ( = 3 ± 9 PgC) at global scale (Fig. 4a). However, models align closer at regional scales, where northern warming has stimulated productivity ( = 11 ± 6 PgC; Supplementary Fig. 5a), and in the tropics, reductions in productivity (due to warming and/or changes in precipitation) lead to losses of vegetation ( = −12 ± 12 PgC; Supplementary Fig. 5b). Overall, these opposing regional trends cancel and explain the small global signal and the model uncertainty at the global scale.

Increased production and loss of biomass enhance soil inputs ( PgC) and drives the growth in global soil carbon stocks (ΔC PgC), albeit with low confidence in the magnitude (Fig. 4b). Our analysis also indicates that enhanced CO2 concentrations increase soil turnover, which is a direct result of enhanced inputs to the fast-turnover litter and surface soil pools. Our methodology simplifies the TRENDY model structures into a single soil pool to represent the entire soil system, and so a relative increase in the faster pools moves the aggregate soil pool turnover (τ) towards that of the litter and surface soil pools—a phenomenon known as false-priming[20]. In addition, changes to soil conditions (e.g. increased soil moisture[13]) following increased atmospheric CO2 can impact τ, although it is difficult to separate the false-priming and actual changes in τ with the set of simulations used in this study. This simplification leads to low confidence in the partition of net soil carbon changes into the input and turnover-driven changes, whereas the net change is well constrained (PgC).

Reduced litterfall (input to soil) due to the conversion of forests to agricultural land is the predominant pathway of soil carbon loss, in line with previous studies[21] ( = −36 ± 76 PgC,  = −25 ± 29 PgC). This is particularly evident in the tropics (Supplementary Fig. 7b) with large-scale deforestation over the last several decades, while deforestation and reforestation are more in balance in the extra tropics. However, there is relatively high uncertainty in the magnitude of LULCC impacts on soil carbon inputs () (Fig. 4b), with models ranging from losses of over 100 PgC to gains of 45 PgC over the period 1959–2020. Land management processes such as crop and wood harvesting (in particular, the frequency of harvest and fraction of biomass removed and respired elsewhere) are leading uncertainties in simulating LULCC impacts on soil carbon (Supplementary Fig. 4)[22,23]. Turnover-driven changes play less of a role in soil C changes in our models. This could be due to offsetting effects such as a decrease in soil C due to the altered quality of the litter input (a higher fraction of faster decomposing material for agriculture as compared to forests) or a reduction in fire-related losses with the transformation of natural ecosystems[21].

Increased heterotrophic respiration rates and subsequent reductions in the turnover time of the soil pool due to global warming entirely offset the increase in soil carbon from climate-enhanced NPP and litterfall (Fig. 4b). Overall, changes in climate caused a net loss ( = −13 ± 12 PgC) of global soil carbon. We find the impact of climate change on carbon storage is accelerating, in particular for northern soil and tropical biomass, where carbon stock changes from increased/decreased productivity (for northern soil and tropical biomass, respectively) have accelerated by 19 TgC yr−2 and 9 TgC yr−2, since the turn of this century (Supplementary Fig. 8). This equates to annual mean carbon gain/losses due to changes in productivity of 1.6 PgC yr−1 and −0.3 PgC yr−1 over 2001–2020 compared to 0.6 PgC yr−1 and −0.2 PgC yr−1 over 1961–2000, for northern soil and tropical biomass, respectively (Supplementary Fig. 8b, e). This acceleration is likely driven by warming in the north (2001–2020 April-August mean 0.9 °C warmer than 1961–2000) and the tropics (2001–2020 annual mean 0.6 °C warmer than 1961–2000)[24,25], as well as potential changes in sensitivity of the biosphere to warming[26].

Discussion

The ‘observed’ GCB global net land sink is captured by the TRENDYv10 DGVM ensemble. We find a robust agreement that rising atmospheric CO2 concentrations and N deposition drive the land carbon sink, and changes in climate along with LULCC lead to a source of carbon in the atmosphere. Further, the DGVMs corroborate that the net land sink is located in northern latitudes, with the tropics more carbon neutral, in line with recent independent top-down estimates[3,8].

However, the models do not entirely agree on the partition of the net land sink between vegetation and soil. Thirteen of 18 models indicate an increase in global vegetation and soil carbon stocks, although the magnitude of change is highly unconstrained (ranging from ~0 to +80 PgC over the past 60 years). The remaining five models simulate a reduction in either vegetation or soil stocks, highlighting the disparity between models regarding internal carbon cycling processes.

There are no observations of global carbon stocks covering the entire study period, but a recent synthesis of observation-based biomass change estimates in global forests suggests a net biomass sink of 0.3–2.1 PgC yr−1 since 2000[27], giving confidence to the models that simulate net gains. However, these observational studies do not constrain the magnitude of the vegetation sink any more than the DGVM mean biomass sink of −0.1 to 1.7 PgC yr−1 (averaged over 2001–2020). Furthermore, we cannot even be certain in the direction of change in global soil carbon as direct observations of large-scale changes in soil carbon stocks do not exist[2,4].

Our analysis indicates that baseline vegetation turnover rates and the turnover response to rising atmospheric CO2 and N deposition are the key uncertainties in modelled ΔC, with a smaller but still significant contribution from changes in NPP. Differences in modelled vegetation turnover are driven by variations in whole-plant mortality rates, simulated tissue lifespan, and the allocation of NPP to plant components with inherently different turnover rates and their responses to changing environmental conditions[14,28]. For example, some models may reduce plant mortality and turnover as atmospheric CO2 concentrations rise. Increased carbon stores could supply maintenance respiration in periods of photosynthetic stress[29,30], and increases in water-use efficiency may alleviate the impact of drought[31,32]. In contrast, shifts in plant community composition or self-thinning dynamics may increase stand-level mortality[33]. These structural differences lead to no agreement on the direction of Δτ in the models presented here.

In addition to turnover-related uncertainties, the increase in NPP and biomass due to rising CO2 is not well constrained. The choice of photosynthesis and stomatal model determines the leaf-level response, with a wide range of implementations and outcomes in DGVMs[34]. Scaling from leaf to canopy level causes further discrepancies due to different structural assumptions between models as to the vertical distribution of light, nitrogen, photosynthetic capacity, and the treatment of sunlit and shaded leaves[34]. Further, NPP and biomass production is constrained by stoichiometric nutrient requirements. Nitrogen availability mediates the biomass response to rising CO2[35], and therefore it is likely that differences in the inclusion and representation of carbon-nitrogen coupling between models contribute to uncertainty in ΔC[36,37]. Overall, these structural uncertainties manifest as a four-fold range in estimates of the CO2 effect on NPP and biomass (gains of 41 PgC (CLASSIC-N) up to 157 PgC (ISAM) over the previous 60 years).

Furthermore, LULCC-driven reductions in NPP and woody carbon stocks have removed a substantial portion of global vegetation carbon. However, the exact magnitude of carbon loss is not well constrained by the DGVMs. Biases in the carbon density of deforested land[38] and model-specific choices of whether to replace forest or grassland for agricultural expansion are leading sources of error when calculating LULCC losses[39]. Our results highlight that, firstly, variation in the representation of relevant processes, in particular, forest management (e.g. treatment of wood harvest) and shifting cultivation, leads to systematic differences between modelled ΔC [7,40]. Second, DGVMs are missing carbon losses associated with degraded forests, which can exceed deforestation losses, indicating simulated LULCC losses could be underestimated[40-42]. Finally, uncertainty exists with LULCC maps used to drive DGVMs as historical land-use is not perfectly known and products differ in the land cover types and transitions included[40].

The models suggest that soil has been a net carbon sink globally over the past six decades, although there is over 100% uncertainty on the magnitude (21 ± 32 PgC). This relatively large range is due to the opposing impacts of rising CO2 and climate/LULCC effects which partly cancel out. Further, the impact of LULCC on soil carbon is also difficult to quantify robustly with the suite of DGVMs used here, as models do not agree on the sign of or following LULCC. In general, changes in soil carbon stocks are difficult to constrain due to a lack of large-scale observations. It is likely historical LULCC caused additional carbon fluxes as a result of land erosion[43,44], degradation of agricultural soils[45,46], and losses from drained peatlands[47]. These management processes are not included in models and most do not simulate peatlands at all, potentially causing errors in simulated soil carbon fluxes.

Overall, the soil carbon sink is caused by increased litter inputs from CO2-driven vegetation growth. There is some evidence additional litter and root exudates enhance soil carbon stocks when nitrogen availability is high[48]. It is important to stress that the model simulations include widespread nitrogen deposition, which may have helped to sustain a strong CO2 response of biomass and soil carbon sinks[5,49].

However, additional soil inputs can accelerate organic matter decomposition to release plant-available nitrogen (via priming effects[50]), which reduces soil carbon storage[51-53]. DGVMs do not simulate priming effects or explicitly account for microbial activity, mineral association, aggregation, or mycorrhizal fungi interactions, all of which regulate soil carbon turnover and accumulation rates[54]. DGVMs represent soil carbon with a set of cascading pools, where decomposition losses are determined by first-order decay rates (although this is an active area of model development[55]). Therefore, the modelled soil sink has to be interpreted in the context of model structure and the limited ability of DGVMs to capture governing soil processes.

In general, it is difficult to attribute processes (soil inputs vs turnover) to changes in soil carbon for two main reasons. First, soil model structure leads to false-priming effects[20]. The apparent reduction in soil turnover is a consequence of increased surface soil carbon relative to deep soil, as opposed to a change in actual turnover rates. Second, our process attribution methodology treats the soil as a single pool. Therefore, the calculation of () may include deep, inactive (on decadal timescales) soil carbon, and underestimate the turnover rate of the active soil. Hence, our process attribution of soil carbon changes is not well constrained. For example, we estimate a 300 PgC increase in soil carbon from increased inputs for CLM5.0, a result which is in part driven by the large carbon stocks in high northern latitudes leading to a large estimate of [56]. We do account for this bias by scaling output to actual changes in modelled carbon stocks (see methodology), but our results highlight the urgent need for evaluation of individual carbon pools and fluxes (leaf, root, wood and each soil pool) separately, and not bulk terms.

In addition to the CO2-driven sink, evidence exists for a substantial impact of forest regrowth on the northern sink over the previous 60 years[18,19]. DGVMs have some capacity to capture regrowth effects (predominantly following land abandonment), however, most do not simulate forest demography or detailed forest management (e.g. enhanced stock densities[57]) and so could underestimate the actual regrowth carbon sink[58]. Further, not all DGVMs simulate disturbance in unmanaged land (e.g. wind and pests), and therefore, natural disturbances beyond fire and the subsequent regrowth will not be captured by the DGVMs here[18]. Moreover, not all models simulate fire mortality, and fire models do not always capture observed spatiotemporal patterns and associated carbon emissions[59], as well as large uncertainties due to issues with forcing data[60].

In addition to direct human-caused carbon losses, alternative drivers of mortality can reduce carbon stocks. For example, climate-induced mortality, combined with a shorter tree lifespan resulting from faster growth[61,62], currently weaken the Amazon forest sink[63]. However, DGVMs show no sign of increased carbon turnover due to changes in climate or faster growth, although this discrepancy is not unexpected as DGVMs do not include detailed mortality processes, e.g. drought-mortality[64], although this is an area of current model development[65,66]. The DGVMs do indicate a decline in tropical productivity due to climate change, in line with remote-sensing and upscaled in-situ observations[67], indicating current temperatures may exceed tropical plant thresholds[24].

Here we have shown that the long-term evolution of the global net land sink is well estimated by an ensemble of state-of-the-art DGVMs. However, a complete process-based understanding is still lacking due to several model shortcomings. Specifically, above/below-ground partitioning of the sink is highly uncertain due to too simplistic representation of internal carbon cycling. Model improvement into plant allocation, tissue lifespan and mortality, as well as the inclusion of process-based soil carbon and nutrient cycling, should be a high priority moving forward.

Methods

TRENDYv10 models

We analyse output from 18 DGVMs that are part of a recent model intercomparison project, TRENDYv10 and the Global Carbon Budget 2021 (GCB2021)[8]. The models included in the analysis here are CABLE-POP, CLASSIC, CLASSIC-N, CLM5.0, DLEM, IBIS, ISAM, ISBA-CTRIP, JSBACH, LPJ-GUESS, LPJ, LPX-Bern, OCN, ORCHIDEE, ORCHIDEEv3, SDGVM, VISIT, and YIBs (see ref. 8 for model descriptions and setup). One model (JULES-ES-1.1) from TRENDYv10 is not included due to incomplete data. Note, CLASSIC-N[68] and ORCHIDEE[69] were not part of GCB2021 due to the inclusion of alternate model versions (CLASSIC and ORCHIDEEv3) but are included in TRENDYv10 and this study.

The models are forced with a merged monthly Climate Research Unit (CRU)[70] and 6-hourly Japanese 55-year Reanalysis (JRA-55)[71] data set. The models are also forced with atmospheric CO2[72], gridded nitrogen deposition[73] and nitrogen fertiliser[74]. DGVMs use the HYDE (v3.3) land-use change data set[75,76], which provides annual pasture and cropland areas at a global scale, and includes improvements in the spatial distribution of agricultural regions[77]. Several models (CABLE-POP, CLM5.0, JSBACH, LPJ-GUESS, LPJ, and VISIT) also use harmonised land-use change data (LUH2-GCB2021[8]), which provides information on sub-grid-scale land-use transitions.

To isolate the response of the land to each driver (CO2, climate, LULCC), each model performs four simulations: S0 (fixed pre-industrial atmospheric CO2 and land-use, recycled 1901–1920 climate), S1 (transient atmospheric CO2, recycled 1901–1920 climate, and fixed pre-industrial land-use), S2 (transient atmospheric CO2, transient climate, and fixed pre-industrial land-use), and S3 (transient atmospheric CO2, transient climate, and transient industrial land-use). Nitrogen deposition varies temporally in simulations S1–S3. Therefore, the transient CO2 (+N deposition) effect on the ‘natural’ land sink is calculated by S1−S0, S2−S1+S0 is the climate effect on the ‘natural’ land sink, S3−S2 is the LULCC effect, and S3 is the net effect. There exists an artefact in the HYDE3.3 data causing a large land-use transition and emission peak around 1960. To correct this, we replace the LULCC estimates for 1959–1961 with the average of 1958 and 1962 in each DGVM.

Net Biome Productivity (NBP) from the S3 simulation represents the net land sink and can therefore be compared to the ‘observed’ Global Carbon Budget land constraint, which is calculated in the GCB as fossil fuel emissions—atmospheric carbon growth rate—ocean carbon sink (see Fig. 1 in the main text and Table 5 in ref. 8 and data taken from 10.18160/gcp-2021). Atmospheric CO2 concentration measurements began in 1959[78], and so this independent constraint on the land sink covers the period 1959–2020, which defines the study period used throughout our analysis.

Data processing

First, we calculate for each model the global and regional (two regions: north of 30°N and south of 30°N) annual mean values for NBP, net primary productivity (NPP), heterotrophic respiration (), vegetation carbon () and soil carbon (), which includes litter () and coarse woody debris () for 1959–2020. We exclude the product pool from our analysis as we do not have carbon emission data from wood products available.

Process attribution framework

To attribute changes in C pools to different processes, we quantify C pool changes attributable to changes in productivity and changes in turnover time. We do this using an analytical approximation of C pool dynamics corrected for non-steady state and more complicated behaviour of the models. First, we approximate the steady-state C pool size given the inputs and turnover time of a given pool in a given year. We then quantify the change in the steady-state pool sizes estimated for 1959 and the year in question and partition this change proportionally among inputs and turnover time. Finally, to adjust to non-steady state conditions and more complicated model structures (i.e. multiple pool dynamics, not a single pool as in this approximation), we adjust the attributed proportions by adding the difference in actual modelled change in carbon stock and the steady-state approximation, scaled by the relative magnitude of productivity and turnover-driven changes.

The change in carbon pools can be expressed as:where we define as the change in carbon pool per year. To quantify the internal processes impacting the output flux, we can define , where is the mean turnover time (the inverse of the mean turnover rate) of the carbon pool, C. Therefore, we can rewrite Eq. 1 as:

We separate the land into two main pools, vegetation () and soil (). is a direct model output and includes leaf, root and wood biomass. Similarly, , are direct model outputs, which we sum together for our analysis. The rest of the methodology, refers to the sum of these three pools. Not all models simulate litter and CWD pools. ISAM does not simulate litter and only CABLE-POP, CLM5.0, DLEM and LPJ-GUESS simulate CWD. Further, for CLM5.0, the litter pool is incorporated in the soil pool output, and so is also excluded from the sum to avoid double counting.

We can now specify Eq. 2 for vegetation and soil as:where is the carbon flux from vegetation to soil, from either litterfall, mortality or direct transfer from roots to soil.

We calculate annual (using Eq. 3), and using model output as follows:where and are the annual changes in vegetation and soil carbon. We want to be able to define in terms of input fluxes and turnover times only, as these are the two aggregated mechanisms driving changes in the pools. To do this, we use the fact that the pools and the input and output fluxes are generally much larger than the change in carbon pools each year to make the assumption that . Rearranging Eqs. 3 and 4, we can estimate the steady-state carbon pools in any particular year as:where and are steady-state approximations of the actual carbon pools and . We can now define the change in carbon pools over our 60-year study period, 1959–2020, in terms of inputs and turnover time:

In order to separate the impacts of changes in inputs and turnover times, we next define , , and . We now substitute these into Eqs. 7 and 8:

The right-hand side of the equation contains three terms corresponding to the change in carbon storage due to changes in inputs (NPP for vegetation and litterfall for soil), outputs (bulk turnover), and an interaction term. Our approach, therefore, extends the factor separation approach by ref. 21, which applied it to attribute simulated soil carbon changes into input-driven, turnover-driven change and a synergy term to also cover vegetation carbon changes.

The change in vegetation carbon due to each term can be written as: , , and  +  + , with equivalent terms for soil carbon.

Finally, we adjust the one-pool steady-state estimates for the difference in the model simulations that arise due to the combination of steady-state assumption (Eqs. 5–8), the grouping of all pools into single vegetation and soil pools (as the distribution between component pools changes over time), and the linearisation not capturing the full dynamics of land carbon cycling (Eqs. 9 and 10). Therefore, the estimates of changes in carbon pools ( and ) do not equal the actual change in C pools ( and ), which are a direct model output (see Supplementary Fig. 9).

To address this issue, we calculate for each model the difference between the actual change in carbon stocks and our approximation ( and ). We make the simple assumption that productivity and turnover are proportionally responsible for the mismatch between steady state and the actual carbon stock, and so adjust each term on the right-hand side of Eqs. 9 and 10 by the relative magnitude of each term multiplied by or for vegetation and soil, respectively. For example, is adjusted by adding:

We calculate Eqs. 9 and 10 for each of the eighteen models, apply the appropriate scaling and use them to produce Fig. 4.

Uncertainty analysis

The spread of model estimates of changes in carbon stocks and fluxes is used to determine the uncertainty of our results. For land flux estimates (NBP) and and , we show ‘mean ± std. dev.’. We decompose the uncertainty (defined as the model spread) in and into four components: baseline input ( or ), change in inputs ( or ), baseline turnover ( or ), and change in turnover ( or ). For each of these terms in Eqs. 9 and 10, we calculate the standard deviation in and using the multi-model mean values of all other terms in the equations and the individual model values for that term.

Supplementary Information File