Climate policies for carbon neutrality should not rely on the uncertain increase of carbon stocks in existing forests

Global climate change puts increasing pressure on societies and ecosystems. The international community, through treaties such as the Paris agreement, has committed to keeping global warming well below 2 °C in order to avoid the worst consequences of climate change. To achieve this ambitious target, a combination of carbon emission reductions and strong carbon sequestration is necessary (Pathak et al 2022). Most of the carbon sequestration currently occurs in the oceans and forests (Friedlingstein et al 2022, Sha et al 2022), and the strategies laid out to reach net carbon neutrality in the coming decades depend, to a large extent, on these natural systems to persist as strong carbon sinks (e.g. European Commission 2020).

On the long term, forests remaining an active carbon sink is by no means guaranteed (Baccini et al 2017). At larger spatial and temporal scales in a stable climate, forests tend towards equilibrium, where carbon sequestration by photosynthesis is balanced out by respiration, tree decay and mortality (Looman 1976, DeAngelis and Waterhouse 1987). A changing climate can alter this balance in several ways. First, tree growth conditions depend on both the temperature regime and the water available to trees. Second, the changing climate also impacts the likelihood of occurrence and severity of some natural disturbances, such as windthrows, insect outbreaks and forest fires (McDowell et al 2020). Third, increased CO₂ levels in the atmosphere directly enhance the rate of tree growth (Chen et al 2022), as CO₂ is a primary driver of photosynthesis. This mechanism, defined as CO₂ fertilization (CFE), has contributed to the increase of the terrestrial sink in the past decades (Friedlingstein et al 2022). However, CO₂ is only one of the environmental factors affecting tree growth. Forests can only benefit from higher CO₂ concentrations as long as other drivers like water, energy and nutrients are not limiting primary productivity and adverse environmental conditions remain absent (Cabon et al 2022, Green and Keenan 2022). At global scale, evidence suggests that primary productivity is currently shifting from a period of growth enhanced by CFE, to a condition where climate change and climate related disturbances are increasingly limiting productivity (Peñuelas et al 2017).

To evaluate the long-term persistence of the forest carbon sink, the changes in the carbon stock in forests at equilibrium under changing climate conditions need to be understood. To model this, dynamic global vegetation models (DGVMs) have been deployed in conjunction with global atmosphere and ocean models within the frame of global Earth System Models (ESMs). However, as recently shown (Terrer et al 2019, Cabon et al 2022), these models likely overestimate the strength of the CO₂ fertilization effect due to an incomplete representation of the other limiting factors (such as the availability of nitrogen and phosphorus for growth) and would consequentially provide an overoptimistic view of the long-term persistence of the forest carbon sink. Proxies for primary productivity derived from satellite data suggest a steeper declining trend of CO₂ fertilization than that predicted by DGVMs (Wang et al 2020). As an alternative to ESMs, attempts have been made to predict future forest carbon storage using data driven statistical models. However, these modelling attempts have so far not included the CO₂ fertilization effect (Bastin et al 2019, Walker et al 2022), which would likely lead to an underestimation of the amount of carbon that can be stored in forests.

In this work we aim to quantify the possible persistence of the forest carbon sink by comparing potential biomass under current and future climate conditions. The objective is to better approximate changes in forest growth potential and disturbance regimes - and their respective uncertainties -, as well as to quantify the effect of diverging estimates of CO₂ fertilization on the global forest carbon pool. This is done with a novel data-driven modelling framework that separately parameterizes and compares the CO₂ fertilization effect based on ESM data, satellite proxies and free-air CO₂ enrichment (FACE) experiments.

