Enhanced growth – more than reforestation – counteracted biomass carbon emissions (1990-2020)

Understanding the carbon (C) balance in global forest biomass is key for climate change mitigation. However, land-use and environmental drivers affecting forest C fluxes remain poorly isolated and quantified. Following a counterfactual modelling approached based on global Forest Resource Assessments, we show that in the hypothetical absence of changes in forest (i) area, (ii) harvest or (iii) burnt area, global forest biomass would reverse from an actual cumulative net C source of c. 2.4 GtC to a net C sink of 23.6, 2.7 and 5.2 GtC, respectively in 1990-2020. In contrast, cumulative emissions would be 12.9 GtC, i.e. 6 times higher, in the absence of forest growth changes. A typology systematically assessing C flux drivers reveals that enhanced growth, more than reforestation, counteracted C emissions. This sink function may, however, be discontinued in the future, thus alerting to the need for safer strategies to effectively preserve or enhance C sequestration.


Introduction
Terrestrial ecosystems globally act as net carbon (C) sinks, thus providing a key ecosystem service for global climate change mitigation 1 . By contrast, global forest biomass has acted as a net C source to the atmosphere over the last three decades, according to the most recent forest resource assessment 2 . This global source arises from the complex interaction of several drivers, operating at different time scales 3 , including forest area change, and the balance of gross primary production, respiration in vegetation and losses through disturbances or extraction 1 . While depletion of C stocks by deforestation and other disturbances is immediate, regrowth is slow and depends on forest age-structure 4 and on the management and climatic conditions affecting growth [5][6][7][8] . In the context of climate change mitigation, it is pivotal to disentangle the roles of these drivers to tap the potentials of different forest-based mitigation strategies 9 . Nevertheless, while several studies isolated the role of land-use and land-cover change on forest biomass C dynamics [10][11][12] , a consistent comparison of the impacts of different drivers at the global scale has never been performed.
Here, we fill that gap by combining the most recent and consistent global forest dataset provided by the Forest Resource Assessment 2 -an authoritative data source 13 -with the parsimonious forest C model CRAFT (CaRbon Accumulation in ForesT 14 ). FRA contains information on biomass Cstocks, harvest and forest area changes, but this information is not sufficient to quantify the relative importance of specific drivers, because forest growth information is required to close the balance between the dynamics of these three parameters. CRAFT builds upon the establishment of a place and time specific relationship between net primary production (NPP) and biomass C stocks (see methods and SI 1) and so provides a suitable tool to systematically integrate the available information. Thus, we could isolate and quantify the relative impact of various drivers on forest change, including, for the first time, changes in forest growth. We calculate the temporal dynamics of production and primary (i.e. unused) forest growth in 152 countries in 1990-2020 and couple counterfactual scenario development with a typology approach (see Method section) in order to answer the following question: Which role do the individual drivers changes in area, harvest, burnt area, and forest growth play for the observed C stock changes in national and global forest biomass dynamics over the last three decades? Answering this question is essential for assessing the efficiency of various forest-based climate change mitigation strategies.

Trends in global biomass C stocks
The CRAFT model reliably reproduces the observed trends in primary and production forest biomass C stocks (including both aboveground and belowground biomass) in 1990-2020 with a relative root mean square error (RMSE) between simulated and observed biomass C stocks of 0.7% at the global level. Thus, the CRAFT simulations corroborated the FRA observations while adding information on annual estimations of forest C stocks, rather than the 5-years interval data provided by the FRA (Fig. 1a), and annual net C emissions ( Fig. 1c and d) from production and primary forests. In line with the FRA data, the main trend is a decline of total biomass C stocks following three phases: stagnation, increase and again stagnation, resulting in net C emissions from forest biomass (Fig 1c) by 2.4 GtC between 1990 and 2020, contrasted by an opposite trend of increasing biomass density from 68 to 72 tC/ha in total forest ( Fig. 1b and d). This figure is also within the range of the estimated sink in forest soil and biomass of 0.1 ± 7.3 GtC/yr in 2001-2019 found by Harris et al 15 . Here we find that the net C emissions mostly arise from primary forests, which undergo area loss, but also biomass thickening ( Fig. 1b and d). By contrast, in spite of area loss, production forests act as C sinks following biomass thickening ( Fig. 1b and d). Increasing biomass density is therefore key to counteract net C emissions from forest biomass in 1990-2020.
While both harvest rate and burnt area increase globally over the period of observation, the increased forest growth that we calculate with CRAFT for both primary and production forest over 1990-2020 emerges here as the only factor explaining increased biomass density at the global level.

