Tracking 21st century anthropogenic and natural carbon fluxes through model-data integration

External datasets

Woody biomass carbon data

The dataset by ref. 16 maps annual global woody biomass carbon densities for 2000–2019 at a spatial resolution of ~10 km. The annual estimates represent averages for the tropical regions and growing-season (April–October) averages for the extra-tropical regions. Ref. 16 analyse global trends of gains and losses in woody biomass carbon for 2000–2019. Overall, they find that grid cells with (significant) net gains of vegetation carbon are by a factor of 1.4 more abundant than grid cells with net losses of vegetation carbon, indicating that there is a global greening trend when only considering the areal extent of biomass gains and not the magnitude of carbon gains. Their regionally distinct analysis of trends shows that almost all regions, except for the tropical moist forests in South America and parts of Southeast Asia, experienced net gains in biomass carbon. On the country scale, the largest net increase in biomass carbon is shown in China, which is mainly attributed to the large-scale afforestation programs in the southern part of the country and increased carbon uptake of established forests. On the other hand, the largest vegetation carbon losses are shown for Brazil and Indonesia, which is partly attributed to deforestation, degradation, and drought events. All of the mentioned trends have been found to be significant¹⁶. The decreasing carbon sink in Brazil is in line with ref. 44, who, considering both natural and anthropogenic fluxes, show that the southeastern Amazon has even turned from a carbon sink to a carbon source, mainly owing to fire emissions from forest clearing. Isolating carbon fluxes in intact, old-growth Amazonian rainforests (i.e., S_LAND,B), ref. 45 also find evidence for a significantly decreasing carbon sink due to the negative effects of increasing temperatures and droughts on carbon uptake since the 1990s.

The dataset was remapped to the BLUE resolution of 0.25^∘ through conservative remapping (i.e., area-weighted averaging).

ERA-5 data

The ERA-5 variables were downloaded from the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/cdsapp#!/home). Monthly air temperature (T_a) at 2 m height was averaged over each year, and annual precipitation was calculated by taking the sum of the monthly total precipitation (P). Both variables were regridded from the original resolution of ~0.1° to 0.25° resp. to the TRENDY resolution of 0.5° through conservative remapping.

TRENDY data

We used the TRENDY model ensemble version 8 (conducted for the 2019 GCB; ref. 8). We used net biome production (NBP) and annual vegetation carbon stocks (cVeg) for 2000–2018 from four different model setups (S2, S3, S5, and S6) and eight resp. 13 DGVMs (depending on the data available). The selection of DGVMs is done as in ref. 19 (Supplementary Tab. 3), but we included one additional model (ISAM) for the S2 simulations. The terrestrial biomass carbon sink (S_LAND,B) was calculated for 13 DGVMs following the GCB 2020 approach, i.e., from the S2 simulation, which is the simulation without LULCC (i.e., fixed pre-industrial land cover) under transient environmental conditions (climate, nitrogen deposition, CO₂ evolution). S_LAND,B is the annual difference of cVeg and makes no statements about the further fate of biomass if cVeg decreases. S_LAND,B, therefore, should not be interpreted as equivalent to the flux to/from the atmosphere, since parts of cVeg may be transferred to litter, dead wood, or soil. The same applies to our BLUE estimates of S_LAND,B, ensuring comparability between our BLUE estimates and the TRENDY estimates. Increases (decreases) of cVeg between two years are a net uptake (release) of carbon from the terrestrial biosphere. The global sums of biomass carbon stocks under transient climate and CO₂ were calculated from the S3 setup (LULCC under historical environmental conditions), whereas the S5 setup provides biomass carbon under constant present-day environmental forcing (closest to the classical bookkeeping approach). In line with the GCB, E_LUC was calculated under historical environmental conditions as the difference in NBP between the S2 and S3 simulations (E_LUC = NBP_S2 – NBP_S3). E_LUC under constant present-day environmental forcing was calculated as the difference in NBP between the S6 (fixed pre-industrial land cover under present-day environmental forcing) and S5 simulations (E_LUC = NBP_S6 – NBP_S5)¹⁹. All datasets were remapped to a common resolution of 0.5^∘ through conservative remapping (area-weighted average) for the data analysis.