We establish the long-term persistence of the global forest carbon sink by estimating the above-ground biomass carrying capacity (BCC) under present and future climatic conditions (Roebroek 2024). BCC is the biomass of a forest 'under prevailing environmental conditions and natural disturbance regimes, but excluding anthropogenic disturbance' (Keith et al 2009, Roebroek et al 2023). Under the assumption that human influences (e.g. harvesting, planting, fire suppression etc) on forests remain equal over time, the expected change in BCC under a given climate projection is equal to the change of the actual biomass under that climate scenario. In doing so, we can separate the biological constraints on biomass trajectories in the future from changes in human pressure on forests such as those related to land use change.

The estimates of BCC are made with a modular machine learning model framework. It is based on a simplified model definition, expanding on previous work (Roebroek et al 2023). It describes changes in BCC as a product of growth potential (B_p , potential biomass), natural disturbance regime ($\left( {1 - D} \right)$, fraction of potential biomass that is expected to be sustainable under given climate conditions) and CO₂ fertilization (${F_{{\text{C}}{{\text{O}}_{\text{2}}}}}$, a factor adapting potential biomass and the disturbance regime to enhanced CO₂ scenarios, systematically larger than 1):

$$\begin{align}{\text{BCC}} = {B_p} \cdot \left( {1 - D} \right) \cdot {F_{{\text{C}}{{\text{O}}_{\text{2}}}}}\end{align} \tag{ 1 }$$

Potential biomass (${B_p}$) is derived from satellite-based estimates of above-ground biomass (Santoro 2018) by estimating the maximum biomass achievable under given environmental conditions (climate, soil and topographical features, see table 1) using a quantile regression forest (QRF) to estimate the 100th percentile (implementation in the python package sklearn-quantile) (Meinshausen et al 2020). In doing so, potential biomass estimates approximate the hypothetical potential, and will thus be referred to as potential realized biomass (PRB). The disturbance regime (D) is modelled as the average relative deviation of biomass in current intact forests (Potapov et al 2017) to the PRB (1—biomass of intact forests/PRB) with random forest regression (as implemented in the scikit-learn python package) under the same environmental conditions as used to model the PRB. Both the QRF and regular random forest regressor have 50 trees and share their parameterization. Hyperparameters were chosen with a grid-based 5-fold cross validation approach, minimizing the pinball loss (with alpha set to 0.995) for the predictions of PRB. This procedure is described in more detail in Roebroek et al (2023) and visualized in figure 1. The quality of the estimates of the natural disturbance regime is dependent on the intact forest data. Excluding the error in mapping untouched forests, uncertainty predominantly comes from the fact that all global forests have indirectly been altered to some extent by human actions, such as by increased likelihood of drought and increasing global temperatures. These indirect effects have not been included in the current analysis.

Zoom In Zoom Out Reset image size