Proximate drivers of net C emissions
We develop six counterfactual scenarios 15,16 in order to investigate how forest biomass density and forest biomass C stocks would evolve in the absence of: (i) change in harvest (CF1); (ii), change in forest growth (CF2); (iii) change in burnt area (CF 3); (iv) change in forest area (CF 4); (v) harvest (CF5); (vi) burnt area (CF 6) (see methods). The comparison of observed and simulated trends allows to isolate and quantify the influence of these four main drivers on global forest C stock changes at national resolution (CF1 to 4) as well as to quantify the overall effects of total wood extraction and burnt area (CF5 and 6).
At the global level, we find that loss of forest area (CF4) is the main driver of the net C emissions from forest biomass (Fig. 2a). In the absence of changes in area, global forest biomass would act as a cumulative net C sink of c. 24.0 GtC in the study period, creating a difference of 26.4 GtC between the actual and the CF4 C budget. This effect, however, is a composite of a gross additional These global trends are the combined result of diverging national forest dynamics ( Fig. 2c to h).
In particular, shifts in forest area (CF4) contribute to global net C emissions only in the Global South, excluding Vietnam, India and Chile (Fig. 2f). The impacts of changes in burnt area and harvest are similarly heterogenous, with considerable effects only in some regions (e.g., Vietnam, Mozambique, Fig. 2c and e). In contrast, changes in forest growth are more ubiquitous, mainly positive (leading to C-sinks) for most countries, with a few notable exceptions, mainly in arid or boreal regions (e.g., India, Spain, Argentina, Canada; Fig 2d). Possible reasons explaining the negative effect of change in forest growth are forest degradation, increasing drought, cloudiness, or insect outbreaks [15][16][17][18][19] . Over the period 1990-2020, the strongest harvest impacts are observed in countries with large area of production forest and high harvest pressure, mostly located in temperate and subtropical areas (CF5; Fig. 2g), while fire impacts are strong in fewer countries in many different biomes (CF6; Fig. 2h).

Typology of forest biomass change
In order to identify spatial and temporal patterns of drivers in forest biomass trends, we establish a typology of the main drivers over the period 1990-2020 (Fig. S3). The typology we established is based on the positive versus negative shift in biomass C stocks, and highlights the most important driver of the shift assessed through the counterfactual assessment, irrespective of the relative importance of the other drivers shown in Fig. 2. By pinpointing the major drivers of forest change at national levels, such an approach enables to identify major levers for forest conservation.
Deforestation was the dominant driver of net C emissions from forest biomass in most countries of South America and Sub-Saharan Africa, corroborating findings from the literature 10,16,17 ( Fig.   3a and b). The net C emissions by countries where deforestation is the most significant driver reach c. 21 GtC, representing c. 90% of the 23 GtC net emissions arising from all countries acting as net C sources (Fig. 3b). Changes in forest growth act as the primary drivers in most countries experiencing a net C sink over the period ( Fig. 3a and b). The net C sinks by countries where changes in forest growth are the main drivers reach c. 16 GtC, representing c. 76% of the 21 GtC net sink created by all countries acting as net C sinks (Fig. 3b).
Forest area expansion from 1990 to 2020 is the main driver of forest biomass net C sink in only a few Northern countries but also some Southern countries, namely Vietnam, India and Chile, in line with findings reported for these countries [18][19][20] . Similarly, changes in harvest as well as changes in burnt areas are the main drivers of net C sink or net C source for a handful of countries in 1990-2020 (Fig. 3a). Finally, declining forest biomass growth is the primary driver of net C emissions only in Canada, which is consistent with other studies highlighting slower growth, higher mortality and insect outbreak events in Canadian forests [21][22][23] .

