Mapping global forest age from forest inventories, biomass and climate data

Forest age can determine the capacity of a forest to uptake carbon from the atmosphere. However, a lack of global diagnostics that reflect the forest stage and associated disturbance regimes hampers the quantification of agerelated differences in forest carbon dynamics. This study provides a new global distribution of forest age circa 2010, estimated using a machine learning approach trained with more than 40,000 plots using forest inventory, biomass and climate data. First, an evaluation against the plot level measurements of forest age reveals that the data-driven method has a relatively good predictive capacity of classifying old-growth vs. non-old-growth (precision = 0.81 and 0.99 for old-growth and non-old-growth, respectively) forests and estimating corresponding forest age estimates (NSE = 0.6 and RMSE = 50 years). However, there are systematic biases with overestimation in young and underestimation in old forest stands, respectively. Globally, we find a large variability of forest age with the old-growth forests in the tropical regions of Amazon and Congo, young forests in China and intermediate stands in Europe. Furthermore, we find that the regions with high rates of deforestation or forest degradation (e.g., the arc of deforestation in the Amazon) are composed mainly of younger stands. Assessment of forest age in the climate-space shows that the old forests are either in cold and dry regions or warm and wet regions, while young-intermediate forests span a large climatic gradient. Finally, comparing the presented forest age estimates with a series of regional products reveals differences rooted in different approaches and different in-situ observations and global-scale products. Despite showing robustness in crossvalidation results, additional methodological insights on further developments should as much as possible harmonize data across the different approaches. The forest age dataset presented here provides additional insights into the global distribution of forest age to better understand the global dynamics in the forest water and carbon cycles. The forest age datasets are openly available at https://doi.org/10.17871/ForestAgeBGI.2021 (Besnard et al., 2021).

Abstract. Forest age can determine the capacity of a forest to uptake carbon from the atmosphere. However, a lack of global diagnostics that reflect the forest stage and associated disturbance regimes hampers the quantification of agerelated differences in forest carbon dynamics. This study provides a new global distribution of forest age circa 2010, estimated using a machine learning approach trained with more than 40,000 plots using forest inventory, biomass and climate data. First, an evaluation against the plot level measurements of forest age reveals that the data-driven method has a relatively good predictive capacity of classifying old-growth vs. non-old-growth (precision = 0.81 and 0.99 for old-growth and non-old-growth, respectively) forests and estimating corresponding forest age estimates (NSE = 0.6 and RMSE = 50 years). However, there are systematic biases with overestimation in young and underestimation in old forest stands, respectively. Globally, we find a large variability of forest age with the old-growth forests in the tropical regions of Amazon and Congo, young forests in China and intermediate stands in Europe. Furthermore, we find that the regions with high rates of deforestation or forest degradation (e.g., the arc of deforestation in the Amazon) are composed mainly of younger stands. Assessment of forest age in the climate-space shows that the old forests are either in cold and dry regions or warm and wet regions, while young-intermediate forests span a large climatic gradient.
Finally, comparing the presented forest age estimates with a series of regional products reveals differences rooted in different approaches and different in-situ observations and global-scale products. Despite showing robustness in crossvalidation results, additional methodological insights on further developments should as much as possible harmonize data across the different approaches. The forest age dataset presented here provides additional insights into the global distribution of forest age to better understand the global dynamics in the forest water and carbon cycles. The forest age datasets are openly available at https://doi.org/10.17871/ForestAgeBGI.2021 (Besnard et al., 2021).