Assimilation of observed woody biomass carbon in BLUE

The observed woody biomass carbon densities by ref. 16 are assimilated in BLUE in several steps.

Carbon transfer in the default setup of BLUE

The BLUE simulation is started in AD 850. Biomass and soil vegetation carbon densities are based on ref. 17 (see ref. 5 for details). These carbon densities are specific for eleven natural PFTs (Supplementary Fig. 3), which are assigned one of four land cover types (primary land, secondary land, cropland or pasture). The LULCC forcing is based on the LUH2 dataset¹⁸, defining the vegetated fractional area of each grid cell that is affected by a land-use transition. Each transition may lead to a change from one land cover type (=source land cover type: j) to another land cover type (=target land cover type: (j^{prime})). In the case of wood harvesting on secondary land, (j=j^{prime}), whereas all other transition types (e.g., clearing for agricultural expansion, abandonment of agricultural lands) induce a change in land cover. The fractional grid cell areas undergoing transitions are further distributed across PFTs proportionally to the temporally constant PFT area fractions (Supplementary Fig. 3). Upon each land-use transition, biomass carbon is transferred between the source land cover type and the target land cover type, whereby the amount of transferred carbon depends on the biomass carbon density of the source and target land cover types (in the respective PFT) and the area affected by the transition (in the respective PFT). The temporal evolution of the biomass carbon pool after any type of land-use transition is approximated by an exponential function with different time constants for decay and regrowth, depending on the type of land-use transition. The time constants are based on linear estimates by ref. 17, which are converted to exponential time constants. A detailed explanation of the exponential model can be found in ref. 5.

While in the default setup, changes are only due to LULCC, our assimilation approach now introduces environmental effects on woody vegetation carbon by assimilating the observed woody biomass carbon densities in BLUE from 2000 onward according to the methodological considerations explained below.

Calculation of woody biomass carbon densities for different land cover types and PFTs

Within each 0.25° cell of the global grid, the (remapped) woody biomass carbon density from ref. 16 must be the sum of woody biomass carbon stored in all woody PFTs of all woody land cover types. The distribution of the woody biomass carbon across PFTs and land cover types is achieved by distributing the observed (i.e., actual) woody biomass carbon densities (ρ_Ba) from ref. 16 across the two land cover types (j) and the eight PFTs (l) that can be woody vegetation (primary land, called virgin, “v” in BLUE and secondary, “s”, land) according to the fraction of total woody biomass carbon (f_B) contained in each land cover type and each PFT (f_B,j,l) as estimated by BLUE. f_B,j,l varies for different PFTs and land cover types, depending on their history of LULCC and their potential for carbon uptake (i.e., the potential carbon densities).

f_B,j,l is extracted from the default simulations for the first year of the time series (i.e., 2000) and calculated for subsequent years from the BLUE simulations using the assimilated woody vegetation carbon densities for that year:

where C_B is the woody biomass carbon stock.

Consequently, the assimilated woody biomass carbon stock per cover type and PFT (C_{B_as,j,l}) at each time step can be calculated as:

with j{v, s}; l{1. . 8}; t{2000. . 2019}. A is the area per grid cell.

Thresholds for excluding inconsistent woody biomass carbon densities