Implications for forest-based solutions
Our results allow to identify major mechanisms behind observed forest biomass C changes that are immediately relevant for forest-based climate change mitigation strategies. We show that deforestation, increasing harvest and burnt area have driven the net C emissions from forest biomass over the last three decades. Deforestation is the dominant driver, corroborating that protection from deforestation is indispensable 1,10,24 . On the other hand, forest growth is identified as the major driver counteracting net emissions ( Fig. 2a and d). In fact, most of the temperate and boreal countries, with the noteworthy exception of Canada, fall under a type in which enhanced forest growth is the major driver of a net C sink (Fig. 3b). Besides, even countries dominated by deforestation in the tropics, show significant increases in growth ( Fig. 2d and S4). These results highlight that enhanced growth, rather than reforestation, is the main driver counteracting biomass C emissions in 1990-2020.
These increases in forest growth may arise from diverse processes, including climatic and landuse drivers. On the one hand, several studies highlight the feedback of environmental variablessuch as warming, atmospheric carbon dioxide (CO2) and nitrogen (N) fertilization 1,5,7,10,25 -on the terrestrial C sink. On the other hand, changes in forest growth can also be driven by effects of shifts in forest management practices, such as tree species selection, forest recovery from past degradation and lesser litter grazing 11,26,27 . Advancing the understanding of the underlying processes of forest growth change is key for forging climate-change mitigation strategies, but it is not straightforward to isolate climatic (e.g. altered CO2 concentration or temperature) from landuse drivers (e.g. non-timber forest uses such as grazing 28 ). Still, a comparison of trajectories in primary and production forest growth change based on our results allows to derive insights into the interplay of these different drivers (Fig. 4). From the fact that only 12% of primary forest carbon stocks show declining growth trends (Fig. 4b) while a relatively larger carbon stock in production forest (25%) is affected by declining growth trends (Fig 4c), we can infer that in overall terms -and assuming primary and production forests of a given country to be similarly affected by climatic drivers -, land use is likely to exert a degrading effect on growth dynamics.
Nevertheless, some countries reveal declining growth in primary forest but increasing growth in production forest, thus suggesting that forest management may have an improving effect on forest growth in those countries (e.g, USA, Fig. 4b and c, see also Fig. S5). In overall term, this result yet suggests that globally a reduction of forest use may have the potential to enhance growth.
However, these interpretations warrant a caveat that the national view is a high aggregate, which does not allow for more detailed analyses. In addition, primary versus production forest growth changes are derived from the FRA data and a state-of-the-art of the literature on changes in primary forest density (see methods and Fig. S8-9), the latter being associated with higher uncertainties (see sensitivity analysis in methods and SI1). Therefore, we conclude that, while increasing forest growth is the dominant driver counteracting the global net C emissions from forest biomass in the past three decades, it is against a precautionary principle to forge climate strategies that rely on a continuous net C sink effect from the same processes in the future. By contrast, our results suggest that reducing wood harvest (Fig.   2g) and halting deforestation (Fig.2c) have the potential to be key strategies to address the challenge of climate change mitigation. By consequence, increasing forest harvest volumes -a strategy often promoted in the course of climate change efforts embraced as the "bioeconomy"might have strong unintended consequences 38 , e.g., increasing harvest might accelerate carbon turnover times resulting in a loss of C sink capacity 39 . The tradeoffs between the wood provisioning and regulating ecosystem services should, however, be considered in all their complexity, because declining wood harvest might occur together with wood product substitution by other materials generating additional emissions 40-42 . Overall, our results plead for a double strategy to enable future forest-based solutions for climate change mitigation: in the Global South, ending deforestation is the main priority to reverse net C source towards net C sink, while in the Global  From these input and calibration data, we optimized the best growth parameters for each country: r (forest growth, yr -1 ) and K (theoretical maximum carrying capacity, tC ha -1 ) within a range of possible values (see below). Because the forest growth parameter r might be variable over time [3][4][5][6] , we assumed that it may linearly increase or decrease from 1990 to 2020. We considered that the value of the r parameter in 2020 could be in a range between 80% to 120% of its estimated value The optimized growth parameters described the relationship between Net Primary Production (NPP) and standing biomass density (B) (Eq. 1) and, together with the other input data, allow to calculate the forest biomass dynamic by recurrence such as:

Principles of the CRAFT model
With, n the time step (yr); NPP the net primary production (tC ha -1 ); B the biomass density (tC ha -1 ); r the forest growth parameter (yr -1 ), K the theoretical maximum carrying capacity (tC ha -1 ); CL the cutting losses (%); H the harvest rate (tC ha -1 ); FL the fire losses (tC ha -1 ) (see below how the fire losses were estimated); m the mortality rate (yr -1 ) and; A the forest area (ha). The initialization was set up by using the calibration data on forest biomass density.

Range of growth parameter values
To derive a range of possible values for the forest growth parameters, we relied on site measurement data on standing biomass and NPP from all ecozones of the world, compiled by Cannell 7 , Luyssaert et al. 8 and Anderson-Teixeira et al. 9 , contributing 372, 503 and 536 data points, respectively (see table S1). The expansion factors, i.e., the allometric coefficients between tree stems and other tree parts, per ecoregion given in table S2, were used when data reported on biomass missed specific tree organs. Eq.1 was fitted to best reproduce the site measurements in each ecoregion. The r and K parameters of the curve with the smallest residues to the data points were identified (using the MATLAB function fmincon), under the constraints of r being between 0.01 and 1.00, and K being between 0 and 1000. The agreement of fit for the fitted function was calculated. To arrive at a range of r and K values, data was only selectively used for the fitting of to 720 tC ha -1 , the lower and upper deciles of the fitted curve. We used these ranges of values for all countries of the world to optimize nationwide growth parameters for both primary and production forests (see section above).