Introduction
Forests cover about 30% of the terrestrial surface of our planet and store a large part of the terrestrial carbon, indicating their fundamental role in the terrestrial carbon cycle ( Bar-On et al., 2018). However, drivers controlling the capacity of the terrestrial biosphere to sequester carbon remain poorly characterized, limiting our understanding of the global land carbon sink's location (Cook-Patton et al., 2020). Such uncertainties on the geographical distribution of the carbon sink have been partly attributed to the fact that forest regrowth and demography are not systematically considered for understanding changes in the forest carbon sink Zscheischler et al., 2017).
While the recent increase in the forest carbon sink is controlled by environmental changes such as carbon dioxide (CO2) fertilization, nitrogen deposition, and climate change (Zhu et al., 2016;Winkler et al., 2021), the forest carbon balance dynamics are also attributed to disturbance history and forest regrowth (Besnard et al., 2018;Pugh et al., 2019;Amiro et al., 2010). Forest disturbances cause physical damages to vegetation properties and changes in forest demography, thereby affecting the balance of terrestrial CO2 exchange with the atmosphere by temporarily increasing respiration and reducing photosynthesis (Birdsey et al., 2006;Johnson and Curtis, 2001;Liu et al., 2011;Williams et al., 2012;Woodbury et al., 2007). The changes in the strength of carbon uptake or release can alter the forest carbon balance by converting forest ecosystems from carbon sinks to sources (Amiro et al., 2010;Bowman et al., 2009;Ciais et al., 2014;Moore et al., 2013). (Odum, 1969) hypothesized the first theory to describe the ecosystem development in the absence of significant disturbance, suggesting that the age of forests and demographic changes drive carbon accumulation.
Nevertheless, stand age distribution can be modified to varying degrees of changes in environmental conditions and disturbance, therefore slowly change along with stand age or successional continuum (Irvine et al., 2005;Piponiot et al., 2018).
Despite the sensitivity of the forest carbon balance to disturbance and regrowth (Buitenwerf et al., 2018;Sulla-Menashe et al., 2018), existing empirical models and current bottom-up spatiotemporal assessment of CO2 fluxes do not explicitly account for these effects (Jung et al., 2011;Tramontana et al., 2016;Jung et al., 2020). The forest carbon balance in regions with newly disturbed and old-growth forests may not be realistically estimated by not explicitly constraining data-driven statistical models with disturbance history or forest demography. For instance, large discrepancies are observed between bottom-up statistical approaches (e.g., FLUXCOM initiatives, http://www.fluxcom.org/) and atmospheric inversions in estimating net ecosystem exchange (NEE), particularly in the tropics where site history plays a substantial role in NEE magnitude . To account for the contribution of disturbance on the land carbon sink, we need information on the geographical distribution of disturbance, albeit such information is somewhat limited at the global scale (Ciais et al., 2014). Forest age, related to time since disturbance, can be seen as a useful surrogate in analyzing the impact of disturbance on the ability of forests to store carbon.
Incorporating forest age into terrestrial biosphere modelling offers a starting point to characterize disturbance history, so getting more insights into the location of the terrestrial carbon sinks . Reliable estimates of forest age at the global-scale are, therefore, an essential and needed source of information. The recent advances in describing the geographical distribution of forest demography globally (Huang et al., 2010;Kennedy et al., 2010;Poulter et al., 2019) have paved the way to consider forest age and disturbance history in carbon cycle studies.
Here, we aim to provide a new gridded global forest age dataset circa 2010 inferred from a compilation of forest inventory, biomass and climate data. More specifically, we introduce the in-situ forest inventory dataset, the modelling framework used in this study and the predictive capacity of the presented model to estimate forest age at the plot level.
We further describe the global and regional patterns of the forest age product and their uncertainties. The presented forest age dataset is finally benchmarked against a series of independent regional and global datasets.

Forest inventory and climate data
The globally gridded forest age dataset was developed by collecting in-situ plot level stand age and aboveground biomass (AGB) estimates from a series of forest inventory databases (Álvarez-Dávila et al., 2017;Anderson-Teixeira et al., 2018;Anderson Teixeira et al., 2016;Baker et al., 2016;Johnson et al., 2016;Lewis et al., 2013;Mitchard et al., -2014;N'Guessan et al., 2019;Poorter et al., 2016;Schepaschenko et al., 2017;Somogyi et al., 2008;Sullivan et al., 2017 Where i is a decadal class and n is the number of observations. The methods used in inventory surveys to estimate stand age relied on expert knowledge, tree diameter measurements, tree rings from cores of selected trees (e.g., Burill et al. 2018), or semi-directive interviews (e.g., (N'Guessan et al., 2019)). Forest inventory plots were classified as old-growth forests when stand age was more than or equal to 300 years. In total, the final dataset had around 25,000 plots and around 44,000 observations. Geographical biases were observed in the compiled dataset ( Fig. 1), with North America, Europe and South East of China being well represented, while Africa, Indonesia, and Australia were either underrepresented or not represented. The Amazon basin and the West part of Eurasia were relatively well represented. Besides, stand age data were generally collected in locations easily accessible; therefore, unmanaged forests in remote areas were very likely less represented than managed forests.