We eliminate unrealistically large values for woody biomass carbon densities that our assimilation framework produces. Woody biomass carbon densities in BLUE that exceed the highest value (~374 t ha⁻¹) of the original dataset indicate inconsistencies between the observed woody biomass carbon estimates and the fractional grid cell areas per PFT and land cover types that BLUE simulates. To account for uncertainties related to the criteria for exclusion of grid cells, multiple threshold approaches are applied and the results are compared. To maintain a temporally and spatially consistent time series of woody biomass carbon, grid cells that are excluded according to the chosen threshold approach are interpolated through linear barycentric interpolation. A first approach relies on a uniform upper threshold of <375 t ha⁻¹ for woody biomass carbon densities. This approach leads to the exclusion of ~3% of all grid cells, but is considered conservative in the sense that it may lead to an overestimation of woody biomass carbon densities of non-forested land, since it is expected that the maximum value of ~374 t ha⁻¹ occurs in heavily forested grid cells only. To account for this potential overestimation, additional threshold approaches are applied by cutting the distribution of grid cells with woody biomass carbon densities smaller than 375 t ha⁻¹ to a range of specific percentiles and choosing the values corresponding to each percentile as upper thresholds for the exclusion of further grid cells. In the first step, we choose the 97th, 98th, and 99th percentiles and evaluate the resulting dynamics of total vegetation carbon in terms of their agreement with the original dataset by ref. 16. The evaluation is done by analyzing the results from each percentile threshold approach in terms of the global dynamics of the biomass carbon stocks in comparison to the estimates from ref. 16 (see Supplementary material). This analysis reveals that the annual dynamics (i.e., increase/decrease) of the woody biomass carbon stocks start to diverge strongly from the original time series for thresholds smaller than the 99th percentile. This is related to an enhanced loss of spatial and temporal variability of the assimilated biomass carbon stocks due to an increased number of interpolated grid cells with smaller percentile thresholds. Consequently, we choose the two approaches with (1) <375 t ha⁻¹ and (2) <99th percentile of 375 t ha⁻¹ as upper limits for the exclusion of inconsistent biomass carbon densities and use their average, unless indicated otherwise. Both threshold approaches are applied to each woody PFT and the two woody land cover types separately over the whole time series (2000–2019). Consequently, it is ensured that differences in carbon storage potential between different PFTs and land cover types are considered within the percentile threshold approach.

Model initialization

In our transient woody biomass carbon approach, we need to initialize the woody biomass pools in BLUE at each time step (i.e., each year) to account for changes in biomass carbon densities due to environmental processes. As we do not assimilate soil carbon densities in our approach, the soil carbon pools are initialized once at the beginning of the BLUE simulations (described below) and subsequently only altered by LULCC. The re-initialization for the woody biomass pools at each time step is necessary, as BLUE only explicitly simulates annual changes in biomass carbon densities due to LULCC. In the default approach, total biomass carbon is partitioned between equilibrium pools (({bar{C}}_{B,j,k,l})) and excess pools (δ_B,j,k,l). The former mark the carbon stock that the biomass pools strive to reach (i.e., the carbon stock that the PFT and land cover type would reach after a sufficiently long time after a land-use disturbance), while the latter indicate whether the current biomass carbon stock is in equilibrium (δ_B,j,k,l = 0) or in excess of equilibrium (δ_B,j,k,l ≠ 0) (Note that k is the (land-use) history type, including clearing (“l”), harvest (“h”), abandonment (“a”), other (“g”)). An in-depth explanation of the different pool types in BLUE can be found in the original documentation by ref. 5. Biomass carbon is assumed to be in equilibrium upon model initialization, i.e., the equilibrium pools contain all biomass carbon and the excess pools are zero. Upon each land-use transition, the equilibrium and excess biomass carbon pools are altered, depending on the transition type.

In our approach, the model initialization is done by distributing the assimilated woody biomass carbon among the equilibrium biomass pools for all woody PFTs and all land cover types (see Supplementary materials for the handling of non-woody land cover types). This means that the equilibrium biomass pools and all excess biomass pools are then re-initialized at each time step (annually) of the simulation from 2000 onward. The excess carbon pools are changed upon each land-use transition, whereby the spatially explicit actual woody biomass carbon densities derived from ref. 16 replace the woody biomass carbon densities based on ref. 17 from 2000 onward. The actual woody biomass carbon densities from ref. 16 are assimilated in BLUE at the beginning of each year X and subsequently altered by the land-use transitions in year X. Consequently, the BLUE output of carbon stocks for year X represents the end of year X and changes in carbon stocks between year X+1 and year X are attributed to year X+1. Legacy fluxes (i.e., carbon fluxes from land-use that do not occur in the same time step as the corresponding land-use event) are tracked according to the approach explained below.