Input data forest area
Data on forest area were extracted from the last Forest Resource Assessment (FRA) (https://fradata.fao.org/), which distinguishes between 'primary forest', 'mangrove', 'planted forest' and, 'naturally regenerated forest'. As mangrove can be both primary and production forest (FRA, 2020), we calculated the primary forest area as the sum of the 'primary forest' area and half of the 'mangrove' area reported by the FRA. As the mangrove area is generally very small in comparison with primary forest area, this assumption does not significantly influence our simulations. The production forests, i.e., the forests exploited for socio-economic activities, are taken as the total forest area reported by the FRA minus our estimations of primary forest areas.

Input data harvest
Wood harvest data were extracted from FAOstat (http://www.fao.org/faostat/en/). We extracted data for 8 wood categories: wood fuel, pulpwood, sawlogs and veneers logs, other industrial roundwood for both coniferous and deciduous wood (provided in m 3 ). To convert wood extraction from m 3 to tC, we used coefficients on wood density (t m -3 ), bark fraction (tbark twood -1 ) and C content (tC twood -1 ). These coefficients were specific to product and countries following Haberl et al. 10 . We then summed all wood product harvest at the country level to have an aggregated value of annual wood harvest (tC yr -1 ). As the FAO stat provided wood extraction data only from 1990 to 2019, we assumed the same harvest figures in 2020 as in 2019.

Input data on burnt area and fire losses
Data on forest area burnt were extracted from the last FRA (https://fra-data.fao.org/). From this data, we calculated that annual wood losses by fire such as: With FL the annual loss by fire (tC/yr); α the fire efficiency coefficient, i.e., the % of burnt wood in an area affected by fire; (Ab/A), the fraction of area burnt over the total forest area and B the standing biomass (tC). We assumed that the fraction of area burnt was the same in production and primary forest as in total country-level forest.
We used the assessment by Ito and Penner 11 , which provides combustion efficiency coefficients for 7 regions of the world, both for woodland (tree cover representing 40-60% of the land cover) and forest (tree cover representing more than 60% of the land cover) (see table S4). We also compared the data reported by Ito and Penner 11 with other data reported at the national level in the literature and found the Ito and Penner data to lie well within the range of other data reported (see SI 1). We estimated fire efficiency coefficients at the national level by deriving the fraction of national forest belonging to one or several macro-regions described by Ito and Penner (2004) as well as the proportion of woodlands and forests within each country. Forest and woodland area data were based on Erb et al. 12 , including wild, used, natural and artificial forests and wooded lands.
Additionally, because the extent of forest area burnt may overlap with deforestation area, we corrected the data provided by the FRA following a 'best guess' approach. We assumed that overlap between forest area burnt and forest area loss were likely to happen in tropical countries 13,14 . By contrast, in boreal and temperate countries fires are rarely purposeful and leading to deforestation 14,15 . Therefore, for tropical forest, we assumed that: if deforestation occurred (i.e., ΔA < 0) and if the extent of fire area was higher than the extent of deforestation (i.e., Ab -ΔA > 0), then the area burnt was corrected by the deforested area, thus assuming that forest area loss was already included in burnt area (i.e., Ab = Ab -ΔA). If deforestation occurred and if the extent of deforestation is higher than the extent of area burnt (i.e., Ab -ΔA =< 0), then the area burnt is corrected by zero, thus assuming that the extent of burnt area was already included in the deforested area (i.e., Ab = 0). If forest area expansion occurred (i.e., ΔA >= 0), then deforestation cannot overlap with burnt area, thus we applied the data on the extent of forest area burnt as provided by the FRA. For all other countries we use the data on the extent of area burnt as provided by the FRA without applying correction.

Calibration data: C stocks in production and primary forest
Calibration data on standing biomass C stocks (including both above-ground and belowground biomass) were extracted from the last FRA (https://fra-data.fao.org/). Conversely to the forest area, the FRA did not provide data on biomass C stocks distinguishing between primary and production forest but only for total forest biomass C stocks at the country level. In the present study we thus calculated the production versus primary forest biomass C stocks by using the benchmark values of primary forest biomass density provided by Erb et al. 12 for the year 2000 at the country level.
From this data, we could calculate the biomass C stocks in both primary and production forest in 2000 such as: With Bprim2000, Bprod2000 and Btot 2000 the C stocks biomass in 2000 respectively in primary, production and total forest. Before and after 2000 forest biomass density of production and primary forest may follow different trends, as empirical studies revealed that primary forests are not in equilibrium and may undergo changes in biomass density 16 . In order to consider these changes, we carried out a literature survey to derive temporal trends in boreal, temperate, paleotropical and neotropical primary forest densities (see table S8 in SI 1). Subsequently, we used these macroregional coefficients to derive biomass C stocks in production and primary forests at the national With BPrim,y, BProd,y and BTot,y the biomass C stocks (MtC) in primary, production and total forest in year y, with y belonging to ; APrim,y the primary forest area in year y (ha); BDPrim,2000 the biomass C stock density (tC ha -1 ) in primary forest in the year 2000 (as provided by Erb et al. 12 ) and; δ the annual change in primary forest density (%/ yr -1 ) derived from the literature survey (see table S9 in SI 1).

Counterfactual scenario analysis
In order to isolate and quantify the relative impact of major proximate drivers on forest biomass C change, we develop six counterfactual scenarios 17,18 . In these counterfactual scenarios we investigate how forest biomass density and forest biomass C stocks would evolve in the absence reported in black, while data extracted from the FAOstat and FRA database are reported in blue (see SI2).

Sensitivity analyses
As the CRAFT model reproduces the observed trends in forest biomass C stocks, its highest uncertainties lie in the estimation of the growth parameters. Here, the two main hypotheses of the model which are likely to affect these calculations are: (i) the coefficient derived to estimate the temporal trend in primary forest biomass density and (ii) the assumption regarding the optimization of temporal changes in growth parameters. To assess the uncertainties associated to those assumptions, we ran 2 sensitivity analyses: (1) We tested one alternative hypothesis in which biomass density in primary forest was taken as constant (using the benchmark provided in 2000 by Erb et al. 12 ) against the default assumption based on an estimation of changes in primary forest biomass density. (2) We tested another alternative hypothesis in which both the temporal trend of the r and K parameters were optimized against the default assumption based on the optimization of the temporal change of the r parameter only.
The differences in the model outputs and counterfactual scenarios assessments between the two sensitivity analyses and the default model assumptions are presented in Fig. S6 to S10. Both sensitivity analyses enabled to reliably reproduce the observed data on production and primary forest biomass C stocks (Fig. S6, S10-11). Unsurprisingly, the sensitivity analysis revealed that different hypotheses would affect the calculation of the r and K parameters as well as the temporal change of the r parameters from 1990 to 2020 but the trends presented in the main manuscript would remain the same (Fig. S6, S8-11). The sensitivity analyses also revealed that the results of the CF assessment would be affected, but not the general trend arising from them (Fig. S7).