Figure 1
Spatial distribution of the forest inventory plots used for the forest age maps. Each dot represents the total number of plots within 10x10 degrees.
A comprehensive meta-analysis of the compiled dataset ( Fig. 2) revealed that the observations covered a large spectrum in the climate-space ( Fig. 2A), although there were few plots in hot and dry regions, likely due to the low presence of forest ecosystems in such regions. We further described the age spectrum covered at the regional scale and found that a large spectrum of forest age was covered in North America (Fig. 2B) and Eurasia (Fig. 2C), while in the tropics, biases were observed (i.e., a significant fraction of tropical old-growth forest and relatively young forests) (Fig. 2D). For each forest inventory plot, we extracted bioclimatic variables from the WorldClim version2 (Fick and Hijmans, 2017). In addition, we extracted all the soil-related variables of the Harmonized World Soil Database v 1.2 dataset.
Finally, we derived a series of proxies for disturbance and management regimes from the Hansen tree cover dataset (Hansen et al., 2013):  The intensity of tree loss from the Hansen tree cover loss layer. This metric was derived by counting the 30m pixels that experienced a tree cover loss for 2000-2019 within a 1km gridcell.
 Last time since tree cover loss from Hansen tree cover loss layer -standard deviation metric. This metric was calculated as the last time from 2019 since a 30m pixel experienced tree cover loss, and we further computed the standard deviation of this last time since tree cover loss within a 1km gridcell. Table S1 summarizes the list of covariates considered in our study. Two datasets were further created. First, we created a dataset that contained the plots with a reported stand age estimates ranging from 1 to 299 years (hereafter non-old-growth forests dataset). Second, we created a binary dataset reporting whether an observation had an age estimate less than 300 years old or whether an observation had an age estimate more than or equal to 300 years old or not reported but considered as old-growth tropical forests (0= non-old-growth forest and 1= old-growth forest) (hereafter old-growth forests dataset).

Feature selection and model training
From the set of predictors related to vegetation and climatic conditions (Table S1) dataset, while the RFclassifier model was used to classify old-growth vs. non-old-growth forests using the old-growth forests dataset. The performances of the two models were assessed using leave-one-cluster-out cross-validation to reduce possible spatial auto-correlation between the training and test sets (Ploton et al., 2020). A cluster of plots contained all the plots within the same pair of latitude and longitude coordinates rounded to the nearest 10th degree (e.g., latitude= 20 degrees and longitude= 110 degrees) (see Fig. 1). For the RFregressor model, the root-mean-square error (RMSE), the normalized root-mean-square error (NRMSE) and Nash-Sutcliffe model efficiency coefficient (NSE) were used for assessing the predictive capacity of the model for predicting forest age. For the RFclassifier model, we reported the precision (i.e., the number of correctly-identified members of a class divided by all the times the model predicted that class), recall (i.e., the number of members of a class that the classifier identified correctly divided by the total number of members in that class) metrics, and F1-score (i.e., the combination of precision and recall).
Additionally, we explored functional relationships between the variables selected by the feature selection procedure and stand age in the RFregressor model by using the Tree SHapley Additive exPlanations (TreeSHAP) algorithm (Lundberg and Lee, 2017;Lundberg et al., 2019). A negative SHAP value for a given variable X translates a negative contribution to the local changes of forest age and vice-versa.