Handling of legacy carbon fluxes

Due to repeated initialization of the (equilibrium and excess) biomass carbon pools at each time step, legacy fluxes are not accounted for and need to be tracked separately. Such legacy fluxes from/to the atmosphere to/from the terrestrial woody biomass occur due to LULCC prior to the current time step, e.g., because the forest regrows slowly or because cleared biomass decomposes slowly on site or in products. We track these legacy fluxes separately for those from the LULCC prior to the assimilation period (2000–2019), and those occurring during the assimilation period, which are caused by the LULCC transitions and other biomass changes. To track the former, we introduce an additional set of excess pools δ_B,leg<2000 that include all excess woody biomass carbon from land-use transitions prior to 2000 upon initialization of the actual woody biomass carbon pools in 2000. Legacy carbon fluxes from land-use transitions prior to 2000 (θ_{B,leg<2000,j,k,l}) are calculated as in the default approach (Note: carbon fluxes from LULCC in time step (t) from/to the land to/from the atmosphere are realized at the beginning of time step (t+1) in BLUE):

with j{v, s, p, c}; l{1. . 8}; t{2000. . 2019}; k{l,h,a,g}. τ is the time constant for relaxation processes, which varies for different pool types (biomass or soil), land cover types, and PFTs.

Excess woody biomass carbon from transitions from 2000 onward is tracked in another set of separate pools (δ_B,leg≥2000) to account for ≥2000 legacy fluxes. δ_B,leg≥2000 is adjusted at the beginning of each time step for all excess woody biomass carbon from the previous time step minus fluxes to/from the atmosphere (Eq. (4a)) (θ_{B,leg≥2000,j,k,l}) from relaxation processes in the respective time step (Eq. (4b)):

Carbon fluxes between the terrestrial woody biomass pool and the atmosphere pool at each time step, including all legacy fluxes (θ_B,j,k,l), can then be calculated as the sum of instantaneous carbon fluxes at the current time step (resulting from LULCC in the previous time step), legacy carbon fluxes prior to 2000 and legacy carbon fluxes from 2000 onward, but prior to the current time step (resulting from LULCC prior to the previous time step).

Derivation of the terrestrial woody biomass carbon sink

To isolate anthropogenic from environmental (=S_LAND,B) carbon fluxes from woody vegetation, we performed two simulation setups based on different approaches for assimilating the observed woody biomass carbon densities. The biomass estimate by ref. 16 includes carbon stored in living woody vegetation (trees and shrubs), whereas carbon stored in dead plant material (litter, harvested wood products) is not included in the estimate. Consequently, the change in woody biomass carbon stocks within a certain time step results from carbon sources and sinks driven by LULCC in the respective time step, from carbon sinks due to regrowth of vegetation driven by past LULCC (i.e., prior to the respective time step) and from all environmental processes (on managed and unmanaged lands) in the respective time step:

where ΔC_source resp. ΔC_sink are sources resp. sinks of biomass carbon due to LULCC in the current time step, ΔC_{reg_leg} are sinks of biomass carbon due to regrowth of vegetation from LULCC prior to the current time step and ΔC_source,env resp. ΔC_sink,env are sources resp. sinks of biomass carbon due to environmental processes in the current time step. We performed additional BLUE simulations with fixed (i.e., stationary in time) woody biomass carbon densities to split the carbon fluxes from woody vegetation into the anthropogenic and environmental terms of Eq. (6). The fixed woody biomass carbon setup is based on the 2000 estimates derived from ref. 16. As in the transient simulations, model initialization is done from the same state of woody biomass carbon in 2000 (i.e., anthropogenic and environmental effects on woody biomass carbon prior to 2000 are implicitly captured), but changes in woody biomass carbon in the subsequent years are only driven by LULCC in the fixed setup. In the fixed woody biomass carbon simulations, there is no need for separately tracking legacy fluxes from 2000 onward, since the biomass carbon pools are only initialized in 2000 and the excess pools are altered subsequently without re-initialization. Legacy carbon fluxes from land-use transitions prior to 2000 are considered following the same approach as in the transient woody biomass carbon density setup. The terms of Eq. (6) that capture environmental changes in biomass carbon stocks are only included in the BLUE simulations with transient biomass carbon, whereas biomass carbon changes driven by land-use change are captured in both the transient and fixed biomass carbon simulations. Consequently, our BLUE simulations allow us to isolate all environmental effects on woody biomass carbon by taking the difference in woody biomass carbon stocks between the two BLUE simulation setups:

This term, the natural carbon sink in terrestrial woody vegetation, represents the net effect of environmental processes on managed and unmanaged lands on the terrestrial woody living vegetation.

Uncertainties

The main sources of uncertainty that affect the results from our biomass assimilation approach are (1) model assumptions regarding the global LULCC dynamics and the rates of vegetation regrowth, (2) potential misattributions of anthropogenic fluxes as natural fluxes owing to incomplete data on LULCC, and (3) uncertainties within the original time series of woody biomass carbon densities by ref. 16. We analyse the different sources of uncertainty as described in the following.

(1) The difference between the observed woody vegetation carbon stocks from ref. 16 and the woody vegetation carbon stocks at the beginning of each time step in the transient BLUE setup (=“assimilated woody biomass carbon”) can be used to evaluate the LULCC forcing and the PFT distribution in BLUE. Since the observed and the assimilated woody vegetation carbon time series are not independent of each other, the comparison solely aims at identifying potential model uncertainties. The observed woody vegetation carbon densities are assimilated into BLUE at each time step according to the spatial distribution of the land cover types and PFTs (see Eq. (2)). Consequently, a larger difference between the observed woody vegetation carbon stocks and the assimilated woody vegetation carbon stocks would indicate that the actual LULCC dynamics and/or the spatial distribution of PFTs are not captured well in BLUE. We call this difference the “biomass assimilation bias”. The average global biomass assimilation bias (±1 SD) amounts to 29 ± 4 PgC between 2000 and 2019. The agreement between the observational dataset by ref. 16 and the assimilated woody biomass carbon in terms of the trends in global biomass carbon is quantified as the number of years that show the same trend in both datasets related to the previous year divided by the total number of years. Following this definition, a temporal agreement of 100% would mean that the observed dataset and the assimilated dataset show the same trend in biomass carbon for all years. The regionally averaged agreement in the estimated biomass carbon trends is generally high (>80%) (Supplementary Fig. 8) but smaller in regions with strong LULCC dynamics (tropics and Europe). Some local hotspots exist (Fig. 3), where differences between the observed dataset and the assimilated dataset are larger. These are mainly located in South- and Southeast Asia, Europe, and Equatorial Africa (Fig. 3a), where clearing and wood harvesting rates of the forest as prescribed in the LULCC forcing (Fig. 3b) are very high, leading to much lower biomass carbon estimates in our assimilated woody biomass carbon estimates than in the observed time series. This suggests that the clearing and/or wood harvesting rates are overestimated in the LULCC forcing and/or the rate of vegetation regrowth is underestimated in BLUE, leading to a high biomass assimilation bias for the mentioned regions, which further affects E_LUC (see Results).

