As enduring carbon sinks, forest ecosystems are vital to the terrestrial carbon cycle and help moderate global warming. However, the long-term dynamics of aboveground carbon (AGC) in forests and their sink-source transitions remain highly uncertain, owing to changing disturbance regimes and inconsistencies in observations, data processing, and analysis methods. Here, we derive reliable, harmonized AGC stocks and fluxes in global forests from 1988 to 2021 at high spatial resolution by integrating multi-source satellite observations with probabilistic deep learning models. Our approach simultaneously estimates AGC and associated uncertainties, showing high reliability across space and time. We find that, although global forests remained an AGC sink of 6.2 PgC over 30 years, moist tropical forests shifted to a substantial AGC source between 2001 and 2010 and, together with boreal forests, transitioned toward a source in the 2011-2021 period. Temperate, dry tropical and subtropical forests generally exhibited increasing AGC stocks, although Europe and Australia became sources after 2011. Regionally, pronounced sink-to-source transitions occurred in tropical forests over the past three decades. The interannual relationship between global atmospheric CO2 growth rates and tropical AGC flux variability became increasingly negative, reaching Pearson's r = -0.63 (p < 0.05) in the most recent decade. In the Brazilian Amazon, the contribution of deforested regions to AGC losses declined from 60% in 1989-2000 to 13% in 2011-2021, while the share from untouched areas increased from 33% to 76%. Our findings suggest a growing role of tropical forest AGC in modulating variability in the terrestrial carbon cycle, with anthropogenic climate change potentially contributing increasingly to AGC changes, particularly in previously untouched areas. In the face of ongoing climate change, sustaining global forest AGC sinks calls for large-scale afforestation and reforestation efforts, targeted ecosystem management, strengthened deforestation control, and, critically, the preservation of intact forests. Forest ecosystems play a key role in the terrestrial carbon cycle, storing approximately 45% of terrestrial carbon and acting as long-term carbon sinks [1,2]. They also support essential ecosystem services, and are an important part of nature-based solutions in the context of anthropogenic climate change, in line with the Paris Agreement and global carbon neutrality goals [3]. Aboveground carbon (AGC) is the carbon stored in living vegetation biomass above the surface. In forests, it accounts for approximately 30-50% of total carbon pools and serves as a crucial indicator of forest function and productivity [4,5]. Recognizing its significance for global carbon budgets and climate stability, AGC is classified by the as an essential climate variable for monitoring biospheric responses to environmental change [6]. However, changes in climatic conditions and disturbance regimes, including natural ones like fire, storms, and insect outbreaks, but also anthropogenic ones like logging, mining, or environmental pollution, can shift forest AGC from being long-term sinks to sources, or vice versa [7,8]. Such long-term, persistent shifts, beyond short-term natural variability, can have profound impacts on the terrestrial carbon budget and cycle [2]. Understanding long-term AGC dynamics and decadal-scale sink-source transitions is thus essential for informing sustainable forest management and climate policy [9,10]. Recent international initiatives, such as ESA's Climate Change Initiative (CCI) and NASA's (GEDI), have improved the ability to map global and regional AGC stocks for recent years through combined field inventories and satellite remote sensing [11,12]. However, existing AGC maps typically only represent temporal snapshots or relatively short periods [13]. This makes it challenging to quantify long-term AGC dynamics, and especially the interannual variability in AGC fluxes [14]. Microwave-based vegetation optical depth (VOD) measurements have emerged as a promising proxy to empirically estimate AGC dynamics over longer periods, given the established approximately linear spatial relationship between VOD, especially at longer wavelengths such as the L-band, and vegetation biomass [15]. This empirical relationship has facilitated several recent assessments of the interannual variability in global and regional AGC fluxes [13,14,[16][17][18][19][20]. However, L-band VOD records are only available from approximately 2010 onwards, limiting the ability to assess AGC dynamics beyond the past decade. Moreover, emerging evidence indicates that existing VOD-derived AGC estimates are inherently confounded by soil moisture variations, potentially overstating biomass fluctuations because of changes in water stress [21]. We empirically find inconsistencies among existing AGC products (see Section 4.8 in Methods) even when derived from similar methods and data sources. Substantial uncertainties remain in understanding long-term AGC dynamics and sink-source shifts driven by anthropogenic climate and environmental change. There is, therefore, a critical need for a robust and internally consistent long-term AGC dataset at global scales. Here, we address these critical knowledge gaps by integrating multi-source, remotely sensed, and ground-based observations with a deep learning approach. We produce harmonized and spatially explicit global AGC estimates at 0.25 deg resolution from 1988 to 2021, along with interval predictions that capture uncertainties arising from both observational data and modeling processes. Relying on such estimates, we conduct a systematic assessment of long-term AGC dynamics, examining spatiotemporal patterns, trends, and decadal-scale shifts between carbon sources and sinks across multiple spatial scales. With a particular emphasis on tropical forests, we further explore the evolving role and interplay of natural and anthropogenic disturbances on AGC fluxes. We first categorize several recently released remotely sensed optical and microwave vegetation indices and environmental variables into dynamic (time-specific) and static (time-averaged or timeindependent) predictors of AGC, based on their temporal characteristics (Methods Section 4.1). Dynamic predictors include growing-season statistics of CXKu-band VOD, normalized difference vegetation index (NDVI), and leaf area index (LAI), as well as plant functional types (PFTs) of trees and forest cover fractions derived from land cover data. Static predictors comprise aggregated Lband VOD and photosynthetically active radiation (PAR), the global digital elevation model (DEM), and geographic coordinates (Table S1). Using probabilistic convolutional neural networks (CNNs), we then model their spatial relationships with the ESA CCI AGC [11], which is available as yearly snapshots for the years 2010 and annually from 2015 to 2021 (Sections 4.4 and 4.5 in Methods). The probabilistic CNNs are integrated with uncertainty quantification techniques to enable the estimation of AGC intervals that account for both aleatoric (data-related) and epistemic (model-related) uncertainties.
Under a space-for-time substitution framework, we applied the trained CNNs to generate continuous time series of global AGC maps for the time period from 1988 to 2021 along with corresponding uncertainty estimates, using the identified predictors as input (Fig. S1). Our harmonized AGC estimates show strong spatial agreement with reference datasets for the limited times where the latter are available, with pixel-wise correlation coefficients (Pearson's r, referred to as correlation or r hereafter) reaching up to 0.92 globally (Fig. S2). The AGC stock time series, derived by spatially aggregating pixel-level AGC density across both regional and global scales, also show broad agreement with reference datasets, with temporal correlations up to 0.70 (Fig. S3a and Table S2). This agreement is dominated by common long-term trends (usually positive due to the global land carbon sink). The interannual variability around this trend differs substantially among the reference datasets themselves, hence correlation values are generally lower and less consistent (Fig. S3b and Table S2). Notably, our AGC estimates exhibit more consistent agreement with individual reference datasets than these references do with each other, in terms of both long-term stock trajectories and interannual variability (Fig. S3 and Table S2), as well as pixel-level AGC density time series (Fig. S4). This indicates that our new dataset is more precise in capturing spatiotemporal forest AGC dynamics than existing datasets. While early studies show a strong spatial correlation between VOD and AGC [13-15, 18, 19], our comparative validation using the held-out ESA CCI AGC test data indicates that relying solely on VOD products as AGC predictors -whether CXKu-band or L-band -with the common approach of using empirical equations achieves low pixel-wise coefficients of determination (r 2  0.2-0.6; see Section 4.7 in Methods for details). In contrast, combining multi-source data as input for a nonlinear deep learning model substantially improves performance with r 2 = 0.97 (Table S3). Moreover, purely VOD-derived AGC exhibits weaker temporal consistency, showing statistically insignificant correlations (r = 0.35, p = 0.13) with aggregated reference stock time series [10], whereas our multi-source AGC estimates align well (r = 0.70, p < 0.001; Table S3). Our CNN model also outperforms other classic machine learning approaches in both spatial and temporal dimensions using the same experimental setup (Section 4.7 in Methods). For example, while Lasso and ridge regression achieve lower spatial accuracy (r 2  0.8) and temporal consistency (r  0.6, p < 0.01), random forest performs comparably in spatial terms (r 2 = 0.98) but yields less reliable temporal dynamics (r = 0.45, p < 0.05). Spatially, our estimates reveal high AGC densities in equatorial forests and lower values at higher latitudes . The absolute uncertainty (i.e., the standard deviation representing both data and model uncertainties; see Section 4.5 in Methods) is similarly elevated in dense tropical forests , possibly reflecting higher observation noise and lower predictive skill of CNNs in these areas. However, higher relative uncertainty values are more prevalent in low-biomass regions, where small deviations in estimated biomass can result in disproportionately high relative uncertainty ratios (Fig. S5). Nonetheless, latitudinal profiles of AGC stocks demonstrate strong consistency with other references (Methods, Table S4). We categorize the long-term AGC averages into four major biomes (Fig. S6b shows masks, Fig. 1d,e shows results). Moist tropical forests account for the largest share (119 PgC, 52% of total stocks) and exhibit the highest median AGC density (113 MgC ha 1 ). Temperate forests comprise 23% (52 PgC, 37 MgC ha 1 median density), dry tropical & subtropical forests contribute 15% (35 PgC, 18 MgC ha 1 median density), and boreal forests contribute 11% (24 PgC, 17 MgC ha 1 median density). These differences align with the latitudinal gradient of biomass and soil carbon distribution, where equatorial forests have a higher proportion of their total carbon in aboveground biomass (AGB) (44%), declining to 15% at higher latitudes [13,22]. Over the period 1988-2021, spatially explicit AGC densities exhibit an overall increasing trend globally, with a median value of 0.07 MgC ha 1 yr 1 . Among biomes, temperate forests show the highest median increase (0.10 MgC ha 1 yr 1 ), followed by boreal forests (0.06 MgC ha 1 yr 1 ), and dry tropical & subtropical forests (0.05 MgC ha 1 yr 1 ). In contrast, moist tropical forests show a median decline in AGC density of 0.04 MgC ha 1 yr 1 . Our more than threedecade-long time series of AGC stocks indicates that global forests have sequestrated 6.20 PgC in 1988-2021 (see Methods; Fig. 2c; Table S7). However, this overall carbon gain masks regional differences. Temperate forests contributed the largest share of the overall gain at 3.10 PgC, followed by 1.96 PgC for dry tropical & subtropical, and 1.25 PgC for boreal forests, whereas moist tropical forests acted as a weak source of -0.11 PgC overall. The widespread AGC sequestration, especially outside the tropics, may be driven by elevated atmospheric CO 2 and nitrogen fertilization [23,24]. Conversely, deforestation, degradation and climate change may counteract such benefits in the moist tropics [9,17,25,26]. We observe a substantial drop in AGC during 1991 (3.16 PgC), with 92% of this decrease occurring in tropical and subtropical biomes (2.92 PgC). Two major factors may have contributed to this decrease: (i ) the compound climate stresses caused by the El Nino/Southern Oscillation (ENSO) and the Mount Pinatubo eruption, as well as the reduction in photosynthetically available radiation, which could have caused widespread vegetation mortality [27][28][29], and (ii ) volcanic aerosols interfering with remote-sensing data, potentially introducing systematic biases into our AGC retrievals [30]. To reduce potential measurement uncertainties from satellites and ensure the robustness of our findings, we exclude the years 1991-1992 from analyses of AGC stock trends, and 1991-1993 from the analysis of interannual variability in AGC fluxes, following ref. [31]. Decadal AGC fluxes derived from our estimates reveal pronounced variability when compared across three subsequent time periods (i.e., 1988-2000, 2001-2010, and 2011-2021; Fig. 3, Fig. S7 and Table S7).
To quantify the change in AGC over a decade, we compute the difference in AGC between the last and the first year, converted to a rate by dividing by the number of years, following ref. [9,10]. To complement this measure, we also compute linear AGC stock trends over each period using the Theil-Sen estimator, providing a more robust view of long-term AGC dynamics (see Section 4.6 in Methods). Overall, global forests have remained a carbon sink: Decadal changes as defined above are 210.1 TgC yr 1 during 1988-2000, 43.0 TgC yr 1 during 2001-2010, and 196.3 TgC yr 1 during 2011-2021, with positive values indicating carbon uptake by forests. Decadal stock trends have continuously increased, with estimated rates of 105.5, 125.5, 265.9 TgC yr 1 for the three respective periods, further confirming the enduring and even strengthening AGC sink of global forests. However, moist tropical and boreal forests exhibit prominent decadal sink-source shifts . During 2001-2010, moist tropical forests shifted from a global carbon sink (90.2 TgC yr 1 in 1988-2000) to a substantial carbon source (190.6 TgC yr 1 ). This transition coincided with a period marked by repeated extreme events, such as droughts [32,33] and fires [34], as well as accelerated deforestation [35], which has likely contributed substantially to the observed carbon losses. Similar sink-to-source shifts were also observed during 2001-2010 in countries with extensive moist tropical forests, including Brazil, Indonesia, and Peru (Table S7). In the following decade (2011-2021), moist tropical forests appeared to transition toward a weak carbon source (4.7 TgC yr 1 ). This apparent partial recovery is likely due to reductions in deforestation rates [36] and regrowth in previously cleared or degraded areas [37], although the specific causal drivers of AGC change remain to be further investigated. Boreal forests transitioned from a carbon sink (9.9 TgC yr 1 and 85.8 TgC yr 1 in the first two decades from 1988 to 2010) to a weak carbon source (2.5TgC yr 1 ) in 2011-2021 (Table S7). Growing disturbance pressures, including fire, insect outbreaks, and logging [8] may have contributed to these substantial AGC losses. Spatial patterns highlight marked declines across parts of eastern Eurasian boreal zones , while some areas, such as western Siberia, show signs of recovery in 2011-2021 (Table S7). Canada's boreal forests exhibited strong decadal variability, alternating between source and sink over the three decades covered by our dataset (Table S7), likely reflecting natural variability, possibly in combination with the influence of regional disturbances [38]. While temperate forests sequestered an average of 94.0 TgC yr 1 and dry tropical & subtropical forests stored about 59.2 TgC yr 1 over 1988-2021, regions such as Europe and Australia show sink-to-source shifts accompanied by high net AGC losses. In Europe, forests have transitioned to a weak source of 10.5 TgC yr 1 in the last decade, despite retaining a small positive stock trend of 2.6 TgC yr 1 (Fig. S8). Previous studies attributed these losses to climate-related storms, pests (e.g., bark beetles), droughts, and fires [39]. Australian forests remained a net carbon source over the full period since 1988 (2.5 TgC yr 1 on average; Fig. S8), as short-term recovery was offset by major drought and fire events, particularly in 2019-2020 [40]. To estimate the contribution of different biomes to overall interannual variability in AGC fluxes, we adopt the partitioning approach developed by ref. [41], which considers both the magnitude and temporal alignment of regional flux anomalies relative to the global signal. We estimate that dry tropical and subtropical forests contributed 37% of the interannual variability in global AGC fluxes over the past 30 years (1989-2021), followed by moist tropical forests (25%), boreal forests (23%), and temperate forests (15%) (Fig. S9). These fractions indicate that dry tropical & subtropical terrestrial forest ecosystems are major drivers of the global interannual variability, consistent with earlier findings [41], and potentially linked to the vulnerability of these ecosystems to climatic impacts. However, these contributions have changed over recent decades (Fig. S9). Temperate and boreal forests collectively dominate in the first (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and third decade (2011-2021), accounting for 71% and 90% of the variability, respectively. In contrast, during the second decade (2001-2010), moist tropical and dry tropical & subtropical forests take the lead, contributing 86%. These temporal shifts may reflect changing disturbance regimes across regions. For instance, warming-related stressors and insect outbreaks may have exerted a stronger influence in high-latitude forests during the first and third decades [8], whereas tropical ecosystems in the second decade appear to have been more strongly affected by deforestation and degradation [26]. Our data reveal more pronounced decadal variability in tropical AGC compared to other regions . We further observe a strengthening negative relationship between detrended interannual AGC fluxes in these regions and the atmospheric CO 2 growth rate over recent decades, reaching r = 0.63 (p < 0.05) in the most recent decade . A 10-year moving-window analysis and uncertainty quantification via bootstrapping corroborate this increasingly negative correlation (Fig. S10). These observations may point to a strengthening role of tropical forest AGC in modulating the terrestrial carbon cycle variability. To better understand spatial patterns within the tropics, we assess the relative contributions of different regions to overall AGC variability and examine the heterogeneity of AGC dynamics. Using the flux partitioning method again [41], we find that the interannual variability in tropical AGC fluxes originates mainly from tropical America and tropical Africa (each contributing  46%), whereas tropical Asia accounts for only 7% (see bars labelled 'All', representing the full period, in Fig. 4b). Within these continents, the Amazon rainforests dominate tropical America's contribution by accounting for 69% of its interannual variability, followed by the Congo rainforests contributing 24% within tropical Africa, and the Indonesian rainforests contributing 6% within tropical Asia. These patterns underscore Amazon's pivotal role in shaping tropical AGC dynamics at the interannual scale.
Although tropical forests gained a total of 1.3 PgC from 1988 to 2021, with a trend of 2.7 TgC yr 1 , this overall balance masks pronounced regional disparities (Table S7). We find that neither tropical American forests nor the Amazon rainforests have fully returned to their 2003 AGC levels, although partial recovery occurred after several subsequent AGC losses . In contrast, AGC in tropical African forests and the Congo rainforest in particular show a relatively stable period between 1998 and 2011, and a substantial drop in 2015/16. Thereafter, their AGC trajectories have diverged: Tropical African forests largely recovered, approaching the 2014 peak, whereas the Congo Basin remained 0.3 PgC below its 2014 level . In tropical Asia, the divergence in AGC dynamics between continental-scale forests and the Indonesian rainforest is particularly pronounced. Tropical Asian forests exhibited an overall upward trend (2.2 TgC yr 1 ), reaching a peak in 2015 followed by a sharp decline in 2016, with partial recovery in subsequent years that has not returned to the peak level. Meanwhile, the Indonesian rainforests followed a long-term declining trajectory (3.6 TgC yr 1 ). The continued increase in deforestation and land-use change in the Congo and Indonesian rainforests likely contributed to their divergence from broader continental AGC trends by altering AGC trajectories and constraining forest recovery [42,43]. Previous studies have reported that certain major tropical drought events were associated with substantial carbon losses, suggesting spatial variability in forest carbon loss and recovery [44,45]. To examine such dynamics in our dataset, we compare regional AGC dynamics with moisture availability conditions quantified as the fraction of pixels with three-month (SPEI3)  1 [46] (see background colouring in Fig. 4c-i). Overall, we do not observe a clear relationship between these drought-affected area fractions and regional AGC dynamics in our results. This may be due to the confounding influence of other concurrent disturbances and recovery processes, which can also strongly affect AGC dynamics and make it challenging to isolate the impacts of individual drought events. However, for the Amazon, we specifically examine four well-documented large-scale drought events: 1997/98 [47], 2004/05 [48], 2009/10 [33], and 2015/16 [49] and observe minima in Amazon AGC stock differences for the latter three drought years. The interannual variability in AGC fluxes in tropical forests results from a complex interplay between human activities, climatic variability, and plant physiological responses to disturbances [50]. Here, we focus on the Brazilian Amazon as a key hotspot [51,52]. We leverage detailed, long-term deforestation records and maps from the Project for Remote Deforestation Monitoring in the Legal Amazon (PRODES) of the Brazilian National Institute for Space Research (INPE) [53] to further investigate the contributions of gross AGC losses driven by both natural and anthropogenic disturbances. We first classified untouched and deforested regions (see Methods; Fig. S11). In untouched areas, gross AGC losses are likely driven by natural factors (e.g., droughts, fires, storms, and tree mortality) and indirect anthropogenic impacts (e.g., deforestation-induced edge effects and precipitation alterations) [26,54]. In contrast, gross AGC losses in deforested regions are more directly associated with land-use changes, including logging, agricultural expansion, and infrastructure development. Our partitioning analysis reveals the shifting contributions of human-and nature-related factors to AGC loss over three decades . During 1988-2000, AGC losses in deforested areas accounted for 60% of the total losses, while losses in untouched regions contributed 33%. In the second decade (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), this balance changed markedly, with the deforested regions dropping to 32% and the untouched areas rising to 59%. This shift might reflect, in part, the influence of stricter deforestation-reduction policies introduced after 2004 , alongside frequent climate anomalies. In the final decade (2011-2021), although deforestation rates began to rise again after 2012 , the average annual deforested area (7,681 km 2 yr 1 ) remained markedly lower than in the preceding two decades (16,959 km 2 yr 1 during 1988-2000 and 16,531 km 2 yr 1 during 2001-2010). This factor likely contributed to a slight decrease in the total magnitude of AGC gross loss compared to the previous decade (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). In 2011-2021, AGC losses in deforested areas fell to a 13% share, whereas AGC losses in untouched regions rose to 76%, indicating a growing role of other disturbances, and potential environmental stress induced by anthropogenic climate change, on untouched forest carbon loss. Spatially explicit maps of decadal gross AGC density loss also corroborate these patterns . From 1988 to 2000, losses strongly coincided with the so-called "Arc of Deforestation" (Fig. S11d), implicating direct human clearing as the principal driver. However, in 2001-2010, losses spread across both deforested and untouched regions, probably reflecting a more widespread effect of disturbances. By 2011-2021, the primary loss occurred in untouched forests, and the deforestation arc signature weakened in the gross-loss maps. These northward-shifting spatial patterns closely correspond to the spatial distribution of degradation drivers identified in previous research [26]. Our results underscore a crucial transition: Despite successful reductions in direct deforestation, both natural factors and indirect anthropogenic impacts appear to play an increasingly prominent role in gross AGC losses. A repeated subset analysis focused on the Brazilian Amazon further supports this emerging pattern, revealing a similarly strengthening negative correlation between interannual AGC fluxes and global CO 2 growth rates (Fig. S12), particularly in untouched forests. Understanding long-term AGC dynamics and their spatiotemporal patterns in natural forests is vital for understanding their carbon sequestration potential under anthropogenic forcing. Although there is no precisely observed "ground truth" for AGC storages and fluxes [19], the availability of stateof-the-art vegetation indices and environmental datasets offers valuable complementary information across different datasets [55][56][57][58][59][60]. Leveraging these multi-sensor data, our probabilistic CNN-based deep learning framework produces long-term, harmonized, and reliable AGC estimates at high spatial resolution, which align closely with available benchmarks in the shorter periods they are available. While the spatial and temporal consistency of our estimates varies by region when assessed against shorter-term reference datasets, our dataset exhibits substantially improved temporal agreement with other references compared to the most recently published product by ref. [61] (see Figs.
S3 and S13, and Table S5). Our global AGC dataset enables a more robust assessment of forest carbon dynamics over the past three decades, extending the temporal duration of global AGC estimates by 10-20 years. Our spatially explicit AGC dataset complements the recent work by ref. [2], which assessed decadal carbon sinks in global natural forests over the past three decades based on compiled forest inventory data. In particular, our AGC estimates enable the identification of critical regions undergoing sink-tosource transitions. Furthermore, the detailed interannual variability in AGC fluxes facilitates a deeper exploration of forest carbon dynamics driven by climatic and anthropogenic disturbances. Conversely, inventory-based carbon estimates in ref. [2], despite their limitations in accuracy for tropical regions due to sparse inventory coverage, provide a comparative benchmark for assessing overall changes in the natural forest carbon sink. Integrating both results can enhance our understanding of global forest carbon dynamics and their underlying drivers. Another recent study reported that the majority of terrestrial carbon accumulation between 1992 and 2019 occurred in nonliving organic matter pools, with gains in living biomass accounting for less than 5% of the total, corresponding to approximately 1 PgC [20]. In contrast, we estimate that forest living biomass increased by approximately 9.5 PgC over the same period. Assuming an AGB-tobelowground biomass ratio from ref. [62], this would account for roughly 25-30% of the total terrestrial carbon gain. Consequently, while our study focuses specifically on forested areas rather than entire terrestrial ecosystems, our findings suggest that although some imbalance likely exists, the contribution of living biomass to the overall carbon sink may be much larger than previously reported. One potential source of discrepancy may stem from the data harmonization approach. In particular, the AGC products in ref. [20] were, following typical procedures, harmonized through linear interpolation and averaging of existing datasets, which may introduce substantial uncertainty. Our analysis highlights that existing AGC reference datasets can differ considerably (see Sections 2.1 and 4.8). Even products based on similar input data and methodologies may exhibit divergent changes or trends; for example, the datasets from ref. [14] and ref. [61], both based on L-VOD and empirical approaches, show substantial inconsistencies in AGC dynamics and interannual variability (Fig. S13 and Table S5). Given these inconsistencies, caution is warranted when interpreting long-term biomass carbon estimates harmonized from established AGC datasets; our modeling approach, based on independent predictors, reduces such uncertainties. The long-term AGC dataset we provide here offers an opportunity to assess global AGC dynamics across biomes, spatial scales, and time periods, particularly for tropical regions. While global forests remain a long-term net AGC sink, moist tropical and boreal forests have shown notable sink-to-source transitions in recent decades. Temperate and dry tropical & subtropical forests generally exhibit steady AGC gains, though regions such as Europe and Australia have shifted toward net carbon sources in the last decade. The long-term AGC flux trends also highlight the strong contribution of dry tropical forests to interannual carbon fluctuations [41]. However, we note that the dominant contributors shifted over time, from high-latitude forests in the first (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) to tropical and subtropical forests in the second decade (2001-2010) and back to high-latitude forests in the third decade (2011-2021). These variations may reflect shifts in disturbance regimes, including changes in deforestation intensity and the heterogeneous impacts of evolving climatic conditions and extreme events across different forest biomes. Interestingly, our analysis reveals an intensifying negative correlation between annual tropical AGC fluxes and atmospheric CO 2 growth rate (Fig. S10). This contrasts the CO 2 fertilization assumption and highlights the growing importance of tropical AGC dynamics in shaping the terrestrial carbon cycle variability. Our results are consistent with an earlier finding showing a strong association between tropical AGC fluxes and CO 2 growth rates during 2011-2017 [14]. Moreover, our results also provide independent support for previous exploratory evidence of a strengthening negative AGC-CO 2 correlation in tropics from 1989-2003 to 2002-2016 [31]. While the authors acknowledged uncertainties in their VOD-derived AGC estimates -particularly the sensitivity of VOD signals to soil moisture [21] -our analysis, based on newly developed AGC estimates, offers continuous decadal-scale evidence from 1989 to 2021, thereby reducing some of the technical limitations in their study and extending the observed relationship. This earlier study [31] further interpreted the intensifying negative AGC-CO 2 correlation as a reflection of strengthened water-carbon coupling in the tropics, likely associated with ENSO teleconnections, based primarily on a comparison between 1960-1989 and 1989-2018. However, our decadal-scale analysis from 1988-2021 does not reveal a consistent relationship between ENSO phases and AGC dynamics (Fig. S14), nor a clear pattern of intensifying negative water-carbon coupling (Fig. S15a). Both the direct correlation between water availability and AGC flux, and the partial correlation after accounting for temperature effects, also show no clear strengthening trend over the 1989-2021 period (Fig. S15b,c). Notably, the decade 2001-2010 exhibits a marked shift in the coupling structures between water availability and both atmospheric CO 2 growth rate and AGC flux. This anomaly in 2001-2010 may be related to more moderate but frequent climatic fluctuations, as well as accelerating deforestation [63]. Our findings suggest that attributing the intensifying negative AGC-CO 2 correlation solely to ENSO-induced water stress is likely insufficient. In addition, improvements in observational data quality over time, such as changes in satellite retrieval algorithms or sensor characteristics, may have introduced additional sources of uncertainty into the observed AGC-CO 2 coupling. These considerations underscore the need for further investigation to disentangle the relative contributions of climatic and anthropogenic drivers. In contrast to previous studies [13,14,45], our results do not show a clear impact of environmental conditions like well-documented droughts on the interannual tropical AGC dynamics, except for the Amazon to some extent. At the pixel level, the relationship between AGC flux and the SPEI is generally weak, with a median correlation coefficient of 0.09 across the pan-tropical forests (Fig. S16). However, stronger positive correlations are observed in dry tropical and subtropical forests, likely due to their higher sensitivity to water availability [41,64].
The overall weak drought fingerprint observed in AGC dynamics may result from several factors: First, interannual variability and long-term stock trajectories in tropical AGC differ across datasets (Fig. S13 and Table S5). Hence, the precision of AGC estimates may still not be high enough to observe drought responses sufficiently well, owing to ongoing uncertainties in existing datasets. In particular, recent studies have shown that purely VOD-derived AGC estimates may be influenced by changes in water stress, which could overstate biomass fluctuations during water scarcity periods [21,31]. Second, AGC changes are influenced by environmental conditions in complex ways, and affected by multiple local disturbance factors, and by tree regrowth and reforestation. For example, drought is associated with lower cloud coverage and hence higher downwelling short-wave radiation, which may in some cases even increase gross primary productivity and hence AGC [65]. Consequently, further research is needed to clarify the extent to which drought events influence AGC storage. Anthropogenic disturbances also modulate the vulnerability of forest AGC sinks. To investigate compound disturbances, we specifically analyzed the Brazilian Amazon, quantifying changes in the contributions of deforested and untouched areas to gross AGC losses over three decades. Our findings reveal a striking transition: Over the three decades, the share of AGC losses from deforested areas has declined from 60% to 32% to 13%, while the share from untouched regions has increased from 33% to 59% to 76%. This pattern may reflect a reduced relative role of direct deforestation and a rising influence of natural disturbances and indirect and spatially remote anthropogenic effects, such as changing climate variability, heat waves, droughts, flooding, reduced moisture recycling across the Amazon basin [66,67], and other impacts of anthropogenic climate change. Our multi-decadal perspective extends the 2010-2019 degradation-dominant pattern reported by [17] and puts it into a longer-term context, revealing that their observed shift constitutes the intensification phase of a three-decade trajectory. The substantial decline in deforestation-related losses (-78%) and concurrent increase in losses from untouched forests (+130%) over three decades may point to a growing influence of climate-related disturbances on AGC losses in the Brazilian Amazon, although this interpretation remains subject to further attribution. While our analysis cannot isolate specific drivers, these findings, together with broader patterns observed at the pantropical scale, are broadly consistent with the hypothesis that anthropogenic climate change plays a growing role in elevating the vulnerability of tropical forests [68]. Resulting shifts from carbon sinks to sources in these ecosystems add to global greenhouse gas emissions and, in turn, accelerate global warming. Our findings regarding the South American Amazon highlight the importance of addressing both human and climate-related drivers in forest carbon management. The region's shift to a substantial AGC source of 180.3 TgC yr 1 in the 2000s illustrates how overlapping deforestation and climate stress can severely weaken carbon sequestration. Notably, we also find that AGC in tropical America and the Amazon has not returned to the 2003 level, which may be driven by the combined effects of climatic anomalies and human-induced forest loss. This pattern is also evident in the dataset from ref. [10] (Fig. S13), although it was not explicitly analyzed or reported in their study. While post-2004 deforestation controls have helped reduce carbon losses specifically associated with forest clearing, the continued rise in AGC losses from untouched forests underscores growing vulnerability to climate extremes and indirect anthropogenic effects. These results suggest that, in addition to reducing deforestation, there is a need for more substantial mitigation efforts and integrated policies that prioritize the protection and management of untouched forests. Targeted adaptation strategiessuch as climate-informed afforestation and restoration -will be essential for sustaining or restoring the carbon sink function of tropical forests under a changing climate. We note that we specifically focus on global forest AGC, not total biomass carbon. Total biomass carbon, including both above-and below-ground components, is commonly estimated using biomespecific root-to-shoot ratios [10,19] or empirical power functions [18]. These approaches introduce additional uncertainty into carbon flux estimates, particularly since the interannual variability in the relationship between AGB and below-ground biomass remains poorly understood, and responses of AGB and below-ground biomass to climate and anthropogenic disturbances can also differ [62,69]. While a detailed attribution of AGC changes is beyond the scope of this study, we interpret observed AGC dynamics in the context of existing literature, and link them to both climatic and anthropogenic influences. Future studies explicitly quantifying the impacts of extreme climatic events and human disturbances using our long-term AGC estimates would provide valuable insights into the evolving drivers of global forest carbon dynamics. 4 Methods We used multi-source vegetation indices and environmental variables to establish spatially explicit aboveground carbon (AGC) estimates from 1988 to 2021. All datasets were harmonized to 0.25 deg spatial resolution (ca. 25 km at the equator) with annual temporal resolution (see Table S1 and Fig. S1 for detailed data information and processing steps). AGC is normally derived from AGB using empirical coefficients (e.g., 0.5) [16]; therefore, the following content will focus on the estimation of AGB. To establish the relationship between AGB and multiple vegetation indices and environmental variables, we classified all predictors as either dynamic or static based on their temporal extent or variability. Dynamic predictors are time-dependent and hence capture interannual fluctuations when extrapolating the derived spatial relationship to other years, whereas static predictors represent baseline environmental properties. As dynamic predictors, we selected the recently released VODCAv2 CXKu-band VOD [55], the PKU GIMMS NDVI [56], and GIMMS LAI 4g [57], and additionally derived PFTs of trees and forest cover fractions from ESA CCI Land Cover data [70]. Passive microwave-based (CXKu-band VOD) and optical-based (NDVI, LAI) indices capture complementary aspects of vegetation dynamics, providing a more comprehensive characterization of AGB changes [71]. Meanwhile, the ESA CCI products provide long-term structural information on spatio-temporal forest variations. Since CXKu-band VOD, NDVI, and LAI can have missing observations over frozen or snow-covered surfaces in northern latitudes during winter, we extracted growing-season values to improve representativeness.
However, we acknowledge that this approach may introduce seasonal bias, particularly in regions where vegetation dynamics outside the growing season play a non-negligible role. To maintain consistency across all three vegetation indices, we adopted fixed-month windows for latitudinal regions following ref. [72], instead of defining growing seasons per grid cell. We calculated the annual mean values of CXKuband VOD, NDVI, and LAI during the growing season, along with the 95% percentile of CXKu-band VOD, as it exhibits greater sensitivity to interannual carbon dynamics [73]. This combination also helps balancing the contribution of passive microwave and optical data in modeling. We employed the user tool of ESA CCI Land Cover to generate PFTs of trees at 0.25 deg resolution [58], and then derived forest cover fractions by summing relevant forest-related land-cover fractions [74]. Because the derived annual LAI means (1988-2020) did not extend to 2021, we used an autoregressive integrated model to project them by one additional year. Specifically, we implemented an ARIMA(1,1,0) model for each spatial grid cell using annual time series data. For each pixel, if the annual LAI series contained no missing values, we fit the ARIMA model to the 1988-2020 data and forecasted a single value for 2021. This process was applied independently to each grid cell. In cases where the time series contained missing data or the model failed to converge, the 2021 value was set to nodata. Similarly, to align PFT and forest cover fraction data (originally from 1992-2021) with our study period of 1988-2021, we applied the same model-based extension for the earlier years. This short temporal extrapolation was used solely to ensure temporal alignment across model inputs and is not intended for independent interpretation. Moreover, since other input variables (e.g., VOD and NDVI) have sufficiently long time series, they provide complementary information that helps constrain the model during training, thereby mitigating potential uncertainties introduced by this extrapolation. As static predictors, we employed the VODCAv2 L-band VOD, the MODIS PAR product (MCD18C2) [59], the ETOPO 2022 Global DEM [60], and geographically encoded coordinates, where longitude was transformed using sine and cosine functions to preserve its cyclical nature. L-band VOD is only available from 2010 to 2021, but it offers stronger penetration capability than higher-frequency microwave products (e.g., CXKu-band VOD) and exhibits a stronger relationship with AGB [15]; hence, we treated it as a static predictor for capturing aspects of vegetation density. Likewise, we used PAR data, available from 2000 to the present, to represent radiation availability and the photosynthetic conditions of global regions. We aggregated the available L-band VOD (growing-season) data to compute total means and 95% quantiles, while we used total means for PAR. The ETOPO 2022 DEM was included as a static predictor due to its near-constant nature; despite being a 2022 product, large-scale topography generally remains unchanged over decadal timescales. Finally, we applied sine and cosine transformations to longitude to avoid discontinuities at the 180  W/180  E boundary. All aboveground biomass (AGB) reference data used in this study are publicly available and either originally provided at, or aggregated to, a spatial resolution of 0.25 deg -matching the resolution of predictors. These AGB reference maps were also converted to AGC by applying a fixed conversion factor of 0.5, assuming that approximately 50% of AGB is composed of carbon [16]. Table S4 details the data sources. We employed the ESA CCI AGB maps (version 5.01) as our primary target in model development because they integrate synthetic-aperture radar (SAR) backscatter and spaceborne light detection and ranging (LiDAR) metrics with globally consistent algorithms, yielding greater reliability in woody biomass estimation than many model-based AGB products [11]. ESA CCI provides eight AGB maps (one for 2011, plus annual maps from 2015 to 2021) accompanied by standard deviation layers. However, these maps vary in quality due to observational differences: the 2017-2021 annual maps are considered equally reliable and supersede any earlier releases [75]. Moreover, according to ESA CCI's quality flags, over 95% of forested regions exhibit changes labeled as improbable or unrealistic (see Fig. S17). Consequently, we used only the individual ESA CCI AGB maps themselves, excluding any temporal changes within them from our modeling framework. Additional AGB references include global annual maps for 1993-2012 by ref. [13], a global forest static map (ca. 2000) by ref. [76], pan-tropical annual maps for 2010-2017 by ref. [14], global annual maps for 2000-2019 by ref. [10], global annual maps for 2011-2023 by ref. [61], and a global static map (ca. 2020) by GEDI L4B [12]. The products from refs. [13,14,61] use VOD-based empirical methods (with Ku&X-band VOD in ref. [13], and L-band VOD in the latter two), whereas the dataset from ref. [76] merges three earlier maps [77][78][79]. The total living biomass maps by ref. [10] are derived from microwave and optical remote-sensing imagery with a machine learning model; we converted them to AGB by applying the pixel-wise AGB-to-belowground biomass ratio from ref. [62]. Additionally, GEDI L4B provides a static, relatively sparse global AGB sampling around 2020 (covering missions from week 19 to week 138) [12]. We used AGC maps converted from them to assess the spatial and temporal reliability of our estimates. Furthermore, the temporally resolved datasets from refs. [10,13,14,61] allow us to validate our findings on AGC changes across scales discussed in the main text, given their independent time-series coverage and prior validation in other studies. To ensure that AGC change analyses are confined to consistently defined long-term forested regions with a spatial resolution of 0.25 deg over the past three decades, we followed the forest-masking method described by ref. [74]. We derived an annual forest cover fraction from the ESA CCI Land Cover data (Section 4.1) and calculated its 5% quantile across all periods. We then retained only those pixels with quantile values exceeding 20% forest cover [74]. The resulting forest mask is shown in Fig. S6a. This approach ensures consistency with other ESA CCI-based predictors, such as PFTs and forest cover fractions, as well as ESA CCI AGB. The biome masks are consistent with ref.
[10], which classified global vegetation biomes (i.e., moist tropical, dry tropical & subtropical, temperate, and boreal biomes) using the MODIS IGBP (International Geosphere-Biosphere Programme) land cover product from 2001. To ensure consistency with our global forest area, we aggregated the biome mask to a 0.25 deg resolution and used it to intersect with our global forest regions, identifying four distinct forest biomes (Fig. S6b). For forest grid cells that did not overlap with any biome, the nearest biome class was assigned. The deforested regions are derived from the INPE PRODES, which has been conducting satellitebased monitoring of clear-cut deforestation at approximately 30-meter resolution in the Legal Amazon since 1988, producing annual deforestation rates for the region. PRODES aggregated all deforestation maps from 1988 to 2007 into a single layer, while maps from 2008 onward are provided annually. To create a long-term deforestation map, we combined deforestation data from 1988 to 2021. The 30meter resolution maps were aggregated to a 0.25 deg grid, and areas with a deforestation ratio of 10% or higher were classified as deforested (Fig. S11c). The untouched region boundary was obtained from the Intact Forest Landscape (IFL) project; we used the 2020 map, extracted the Brazilian Amazon region, and converted it to a 0.25 deg grid (Fig. S11b). Other masks include the regions of countries with the largest forest areas, as reported by the Food and Agriculture Organization of the United Nations (FAO), which were obtained from the Database of Global Administrative Areas (GADM). The pan-tropical region is defined as the area within approximately 23.5 deg north and south of the equator. The boundaries of the Amazon Basin and Congo Basin were sourced from ArcGIS Online ( html?id=f2c5f8762d1847fdbcc321716fb79e5a for Amazon Basin and item.html?id=23ed36037ac14b479dd787d7fd418b9e for Congo Basin), while the boundary of Indonesia was obtained from GADM. Rainforest areas within these regions were extracted by intersecting them with the moist tropical forest mask. All masks were harmonized to a 0.25 deg grid (Fig. S6c-e). We adopted a space-for-time substitution approach (see Fig. S1), given the scarcity of fine temporal resolution AGC data [16]. Rather than relying on empirical VOD-AGB relationships [13,14,[16][17][18][19] or classic machine-learning approaches (e.g., random forest) that capture nonlinear relations on a per grid cell basis [10], we developed a probabilistic CNN. This CNN simultaneously learns the relationship between predictors and AGB targets, as well as the spatial interdependence among neighboring pixels (i.e., texture features). Our network architecture follows ref. [80], which was successfully applied in global canopy-height regression, featuring residual blocks with separable convolutions to stabilize training and enhance computational efficiency. We employed learnable 3x3 convolutional filters without pooling to retain critical information from limited input datasets. We extracted 15x15-pixel patches from all predictor variables, and the corresponding ESA CCI AGB maps. At a spatial resolution of 0.25 deg , each patch spans approximately 3.75 deg x3.75 deg , enabling the CNN to capture regional interconnections across space. Because some pixels in these patches may lack data (e.g., from non-forest areas or missing observations), we excluded these pixels from the loss function computation. Given inaccuracies in ESA CCI's annual AGB changes (Fig. S17), we trained separate CNNs for each year from 2015 to 2020 (using that year's predictors and ESA CCI AGB), then averaged their outputs to form an ensemble prediction. We reserved the 2010 and 2021 data for testing. Further hyperparameter details are provided in Table S6. Finally, we applied the ensemble models to predictor patches spanning 1988-2021 to generate spatially explicit AGB estimates, which we converted to AGC via a fixed 0.5 factor. It is also important to note that the observational data used in our modeling contain spatially and temporally heterogeneous gaps (missing data), primarily due to limitations of optical and passive microwave remote sensing, such as cloud cover, data availability, and sensor-specific constraints. Although aggregating observations over the growing season helps mitigate much of the missing-data issue, some gaps still persist, and the AGC estimates presented exhibit similar patterns of missing values. Such data gaps are common in existing AGC datasets derived from remote sensing [14,19,61]. This may lead to a relative underestimation of the total AGC stock in corresponding regions when conducting spatially explicit analyses. However, since the primary focus of this study is on AGC changes, trends, and fluxes, such limitations do not affect the main findings. To ensure the reliability of temporal analyses, we excluded grid cells with missing values across time. This masking approach guarantees that AGC change and flux estimates are based on consistent spatial coverage without the influence of data gaps. Predictive uncertainty in deep learning can be categorized into aleatoric (data-related) and epistemic (model-related) uncertainties [81]. In our context, aleatoric uncertainty stems from noise in both predictor variables and the AGB target, whereas epistemic uncertainty reflects the CNN's limitations in modeling, and the limited information content of the input data. We employ deep ensembles with a weighted Gaussian negative log-likelihood (NLL) loss (Equation 1) to jointly capture these two sources of uncertainty [80]. Specifically, at each pixel i, we represent the output of model  as a conditional Gaussian distribution over AGB target y i , given the input data x i . We estimate both the mean  (x i ) and variance 2  (x i ) of this distribution [81,82]. The weight w i corresponds to the ratio between the ESA CCI AGB value and its associated standard deviation at pixel i, thereby emphasizing pixels with lower relative uncertainties (i.e., higher reliability in ESA CCI AGB). To promote stable convergence, we begin training with a weighted mean-square-error loss and then fine-tune with the weighted Gaussian NLL loss: We trained an ensemble of five CNNs for each year between 2015 and 2020 (each initialized randomly). At inference time, we thus ultimately obtained M = 30 CNN realizations of AGB means and variances. Following ref. [82], we combined these outputs as a mixture of Gaussians to produce the final mean  * and variance 2 * for each pixel i (Equations 2 and 3).
For the temporal analysis of annual AGC stock trends and changes, we quantified uncertainties via a Monte Carlo approach. For each pixel and each year, we sampled 1,000 times from the final per-pixel Gaussian distribution ( * and 2 * ) to generate an ensemble of AGC estimates. Summing over all pixels for each sampled realization yields a distribution of total (or regional) AGC from which we derive the mean and the 95% uncertainty interval (the 2.5% quantile as the lower bound and the 97.5% quantile as the upper bound). To determine whether global or regional vegetation acts as a sink or source, we primarily rely on the sign of the AGC flux estimate. We consider the source or sink to be insignificant in two cases: (1) when the 95% uncertainty interval of the flux includes zero, or (2) when the sign of the long-term AGC trend differs from that of the flux. To quantify AGC dynamics, we computed a set of complementary metrics at pixel, regional, biome, and global scales. Definitions and calculation methods are as follows: 1. AGC density. Each pixel in our estimates represents AGC density, expressed in MgC ha 1 . 2. AGC stocks. To obtain AGC stocks at larger spatial scales (e.g., region, biome, global), we aggregated pixel-level AGC density by multiplying each pixel's density by its corresponding area (accounting for latitudinal weighting) and summing the total AGC across all relevant grid cells, expressed in PgC (1 PgC = 10 3 TgC = 10 6 GgC = 10 9 MgC). 3. AGC net change. The net change over a specific period was defined as the difference between AGC in the first and last years of the time window. To ensure comparability across different time spans, we reported net change as the annual mean value. Spatially explicit net change maps were calculated from AGC density maps and expressed in MgC ha 1 yr 1 , while regional and global net changes were obtained by aggregating spatial maps, with results reported in TgC yr 1 . In a given region and time period, a positive AGC net change is interpreted as a carbon sink, while a negative net change indicates a carbon source. 4. AGC gross loss. If a pixel showed a decrease in AGC over the period, the magnitude of the loss was recorded as gross loss. Gross loss at the regional or global scale was computed by summing the negative components of the net change map, and is expressed in TgC yr 1 . A spatially explicit gross loss map was also derived, representing the average loss per pixel over the period, with units of MgC ha 1 yr 1 . Gross gain was calculated using the same approach but retaining only the positive components. In this study, we primarily focused on gross loss, as it is more directly related to the impacts of natural and anthropogenic disturbances. 5. AGC trend. We estimated long-term trends in AGC density or stocks using the Theil-Sen slope estimator and assessed their statistical significance with the Mann-Kendall test. Unlike net change, which only considers the first and last years, the trend reflects continuous temporal patterns across the full period. Pixel-level trend maps are reported in MgC ha 1 yr 1 , while aggregated trends are expressed in TgC yr 1 . 6. AGC flux. AGC flux represents the annual difference in AGC stocks at regional or global scales. Positive values indicate a gain in AGC stock, while negative values indicate a loss. Fluxes are reported in TgC yr 1 or PgC yr 1 . This metric captures the year-to-year variability in net carbon balance and is particularly relevant for assessing interannual responses to climatic or disturbance events. To evaluate our model's predictive performance, we compared it against both VOD-AGB empirical equations (arctan or logistic formulas derived from earlier studies [13,14,[16][17][18][19]) and classic machinelearning methods (ridge regression, lasso regression, and random forest). All models were trained using the same procedure as our CNN (i.e., trained annually and then ensemble-averaged). We validated their pixel-wise accuracy against held-out ESA CCI AGB maps in 2010 and 2021 and assessed the interannual plausibility of global AGC stock trends using the dataset from ref. [10]-the longestrunning records available. Table S3 and Fig. S18 illustrate that purely empirical VOD-AGB methods perform sub-optimally (r 2  0.2 and RMSE  70 MgC ha 1 for CXKu-band), whereas L-band VOD-based methods yield moderate improvements (r 2  0.64 and RMSE  47 MgC ha 1 ), due to greater signal penetration of woody vegetation. We also observe that VOD-derived AGC tends to be highly sensitive to climate-driven hydrological anomalies, showing abrupt and large-amplitude changes (Fig. S18), suggesting that the VOD-derived AGC products are likely to be also dominated by water stress [21]. Incorporating multiple predictors provides notable gains, even for linear regressions (r 2  0.8; RMSE  35 MgC ha 1 ). Nonlinear machine learning models further improve prediction: random forest achieves r 2 = 0.98 with an RMSE = 10.44 MgC ha 1 , and our CNN attains r 2 = 0.97 and RMSE = 12.64 MgC ha 1 . While the random forest model slightly outperforms the CNN in terms of pixel-wise accuracy when fitting spatial relationships on the ESA CCI AGB test sets, the CNNderived global AGC trend exhibits stronger temporal consistency with the reference from ref. [10] (r = 0.70, p < 0.001), whereas random forest shows a weaker association (r = 0.45, p < 0.05). Crucially, the CNN-derived trend reflects the established persistence of global forest carbon sinks [2], whereas random forest produced a probably implausible decline in the 2000s (see Fig. S18). We further investigated predictor importance using integrated gradients, an explainable artificial intelligence technique. Fig. S19 shows that dynamic variables collectively contribute more than 50% to the final AGC estimates, with LAI having the largest single impact (16.8%). Among plant functional types (PFTs), needle-leaved evergreen trees contribute 9.51%, while other vegetation indices and forest cover fractions each account for around 5%.
Notably, photosynthetically active radiation (PAR), an otherwise static predictor, captures 10.84% importance, reflecting the influence of available light on forest biomass. Both CXKu-band and L-band VOD collectively contribute a moderate share to the final predictions (7.39% and 11.05%, respectively), indicating that while VOD provides valuable information, relying on it alone would overlook other essential vegetation and environmental signals. Overall, these findings underscore the importance of integrating multiple sensor-based vegetation indices and environmental factors to achieve robust, large-scale forest AGC modeling. We evaluated both the spatial and temporal reliability of our AGC estimates by comparing them against multiple reference datasets (Section 4.2, Figs. S2, S20 and S3). To assess spatial consistency, we analyzed pixel-wise Pearson correlations between our AGC estimates and various reference products on overlapping time periods, at global and biome scales (Figs. S2 and S20). Across all global grids, correlation ranges from 0.70 to 0.87 and reaches 0.92 in tropical forests. When startified into four biomes, correlations decrease but remain high, with maximum values of r = 0.78 in moist tropical forests, r = 0.67 in dry tropical & subtropical forests, r = 0.77 in temperate forests, and r = 0.78 in boreal forests. All correlations reported here are statistically significant with p  0. Lower correlations occur with the Liu AGC and GEDI L4B AGC datasets, likely reflecting biases from high-frequency Ku&X-band VOD and empirical methods in Liu AGC (Section 4.7), as well as the sparse distribution of GEDI sampling, which can introduce systematic errors when aggregated to 0.25 deg grids. We also compared the reference datasets against each other (Fig. S21, finding that correlations generally decline from the global scale to individual biomes, which underscores inconsistencies in data acquisition and processing methods across references. We then examined temporal agreement by comparing our AGC time series (carbon trajectories and annual fluxes) to the reference products at global and biome scales (Fig. S3). Globally, our AGC stock trends closely match Xu (2000-2019; r = 0.70, p < 0.001), yet diverge from Boitard (2011-2021; r = 0.01, p > 0.05) and Liu (1993Liu ( -2012; r = 0.36, p > 0.05), which shows a biologically implausible decline in global forests [2,19]. Biome-specific comparisons also exhibit good alignment with reference datasets, with statistically significant high correlation coefficients (see Fig. S3a): for instance, r = 0.64, p < 0.01 in moist tropical forests (1993-2012), r = 0.85, p < 0.001 and r = 0.84, p < 0.001 in dry tropical & subtropical forests (2000-2019 and 2011-2021), r = 0.70, p < 0.001 and r = 0.92, p < 0.001 in temperate forests (1993-2012 and 2000-2019), and r = 0.61, p < 0.01 and r = 0.59, p < 0.01 in boreal forests (1993-2012 and 2000-2019). The comparisons of annual fluxes show lower correlation and less consistency across regions, with r > 0.5 in most biomes except boreal forests (Fig. S3b). We note that correlations among the time series of various references themselves reveal large discrepancies in stock trajectories and annual fluxes. Our AGC estimates exhibit more consistent agreement with individual reference datasets than the references do with each other, both in terms of long-term stock trajectories and interannual variability (Fig. S3), as well as pixel-level AGC density time series (Fig. S4). The space-for-time substitution strategy employed -training CNNs on ESA CCI AGC data from 2015 to 2020 and applying them retrospectively -assumes the spatial relationships learned by the model remain consistent over time. However, changes in climate or land-use patterns could potentially weaken this assumption, particularly for earlier periods. Detecting such temporal mismatches remains inherently challenging. Fortunately, our CNN framework integrates uncertainty quantification, enabling explicit measurement of out-of-distribution predictive uncertainty [82]. If substantial changes in the underlying spatial-temporal relationships had occurred, we would expect significantly wider uncertainty intervals in earlier periods. However, our analyses show no notable temporal variation in uncertainty estimates, indicating that the spatial relationships captured by the CNN have remained relatively stable over the past three decades. Overall, the spatiotemporal validation indicates the high reliability of our AGC estimates. Our comparative analysis suggests that systematic inconsistencies among existing reference datasets can introduce considerable uncertainty when these products are merged for long-term AGC analyses. This underscores the importance of cross-comparison and underscores the need for methodological harmonization in global and biome-level biomass assessments. Acknowledgements. We thank and Sung-Ching Lee for their suggestions and comments. Z. Acknowledgements. We thank and Sung-Ching Lee for their suggestions and comments. Z. Author contribution. Z.Q, S.B. and N.B. conceived and designed the study. Z.Q. performed the data collection and computations and analyzed the results. All authors discussed the results. Z.Q. wrote the paper with contributions from all authors. Competing interests. The authors declare no competing interests. Supplementary information. Supplementary information is available for this paper. Data availability. All predictors and AGC reference data are publicly available from earlier studies, the sources of which are described in Methods and Supplementary Information. The sources of all masks used in this study are provided in Methods. (SPEI) data are publicly available at The annual mean global carbon dioxide growth rates are publicly available at gr.html. The long-term data on deforestation areas in the Brazilian Amazon are publicly available at The temperature and precipitation data are from the Climatic Research Unit (CRU), publicly available at uk/groups-and-centres/climatic-research-unit. The AGC estimates produced in this study will be publicly available upon acceptance of this manuscript for publication. Code availability. The code for modeling and analysis will be publicly available upon acceptance of this manuscript for publication. Author contribution. Z.Q, S.B. and N.B. conceived and designed the study. Z.Q. performed the data collection and computations and analyzed the results. All authors discussed the results. Z.Q. wrote the paper with contributions from all authors. Competing interests. The authors declare no competing interests. Supplementary information. Supplementary information is available for this paper. Data availability. All predictors and AGC reference data are publicly available from earlier studies, the sources of which are described in Methods and Supplementary Information.
The sources of all masks used in this study are provided in Methods. (SPEI) data are publicly available at The annual mean global carbon dioxide growth rates are publicly available at gr.html. The long-term data on deforestation areas in the Brazilian Amazon are publicly available at The temperature and precipitation data are from the Climatic Research Unit (CRU), publicly available at uk/groups-and-centres/climatic-research-unit. The AGC estimates produced in this study will be publicly available upon acceptance of this manuscript for publication. Code availability. The code for modeling and analysis will be publicly available upon acceptance of this manuscript for publication.