Upscaling procedure
To upscale the two trained models (i.e., RFclassifier and RFregressor models) from plot-level data to the global scale, we collected climate grids from the WorldClim dataset (Fick and Hijmans, 2017) (Fick and Hijmans, 2017) and a series of AGB grids circa 2010 (i.e., corrected for tree cover with thresholds of 0%, 10%, 20% and 30%) from the Globbiomass project (http://globbiomass.org/) (Santoro et al., 2021). The tree cover correction was done by masking out the 100-meter pixels in the original AGB product (i.e., 100m resolution) having tree cover estimates (Hansen et al., 2013) below one of the tree cover thresholds mentioned above within a 1km extent. The original filtered AGB maps were aggregated from 100m to 1km spatial resolution with a bilinear resampling method.
The upscaling procedure was done in two steps. First, each 1km pixel was classified as old-growth or non-old growth forests using the trained RFclassifier model. Second, the 1km pixels classified as non-old growth were assigned with an age estimate ranging from 0-299 years inferred from the RFregression model, while the pixels classified as old-growth forest were assigned a default age value of 300 years. In total, four gridded forest age maps with a 1km spatial resolution were obtained using the different AGB maps derived from the different tree cover thresholds as mentioned above (hereafter MPI-BGC forest age datasets). We also created maps from the 1km resolution forest age maps that reflected the fraction of several age classes (0-300+ with decadal resolution) within each 0.5-degree grid cell resolution.

Model development and evaluation
We used the ten most important variables from the set presented in Table S1 identified by the RFE algorithm procedure for the RFregression and the RFclassifier models (Table 1). This set of selected variables was further used to train the two models in the cross-validation analysis and the global upscaling procedure. By assessing the cross-validation results, we found that the RFclassifier model could accurately partition old-growth and non-old-growth forests with precision estimates of 0.81 and 0.99 for old-growth and non-old-growth forests, respectively (Fig. 3A). Furthermore, we found recall values of 0.94 and 0.98 for old-growth and non-old-growth forests, respectively, while we found F1-scores of 0.87 and 0.99 for old-growth and non-old-growth forests, respectively (Fig.   3A). The performance of the RFregression model was relatively high (NSE= 0.60, RMSE= 47.63 years and NRMSE = 51.52%) (Fig. 3B), while the model residuals across 10-degree latitudinal averages were relatively low (Fig. 3C).
However, the quantile-quantile plot depicted biases in both very young and old forests (Fig. 3D). More precisely, the observed forest age estimates from the regression model (B). In C, the average model residuals standard deviation ∓ within 10-degree latitudinal beans is shown. The quantile-quantile plot (D) is also shown.
We further investigated the variable importance of the selected variables and the functional relationships learned by the RFregression model between forest age and these selected variables. For this, we computed the SHAP values for each predictor to show how each predictor contributes, either positively or negatively, to the forest age estimates. First of all, we observed that vapr was the most important variable, followed by agb and MeanTemperatureWettestQuarter (Fig. 4).
The importance of atmospheric water demand in explaining stand age variability could indicate how biomass is 9 206 207 208 209 210 211 212 213 associated with stand age across different climate regimes. More precisely, such observations could imply that high atmospheric water demand limits growth rates and maximum biomass, thereby indirectly controlling how biomass relates to age. In addition, high atmospheric water demand might influence fire frequency (Mueller et al., 2020) and indirectly control forest age distribution through the effect of fire on biomass. Biomass estimates contain information about the current state of the forest, integrating the cumulative effect of land-use change, management and disturbance history. Therefore, having biomass (i.e., agb) as an important variable in predicting forest age confirmed strong management and disturbance regime controls on the forest age distribution (Amiro et al., 2010). The emergent relationships revealed that an increase in AGB was associated with an increase in the forest age estimates (Fig. 5B). This relationship was expected as older trees have more carbon stored in the aboveground components than younger forests. The modelled forest age estimates appeared to be also relatively sensitive to the climatic conditions. For instance, we observed that climatic conditions with low water atmospheric demand (i.e., low vapr) (Fig. 5A) or conditions with high solar radiation (Fig. 5E) increased forest age. Similarly, we observed that forest age variability was also associated with air temperature conditions (Fig. 5C, F, G and H) and precipitation regimes ( Fig. 5D and I).

Global forest age patterns and regional overview
The MPI-BGC forest age product shows an extensive range of forest age across the globe (Fig. 6). We observed that the most represented age class was the old-growth forests with around 1,1 billion hectares, while a limited fraction of very young forest was observed (i.e., < 10 years old) (Fig. 6A). Not surprisingly, most of the old-growth/undisturbed forests (+300 years old) can be found in the Amazon basin (Fig. 6B), the Congo basin (Fig. 6C) and part of the Indonesian peninsula (Fig. 6H), where the minimal human disturbance occurred. A large area occupied by very young forests was found in the Southeast part of China (Fig. 6H), probably due to afforestation/reforestation policies and natural disturbances (Zhang et al., 2017). Similarly, young and intermediate forests were found in the African tropical dry forests (i.e., Sahel and Miombo regions) (Fig. 6C), where the frequency of the fire regimes is very high, resulting in a relatively young age-class structure . Large scale fires in the North American boreal region also resulted in widespread patches of younger forests and a mosaic of stands of different ages since they last burned (Fig.   6G).
Furthermore, the unmanaged part of the North American boreal region near the ecotone, where fires are infrequent, revealed older stands (Fig. 6G). Forests in British Columbia were generally old, although patches of younger forests were probably in the early stages of disturbance recovery. European forests were in young/intermediate stages of forest succession (Fig. 6E). The increased harvested forest area and considerable afforestation practices (Naudts et al., 2016)  probably explained a relatively young to intermediate forest demography and a mosaic of different age classes in the European region. The region of Siberia revealed a gradient of younger to older forests going from the South to the North part of the Siberian region (Fig. 6F). Such an observation could suggest different fire regimes between Southern and Northern Siberia (Shorohova et al., 2011) and confirm harvesting practices identified in Southern Siberia (Curtis et al., 2018). Finally, Australian forests were relatively young in the North part of the country while a mosaic of age class dominated the Southern part of Australia (Fig. 6D). The age patterns observed in the Northern part of Australia somehow correspond to the fact that forests are regrowing in this region . However, it is essential to note that the few forest inventory plots in regions such as Australia (Fig. 1) could limit our certainty on the forest age estimates attributed by the statistical approaches due to, for instance, extrapolation issues. Another limitation is that we assumed forest homogeneity within a 1 km grid-cell, which would reduce the extremes of low and high biomass estimates in the gridded global products that the models have learned in the plot-level training data. This limitation might, for instance, explain the relative dearth of very young stands (1-10 years old) in the MPI-BGC global age product (Figs.6A).

Global forest age relationships with atmosphere, hydrosphere and vegetation conditions
We further investigated the distribution of the forest demography in the climate and vegetation spaces (Fig. 7).
Generally, we observed that with warmer (i.e., air temperature) and drier (i.e., VPD) conditions, forests appeared to be younger with the expectation of old-growth tropical forests located in relatively warm climatic conditions (Fig. 7A). Not surprisingly, we found that most of the old-growth tropical forests were located in regions with high productivity (i.e., high GPP and high biomass) (Fig. 7B), which coincides with our previous results investigating the structure of the statistical model showing that an increase in forest biomass was coupled with an increase in forest age (Fig. 5A). On the other hand, we observed that younger-intermediate forests were more productive than older forests outside the tropical old-growth forest envelope. More precisely, we found that forests being less productive will belong to an older age class for similar carbon stocks. Mature forests were found in cool temperatures and moderately low precipitation conditions (Fig. 7C), where rates of fast growth but slow decomposition generally drive forest dynamics. Younger stands were found in relatively warm conditions but a wide range of precipitation regimes (Fig. 7C). Finally, while a significant fraction of young forests were located in regions with low water availability and high water atmospheric demands, we also observed that above a certain threshold of water availability (i.e., > 0.4-0.5), the amount of water available for trees (i.e., IWA) was not directly associated with changes in forest age unlike VPD (Fig. 7D).

Figure 7
Forest age distribution in the climate, hydrological and productivity spaces defined by air temperature, vapour pressure deficit, total precipitation, soil water availability, GPP and above-ground biomass. The forest age map used here corresponds to a tree cover threshold of 10% aggregated to 0.25 degree using a weighted average of all non-NODATA contributing pixels. GPP is gross primary productivity derived from the FLUXCOM RS+meteo product (Tramontana et al., 2016;Jung et al., 2011Jung et al., , 2020, and IWA is an index for soil water availability (Tramontana et al., 2016).
The climatic variables were retrieved from the ERA5-reanalysis data (https://apps.ecmwf.int/datasets/licences/copernicus/). For all the climatic variables, we computed an annual mean for the year 2010.

Sensitivity analysis, uncertainties and comparison with previous products
We performed a sensitivity analysis using a series of AGB gridded products filtered with different tree cover thresholds to produce different global age products (see method section) (Fig. 8). This analysis showed that in South America, mainly the dry regions were sensitive to the tree cover threshold being applied, with forest age estimates being lower when no tree cover threshold was applied compared to a 30% tree cover correction (Fig. 8A). Similarly, we observed that the dry parts of the Congo basin depicted a sensitivity to the applied tree cover thresholds (Fig. 8B). In Europe, we observed widespread differences between the forest age estimated without a tree cover correction and with a tree cover correction (Fig. 8C). Generally, forest age estimates were higher when the 30% tree cover correction was applied. In Siberia (Fig. 8D), North America (Fig. 8E) and Southeast Asia (Fig. 8F), there were also large patches of forest where correcting the biomass maps with a tree cover threshold led to substantial differences in the age estimates. Overall, such observations were expected because of management practices or disturbance regimes, resulting in mosaic vegetation within a 1km grid cell. Such mosaic vegetation in regions such as the dry tropics (forest/grassland/shrubland), Europe (forests/croplands) and Northeast of the United States (forests/croplands) could explain the sensitivity of the forest age estimates to tree cover thresholds in these regions.

Figure 8
Sensitivity of the presented age product using 30% tree cover correction thresholds or no tree cover correction.
The differences between the age estimates derived from a forest biomass product using a 30% tree cover correction and the age estimates derived from a forest biomass product not using a tree cover correction are shown. Blue colour means that the age estimates are higher with the 30% tree cover correction than without correction, while the red colour means that the age estimates are lower with the 30% tree cover correction than without correction.
Besides, we explored uncertainties associated with the two statistical models used for the upscaling procedure (Fig. S1, Fig. S2 and Fig. S3). First, we observed that the RFclassifier model had very high probabilities of classifying either a non-old-growth or an old-growth forest at pixel level as the fraction of the random forest ensemble to classify the two forest classes was generally close to one. (Fig. S1 and Fig. S2), suggesting relatively high confidence in the partitioning between old-growth and non-old-growth forests in the MPI-BGC forest age product. The regions at the edge of the Amazon and the Congo basins appeared to have the lowest confidence in classifying old-growth vs. non-old-growth forests (Fig. S1A, Fig. S1B and Fig. S2) with a probability close to 0.5. On the other hand, we observed relatively high probabilities for classifying non-old-growth forests in Europe (Fig. S1C), Siberia (Fig. S1D), North America (Fig. S1E) and Southeast Asia (Fig. S1F). We also provided uncertainties in predicting forest age estimates by retrieving the 25%, 50%, and 75% quantile predictions from the RFregressor model for computing the inter-quantile range (IQR, quantile 75% -quantile 25%) divided by the median (i.e., quantile 50%) of the forest age estimates (IQR/median) (Fig. S3).
While in Europe (Fig. S3C), China (Fig. S3F) and the Eastern United States (Fig. S3E), the IQR/median estimates were relatively low, we observed high IQR/median estimates in Northern North American regions (Fig. S3E) as well as in large patches of Siberia (Fig. S3D) and the dry tropics ( Fig. S3A and Fig. S3F).
We further compared the spatial patterns of the MPI-BGC forest age dataset with a series of independent regional and global forest age products Pan et al., 2011;Poulter et al., 2019;Zhang et al., 2017) (Fig. 9, Fig.   10, and Fig. S5). In the Amazon basin, we found that the MPI-BGC forest age product depicted widespread higher forest age estimates (i.e., blue colour) than the Chazdon et al. (2016) dataset (Fig. 9A), resulting in a more extensive area of tropical old-growth forest in the MPI-BGC forest age product (Fig. 9B). On the other hand, we observed lower forest age estimates in the regions of Rio Grande Do Norte and Paraiba in the MPI-BGC forest age product (i.e., red colour). Such disagreement between the two products could be related not only to the different methods used to infer forest age (i.e., statistical method vs. age-AGB chronosequence approach for the MPI-BGC forest age and the Chazdon products, respectively) but also to the uncertainties of the RFclassifier for classifying old-growth vs. non-old-growth forests in this region ( Fig. S1 and Fig. S2). Similarly, the presented product and the Pan dataset revealed widespread discrepancies in the North American region, particularly in the Western part of the United States and the North American boreal forests (Fig. 8E). More precisely, the Pan dataset had a higher fraction of young forest patches than the MPI-BGC forest age product (Fig. 9F). Methodological differences between the Pan and the MPI-BGC forest age datasets could explain such differences. While the Pan dataset integrates forest inventories, disturbance datasets, and land-use/land cover change data to retrieve forest age estimates in the Pan dataset, the MPI-BGC forest age product relied on forest inventory, climate data and statistical methods. Additionally, forest inventory plots used to derive the MPI-BGC forest age product were relatively sparse in Canada (Fig. 1), which might limit the statistical methods used for the MPI-BGC forest age product to predict realistic forest age estimates (i.e., extrapolation issues). Finally, the forest age estimates of China's MPI-BGC forest age product were consistent with the Zhang dataset (Fig. 9C). The area distribution across age classes of the two products appeared to have a relatively good agreement in China (Fig. 9D).

Figure 9
Comparison between the forest age dataset from this study and independent forest age dataset: Amazon basin (A and B), China (C and D) and North America (E and F). For a fair comparison with the independent age datasets, the MPI-BGC forest age map used here is the one without tree cover correction applied to the AGB dataset. Differences were computed using weighted age estimates from the fraction of the decadal age classes within each 0.5-degree grid cell resolution.
We also found significant and widespread discrepancies between the MPI-BGC forest age dataset and the global forest age dataset (GFAD)   (Fig. 10). Overall, the GFAD product had higher fractions of very young and old forests ( Fig. 10B and C). Because the GFAD used a different AGB product for the pan-tropical region and mainly relied on statistics from national forest inventories for the Northern hemisphere, widespread differences were expected between the GFAD and the MPI-BGC forest age maps. The MPI-BGC forest age dataset depicted older forests in the Western part of the United States (i.e., blue colour), while it showed younger forests across Europe than the GFAD product ( Fig. 10A). Differences were also apparent in the dry tropics, where the GFAD product showed younger forests than the MPI-BGC forest age dataset, particularly in the Miombo regions. Such discrepancies could be explained either by the use of a biomass-age approach in this region or by integrating MODIS fire information in the GFAD forest age dataset. We adjusted the MPI-BGC forest age dataset with the forest age product inferred from the MCD45A1 MODerate-resolution Imaging Spectroradiometer (MODIS) fire product at 1 km resolution (Giglio et al., 2018;Poulter et al., 2019), which was used in the GFAD product. In this MODIS-age product, forest age was determined as the last time since a fire event occurred within a grid cell for 2000-2015, thereby assuming that the entire pixel was burned down. For instance, forest age within a 1 km grid cell was five years old if the last time a fire occurred within this grid cell was in 2010. The latter took precedence over the former dataset when adjusting the MPI-BGC forest age dataset with the MODIS-age product. As expected, we observed a higher fraction of younger forest in the adjusted MPI-BGC forest age dataset (Fig. S4B), particularly in regions relying on the biomass-age approach in the GFAD product ( Fig.   S4A and C). However, significant discrepancies between the two products remained when comparing the weighted average forest age estimates at the pixel level, particularly in European forests (Fig. S4A). Yet, we acknowledge that a comparison between the GFAD and the MPI-BGC forest age maps has to be taken with caution when evaluating the MPI-BGC product as substantial methodological differences exist between the two products.

Figure 10
Difference map between the forest age estimates derived from the MPI-BGC and GFAD products (A). Areas per age class were also compared between the two products for regions relying on national inventories (B) and relying on biomass-age approach (C) in the GFAD product. Differences were computed using weighted age estimates from the fraction of the decadal age classes within each 0.5-degree grid cell resolution.

Data availability
The dataset of the different forest age products presented in this study can be downloaded from the Data Portal of the Max Planck Institute for Biogeochemistry at https://doi.org/10.17871/ ForestAgeBGI.2021 (Besnard et al., 2021). For anonymous access during review, the data are available at https://nextcloud.bgc-jena.mpg.de/s/Xwt8AdkHkgL3TTc.

Conclusion
We presented a new forest age dataset derived from forest inventory, biomass, climate and remote sensing data.
Generally, the statistical model used to create the gridded age datasets had a relatively good capacity to predict forest age estimates at plot level (precision of 0.81 and 0.99 for classifying old-growth and non-old-growth, respectively and NSE of 0.6 for predicting non-old-growth forests). At the same time, biases were observed, mainly when predicting older forests (i.e., > 150 years old). The functional relationships between biomass and forest age learned by the statistical models appeared to agree with forest age theory and the role of environmental/climate in modulating the relationship. The proposed gridded datasets allowed us to assess the global patterns of forest age and provided insights into regional forest demography. For instance, relatively young-intermediate forests were observed in Europe and China, where predominant management practices and afforestation/reforestation activities. We could also demonstrate that old forests are primarily represented in very wet, warm, and cold regions. However, comparing the MPI-BGC forest age product with independent forest age datasets revealed large discrepancies, suggesting high uncertainties in mapping forest demography globally. Overall, this forest age product provides a new source of information related to disturbance history and forest regrowth, which is crucial to better understanding the location of the forest carbon sinks and sources.