We further assessed the validity of our S_LAND,B estimates in terms of the high IAV shown in Canada, Russia, Brazil, and Europe. The comparison of our time series of assimilated woody biomass carbon to the original time series by ref. 16 shows that our assimilated dataset is very close to the original dataset in the respective regions and that the high IAV is also shown in the original time series¹⁶. Consequently, we conclude that the high IAV is not introduced by uncertainties in our model-data integration. Nevertheless, we acknowledge that the estimated IAV in NAM (especially Canada) may seem high compared to the IAV of carbon fluxes due to natural disturbances estimated by the National Inventory Report of Canada⁴⁶ or to specific disturbance events, such as the mountain pine beetle outbreak in British Columbia in the early 2000s⁴⁷. However, our estimated IAV of S_LAND,B of up to 2 PgC yr⁻¹ in NAM is not inconsistent with previous studies estimating annual changes of the total land sink⁴⁸ resp. the natural land sink^49,50 of up to 3 PgC yr⁻¹, including a switch in sign of the flux. Furthermore, ref. 51 combine atmospheric CO₂ measurements with inverse modelling and show that the North American net ecosystem exchange (NEE) between 2007 and 2015 and its variability was strongly driven by el Niño (more than average carbon uptake) and la Niña (less than average carbon uptake) conditions, with monthly anomalies (related to the mean for 2007–2015) of up to ±1.5 PgC yr⁻¹. They further find that the boreal coniferous forest is among the ecosystems with the largest difference in NEE anomalies between el Niño and la Niña periods. Our S_LAND,B estimates broadly follow the dynamics described by ref. 51, with higher than average carbon uptakes during el Niño conditions in 2010 and 2015 and lower than average carbon uptakes resp. carbon releases during la Niña conditions in 2011. Furthermore, we wish to clarify that we assume that the magnitude of fluxes to/from the atmosphere from/to the biomass is smaller than our estimated changes in S_LAND,B, since we include depositions to dead matter and soil in our estimates. According to the global carbon budget decomposition by ref. 49, around 50% of annual gross anthropogenic emissions (including natural fluxes on managed lands, as described in ref. 52) between 2000 and 2015 was due to direct emissions (i.e. in the year of the respective disturbance), while the rest was attributable to legacy fluxes (onsite decomposition and wood product degradation). Assuming a similar dynamic for (solely) natural disturbances, we expect that the magnitude of annual carbon fluxes to/from the atmosphere from/to the biosphere is lower – depending on the degree and type of disturbance – than the annual changes in S_LAND,B. The mentioned considerations highlight the need for future, independent estimates – especially on a regional scale – to foster our understanding of the IAV of terrestrial carbon fluxes. Furthermore, ref. 16 compared their estimated IAV of woody biomass carbon with FAO estimates of forest carbon, showing that there is no systematic overestimation of IAV in the mentioned regions.

(2) A general shortcoming of all accounting approaches based on observational datasets is their inability to capture all anthropogenic activities related to LULCC. Consequently, there might be anthropogenic carbon fluxes that are classified as environmental fluxes in our approach, simply because the LULCC forcing is not capturing the underlying anthropogenic activities completely. This caveat applies foremost to certain types of anthropogenic degradation: while the LULCC forcing covers logging (implemented as wood harvesting in BLUE) and rangeland degradation, it does not account for degradation caused by anthropogenic fires, which might lead to misattributions of the related fluxes towards S_LAND,B. However, there is currently no (global) dataset available that separates anthropogenic and natural degradation fires and that would allow us to provide an uncertainty estimate for the misattributed fluxes.

(3) Reference 16 defines errors in the order of ±0.5% (2 PgC) related to the 2000–2019 average global sum of carbon contained in woody vegetation (381 PgC). The error estimate includes pixel-level uncertainty and modeling uncertainty from parameter estimation. The global uncertainty range of ±0.5% is considered in all of our aggregated global estimates of woody biomass carbon. This means that in Table 1, the uncertainty range of ±2 PgC for the global living biomass (i.e., woody plus herbaceous vegetation) carbon stocks refers to the woody vegetation estimate only (357 PgC).

Global maps of absolute and relative pixel-level uncertainty (Supplementary Fig. 9) are provided and can be used as a reference to evaluate the accuracy of our estimates for different regions.

Source: Ecology - nature.com