The original model description in Roebroek et al (2023) was purely based on current environmental conditions, for which various predictors are available. In order to transpose the model framework to future climate scenarios, the set of environmental predictors had to be adapted. Soil water availability and some soil characteristics deemed most dependent on the state of the climate (i.e. nitrogen, pH, organic carbon and cation exchange capacity) were excluded from the analyses as no comparable datasets exist for future climate conditions. The remaining soil characteristics in the SoilGrids V2 (Poggio et al 2021) describe soil physical properties (fraction sand/silt/clay and bulk density) and are assumed to remain constant to the end of this century. Current and future climate features (bioclimatic variables) were obtained from the WorldClim V2.1 database (Fick and Hijmans 2017), which are generated from the CMIP6 model ensemble outputs and bias corrected using current weather station observations. From this data we have selected all (25) climate models, 3 future climate scenarios (sustainability (SSP1-2.6), middle of the road (SSP2-4.5), regional rivalry (SSP3-7.0) for the end-of-century time frame (period 2081–2100) as well as data on current (1970–2000) climate conditions. Each of the three future climate scenarios is available from all climate models, with few exceptions (see table S1 in the supplementary information) (Fick and Hijmans 2017).

This modelling approach is a space-for-time analogy approach in which future predictions of potential biomass and disturbance regimes are based on currently observed values found in locations with similar environmental characteristics (otherwise known as climate analogues). We excluded the strongest climate change scenario (fossil-fueled development (SSP5-8.5)) available in the WorldClim V2.1 database, as climate analogues for this scenario are unlikely to be found. Additionally, we focus on the end-of-century time frame as this model does not have any explicit time dimension. Results presented here need to be interpreted as if the climate seen in the bioclimatic features is stable at the end of the century and forests have had the time to adjust to these new conditions. The analysis is performed at a 10 km spatial resolution using an equal-area Eckert IV projection for the entire global domain.

The last component of the modelling setup, CO₂ fertilization, is calculated separately with three different methodologies: (1) ESMs, (2) time series of satellite-derived proxies for tree growth and (3) experiments of CO₂ enrichment (free-air concentration enrichment experiments). First, to derive a global CO₂ fertilization factor from ESMs we used the 'Biogeochemical' (BGC) experiment of the C4MIP framework (Jones et al 2016) (containing a subset of the CMIP6 ensemble). The C4MIP framework is an experimental model setup in which the CO₂ concentration is increased consistently with 1% per year, starting from preindustrial CO₂ concentrations. It is subdivided in three runs: (1) only the biogeochemical cycles 'see' the increase in CO₂ concentrations (C4MIP: 1% BGC), while the radiative forcing protocol keeps running with preindustrial CO₂ concentrations; (2) only the radiative forcing 'sees' the increase in CO₂ concentrations, while the biogeochemical cycles do not; and (3) a control run, where both biogeochemical and radiative forcing respond to enhanced CO2 concentrations.

To derive a CO₂ fertilization factor, we used the results of the first experiment of all 11 models used the C4MIP runs by calculating the fractional increase of total live carbon (C_LIVE, an approximation of biomass) for future predicted CO₂ concentrations under the socio-economic pathways as described in Meinshausen et al (2020), compared to biomass under current CO₂ concentrations (390 ppm corresponding to the level at the time of the original biomass data (Santoro 2018)). Note that this factor is different from the β-parameter (CO₂ effect on the carbon cycle) used in IPCC literature.

In the second approach, we integrated satellite-based proxy data of CO₂ fertilization of primary productivity (GPP, instead of biomass stocks) with the approach outlined above. This is because deriving a global CO₂ fertilization factor from satellite observations alone is currently not feasible. Calculation methods are identical, but instead of using future atmospheric CO₂ given for the different climate scenarios, we use the concentration where the trend in CO₂ fertilization of GPP observed by Wang et al (2020) reaches zero (approximately 405 ppm).

The third approach relies on results from FACE experiments. Terrer et al (2019) derived a global ${F_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ factor of 1.12 ± 0.03 by extrapolating the results of experimental CO₂ enrichment studies, often called free-air concentration enrichment (FACE) experiments. This factor has only been derived for the SSP2–4.5 scenario (at 625 ppm) due to the setup of the FACE experiments, and we will use this value directly in this context.

To interpret the results, some considerations need be made about the CO₂ fertilization factor. All derivations have strong caveats. The factor derived from most ESMs is overoptimistic, as comprehensively described in Terrer et al (2019). This is mainly due to the incomplete representation of other limiting factors to tree growth within the model implementations, such as nitrogen (some models) and phosphorus availability (most models). We included the factor in this analysis as a reference, as it underpins the current state of the art future carbon projections, the IPCC literature, and the international climate change mitigation policies. The second, satellite derived factor, is based on proxies of GPP instead of carbon stocks. To extract a factor, we assumed that when the sensitivity of GPP on the atmospheric CO₂ concentration is reduced to zero, biomass will reach equilibrium. In practice, biomass can still increase even after GPP becomes insensitive to further CO₂ increase, although this effect will likely be minor relative to the current forest carbon stock (Luyssaert et al 2008). The third, FACE derived factor brings high uncertainty because the global distribution of the FACE experiments is strongly concentrated in temperate regions, with no or few experimental sites in arid and tropical regions.

Additionally, it needs to be noted that all three schemes boil down to a single global constant fertilization factor. As the carbon fertilization effect is strongly dependent on microclimate, soil water availability and specific species composition, this choice likely leads to an underestimation of the CO₂ fertilization effect in water-limited areas, and an overestimation in others (Meinshausen et al 2020). Although the CO₂ fertilization component brings a lot of caveats, its inclusion is fundamental to this approach as the climate parameterization from WorldClim only includes the indirect effect of CO₂. The model framework presented here is modular and can deal with various sources of information in a homogeneous way. It expands on the work of Terrer et al (2019) by including a satellite derived CFE factor and uses climate scenarios to study various future pathways of the possibility of the carbon sink to persist over the coming decades.

BCC changes substantially under all climate change scenarios (figure 2). In the high-latitude regions, BCC will increase under all scenarios, independently from the CO₂ fertilization calculation method. The tropical regions, however, show a large divergence, with the C4MIP CO₂ fertilization causing the BCC to increase under stronger climate change scenarios, while a net loss is predicted when following satellite derived CFE. The strength of the CFE from the FACE experiments falls in the middle between the previous two values.

Zoom In Zoom Out Reset image size

Global aggregate values allow us to study the shifts in figure 2 in more detail (figure 3 and table 2). Potential biomass—or the potential growth of forests—increases almost everywhere in the world, with global values ranging from 120 Pg to 155 Pg of AGB. This is partly offset by an increase in natural disturbances, particularly notable in the tropical rainforests (supplementary figure 1), which increases under more severe climate scenarios. Combined, but without including CFE, they indicate an increase between 49 and 58 Pg of AGB. These changes follow major trends described in the literature. First, the likelihood of disturbance is reduced in the high-latitude areas, while increasing in the temporal and tropical zones. Second, the growth potential in the high-latitude regions increases, with reduced temperature variability, higher mean annual temperatures and a consequential longer growing season. Lastly, other regions have a more variable shift in growth potential, based on local energy- and water availability/limitations (Keenan 2015, Forzieri et al 2022).

Zoom In Zoom Out Reset image size

The three CO₂ fertilization schemes show much more diverging results, suggesting a high degree of uncertainty in this component. Following the C4MIP derived CO₂ fertilization scheme, BCC increases with the severity of the climate change scenario, with the maximum value of 383 Pg of AGB, achieved at the end of the century under SSP3-7.0. When looking at the satellite derived CO₂ fertilization, conclusions are starkly different. Under all scenarios, CO₂ fertilization equates to approximately 14 Pg AGB, with the resulting change in BCC falling between 63 and 73 Pg AGB. Under the SSP3-7.0 scenario, this is a difference of 310 Pg of AGB between the carbon fertilization schemes, equivalent to about 155 PgC. The CFE derived from the FACE experiments show intermediate values at 140 Pg of AGB increase (SSP2–4.5). As a reference, the last row in figure 3 shows the difference in the live carbon pool in the C4MIP modelling intercomparison product. The estimate lies close to the change in BCC, given that the calculations are developed largely independently from C4MIP and only share their relationship between atmospheric CO₂ concentrations and its direct effect on vegetation.

3.1. Changes in water and energy availability

To better understand the main drivers of the changes presented above, we assessed the change in the main climate parameters of the modelling framework in regions with future gains or losses in BCC compared to current values (SSP2-4.5 at the end of the century, without considering CFE). We used four major climate determinants of the spatial distribution of forest growth: (1) annual mean temperature, (2) temperature seasonality, (3) annual precipitation and (4) precipitation of the driest quarter. The change in BCC yielded by increases in annual temperature is largely positive but does not substantially capture the effective change in BCC (figure 4). Temperature seasonality instead shows a much clearer response. The stronger the decrease in temperature seasonality (or the lower the increase), the more beneficial it is for forest growth, likely because of the positive impact on the length of the growing season. This effect is especially pronounced in the high-latitude regions, where currently the strong seasonal cycle limits the growing season.

Zoom In Zoom Out Reset image size

A stronger increase in annual precipitation almost everywhere correlates with a higher chance in increased expected biomass. An exception to this pattern is the high-latitude regions. This could potentially be explained by a reduction in energy availability (by increased cloud cover) rather than an actual reduction in BCC caused by higher water availability, although this could also occur in regions where waterlogging is currently already hindering the growth of trees. Precipitation in the driest quarter, a fundamental parameter for understanding the establishment and survival of trees, shows similar (and even stronger) patterns as annual precipitation. In the low-latitude regions an increase in expected biomass overlaps with an increase in water availability in the driest quarter, while a loss of BCC corresponds to a decrease. Again, at the highest latitudes the relation seems to reverse, but the reasoning applied above should equally hold in this case: a loss in energy-availability might be the culprit in the loss of BCC in such energy-limited systems.

For the interpretation of the results, it needs to be considered that the model framework presented here is an equilibrium forest model and does not consider any temporal dimension other than the one coming with the changes in climate predicted by the CMIP6 ensemble. This means that evolution of forests is not directly modelled but is assessed at quasi-equilibrium with the changing climate conditions. This likely implies that actual carbon stocks would lie below the one predicted under the equilibrium conditions considered by the presented modelling framework. Reasons are that (A) increases in biomass occur relatively slowly, while a more severe disturbance regimes can show its effect on much shorter time scales and (B) the species composition at optimal (at least for the point of view of carbon storage) equilibrium conditions might differ drastically from currently present composition. For example, in the boreal zones the evergreen needle-leaf forests that are currently present have been found to struggle under changing conditions because of the increased growing season length and the consequential increase of transpiration (and thus the lower water availability during the growing season) (Buermann et al 2018). Changes in species composition, especially in the world's coldest regions, are a relatively slow process which might require decades to centuries before forests are again in equilibrium with local climate.

Our results show a high level of uncertainty in predicting future carbon storage in forests, with the estimates strongly depending on the method, data sources and modelling strategy. We designed the machine learning model framework to be flexible and allow for the ingestion of different estimates and scenarios of a changing climate, as well as (partially) independent estimates of the CO₂ fertilization effect. In doing so, however, we propagate the assumptions and uncertainties contained in the data sources. From our results we conclude that the most uncertain component in estimating future carbon stocks in forests, other than the climate scenarios themselves, is the CO2 fertilization effect. Following the evidence from the climate models yields much higher values than by constraining the fertilization with satellite observations, with the results from the FACE experiments falling in between. These three methods all have their own shortcomings: climate models may overestimate the sensitivity of photosynthesis to the atmospheric CO₂ concentration because not all constraints of primary productivity are included. Trends derived from satellite records are based on proxies of photosynthesis, which then need to be translated to effective changes in carbon fluxes, and the separation of the effects of the multiple co-drivers is challenging. FACE experiments might yield the most ecologically accurate representation of how vegetation would react to an increase in atmospheric CO₂ concentration, but they are limited in spatial extent. To extrapolate these plot experiments to realistic global estimates, however, would require a vast network of such sites distributed over ecosystem types and age classes. Currently, most FACE experiments are in temperate climates and focus on young forests and therefore the extrapolation to a global value necessarily generates large uncertainties. Adding more data sources or using different modeling strategies that include a disturbance regime and various estimates of CO₂ fertilization would likely not reduce uncertainty in the estimates and would likely leave our conclusions unchanged.

4.1. Human interventions

4.2. Implications

The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.6084/m9.figshare.25002353.

Conceptualization was done by all co-authors. C T J R designed the methodology with the contribution of A C. C T J R performed the formal analysis. The text writing of the original draft was done by C T J R, A C, G D, L C, R A. E L D and S I S supervised the work.