Large live biomass carbon losses from droughts in the northern temperate ecosystems during 2016-2022

Domain: earth_science
Published: 2025-06-02
Editorial outcome: minor_revision

Abstract
Northern ecosystems (≥30°N) have been accumulating vegetation biomass carbon in recent decades, but increasing droughts and wildfires threaten this carbon sink. Here, we analyse annual changes in live vegetation biomass in northern ecosystems using low-frequency microwave satellite observations at 25 km spatial resolution from 2010 to 2022. We find that live biomass carbon stocks have undergone a reversal from a positive to a negative trend during the study period with 2016 marking the turning point. During 2016–2022, live biomass carbon stocks decreased at a rate of −0.20 ± 0.11 PgC yr−1 across northern ecosystems, primarily in temperate biomes (−0.26 ± 0.17 PgC yr−1). The annual mean gross loss of 4% of live biomass carbon in this region during 2016-2022 reflects high interannual variability, with significant losses associated with droughts and a further drop of −0.60 ± 0.47 PgC in the very dry year of 2022. Our findings highlight the vulnerability of live biomass carbon stocks to emerging climate-induced disturbances in northern ecosystems, challenging the sustainability of the current large terrestrial carbon sink in this key region for the global carbon balance.

Introduction
Northern ecosystems (north of 30°N) hold about 41% of the world’s forest area1, and are critical for mitigating global warming2,3, contributing an estimated 1.4-2.0 PgC annually to the global terrestrial carbon sink4. In recent decades, rising CO2 levels and global warming have been reported to increase mid- and high-latitude vegetation productivity by lengthening the growing season5, which has resulted in a ‘greening’ trend as observed by optical satellite systems6. Concordantly, multiple lines of evidence from forest inventories7,8 and satellite-based studies9–11 suggest that biomass carbon has increased in northern ecosystems in recent decades.

Climate hazards, however, may have affected the biomass carbon budget of northern ecosystems as a result of increased wildfire frequency12, extensive drought-induced tree mortality13,14 and unprecedented scale of insect outbreaks15. Previous studies have indicated that threats to forest carbon sinks have increased in North America16–18 and Europe19–21 over recent decades due to tree mortality caused by very high temperatures and droughts. Wildfires, a main driver of Siberian and Canadian boreal forest disturbances, have been reported to cause large carbon emissions during recent extreme fire years12,22 and have even led to parts of Siberia turning into a carbon source13,23. Compound droughts-heat waves (typically characterized by low soil moisture and high air temperatures) and wildfires have affected the northern Hemisphere in recent years, such as in 201519,24, 201820, 202112,25 and 202221,26.

Improving our knowledge of the interannual variability and trends in biomass carbon is crucial for a timely understanding of how climate change and human management affect the terrestrial carbon balance. Even in regions and countries with dense forest inventory sampling, such as the United States and Europe, the spatial sampling design remains insufficient to assess regional scale changes27. Further, the practice of conducting inventory revisits in multi-year cycles omits unmanaged forests and relies on simple methods such as static regrowth curves in countries using the gain-loss inventory method, which may ultimately not capture recent climate impacts28. Large uncertainties have also been reported in the quantification of carbon fluxes and sinks, for which biomass carbon change represents a key component, whether estimated from inventories8, dynamic global vegetation models (DGVMs)29, or atmospheric inversion models30. Alternatively, satellite data from light detection and ranging (LIDAR), optical and synthetic aperture radar (SAR) have been combined to map biomass carbon31,32 or track its annual changes over large regions11,27. However, due to the lengthy process of acquiring, processing, and formatting new data, the use of methods based on multi-source remote sensing technologies poses challenges in achieving consistent and continuous assessment, particularly under the recent climate change impacts in northern ecosystems. The net decrease of −0.20 ± 0.11 PgC yr−1 is the result of gross changes in biomass carbon as a consequence of northern ecosystems in the 2016-2022 period, which is reflected by the imbalance between gross gains and gross losses along the latitudinal band from 30° to 52°N (Fig. 2b, c). The majority of net losses occurred in temperate ecosystems (−0.26 ± 0.17 PgC yr−1) including North America (−0.15 ± 0.10 PgC yr−1) and Eurasia (−0.11 ± 0.07 PgC yr−1), while boreal biomes in these two regions experienced a small net gain (+0.07 ± 0.08 PgC yr−1) and loss (−0.03 ± 0.02 PgC yr−1) of carbon, respectively (Fig. 2b). Notably, a gross live biomass carbon loss of 1.50 ± 0.77 PgC yr−1 occurred in temperate biomes, which was partially offset by a gross gain of +1.24 ± 0.49 PgC yr−1. This gross annual loss represents about 4% of the total live biomass carbon in northern temperate regions (35.05 ± 2.76 PgC in 2016) (Supplementary Fig. 5), with hotspots in the south-central United States, southwestern Europe, and western Russia (Supplementary Fig. 4).

Annual gross losses in boreal biomes amounted to about 2% (1.29 ± 0.27 PgC yr−1) of the total live biomass carbon stock (54.26 ± 1.68 PgC in 2016), with higher losses in the Siberian region, where the spatial loss patterns are quite consistent with burned areas (Fig. 2b, d). Fan et al.13 reported that forest disturbances, attributed to fires and droughts in recent years, have reduced the biomass accumulation in Siberia, bringing the biomass carbon balance of the region close to neutrality. Gross gains in boreal regions were +1.32 ± 0.46 PgC yr−1, mainly located in northwestern Canada, Alaska, and northern Finland (Supplementary Fig. 4a). The fact that the reported values of gross fluxes were much larger than those of net fluxes suggests that changes in live biomass carbon stocks were highly dynamic during the study period, especially for the temperate biomes11.

Here, we estimated the wall-to-wall annual changes in total (the sum of above- and below-ground) live biomass carbon in northern ecosystems from 2010 to 2022 using L-band (1.4 GHz) vegetation optical depth (L-VOD) data from the Soil Moisture and Ocean Salinity (SMOS) satellite. L-VOD relies only on passive microwave observations and therefore it can be rapidly updated to probe the biomass components of living woody vegetation33,34. Moreover, the high penetration of low-frequency radiation through the canopy makes L-VOD less sensitive to changes in vegetation structure and more sensitive to biomass carbon stocks, which has been validated in previous studies35–37. The above-ground biomass carbon (AGC) was first retrieved from L-VOD (Methods). Total live biomass carbon was then derived from AGC using the root-to-shoot ratio of below- and above-ground biomass from the IPCC Guidelines for National Greenhouse Gas Inventories38. From the L-VOD-derived biomass carbon estimates, we analyse the changes in live biomass carbon at the regional scale and identify factors that drive the spatial variability in biomass carbon density changes across different biomes. The L-VOD-derived biomass carbon density was also compared with the DGVM simulations from TRENDY v12 project39, particularly focusing on recent trends.

Results
Changes in total live biomass carbon
For the years 2010–2022, the L-VOD-derived live biomass carbon stocks showed an average net gain in the northern ecosystems at a rate of +0.25 ± 0.26 PgC yr−1 (the range represents the uncertainties, derived from 12 calibrations of L-VOD against three benchmark maps; Methods). Among ecosystems, boreal regions had a net increase of +0.15 ± 0.16 PgC yr−1 in live biomass carbon, while northern temperate regions remained nearly neutral with a net change of +0.03 ± 0.04 PgC yr−1 (Fig. 1d, e). Our estimates, updated to 2022, are in line with previous global analyses of L-VOD until 201911, when comparing changes in biomass carbon over the same time period and under the same biome delineations (Supplementary Fig. 1). In addition, a comparison with the national inventories from the Food and Agriculture Organization (FAO) Global Forest Resources Assessment indicates that the L-VOD estimated biomass carbon stocks closely match those reported by countries (Supplementary Fig. 2). Moreover, the biomass carbon changes in managed forests reported by L-VOD show a strong spatial correlation with estimates from the national forest inventories across countries (R2 = 0.58, P < 0.01).

The temporal variations in annual biomass carbon shows a shift from a positive to a negative trend after 2016 (Fig. 1c). The turning point occurring around 2016 (P < 0.05) detected by a piecewise linear regression (Methods) is confirmed by the analysis of two independent datasets derived from multiple satellites available up to 2019: one from a high-frequency C-band VOD product40 and the other from global annual biomass maps11 with native spatial resolutions of 25 km and 10 km, respectively. Both sources indicated that live biomass carbon in northern ecosystems peaked in 2016, representing the highest levels observed since 2000 (Supplementary Fig. 3). Spatially, the transition to biomass loss revealed by L-VOD predominantly occurred in the south-central United States, Kazakhstan, and western Russia, while both loss and gain were observed in Europe (Fig. 1a, b). L-VOD-derived results indicated that the largest decline in biomass carbon during the 2010-2022 period occurred in 2022, with a net change of −0.67 ± 0.41 PgC across northern ecosystems, of which −0.60 ± 0.47 PgC occurred in temperate regions and −0.14 ± 0.12 PgC in boreal regions (Fig. 1c, d, e).

Between 2016 and 2022, net L-VOD-based biomass carbon estimates showed a significant decrease of −0.20 ± 0.11 PgC yr−1 across northern ecosystems, which is the result of gross changes in biomass carbon as a consequence of northern ecosystems in the 2016-2022 period, which is reflected by the imbalance between gross gains and gross losses along the latitudinal band from 30° to 52°N (Fig. 2b, c). The majority of net losses occurred in temperate ecosystems (−0.26 ± 0.17 PgC yr−1) including North America (−0.15 ± 0.10 PgC yr−1) and Eurasia (−0.11 ± 0.07 PgC yr−1), while boreal biomes in these two regions experienced a small net gain (+0.07 ± 0.08 PgC yr−1) and loss (−0.03 ± 0.02 PgC yr−1) of carbon, respectively (Fig. 2b). Notably, a gross live biomass carbon loss of 1.50 ± 0.77 PgC yr−1 occurred in temperate biomes, which was partially offset by a gross gain of +1.24 ± 0.49 PgC yr−1. This gross annual loss represents about 4% of the total live biomass carbon in northern temperate regions (35.05 ± 2.76 PgC in 2016) (Supplementary Fig. 5), with hotspots in the south-central United States, southwestern Europe, and western Russia (Supplementary Fig. 4).

Annual gross losses in boreal biomes amounted to about 2% (1.29 ± 0.27 PgC yr−1) of the total live biomass carbon stock (54.26 ± 1.68 PgC in 2016), with higher losses in the Siberian region, where the spatial loss patterns are quite consistent with burned areas (Fig. 2b, d). Fan et al.13 reported that forest disturbances, attributed to fires and droughts in recent years, have reduced the biomass accumulation in Siberia, bringing the biomass carbon balance of the region close to neutrality. Gross gains in boreal regions were +1.32 ± 0.46 PgC yr−1, mainly located in northwestern Canada, Alaska, and northern Finland (Supplementary Fig. 4a). The fact that the reported values of gross fluxes were much larger than those of net fluxes suggests that changes in live biomass carbon stocks were highly dynamic during the study period, especially for the temperate biomes11.

Drivers influencing the recent decline in live biomass carbon over northern ecosystems
The determinants of the spatio-temporal variations in live biomass carbon changes are multifaceted and can be related to drought, mortality and fire events in northern ecosystems. Spatially, regions showing a decrease in live biomass from 2016 to 2022 closely matched those of forest area loss reported by the Global Forest Watch (GFW)41 and the areas of fire disturbance monitored by active fire data42 (Fig. 2a, b, d). We separately calculated the total biomass losses in relation to forest area loss and burned area for each subregion of the northern ecosystems defined in the IPCC Sixth Assessment Report (Supplementary Fig. 6). The extent of biomass losses increases with burned area (R2 = 0.86, P < 0.01) and forest loss area (R2 = 0.49, P < 0.01) within each 25 km L-VOD grid cell. The high R2 for the relationship between biomass carbon losses and burned area indicates that wildfires are a significant disturbance factor reducing live biomass accumulation13,43. To further explore the factors behind the spatial variability of net changes in biomass carbon across different ecosystems during 2016-2022, we employed an interpretable machine learning (RF-SHAP) framework to quantify the relative importance of more than 15 factors separately for temperate and boreal biomes (Methods). The factors we considered include climate, vegetation, edaphic, and anthropogenic disturbances (Supplementary Table 2). Altogether, these factors explained 51% and 61% of the spatial variation in live biomass carbon changes observed in temperate and boreal biomes, respectively (Fig. 3). Moreover, the RF models presented good predictive skills in reproducing the observed spatial patterns of biomass carbon changes (Supplementary Fig. 7).

The most important climate-related variable affecting live biomass carbon change (Δbiomass) in both temperate and boreal biomes was the change in surface solar radiation (ΔRadiation), accounting for 13% and 19% of the explained variability, respectively (Fig. 3a, f). However, partial dependence plots revealed that the effect of ΔRadiation on Δbiomass differed between the two biomes, with Δbiomass decreasing as ΔRadiation increased in temperate areas, while a positive relationship was found in boreal areas (Fig. 3c, g). The negative relationship between Δbiomass and Δradiation in temperate regions may be associated with reduced soil moisture (SM) and increased vapor pressure deficit (VPD) as radiation increases (Supplementary Fig. 8a, b). The Spearman rank correlation between L-VOD biomass estimates and SM at low-mid latitudes highlight the negative sensitivity of live biomass carbon stock changes to soil moisture anomalies in this area (Supplementary Fig. 9f). ΔSM contributed 10% of the explained variability, ranking among the higher climate factors, with soil moisture deficits (i.e., negative ΔSM) linked to negative biomass carbon changes (Fig. 3a, e). The most significant factor positively influencing Δbiomass in temperate biome was tree cover with an average contribution of 18%, followed by soil nitrogen content at 11% (Fig. 3a, b, d). For the boreal biome, ΔRadiation is ranked as the most important factor affecting Δbiomass, which may reflect the limitation of solar energy on biomass accumulation, consistent with its role in constraining vegetation photosynthesis (gross primary production) for this region22. The U-shaped relationship, with the lowest point around ΔRadiation ≈ 0 J/m², is likely influenced by interactions with other factors, particularly disturbances related to forest area (Supplementary Fig. 8c, d). The second-most influential climatic factor was ΔVPD, contributing 8%, with increasing VPD associated with decreasing biomass (Fig. 3f, j), reflecting the fact that rising VPD is associated with tree mortality, reduced growth44, and fire activity27. Notably, while forest loss due to stand-replacing wildfires ranked high (9%), non-fire forest disturbances such as logging and drought-induced mortality had a comparatively larger impact on Δbiomass in the boreal areas, accounting for 12% (Fig. 3f, h).

For the years 2016-2022, the northern lands experienced extensive drought events, as indicated by negative anomalies in soil moisture and SPEI values alongside elevated VPD, especially in recent years45, such as the 2022 global widespread heatwaves (Fig. 4a, e, f). Over this period, the net decrease in L-VOD live biomass carbon in the low-mid latitudes (30°N~55°N) coincided with the negative trend in soil moisture (Fig. 4d). Moreover, these conditions of extreme drought were associated with large wildfires, which were particularly pronounced at high latitudes (>60°N) (Fig. 4c). Correspondingly, the area of forest loss per year in the boreal biome reported by GFW showed a clear positive trend at a rate of 1.01 × 106 ha/year (Supplementary Fig. 10). Notably, the live biomass carbon changes revealed by L-VOD showed a decreasing trend around 60°N–65°N during 2016-2022, which coincided with the highest average annual burned area recorded in this region over the same period (Fig. 4h). However, the TRENDY DGVMs with rapid update capabilities did not mirror this trend, partly due to limitations in their prognostic fire modules in realistically reflecting wildfire occurrence and intensity as well as the associated loss of live biomass. Yet, the three DGVMs models, JULES, OCN, and ORCHIDEE, which used prescribed burned areas from the ESA Climate Change Initiative (CCI-Fires) (Methods), revealed better consistency with L-VOD estimates compared to their prognostic fire versions (Supplementary Text 4).

During the summer drought in Siberia in 2021, when record wildfire emissions were reported12, the area of forest loss due to wildfires was 7.58 × 106 ha as monitored by stand-replacing wildfires data46. This event is consistent with a decline in live biomass carbon of −0.12 ± 0.11 PgC that year, which suggests that disturbances in recent years have further reduced live biomass carbon stocks in Siberian forests, consistent with Fan et al.13. The severe drought of summer 2022 had swept widely across northern temperate regions (Supplementary Fig. 11), such as Europe and the contiguous United States45, coinciding with the pronounced decrease in biomass carbon of −0.25 ± 0.25 PgC and −0.18 ± 0.11 PgC observed in these two regions during that year. Latitudinally, the decrease was particularly strong at low to mid-latitudes (30°-50°N), where the biomass carbon anomalies showed a clear congruence with the ERA547 soil moisture anomalies and SPEI (Fig. 4a, b, e). The L-VOD-derived live biomass carbon decrease in this latitudinal band was generally reproduced by the TRENDY DGVMs (Fig. 4g) and the high frequency X-VOD (Supplementary Fig. 9a), suggesting that the compound effects of droughts and heatwaves in 2022 severely suppressed the vegetation carbon uptake in northern temperate biomes21,26.

Live biomass carbon losses from forest loss and degradation during 2016-2022
Live biomass carbon losses in northern ecosystems have been attributed to increased forest fires13,43, insect outbreaks48, logging27 and tree mortality due to climatic disturbances such as windstorms and drought18. Here, we separated the contributions of forest area loss and ‘degradation’ to live biomass carbon losses using the method of Harris et al.49 and Qin et al.36. Since forest loss and degradation cannot be explicitly separated as two distinct processes within a 25 km grid cell, we simplified our definition by categorizing all mechanisms that do not result in forest area loss as degradation13. To do so, we subtracted live biomass carbon losses due to forest loss events, identified by Landsat data41 at a spatial resolution of 30 m, from the total biomass carbon losses. This means that smaller scale (<30 m) forest area loss is implicitly considered as degradation (Methods).

For the period 2016-2022, the gross live biomass carbon losses were −2.96 ± 2.36 PgC yr−1, of which ~38% (−1.12 ± 1.12 PgC yr−1) was estimated to result from forest area loss and ~62% (−1.84 ± 1.26 PgC yr−1) from degradation. Here, the degradation occurred mainly in temperate regions, accounting for 83% of the total biomass carbon loss in this biome, highlighting this as an important process driving biomass carbon loss in temperate biomes. In contrast, the forest area loss (−0.85 ± 0.97 PgC yr−1) was the major factor (~66%) in reducing the live biomass carbon stocks in boreal biomes (Supplementary Fig. 12). The impacts of stand-replacing wildfires and other stand-replacing processes were further separated in each 25 km grid cell in forest loss areas (Methods). The latter, including processes such as timber harvest and insect outbreaks, resulted in ~26% (−0.76 ± 0.68 PgC yr−1) of the gross live biomass losses across the northern ecosystems, with ~70% (−0.53 ± 0.48 PgC yr−1) and ~29% (−0.22 ± 0.20 PgC yr−1) occurring in boreal and temperate forests, respectively. The average annual carbon removed through timber harvest reported by FAO is less than half of the losses from non-fire stand-replacing disturbances (Supplementary Fig. 13), suggesting that harvesting may not be the main factor of these losses but rather other types of disturbances18,27, warranting further investigation. Stand-replacing fires (4.51 × 107 ha) corresponding to ~33% of the forest loss area caused a gross live biomass carbon loss of −0.37 ± 0.34 PgC yr−1 in northern ecosystems, of which ~87% (−0.32 ± 0.23 PgC yr−1) occurred in boreal forests (Supplementary Fig. 12). While this result confirms that wildfires are a key source of forest disturbance in boreal forests relative to other ecosystems50, other types of forest area disturbances within this biome were ~1.7 times larger than those caused by stand-replacing fires.

In addition to stand-replacing events associated with forest area loss, we further investigated the impacts of all wildfires on live biomass carbon stock changes in northern ecosystems. We considered ‘unburned’ areas defined as areas that include pixels with an annual fire area fraction < 1% within a 25 km grid between 2005 and 2022, while all other pixels were considered as partly burned areas, as in Fan et al.13. During 2016–2022, we found a net biomass carbon change of −0.17 ± 0.12 PgC yr−1 in partly burned areas covering ~55% of the study area (Supplementary Fig. 14), resulting in a reduction in live biomass carbon stocks to levels close to those observed a decade ago (Fig. 1c). In contrast, live biomass carbon stocks over unburned areas showed a smaller decline at a rate of −0.03 ± 0.05 PgC yr−1. Over temperate regions, partly burned and unburned areas showed similar net losses of −0.15 ± 0.12 PgC yr−1 and −0.10 ± 0.10 PgC yr−1, respectively. This similarity suggests that biomass carbon loss in this region is not directly driven by fire but by other factors (e.g., droughts), which is in line with the RF-based results (Fig. 3a). In contrast, there is a distinct pattern in the time variations in biomass carbon between fire-disturbed and non-disturbed areas in boreal regions (Fig. 1e). In the latter biome, the gross losses (−0.81 ± 0.76 PgC yr−1) and gains (+0.79 ± 0.86 PgC yr−1) in the partly burned areas were ~40 times larger than the net change (−0.02 ± 0.03 PgC yr−1), suggesting that wildfires were an important driver of biomass dynamics in this region. L-VOD-derived gross live biomass carbon losses in partly burned areas of the boreal regions were 2.7 times higher than those estimated solely from stand-replacing fires (−0.30 ± 0.32 PgC yr−1), implying that Landsat forest loss data may not capture all small-scale fire-burned forest loss patches and/or that disturbances such as drought and low-intensity understory wildfires may not always result in stand replacement but cause biomass carbon losses, consistent with the studies of Fan et al.13 and Rogers et al.51.

Methods
Northern Ecosystems
Our study area is located above 30°N, excluding southern China and southeastern Europe, where L-VOD is susceptible to radio frequency interference (RFI) from human activities33. There are two main ecosystems included in our study area, namely the northern temperate biomes and the boreal biomes, and their boundaries are defined according to Dinerstein et al.62. Due to the potential underestimation of L-VOD in areas where the SMOS footprint contains significant open water bodies10, we have excluded 25 km pixels predominantly characterized by ‘wetland’ land cover from our analysis. Our analysis ultimately focused on pixels (Fig. 1a) that consistently present values for each year from 2010 to 2022, as discontinuous observations could affect the calculation of trends.

SMOS-IC L-VOD
We investigated changes in live biomass carbon over northern ecosystems using the L-VOD product derived from the Soil Moisture and Ocean Salinity (SMOS) satellite, the first L-band passive microwave satellite mission since 2010. L-VOD, which probes the biomass components of live woody vegetation33,34. Moreover, the high penetration of low-frequency radiation through the canopy makes L-VOD less sensitive to changes in vegetation structure and more sensitive to biomass carbon stocks, which has been validated in previous studies35–37. The above-ground biomass carbon (AGC) was first retrieved from L-VOD (Methods). Total live biomass carbon was then derived from AGC using the root-to-shoot ratio of below- and above-ground biomass from the IPCC Guidelines for National Greenhouse Gas Inventories38. From the L-VOD-derived biomass carbon estimates, we analyse the changes in live biomass carbon at the regional scale and identify factors that drive the spatial variability in biomass carbon density changes across different biomes. The L-VOD-derived biomass carbon density was also compared with the DGVM simulations from TRENDY v12 project39, particularly focusing on recent trends.

As in refs. 10, 13, 64, 65, strict filtering of the SMOS data is necessary due to the effects of RFI, which disrupt the natural microwave signals from Earth’s surface as measured by passive microwave radiometers. Here, a triple-filtering method is used to filter the daily SMOS L-VOD data for both ascending and descending orbits by taking into account inter-orbit differences and outliers (Supplementary Text 1). The efficiency of this method has been demonstrated in a recent study of global biomass carbon estimation based on L-VOD10. As L-VOD is directly proportional to the total vegetation water content (VWC), daily values of L-VOD reflect not only changes in water mass due to changes in aboveground biomass (AGB) of live vegetation, but also to changes in water mass per unit of biomass (relative water content, RWC)66. Previous studies have assumed that RWC is relatively stable on an annual scale, or during the wet season, and does not show long-term (typically within a decade) increasing or decreasing trends, so that long-term changes in annual L-VOD can be directly related to changes in AGB10,34–36,58,66. A detailed analysis of these assumptions and a review of validation studies on the use of L-VOD for biomass monitoring is provided in Wigneron et al.34. In this study, to further eliminate a possible contribution of the RWC variations to changes in annual L-VOD, we developed multiple regression models based on the multiplicative relationship considering L-VOD scales with both dry biomass and RWC67 (Supplementary Text 2). Correction for water information in L-VOD showed that RWC had a minimal effect on annual L-VOD changes, changing the trend by less than 5% between corrected and uncorrected data (Supplementary Figs. 15-17). This is consistent with the conclusions of Yang et al.57 in a study covering the tropics.

Benchmark maps of vegetation carbon density
SMOS L-VOD was calibrated according to aboveground carbon density (AGC; Mg C ha−1) using previously published live biomass maps as a benchmark10,13. However, considering that none of the available live biomass maps can be considered completely reliable (because the true AGC map is unknown) and that they all contain uncertainties and biases, we focused on utilizing the spatial gradients of AGC derived from several existing maps. Any systematic errors in the AGC reference maps are also a major external source of uncertainty in the L-VOD-derived live biomass maps. While individual reference maps may have biases, using an averaged calibration from multiple maps across various regions as the final biomass carbon density mitigates potential biases associated with any single map, as outlined in the L-VOD biomass data production and uncertainty analysis. In a second step, the temporal dynamics of AGC was derived from L-VOD. Here, we used three of the most recent global static AGC benchmarks, including those provided by Saatchi et al.31 updated in 2015, Santoro et al.68 in 2010, as well as ESA CCI biomass map52 (https://climate.esa.int/en/projects/biomass/) in 2017, hereafter referred to as the “Saatchi”, “GlobBiomass” and “ESA CCI”, respectively. All these maps were aggregated by averaging to match the spatial resolution of SMOS (i.e., 25 km grid cell).

Total live biomass carbon derived from L-VOD
The method used here to obtain the calibration coefficients between the annual L-VOD and the reference AGC is the same as that used in ref. 13, where it is documented in detail. The AGC was calculated from the annual L-VOD based on the following empirical spatial calibration:

AGC = a × L-VOD^b (1)

where a and b are best fit coefficients. Here, we used the annual L-VOD in 2011 to calibrate Eq. 1, and as shown by ref. 13, the year of L-VOD used for the calibration has a very little impact on the final calibration result.

In contrast to previous studies, instead of creating a single spatial calibration function for the study area, where possible outliers could affect the fitting parameters, we roughly (not according to administrative regions) divided the northern ecosystems (>30°N) into three parts: Russia (40°E-180°E), Europe (20°W-40°E), and North America (20°W-180°W) and then computed the optimal coefficients of the calibration equations for the three parts of the study area separately. Twelve calibrations of Eq. 1 were thereby obtained (Supplementary Table 1). They were generally consistent across subregions and throughout the entire northern ecosystems (Supplementary Fig. 18), suggesting that the statistical relationship between L-VOD and vegetation carbon stocks is robust across regions. To estimate the total biomass carbon (BC) from aboveground biomass, we estimated belowground biomass using the root-to-shoot ratio38 suggested by the IPCC Guidelines for National Greenhouse Gas Inventories (Supplementary Table 4). We produced twelve maps of BC stocks using all twelve calibrated L-VOD AGCs. As in refs. 10, 35, the median of these twelve maps was finally used to calculate annual BC maps for the period 2010–2022.

Uncertainties and evaluation related to the L-VOD-derived live biomass carbon product
There are two main sources of uncertainty associated with the L-VOD-derived live BC estimates, including: i) external uncertainties, which are caused by the systematic errors in the reference biomass maps; ii) internal uncertainties corresponding to calibration and sampling strategy errors. Previous analyses by Fan et al.13 found a strong dominance of external errors over internal uncertainties (by an order of magnitude), and the total relative uncertainties relative for the carbon stock and carbon stock changes are on the order of 20–30%. Our bootstrap cross-validation results showed that all twelve BC maps exhibited high spatial correlation values and lower RMSE values compared to the corresponding reference dataset, suggesting that no obvious regional bias could be found between the reference biomass maps and the L-VOD-derived live BC over the different continents (Supplementary Table 5). In addition, the 95% bootstrap confidence intervals for the BC estimates obtained using each set of calibrations are small, suggesting that the internal uncertainty due to sampling/calibration errors is small. In consideration of the above, only external errors, as indicated by the range (or spread) of the twelve maps, are reported in the main text.

Wigneron et al.34 provided a review of evaluation studies using L-VOD to calculate BC. Here, to further evaluate the L-VOD-derived live BC data, we used biomass carbon data from the National Forest Inventory (NFI) reported by 43 countries to the Food and Agriculture Organization of the United Nations (FAO) and by 34 Annex I countries to the United Nations Framework Convention on Climate Change (UNFCCC), respectively, which are relatively accurate and independent sources. Spatially, the L-VOD-based live BC stocks were not only found to be similar to the carbon stock values reported by the FAO, but also to circumvent the underestimation problem of using a single reference dataset in some countries (Supplementary Text 3). Temporally, the total net carbon change in managed forests estimated by L-VOD for Annex I countries (Supplementary Fig. 19) between 2010 and 2021 correlated well with that provided by the UNFCCC (R2 = 0.58, P < 0.01) (Supplementary Fig. 2). This long-term comparison added confidence to the L-VOD estimates of biomass carbon changes. It should be noted that the UNFCCC NFI could not capture updated interannual biomass variations, as the NFI was limited by its revisit time rates (in most cases, in each country, several years were required before all plots of the NFI network can be revisited) (Supplementary Fig. 20).

Outputs of dynamic global vegetation models
We utilized the ‘cVeg’ outputs, which represent total live biomass as the sum of above- and below-ground biomass, from 17 DGVMs provided by the TRENDY v12 project39: CABLE-POP, CARDAMOM, CLASS, CLM5.0, DLEM, EDv3, IBIS, ISAM, JSBACH, JULES, LPJ-GUESS, JPLml, LPJwsl, LPX-Bern, OCN, ORCHIDEE and SDGVM. The TRENDY DGVMs performed three factorial simulations. S1 was driven only by CO2, using a pre-industrial climate and land-use mask. S2 included both CO2 and climate while retaining a pre-industrial land-use mask. S3 incorporated CO2, climate, and land use, all varying over time. We selected the S3 simulations, which were spun-up by recycling climate data from 1960 to 1969 until carbon pools were in equilibrium, and then applied varying CO2, climate and land-use conditions over the period 1960 to 2022. We used the mean of these DGVMs in the main text, marking uncertainty by the range between maximum and minimum values. It should be noted that these DGVMs, which have rapid update capabilities, handle emissions using a prognostic fire module. In addition, we analyzed three DGVMs (JULES, OCN, and ORCHIDEE), which were the only models available through 2020 at the time of analysis and used prescribed burned areas from the ESA Climate Change Initiative (CCI-Fires) for calculating emissions69. This allowed us to focus on their consistency with L-VOD estimates before and after incorporating fire diagnostics (Supplementary Text 4).

The RF-SHAP interpretation framework
We applied the Random Forest models-Shapley additive explanations (RF-SHAP) framework to quantify the relative importance of the different climatic, environmental and human-induced variables on the live BC density changes. Comprehensive details of these variables are cataloged in Supplementary Table 2. Here, we aim to understand the underlying processes driving net changes in live biomass carbon density (Δbiomass) in the northern region during 2016–2022, considering both recovery and disturbance factors. While this study focuses on net changes, a separate analysis of biomass carbon losses57 and gains58 may offer additional insights and merits more detailed investigation in future work. The RF model is a widely used supervised learning algorithm, known for reducing overfitting through ensemble learning and random feature selection57. It builds multiple decision trees on random subsets of the training data, and aggregates their predictions for the final result. Our modeling approach employed a hold-out validation strategy, using 75% of the dataset for training the RF model and reserving the remaining 25% for validation. To determine the optimal number of trees for the RF model, we used a trial-and-error approach, incrementally varying the ‘number of trees’ parameter from 25 to 1500 in steps of 100. We ultimately set the number of trees to 400, as the model performance on the test set, measured by R², improved with increasing number of trees and plateaued beyond 400 (Supplementary Fig. 21). We built separate RF models for temperate and boreal biomes, and compared the spatial distribution of predictions from these models with observed values (Supplementary Fig. 7). Finally, we used the state-of-the-art TreeExplainer-based SHAP framework to explain each RF model. SHAP, based on cooperative game theory, quantifies marginal contributions through SHAP values and captures complex interactions, offering a more nuanced view of variable importance than traditional metrics70.

The SHAP value represents each variable’s contribution to the model’s prediction, with positive values associated with higher predicted Δbiomass, and negative values with lower predictions. The importance of each feature was determined by the mean absolute SHAP values, which were scaled to percentage contributions for easier interpretation. We applied SHAP to each RF model 30 times and calculated the standard deviation to quantify uncertainty. To elucidate the nonlinear relationships between Δbiomass and its four key explanatory variables across the two ecosystems, SHAP-based partial dependence plots (PDPs) were generated, showing how SHAP values change with each variable while accounting for both its direct effects and interactions with other predictors70. By combining SHAP rankings with SHAP-based PDPs, we not only identify variable importance but also gain insights into how the top variables shape Δbiomass patterns.

To ensure robustness and reliability, the scatter plot highlight points explained by SHAP at least 10 times while rendering others more transparent. The variable range was then divided into 50 bins based on its maximum and minimum values, and the median was computed only for bins containing more than 20 points to clearly illustrate the variable’s contribution pattern to biomass change (Fig. 3).

In addition to the RF model, we also employed a conventional least-squares-based multiple linear regression (MLR) to validate the robustness of the machine learning model. We calculated the variance inflation factor (VIF) for each variable and found that most predictors had a VIF below 3, indicating minimal multicollinearity. Although removing variables with high VIFs (e.g., soil organic matter (SOC) with a VIF ~4 in the northern temperate model) reduced the VIFs for the remaining variables, the Akaike Information Criterion (AIC) increased, indicating a reduction in the model’s explanatory power for the target variable. Consequently, the optimal MLR model retained all of the predictors listed in Supplementary Table 2. However, the MLR models could only explain 13% and 15% of biomass variation in northern temperate and boreal biomes, respectively, which is significantly lower than the 51% and 61% explained by our RF models (Supplementary Fig. 22). This suggests that the linear model is insufficient to explain variations in Δbiomass compared to the RF model.

Vegetation and climate variables
Forest loss data. The global forest change product generated from Landsat imagery by Hansen et al.41 was used to calculate the rate of forest loss. Forest loss provided in this product, called ‘loss year’, is defined as a stand-replacing disturbance, or a change from forest to non-forest. Each 30 m pixel of the ‘loss year’ data is labeled with a number between 0 and 22, where ‘0’ represents no loss and the others correspond to the loss of forest cover (defined as canopy closure for all vegetation above 5 m) in that year during 2001–2022. Here, the percentage forest loss rate for 2016–2022 at the SMOS resolution is calculated as the proportion of the summed area of forest loss for this period within each 25 km grid cell, as done by Fan et al.13.

Stand-replacing wildfires product
This product, updated from van Wees et al.46, estimates the areas where burned areas or active fires overlap with forest losses at a spatial resolution of 500 m. It defines fire-induced forest losses as encompassing various sequences of events involving wildfire and forest loss. This includes cases where fire (e.g., wildfires) and forest loss occur simultaneously, post-fire forest loss (e.g., tree mortality due to fire damage), and fire following forest loss (e.g., burning of felled trees after logging, which mostly occurs in the same year as logging). We ultimately computed the annual stand-replacing fire fraction at the SMOS resolution using the same approach as above for processing forest loss.

Non-fire induced forest loss
Non-fire related forest losses are defined as forest loss outside of the various sequences of events involving wildfires and forest loss described above. Similarly, the annual fraction of non-fire induced forest loss for the period 2010–2022 was calculated as the proportion of the summed area of non-fire related forest loss within each 25 km grid cell.

Annual burned fraction
The annual active fire maps were generated using active fire products from the latest version of the Terra MODIS Level 3 Thermal Anomalies/Fire 8-Day Global 1 km SIN Grid product (MOD14A2, V6.1). The MOD14A2 provides 8-day fire mask composites at 1 km resolution, which contains the maximum values for each pixel class during the composite period42. This widely used dataset provides global coverage, making it ideal for monitoring fire activity across diverse ecosystems, while its high temporal resolution enables timely detection over large areas12,13. The MOD14A2 active fire observations with a high confidence level were first identified using the associated quality information. We then calculated the annual burned fraction at the SMOS resolution for the study period as the proportion of the summed area of active fires assuming that each 1 km pixel was fully burned.

SPEI-12
The Standardized Precipitation-Evapotranspiration Index (SPEI)-12 refers to drought conditions measured over a 12-month timescale, obtained from the global SPEI database (https://spei.csic.es/). SPEI-12 is a multi-scalar index commonly used to quantify drought. Unlike some other drought indices, SPEI accounts for multiple climatological factors, including both precipitation and temperature, making it crucial for evaluating the impact of climate change on drought conditions71.

Additional data
Soil moisture and 2-m air temperature data resampled to the SMOS resolution were obtained from the European Centre for Medium-Range Weather Forecasts Reanalysis v.5 (ECMWF ERA5), which is an advanced and widely used climate reanalysis47. NDVI data from the MODIS Vegetation Indices 16-Day (MOD13C1) V6.1 product were used as a proxy for leaf biomass. Other high-frequency VOD products at the C-band for 2003–2018 (C-VOD) and X-band (X-VOD) for 2012–2022 are provided by the global long-term microwave VOD Climate Archive (VODCA)40 and the global land parameter data record (LPDR)63, respectively.

Net and gross gains/losses in BC
Mean annual net carbon changes were estimated from the difference in BC between the end year and the start year divided by the year interval. A positive value indicates a net gain (sink) in BC, while a negative value indicates a net loss (source) of carbon. Gross losses (gains) were estimated by accumulating negative (positive) differences between the consecutive years over each pixel, following the method of Brandt et al.55.

Identification of the turning point in live biomass carbon
We used the piecewise linear regression method to determine the turning point of total live BC across the northern ecosystems72. The Mann-Kendall nonparametric trend test was then used to calculate the temporal trends in total live BC after the turning point, while the Sen’s slope statistics were used to determine confidence intervals.

Contributions of forest loss and forest degradation to gross live biomass carbon losses
The gross losses of live BC in a SMOS grid cell can be attributed to loss of forest area, forest degradation and other mechanisms such as non-forest carbon loss. Here, we performed a simple separation of the contribution of forest area loss and degradation to live BC loss based on the method proposed by Harris et al.49 and Qin et al.36. It should be noted that forest loss and forest degradation are not independent processes and cannot be explicitly separated on a 25 km grid cell, and for simplicity we defined all mechanisms that do not result in loss of forest area as degradation. As in ref. 36, the contribution of forest loss to gross live BC loss (Gross BC loss_forestloss) was first isolated by multiplying the total area of forest loss (calculated from the 30 m Landsat forest loss map41) during 2017-2022 by the BC density in 2016 (Eq. 2). Degradation was then roughly estimated as subtracting this forest loss from the gross L-VOD-derived BC losses (Eq. 3).

Gross BC loss_forestloss = (Gross forest area loss) × BC density (2)

Gross BC loss_degradation = Gross BC loss - Gross BC loss_forestloss (3)

The contributions of stand-replacing wildfires (Gross BC loss_stand-replacingfire) and other stand-replacing processes (e.g., clearcutting and widespread insect outbreaks) were further separated from the gross live BC losses due to forest loss (Gross BC loss_forestloss). To do so, we multiplied the total area of forest loss due to stand-replacing fires (i.e., the sum of the area of forest loss overlapping with fire in each year during 2017-2022) by the live BC density in 2016.

Gross BC loss_stand-replacingfires = (Gross forest area loss_fire) × BC density (4)

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of “Large biomass carbon losses from droughts in the northern temperate ecosystems during 2016-2022 ” by Li et al. This manuscript analyzes time series of changes in above-ground biomass during the period 2010-2022 over ecosystems north of 30N, with the goal to derive trends in carbon storage. The dataset has annual resolution which allows temporal trends, and extreme years to be analyzed. This leads to the most important result reported here: an apparent reversal of a positive trend in biomass storage into a negative trend, which signifies loss of carbon from forests in the NH since 2016. Reduced carbon uptake in these regions was seen before in this dataset, at a time when it only reached to 2019, while 3-years of additional data now allow significant trends over large subsets of the pixels of the NH. Moreover, earlier speculations of this reversal being in progress, and driven by large-scale disturbances from fires and droughts, can be confirmed by the inclusion of post-2019 large climate extremes in Siberia, Europe, and North America. Since these events were separately studied and reported (some also from the perspective of ΔAGB) in previous publications, the novelty of this manuscript is in the analysis of the longer time series, and the focus on the decline in forest biomass since 2016. In the Methodology, the authors properly benchmark the new dataset against other biomass products, and make a credible case for its agreement with bottom-up derived changes in forest carbon storage at country level. This leaves little doubt that the observed trend is a real signal, and not an artefact of the data used, nor of its known sensitivity to water that partially absorbs the microwave signals from SMOS. Statistical analyses are mostly straightforward and properly tested for significance. Error estimates and propagation of errors are present in the numerical results in the text and tables, but should be added in some of the Figures. The extended analyses of the dataset focus on the possible drivers of the biomass decline, and this part does not fully satisfy me. Quite some words and figures are spent on a machine-learning effort on the spatial patterns of ΔAGB, which ends up with a large number of contributing proxy variables, with unknown and likely non-linear relations to the data trained on. The purpose of this effort seems to have been to explain the variations in the dataset with certain proxies, rather than to find a predictive model for untrained data (or events yet to come). If indeed the purpose is to understand the role of various drivers, then a regular multilinear regression, or a principal component analysis seems more logical to me. These allow you to optimize and quantify the variance explained by each proxy, instead of letting a ML-method combine them in a way that is afterwards untraceable, and interpretable only from their rank in the prediction. Any reply to this concern would at the minimum warrant a justification and explanation of the added value of the use of the ML-model in this manuscript. About the outcomes of the ML-approach: I see that highly informative variables end up to be partly related to forest structure (tree cover, soil nitrogen load, age), and partly to changes in hydrological conditions (Δsoil moisture, ΔVPD, Δradiation, ΔT). The latter are all likely to be strongly covarying as they are different representations of drought conditions, or known consequences (Δfire). The outcomes show some interesting relations (e.g. a sharp gradient in the relation with ΔVPD) which feel like they need to be further discussed or analyzed. It also seems to me that the machine-learning predominantly confirms -qualitatively- the conditions that were described in previously published analyses of these events, rather than new insights. But I am open to arguments from the authors. Similarly, the authors introduce in Figure 4 an analysis that includes three DGVMs. The basis for selecting these three remains unclear, and the chosen figures clarify very little about the actual processes captured by these models that would make them reproduce the observed ΔAGB. In fact, besides the reduced BC density coinciding with low SPEI or low soil moisture, the models reproduce very little of the observed signal. They of course know about hot and dry locations through their driver data (reanalyzed temperature and precipitation), but this has been shown in many studies. Here it matters that they explain changes in AGB. But the trend in ΔAGB per latitude is captured only partly by only one model (Fig 4f), while the aggregate changes even over large NH land areas is not reproduced by any (Table S5). It is therefore not surprising that they are not used for further insights on the processes driving the trends. Without this, I miss a justification to use these process-based models and the paper would be stronger without them. I would like to read a strong justification of their use, and the specific selection of models and co-authors that was made. My suggestion for further improvement of the manuscript would be to keep the opening two sections and Figures of the current work, but to focus the manuscript after that on the spatiotemporal analysis and contributing processes (such as in the strong final section on forest loss and degradation) instead of the proxies. For example, one could discuss specific geographic areas of the NH in relation to known disturbances, and go beyond drought and fire but also include harvest, regrowth, or shifting land-use. The mentioned areas in western Canada, Alaska and Scandinavia can be a start, but also western/central USA, Kazakhstan and its neighboring countries could be highlighted. Another angle to provide more context to the finding of a trend reversal is the fate of the carbon lost: how much of it has returned to the atmosphere and constitutes a true loss of sink, and how much has ended up in semi-permanent stores (such as buildings and furniture), possibly after salvage logging? Overall, I find the longer dataset and its dual-period trend analysis sufficiently interesting to warrant publication. But I find the current manuscript to lack innovation in the analysis of drivers, processes, and regional differences in them due to the use of the uninformative ML method and the uninformative trendy models. The re-discovery of various known recent anomalies moreover asks for a more quantitative comparison to the previous papers published on them. With these elements in place, the work would likely also be more appealing to the wider users of this journal, and establish itself as the go-to dataset for analysis of recent carbon cycle changes. Below, I subsequently add some textual suggestions to consider for a possible resubmission of the manuscript: - The methodology section ‘Benchmark maps of vegetation carbon density’ discusses the importance of calibration against existing biomass maps. To overcome the inevitable uncertainties and biases in existing maps, the authors use spatial gradients of the existing maps. Since this seems to be an important step in the methodology, I suggest adding a brief explanation on the use of spatial gradients and its effect on the conversion from L-VOD to AGC. - Figure 4d shows a combination of several datasets: L-VOD, SPEI and soil moisture. Although I appreciate the attempt to bring these datasets together in one plot, it took me some time to grasp it. Specifically, I think that overlying three distinct line plots, with an inverted y-axis for SPEI, could be confusing to some readers. Hence, I suggest to plot the SPEI in a small separate graph just above the lines representing L-VOD. - Lines 125-126, pixels could be written in plural. - Line 175: ‘relatively’ could be written as ‘relative’. - Line 223: could ‘temperature’ be a mistaken for ‘temperate’ here?", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "Title: Large biomass carbon losses from droughts in the northern temperate ecosystems during 2016-2022 This is a remote sensing-based vegetation biomass study which is critical for current carbon cycle science. The study provides quantifications of the impacts of drought and wildfire of biomass change. If possible, I would suggest the authors to provide some forest area change information. Abstract Nice and clear. But please add the spatial resolution 25km in the Abstract. Also, I noticed that the 4% biomass drop during 2016-2022 in the main text is a very important information. I would suggest adding that 4% in the Abstract if possible, and hopefully you can also point out where it happened (e.g., Russia, U.S., Amazon, etc.). Introduction Line 93: \"... exploring their causes in relation to recent drought and wildfire\" Seems later you also mentioned forest cutting in several places (Line 299, 656)? If cutting impact on biomass loss is smaller than fire, then it is OK to omit it, otherwise, maybe it is good to mention why cutting is omitted (e.g. lack of data, not a focus of this paper, cutting is too small, or cutting can not be captured by remote sensing while fire can be captured?). Line 113-114: \"Both sources indicated that live biomass carbon in northern ecosystems peaked in 2016,\" Suggest providing spatial resolution of those two data products. Line 125-126: \"The partly burned areas are pixel with an annual fire area fraction > 1% within a 25 km grid from 2005 to 2022\" Looks like this gives about 27.5/49.9=55% of the global vegetation area being partially burned. It is far from my thinking. Can you try to find a realistic total area of fire disturbance? The \"partly burned\" may be too loose. Effects of other changes within the partial burned areas are counted as fire effect, I guess. Suppose a fire can kill 50% of live biomass, then 1% area burn could only kill 0.5% total biomass. If bigger C loss was detected, say 1%, could it because of some other factors as well? Line 141-142:\"Notably, temperate ecosystems sustained a gross carbon loss of −1.50−1.77−1.26 PgC yr-1 which represented about 4% of the total biomass carbon\" Are you saying 4% biomass carbon loss in one single year? Please double check. I guess it should be 2016-2022 total biomass carbon loss? Figure 2: Are the gain and loss were calculated based on land cells? Line 176-179: \"These factors include climatic, vegetation, edaphic, and anthropogenic disturbances (Table S4). ...... (Fig. 3)\" This section talks about drivers. if possible, I would suggest adding one sentence to give the relative contribution (X%) of each of the four driving factors to biomass loss. Seems to me Fig.3 is saying that Fire is the biggest influencing factor, and human activity is the smallest influencing factor. But it is for something about 'spatial variability', not the carbon loss quantity. Line 254: \"For the period 2016-2020, biomass carbon losses were −2.96−3.38−2.65 PgC yr-1,\" This carbon loss amount is significant. in Line 143, you mentioned temperate region has a total biomass of 35.05 PgC in2016. I didn't find the number for boreal regionas. Try to add boreal total biomass carbon if you have. Also I believe it shouldn't be a per year value. (OK, coming back to this comment again, I realized maybe your carbon loss is not a net loss?) Line 254-268: The whole paragraph gives me an impression that fire is not the #1 factor of biomass carbon loss compared to cutting and insect. Please double check the PgC numbers and X% numbers to make sure that fire effect is bigger than forest area loss, cutting loss, and insect loss. Line 284-287: \"... indicating that Landsat forest loss data cannot capture small-scale fire-burned forest loss patches and/or that disturbances such as drought and low-intensity understory wildfires may not always result in stand replacement but cause biomass carbon losses, consistent with the studies of Fan et al. (2023)14 and Rogers et al. (2015)50.\" I agree that most wildfires are not stand-replacing. But I am curious that if 30m data cannot capture such fire-burned forest loss, how can 25km data capture the fire-burned forest loss? OK, I see you are using MODIS fire data (Ref # 42, 250m?), then I suggest add one sentence in main text saying that MODIS has a better way to capture fire events. (And it is multiple year data with global coverage, better than Landsat single year data) L295-297: \"The downward trend from 2016 to 2022 is reflected by a net loss of biomass carbon over the northern temperate biome (−0.26−0.33−0.17 PgC yr-1) and boreal Eurasia (−0.03−0.04−0.02 PgC yr-1).\" Here it is obvious that temperate biome has bigger biomass carbon loss than boreal, right? But in Line 254 you mentioned -2.96 PgC for total loss, in Line 142 you mentioned -1.5 PgC. That makes me thinking boreal might lost ~1.46 PgC, almost equal to temperate region. Maybe I read the data wrong, but please double check. Line 300-302: \"Notably, 2022, the year most affected 300 by drought, showed the sharpest decline in biomass carbon stocks over the 2010-2022 period, with a net change of −0.67−0.69−0.41 PgC.\" ~2% biomass loss in the year 2022 is still a huge number at global scale.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Thanks for the positive feedback from the Reviewers. The requests from the reviewers are reiterated below in italics, with our responses and revisions presented beneath each suggestion in blue and red text, respectively. We have fully addressed each request. [Reviewer #1 Comment 1] Review of “Large biomass carbon losses from droughts in the northern temperate ecosystems during 2016-2022” by Li et al. This manuscript analyzes time series of changes in above-ground biomass during the period 2010-2022 over ecosystems north of 30N, with the goal to derive trends in carbon storage. The dataset has annual resolution which allows temporal trends, and extreme years to be analyzed. This leads to the most important result reported here: an apparent reversal of a positive trend in biomass storage into a negative trend, which signifies loss of carbon from forests in the NH since 2016. Reduced carbon uptake in these regions was seen before in this dataset, at a time when it only reached to 2019, while 3-years of additional data now allow significant trends over large subsets of the pixels of the NH. Moreover, earlier speculations of this reversal being in progress, and driven by large-scale disturbances from fires and droughts, can be confirmed by the inclusion of post-2019 large climate extremes in Siberia, Europe, and North America. Since these events were separately studied and reported (some also from the perspective of ΔAGB) in previous publications, the novelty of this manuscript is in the analysis of the longer time series, and the focus on the decline in forest biomass since 2016. [Response]: Thank you for taking the time to review our manuscript. We greatly appreciate your valuable feedback and insights, which have helped us to refine the manuscript. All of your comments have been addressed in the revised manuscript, and please find our point-by-point responses below. [Reviewer #1 Comment 2] In the Methodology, the authors properly benchmark the new dataset against other biomass products, and make a credible case for its agreement with bottom-up derived changes in forest carbon storage at country level. This leaves little doubt that the observed trend is a real signal, and not an artefact of the data used, nor of its known sensitivity to water that partially absorbs the microwave signals from SMOS. Statistical analyses are mostly straightforward and properly tested for significance. Error estimates and propagation of errors are present in the numerical results in the text and tables, but should be added in some of the Figures. [Response]: Thank you for your kind approval and suggestions. We have incorporated the estimation uncertainties into the revised figures, including Fig. 2c, Fig. 3, Fig. 4h, Extended Data Fig.1 and Extended Data Fig.3. For your reference, we have provided the updated Fig. 2c, Extended Data Fig.1 and Extended Data Fig.3, while the other figures are included in the responses to your comments below. [Reviewer #1 Comment 3] The extended analyses of the dataset focus on the possible drivers of the biomass decline, and this part does not fully satisfy me. Quite some words and figures are spent on a machine-learning effort on the spatial patterns of ΔAGB, which ends up with a large number of contributing proxy variables, with unknown and likely non-linear relations to the data trained on. The purpose of this effort seems to have been to explain the variations in the dataset with certain proxies, rather than to find a predictive model for untrained data (or events yet to come). If indeed the purpose is to understand the role of various drivers, then a regular multilinear regression, or a principal component analysis seems more logical to me. These allow you to optimize and quantify the variance explained by each proxy, instead of letting a ML-method combine them in a way that is afterwards untraceable, and interpretable only from their rank in the prediction. Any reply to this concern would at the minimum warrant a justification and explanation of the added value of the use of the ML-model in this manuscript. [Response]: Thank you for your valuable feedback. Indeed, our purpose in using the ML model was to quantify the relative importance of different climatic, environmental, and anthropogenic variables on biomass change. Furthermore, we fully agree with you that multiple linear regression (MLR) is a conventional approach for quantifying variable importance. However, the influence of these factors on biomass carbon change is not a simple linear relationship, as revealed by our dependency plots (as shown in Fig.3 in Comment 4). This suggests that the linear assumptions may be problematic. Moreover, we conducted an MLR analysis using the ordinary least squares method. After accounting for the variance inflation factor (VIF) and Akaike Information Criterion (AIC), we identified such optimal regression models. These models could only explain 13% and 15% of biomass variation in northern temperate and boreal ecosystems, respectively (see newly added Fig. S17 below), which is significantly lower than the 51% and 61% explained by our RF model, further reinforcing the issues with linear models. Instead, the combination of the Random Forest (RF) model and SHapley Additive exPlanation (SHAP) could effectively elucidate the complex nonlinear relationships between variables and biomass change while providing a ranking of variable importance (Lundberg et al., 2020). In response to your comment, we have incorporated the MLR results into the methodology. Additionally, considering that readers may find SHAP values challenging to interpret and may perceive them as a qualitative method, we have normalized them to percentage contributions and refined the description of the RF-SHAP framework in the revised manuscript: “ The RF- SHAP interpretation framework We applied the Random Forest models-Shapley additive explanations (RF-SHAP) framework to quantify the relative importance of the different climatic, environmental and human-induced variables on the BC changes. Comprehensive details of these variables are cataloged in Table S4. The RF model is a widely used supervised learning algorithm, known for reducing overfitting through ensemble learning and random feature selection65. It builds multiple decision trees on random subsets of the training data, and aggregates their predictions for the final result. Our modeling approach employed a hold-out validation strategy, using 75% of the dataset for training the RF model and reserving the remaining 25% for validation. To determine the optimal number of trees for the RF model, we used a trial-and-error approach, incrementally varying the ‘number of trees’ parameter from 25 to 1,500 in steps of 100. We ultimately set the number of trees to 400, as the model performance on the test set, measured by R², improved with increasing number of trees and plateaued beyond 400 (Fig. S16). We built separate RF models for temperate and boreal biomes, and compared the spatial distribution of predictions from these models with observed values (Fig. S4). Finally, we used the state-of-the-art TreeExplainer-based SHAP framework to explain each RF model. SHAP, based on cooperative game theory, quantifies marginal contributions through SHAP values and captures complex interactions, offering a more nuanced view of variable importance than traditional metrics69. The importance of each feature was determined by the mean absolute SHAP values, which were scaled to percentage contributions for easier interpretation. We applied SHAP to each RF model 30 times and calculated the standard deviation to quantify uncertainty. To elucidate the nonlinear relationships between biomass density changes in the two ecosystems and their four key explanatory variables, partial dependencies of the target variables on these factors were calculated, averaging out the effects of other predictor variables70. A bootstrapping approach was also employed to capture uncertainty. Specifically, we systematically repeated the process 1,000 times, randomly selecting 70% of the data for each explanatory variable to generate partial dependence curves. The resulting means and individual estimates were then used for further analysis (Fig. 3). In addition to the RF model, we also employed a conventional least-squares-based multiple linear regression (MLR) to validate the robustness of the machine learning model. We calculated the variance inflation factor (VIF) for each variable and found that most predictors had a VIF below 3, indicating minimal multicollinearity. Although removing variables with high VIFs (e.g., soil organic matter (SOC) with a VIF ~ 4 in the northern temperate model) reduced the VIFs for the remaining variables, the Akaike Information Criterion (AIC) increased, indicating a reduction in the model’s explanatory power for the target variable. Consequently, the optimal MLR model retained all of the predictors listed in Table S4. However, the MLR models could only explain 13% and 15% of biomass variation in northern temperate and boreal biomes, respectively, which is significantly lower than the 51% and 61% explained by our RF models (Fig. S17). This suggests that the linear model is insufficient to explain variations in biomass change compared to the RF model.” [Reviewer #1 Comment 4] About the outcomes of the ML-approach: I see that highly informative variables end up to be partly related to forest structure (tree cover, soil nitrogen load, age), and partly to changes in hydrological conditions (Δsoil moisture, ΔVPD, Δradiation, ΔT). The latter are all likely to be strongly covarying as they are different representations of drought conditions, or known consequences (Δfire). The outcomes show some interesting relations (e.g. a sharp gradient in the relation with ΔVPD) which feel like they need to be further discussed or analyzed. It also seems to me that the machine-learning predominantly confirms -qualitatively- the conditions that were described in previously published analyses of these events, rather than new insights. But I am open to arguments from the authors. [Response]: Thank you for your valuable comments. Firstly, regarding your Comment 1, we fully agree that Δsoil moisture, ΔVPD, Δradiation, and ΔT exhibit strong covariation. However, as noted by the reviewer, they capture drought severity from different perspectives. For instance, soil moisture reflects soil drought, directly influencing plant water uptake and soil-plant-atmosphere interactions. In contrast, VPD represents atmospheric drought, which reduced plant transpiration and CO uptake across a large range of soil moisture conditions (Fu et al., 2022). Decoupling and evaluating the impacts of these different drought types on terrestrial vegetation productivity has been a key research focus (Verma et al., 2024; Zhong et al., 2023; Liu et al., 2020). Furthermore, through the RF-SHAP framework, which handle complex interactions among variables even in the presence of strong covariation (Wang et al., 2023), we identified the distinct contributions of these variables across different biomes, as shown in the revised Fig. 3. Therefore, we believe including these variables in our ML model to analyze their contributions to vegetation biomass changes provides meaningful insights. Most importantly, as shown in the variance inflation factor (VIF) results in Fig. S17, the collinearity among these variables is within a tolerable range, supporting their simultaneous use in the model. To address your second comment, we have added a discussion on VPD based on your suggestion, as follows: “The second-most influential climatic factor was ΔVPD, contributing 8%, with increasing VPD associated with decreasing biomass carbon (Fig. 3f,j), reflecting the fact that rising VPD is associated with tree mortality, reduced growth44, and fire activity28” (Page 10 lines 205–207 in the revised manuscript). For your last comment, we are sorry if the original manuscript gave the impression that our Random Forest-SHAP framework is qualitative. As noted in our response to Comment 3, SHAP is a game theory derived interpretability index that quantitatively explain how the dependent variables of machine learning models depend on predictors (Lundberg et al., 2020). This method has been widely adopted in recent studies published in Nature journals, such as Fu et al. (2024) and Li et al. (2022). However, since SHAP values do not directly provide the percentage contributions of the predictor variables, we applied L1 normalization to scale them into percentage contributions for easier understanding. Additionally, instead of using SHAP's built-in dependence plots to illustrate the effects of key explanatory variables, we chose to plot the partial dependence relationships between the target variable (i.e., biomass density change, Δbiomass) and these factors, as shown in the revised Fig. 3 above. We also revised the relevant descriptions in the revised main text to provide the specific percentages explained by these variables: “The most important climate-related variable affecting biomass carbon change (Δbiomass) in both temperate and boreal biomes was the change in surface solar radiation (ΔRadiation), accounting for 13% and 19% of the explained variability, respectively (Fig. 3a,f). However, partial dependence plots revealed that the effect of ΔRadiation on Δbiomass differed between the two biomes, with Δbiomass decreasing as ΔRadiation increased in temperate areas, while a positive relationship was found in boreal areas (Fig. 3c, g). The negative relationship between biomass and radiation in temperate regions may be associated with reduced soil moisture (SM) and increased vapor pressure deficit (VPD) (Fig. S5a,b). The Spearman rank correlation between L-VOD biomass estimates and SM at low-mid latitudes highlights the negative sensitivity of biomass carbon stock changes to soil moisture anomalies in this area (Extended Data Fig. 4). ΔSM contributed 10% of the explained variability, ranking among the higher climate factors, with soil moisture deficits (i.e., negative ΔSM) linked to negative biomass carbon changes (Fig. 3 a,d). The most significant factor positively influencing Δbiomass in temperate biome was tree cover with an average contribution of 18%, followed by soil nitrogen content at 11% (Fig. 3a, b, d). For the boreal biome, ΔRadiation is ranked as the most important factor affecting Δbiomass, which may reflect the limitation of solar energy on biomass accumulation, consistent with its role in constraining vegetation photosynthesis (gross primary production) for this region23. The U-shaped relationship, with the lowest point around ΔRadiation ≈ 0 J/m², is likely influenced by interactions with other factors, particularly disturbances related to forest area (Fig. S5c,d). The second-most influential climatic factor was ΔVPD, contributing 8%, with increasing VPD associated with decreasing biomass (Fig. 3f,j), reflecting the fact that rising VPD is associated with tree mortality, reduced growth44, and fire activity28. Notably, while forest loss due to stand-replacing wildfires ranked high (9%), non-fire forest disturbances such as logging and drought-induced mortality had a comparatively larger impact on Δbiomass in the boreal areas, accounting for 12% (Fig. 3f, h).” (Pages 9–10 lines 188–210 in the revised manuscript). [Reviewer #1 Comment 5] Similarly, the authors introduce in Figure 4 an analysis that includes three DGVMs. The basis for selecting these three remains unclear, and the chosen figures clarify very little about the actual processes captured by these models that would make them reproduce the observed ΔAGB. In fact, besides the reduced BC density coinciding with low SPEI or low soil moisture, the models reproduce very little of the observed signal. They of course know about hot and dry locations through their driver data (reanalyzed temperature and precipitation), but this has been shown in many studies. Here it matters that they explain changes in AGB. But the trend in ΔAGB per latitude is captured only partly by only one model (Fig 4f), while the aggregate changes even over large NH land areas is not reproduced by any (Table S5). It is therefore not surprising that they are not used for further insights on the processes driving the trends. Without this, I miss a justification to use these process-based models and the paper would be stronger without them. I would like to read a strong justification of their use, and the specific selection of models and co-authors that was made. [Response]: We greatly appreciate the reviewer’s professional feedback. On the first point of the comment, the reason why only these three DGVM models were chosen is that they were the only ones available at the time of the initial manuscript submission and analysis. However, we now have access to additional DGVMs provided by the TRENDY v12 project. To avoid any similar confusion from readers, we have incorporated 17 DGVMs with rapid update capabilities, including CABLE-POP, CARDAMOM, CLASS, CLM5.0, DLEM, EDv3, IBIS, ISAM, JSBACH, JULES, LPJ-GUESS, JPLml, LPJwsl, LPX-Bern, OCN, ORCHIDEE, and SDGVM, and have analyzed them in the main text as a mean ensemble. We have now revised the relevant descriptions in the introduction and updated Figure 4 to: “The L-VOD-derived biomass carbon density was also compared with the DGVM simulations from TRENDY v12 project, particularly focusing on recent trends.” (Pages 3-4 lines 95–97). On the second point of the comment, we fully agree with the reviewer’s point that DGVMs are predictably able to identify hot and dry conditions, as they integrate climate forcing data. This capability is particularly evident in areas where biomass loss is driven by extreme climate events, such as the biomass loss at low to mid-latitudes (30°-50°N) due to the extreme drought of 2022, as identified in this study. However, our focus is to highlight that the TRENDY DGVMs were unable to reproduce the negative biomass carbon trend in boreal regions from 2016 to 2022, instead showing a positive trend in regions where a high average annual burned area was recorded (Fig. 4h). This discrepancy is partly due to the fact that DGVMs do not simulate accurately actual wildfire occurrences and miss extreme fires. In contrast, DGVMs that use prescribed burned area data from CCI observed burned area to calculate emissions, such as ORCHIDEE, JULES, and OCN, show higher consistency with L-VOD estimates. As we considered this discrepancy between observations and models to be an important result, we have kept the DGVMs in the main text, and we hope that this justification is sufficient for your opinion. In addition, no DGVM simulates tree mortality due to drought and biotic disturbances. To emphasize our findings, we have removed Table S5, revised the descriptions of DGVM usage, and enhanced the relevant discussion in the main text, as follows: “However, the TRENDY DGVMs with rapid update capabilities did not mirror this trend, partly due to limitations in their prognostic fire modules in realistically reflecting wildfire occurrence and emissions. Yet, the three DGVMs models, JULES, OCN, and ORCHIDEE, which used prescribed burned areas from the ESA Climate Change Initiative (CCI-Fires) (Methods), revealed better consistency with L-VOD estimates compared to their prognostic fire versions (Text S4).”(Pages 12 lines 234–239). In method of the revised manuscript, we have refined the description of the usage of DGVMs: “Outputs of Dynamic Global Vegetation Models. We utilized the ‘cVeg’ outputs, which represent total live biomass as the sum of above- and below-ground biomass, from 17 DGVMs provided by the TRENDY v12 project67: CABLE-POP, CARDAMOM, CLASS, CLM5.0, DLEM, EDv3, IBIS, ISAM, JSBACH, JULES, LPJ-GUESS, JPLml, LPJwsl, LPX-Bern, OCN, ORCHIDEE and SDGVM. The TRENDY DGVMs performed three factorial simulations. S1 was driven only by CO , using a pre-industrial climate and land-use mask. S2 included both CO and climate while retaining a pre-industrial land-use mask. S3 incorporated CO , climate, and land use, all varying over time. We selected the S3 simulations, which were spun-up by recycling climate data from 1960 to 1969 until carbon pools were in equilibrium, and then applied varying CO , climate and land-use conditions over the period 1960 to 2022. We used the mean of these DGVMs in the main text, marking uncertainty by the range between maximum and minimum values. It should be noted that these DGVMs, which have rapid update capabilities, handle emissions using a prognostic fire module. In addition, we analyzed three DGVMs (JULES, OCN, and ORCHIDEE), which were the only models available through 2020 at the time of analysis and used prescribed burned areas from the ESA Climate Change Initiative (CCI-Fires) for calculating emissions68. This allowed us to focus on their consistency with L-VOD estimates before and after incorporating fire diagnostics (Text S4).” [Reviewer #1 Comment 6] My suggestion for further improvement of the manuscript would be to keep the opening two sections and Figures of the current work, but to focus the manuscript after that on the spatiotemporal analysis and contributing processes (such as in the strong final section on forest loss and degradation) instead of the proxies. For example, one could discuss specific geographic areas of the NH in relation to known disturbances, and go beyond drought and fire but also include harvest, regrowth, or shifting land-use. The mentioned areas in western Canada, Alaska and Scandinavia can be a start, but also western/central USA, Kazakhstan and its neighboring countries could be highlighted. Another angle to provide more context to the finding of a trend reversal is the fate of the carbon lost: how much of it has returned to the atmosphere and constitutes a true loss of sink, and how much has ended up in semi-permanent stores (such as buildings and furniture), possibly after salvage logging? [Response]: Thank you for your constructive suggestions. As addressed in our response to Comments 3-4, the Shapley analysis framework we used is a quantitative method for analyzing the driving factors behind changes in biomass carbon density. In this context, we would like to retain this section with the changes outlined above. Additionally, we agree with you that incorporating more regional analyses would enhance the manuscript, and to this end, we have added the following description: “Spatially, the transition to biomass loss revealed by L-VOD predominantly occurred in the south-central United States, Kazakhstan, and western Russia, while both loss and gain were observed in Europe (Fig. 1a, b).” (Page 4 lines 117-119) “This gross annual loss represents about 4% of the total biomass carbon in northern temperate regions (35.0542.79 PgC in 2016) (Fig. S3), with hotspots in the south-central United States, 28.65 southwestern Europe, and western Russia (Extended Data Fig. 2).” (Page 6 lines 146–149) We also sincerely thank the reviewer’s interesting suggestion to include the perspective on carbon fate. However, as the current manuscript is already quite lengthy, this aspect may be more suitable for future in-depth analysis. To acknowledge this, we have added the following statement to the discussion section of the revised manuscript: “Moreover, carbon losses from non-fire disturbances may not be immediately emitted into the atmosphere in the year they occur but could be delayed (e.g., timber products) or transferred to other carbon pools (e.g., from live biomass to dead wood and litter). In addition, forest fires have a combustion completeness that can be rather low, e.g. in Siberian larch and Siberian fir, or temperate broadleaved forest fires. Therefore, future studies should further quantify the fate of carbon to better contextualize net carbon stock changes within the framework of the global carbon budget and terrestrial carbon emissions and removals12,14.” (Pages 17-18 lines 361–367) [Reviewer #1 Comment 8] Overall, I find the longer dataset and its dual-period trend analysis sufficiently interesting to warrant publication. But I find the current manuscript to lack innovation in the analysis of drivers, processes, and regional differences in them due to the use of the uninformative ML method and the uninformative trendy models. The re-discovery of various known recent anomalies moreover asks for a more quantitative comparison to the previous papers published on them. With these elements in place, the work would likely also be more appealing to the wider users of this journal, and establish itself as the go-to dataset for analysis of recent carbon cycle changes. [Response]: Thank you for recognizing the value of our manuscript. We have addressed your comments regarding ML and Trendy DGVMs in the revised manuscript, and we hope these changes meet with your approval. Our study presents a novel finding: the shift in biomass carbon stocks in northern ecosystems from a positive to a negative trend during 2016–2022. Although a decrease could be detected in a previously published satellite-based model using optical and shortwave microwave data up to 2019 (Xu et al., 2021), its magnitude (0.63 PgC y-1) appears too large to be reconciled with high resolution atmospheric inversions results (Ke et al., 2024). Our study shows a carbon loss in biomass of −0.20−0.11 PgC y-1, which is consistent with inversion −0.26 data showing a sink reduction in biomass, deadwood and soils of 0.3 PgC y-1 between 2016 and 2023. Furthermore, here we used an interpretable machine learning regression model to analyze the underlying drivers. We found that the sign of net biomass carbon loss over northern ecosystems revealed by L-VOD is comparable to other satellite datasets but differs from the inventory data analyzed by Pan et al. (2024). However, while the sign of change is generally consistent across different datasets, only L-VOD, with its rapid update capability, was able to detect the recent reversal in biomass carbon stock trends in northern ecosystems during 2016-2022. The inventory-based integrated data differ from the remote sensing products due to the exclusion of unmanaged forests and non-forest categories, as well as the lack of consistent inventory data in some countries, such as Russia, since the early 2000s. The CCI biomass maps have inter-annual biases due to the use of different satellite sensors in different years (Santoro et al., 2024)53, and the dataset from Xu et al. (2021)12 presents saturation in dense forests due to the use of optical and high-frequency observations, limiting their performance for consistent multi-year change detection (Wigneron et al., 2024). To address your comments, we have also included the following description in the revised manuscript: “While a fair comparison with most previous estimates is not straightforward due to differences in the time period studied, we compiled data from several studies providing biomass carbon estimates over northern ecosystems, overlapping where possible, with the post-2016 decline period identified by L-VOD (Extended Table 1). The recent signs of net biomass carbon loss observed by L-VOD were comparable to those reported by other remote sensing products, particularly for the northern temperate biome, where L-VOD estimates ranged from -0.33 PgC yr-1 to -0.06 PgC yr-1, as estimated using the biomass datasets from Xu et al. (2021)12 and the ESA CCI biomass project54. For the boreal biome, L-VOD and CCI showed biomass carbon neutrality, whereas data from Xu et al. (2021)12 indicated a significant net loss. The inventory-based integrated data55 also indicated a reduction in the biomass carbon sink of northern ecosystems during the 2010s compared to the 1990s and 2000s (Extended Table 1). However, while the overall results are generally consistent between the different data sets, L-VOD, with its rapid update capability, was able to capture the latest reversal in biomass carbon stock trends in northern ecosystems during 2016-2022. The inventory-based integrated data differ from the remote sensing products due to the exclusion of unmanaged forests and non-forest categories, as well as the lack of consistent inventory data in some countries, such as Russia, since the early 2000s. The CCI biomass maps have inter-annual biases due to the use of different satellite sensors in different years53 and the dataset from Xu et al. (2021)12 presents saturation in dense forests due to the use of optical and high frequency observations, limiting their applicability for consistent multi-year change detection35.” (Pages 16–17 lines 324–341) [Reviewer #1 Comment 9] Below, I subsequently add some textual suggestions to consider for a possible resubmission of the manuscript: [Response]: Thank you for your valuable suggestions. They have all been incorporated into the revised manuscript as below. [Reviewer #1 Comment 10] - The methodology section ‘Benchmark maps of vegetation carbon density’ discusses the importance of calibration against existing biomass maps. To overcome the inevitable uncertainties and biases in existing maps, the authors use spatial gradients of the existing maps. Since this seems to be an important step in the methodology, I suggest adding a brief explanation on the use of spatial gradients and its effect on the conversion from L-VOD to AGC. [Response]: Thank you for your advice. We have added a sentence in this section to clarify this. “Any systematic errors in the AGC reference maps are also a major external source of uncertainty in the L-VOD-derived biomass maps. While individual reference maps may have biases, using an averaged calibration from multiple maps across various regions as the final biomass carbon density mitigates potential biases associated with any single map, as outlined in the L-VOD biomass data production and uncertainty analysis.” (Page 24 lines 574–578) [Reviewer #1 Comment 11] - Figure 4d shows a combination of several datasets: L-VOD, SPEI and soil moisture. Although I appreciate the attempt to bring these datasets together in one plot, it took me some time to grasp it. Specifically, I think that overlying three distinct line plots, with an inverted y-axis for SPEI, could be confusing to some readers. Hence, I suggest to plot the SPEI in a small separate graph just above the lines representing L-VOD. [Response]: Thank you for your suggestion. We agree that an inverted y-axis for SPEI could be confusing for some readers. To enhance the readability of this figure, we have reoriented SPEI to display the standard direction and adjusted the transparency to ensure clarity. The revised figure is presented in our response to your Comment 5. [Reviewer #1 Comment 12] - Lines 125-126, pixels could be written in plural. [Response]: Thanks for the suggestion. We have revised the wording accordingly. [Reviewer #1 Comment 13] - Line 175: ‘relatively’ could be written as ‘relative’. [Response]: Thanks for the suggestion. We have revised the wording accordingly. [Reviewer #1 Comment 14] - Line 223: could ‘temperature’ be a mistaken for ‘temperate’ here? [Response]: Thank you for pointing that out. Indeed, it should be 'temperate' in this context. [Reviewer #2 Comment 1] I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts. [Response]: We thank reviewer 2 for taking the time to co-review the manuscript. [Reviewer #3 Comment 1] Title: Large biomass carbon losses from droughts in the northern temperate ecosystems during 2016-2022. This is a remote sensing-based vegetation biomass study which is critical for current carbon cycle science. The study provides quantifications of the impacts of drought and wildfire of biomass change. If possible, I would suggest the authors to provide some forest area change information. [Response]: We sincerely thank the reviewer for taking the time to provide valuable and constructive feedback, which has greatly helped improve our manuscript. All of your comments have been carefully considered in the revised manuscript. Please find our detailed responses to each comment below. Following your suggestion, we have added information on forest area changes. Unfortunately, there is currently no suitable forest area dataset updated to 2022. Even the MODIS tree cover data (i.e., MOD44B) has many blank pixels in northern ecosystems after quality control. Instead, we used global forest loss data from the Global Forest Watch platform, which has been incorporated into Extended Data Figure 4c in the revised manuscript. [Reviewer #3 Comment 2] Abstract Nice and clear. But please add the spatial resolution 25km in the Abstract. Also, I noticed that the 4% biomass drop during 2016-2022 in the main text is a very important information. I would suggest adding that 4% in the Abstract if possible, and hopefully you can also point out where it happened (e.g., Russia, U.S., Amazon, etc.). [Response]: Thank you for your positive feedback on our abstract and for your thoughtful suggestions. Regarding your first suggestion, we have incorporated the resolution information as recommended. For the second suggestion, the 4% figure actually reflects the yearly gross carbon losses as a proportion of the northern temperate biome’s total biomass carbon stocks, highlighting the substantial inter-annual variability in biomass carbon within this biome. To address your comment, we have clarified this information in the abstract to convey its significance and have specified the exact locations in the main text. In abstract, we have added: “The annual mean gross loss of 4% of biomass carbon in this region during 2016-2022 reflects high interannual variability, with significant losses associated with droughts and a further drop of −0.60−0.47 PgC in the very dry year of 2022.” −0.75 In main text, we have added: “Notably, a gross carbon loss of −1.50−1.26 PgC yr-1 occurred in temperate biomes, −1.77 partially offset by a gross gain of +1.24+1.44 PgC yr-1. This gross annual loss +1.09 represented about 4% of the total biomass carbon in northern temperate regions (35.0542.79 PgC 28.65 in 2016) (Fig. S3), with notable losses in the south-central United States, southwestern Europe, and western Russia (Extended Data Fig. 2).” (Page 6 lines 145–149 in the revised manuscript). [Reviewer #3 Comment 3] Introduction Line 93: \"... exploring their causes in relation to recent drought and wildfire\" Seems later you also mentioned forest cutting in several places (Line 299, 656)? If cutting impact on biomass loss is smaller than fire, then it is OK to omit it, otherwise, maybe it is good to mention why cutting is omitted (e.g. lack of data, not a focus of this paper, cutting is too small, or cutting can not be captured by remote sensing while fire can be captured?). [Response]: Thank you for your valuable feedback. While fire significantly impacts biomass changes in boreal biomes, both our machine learning and attribution analyses showed that other disturbances, such as forest cutting and insect outbreaks, have relatively greater impacts. However, the lack of spatial data defining forest harvest disturbances currently limits our ability to further attribute biomass carbon loss to increased harvests versus increased natural disturbances. Both are also related with salvage harvest following beetle attacks in Europe. To address your comments, we have revised this sentence and included the discussion in the manuscript. In Lines 95-96, we have changed the original sentence as: “From the L-VOD-derived biomass carbon estimates, we analyse the changes in live biomass carbon at the regional scale and identify factors that drive the spatial variability in biomass carbon density changes across different biomes.” In discussion, we have added: “Future studies could apply the approach of this study to isolate biomass carbon loss from forest cutting, provided spatial data defining cutting boundaries or reliable proxies (e.g., forest management boundaries) become available58.” (Page 17 lines 359–361). [Reviewer #3 Comment 4] Introduction Line 113-114: \"Both sources indicated that live biomass carbon in northern ecosystems peaked in 2016,\" Suggest providing spatial resolution of those two data products. [Response]: Thank you for your suggestion. As recommended, we have included information on the spatial resolution of these two data. “One from a high-frequency C-band VOD product40 and the other from global annual biomass maps12 with native spatial resolutions of 25 km and 10 km, respectively.” [Reviewer #3 Comment 5] Line 125-126: \"The partly burned areas are pixel with an annual fire area fraction > 1% within a 25 km grid from 2005 to 2022\" Looks like this gives about 27.5/49.9=55% of the global vegetation area being partially burned. It is far from my thinking. Can you try to find a realistic total area of fire disturbance? The \"partly burned\" may be too loose. Effects of other changes within the partial burned areas are counted as fire effect, I guess. Suppose a fire can kill 50% of live biomass, then 1% area burn could only kill 0.5% total biomass. If bigger C loss was detected, say 1%, could it because of some other factors as well? [Response]: Thank you for your valuable comments. The purpose of defining \"partly burned areas\" was to further investigate the impact of all fires, not only stand-replacing ones associated with forest area loss. Through the reference to \"unburned areas,\" defined as 25 km SMOS pixels with an annual burned area fraction of <1% (driven by factors other than fire), we found that fire indeed amplified the biomass carbon dynamics in partly burned areas. In this context, increasing the threshold for defining partly burned areas would inevitably result in unburned areas also encompassing fire-affected regions. We apologize for not clearly explaining this in the original manuscript. To enhance clarity, we have revised two expressions in the original manuscript as follows: For the title of Fig. 1, we have revised the original statement to: “The 'unburned' areas are pixels where the annual fire area fraction was consistently < 1% each year from 2005 to 2022 within a 25 km grid, while all other pixels are considered partly burned areas14.” Additionally, we have revised the following sentences in the main text for better clarity: “In addition to stand-replacing events associated with forest area loss, we further investigated the impacts of all wildfires on biomass carbon stock changes in northern ecosystems. We considered ‘unburned’ areas defined as areas that include pixels with an annual fire area fraction < 1% within a 25 km grid between 2005 and 2022, while all other pixels were considered as partly burned areas, as in Fan et al. (2023)14.” (Page 15 lines 293–297). [Reviewer #3 Comment 6] Line 141-142:\"Notably, temperate ecosystems sustained a gross carbon loss of −1.50−1.77−1.26 PgC yr-1 which represented about 4% of the total biomass carbon\" Are you saying 4% biomass carbon loss in one single year? Please double check. I guess it should be 2016-2022 total biomass carbon loss? [Response]: Thanks for the comment. As explained in response to your Comment 2, the 4% figure actually represents the proportion of annual gross carbon loss relative to the total biomass carbon stock in the northern temperate biome from 2016 to 2022, highlighting the substantial interannual variability of biomass carbon in this biome. To improve clarity, we have revised the sentence as follows: “Notably, a gross carbon loss of −1.50−1.26 PgC yr-1 occurred in temperate biomes, −1.77 partially offset by a gross gain of +1.24+1.44 PgC yr-1. This gross annual loss represents about 4% +1.09 of the total biomass carbon in northern temperate regions (35.0542.79 PgC in 2016) (Fig. S3), with 28.65 hotspots in the south-central United States, southwestern Europe, and western Russia (Extended Data Fig. 2)” (Page 6 lines 145–149). [Reviewer #3 Comment 7] Figure 2: Are the gain and loss were calculated based on land cells? [Response]: Thank you for pointing this out. Indeed, the gain (green line) and loss (red line) in Fig.2c are calculated on a per-pixel basis and then summed across longitude within latitude bands. We have clarified this point further in the caption of Figure 2 to enhance clarity. “c, Latitudinal sums of gross gains and losses in L-VOD biomass estimates from 2016 to 2022, derived from per-pixel gain and loss maps by cumulating positive and negative changes for consecutive years.” (Page 8 lines 163–165). [Reviewer #3 Comment 8] Line 176-179: \"These factors include climatic, vegetation, edaphic, and anthropogenic disturbances (Table S4). ...... (Fig. 3)\" This section talks about drivers. if possible, I would suggest adding one sentence to give the relative contribution (X%) of each of the four driving factors to biomass loss. Seems to me Fig.3 is saying that Fire is the biggest influencing factor, and human activity is the smallest influencing factor. But it is for something about 'spatial variability', not the carbon loss quantity. [Response]: Thank you for your suggestion. Indeed, we used an interpretable Shapley framework to attribute the spatial patterns of net biomass carbon changes in northern ecosystems to 15 driving factors, allowing us to estimate their contributions. Our results show that while forest loss from stand-replacing wildfires ranked high, non-fire forest disturbances such as logging and drought-induced mortality had a comparatively greater impact on biomass carbon changes in the boreal biome, accounting for 12%. Following your suggestion, we have now added the relative contribution percentages in this section: “The most important climate-related variable affecting biomass carbon change (Δbiomass) in both temperate and boreal biomes was the change in surface solar radiation (ΔRadiation), accounting for 13% and 19% of the explained variability, respectively (Fig. 3a,f). However, partial dependence plots revealed that the effect of ΔRadiation on Δbiomass differed between the two biomes, with Δbiomass decreasing as ΔRadiation increased in temperate areas, while a positive relationship was found in boreal areas (Fig. 3c, g). The negative relationship between biomass and radiation in temperate regions may be associated with reduced soil moisture (SM) and increased vapor pressure deficit (VPD) (Fig. S5a,b). The Spearman rank correlation between L-VOD biomass estimates and SM at low-mid latitudes highlights the negative sensitivity of biomass carbon stock changes to soil moisture anomalies in this area (Extended Data Fig. 4). ΔSM contributed 10% of the explained variability, ranking among the higher climate factors, with soil moisture deficits (i.e., negative ΔSM) linked to negative biomass carbon changes (Fig. 3 a,d). The most significant factor positively influencing Δbiomass in temperate biome was tree cover with an average contribution of 18%, followed by soil nitrogen content at 11% (Fig. 3a, b, d). For the boreal biome, ΔRadiation is ranked as the most important factor affecting Δbiomass, which may reflect the limitation of solar energy on biomass accumulation, consistent with its role in constraining vegetation photosynthesis (gross primary production) for this region23. The U-shaped relationship, with the lowest point around ΔRadiation ≈ 0 J/m², is likely influenced by interactions with other factors, particularly disturbances related to forest area (Fig. S5c,d). The second-most influential climatic factor was ΔVPD, contributing 8%, with increasing VPD associated with decreasing biomass (Fig. 3f,j), reflecting the fact that rising VPD is associated with tree mortality, reduced growth44, and fire activity28. Notably, while forest loss due to stand-replacing wildfires ranked high (9%), non-fire forest disturbances such as logging and drought-induced mortality had a comparatively larger impact on Δbiomass in the boreal areas, accounting for 12% (Fig. 3f, h).” (Pages 9–10 lines 188–210 in the revised manuscript). [Reviewer #3 Comment 9] Line 254: \"For the period 2016-2020, biomass carbon losses were −2.96−3.38−2.65 PgC yr-1,\" This carbon loss amount is significant. in Line 143, you mentioned temperate region has a total biomass of 35.05 PgC in2016. I didn't find the number for boreal regionas. Try to add boreal total biomass carbon if you have. Also I believe it shouldn't be a per year value. (OK, coming back to this comment again, I realized maybe your carbon loss is not a net loss?) [Response]: We agree with the reviewer that the loss of (−2.96−2.65 PgC yr-1) is indeed a −3.38 substantial amount. As you understood after coming back to this comment again, this figure represents the average annual gross biomass carbon loss (not net loss) across the entire northern land from 2016 to 2022, calculated by summing the negative differences between the consecutive years and dividing by six. To avoid any confusion, we have clarified the definitions of gross gains and total losses in the methodology section. Additionally, following your suggestion, we have also added the total biomass carbon stock for the boreal region as follows: In Methods part, we have added: “Gross losses (gains) were estimated by accumulating negative (positive) differences between the consecutive years over each pixel, following the method of Brandt et al. (2018)53.” (Pages 29 lines 734–736). In main text, we have added: “Annual gross losses in boreal biomes amounted to about 2% (−1.29−1.22 PgC yr-1) of the total biomass carbon stock (54.2661.69 PgC in 2016), with higher −1.47 48.88 losses in the Siberian region” (Page 6 lines 150–151). [Reviewer #3 Comment 10] Line 254-268: The whole paragraph gives me an impression that fire is not the #1 factor of biomass carbon loss compared to cutting and insect. Please double check the PgC numbers and X% numbers to make sure that fire effect is bigger than forest area loss, cutting loss, and insect loss. [Response]: Our machine learning interpretability analysis indicates that while forest loss due to stand-replacing wildfires ranked high, non-fire forest disturbances such as logging and drought-induced mortality had a comparatively greater impact on biomass carbon changes in boreal areas. Therefore, your observation is correct. To avoid any confusion, we have reorganized this paragraph as follows: “The latter, including processes such as clearcutting and insect outbreaks, resulted in ~ 26% (−0.76−0.68 PgC yr-1) of the gross live biomass losses across the northern ecosystems, with −0.85 ~ 70% (−0.53−0.48 PgC yr-1) and ~ 29% (−0.22−0.20 PgC yr-1) occurring in boreal and temperate −0.60 −0.25 forests, respectively. In contrast, stand-replacing fires (4.51×107 ha) corresponding to ~ 33% of the forest loss area caused a gross biomass carbon loss of −0.37−0.33 PgC yr-1 in northern ecosystems, −0.41 of which ~ 87% (−0.32−0.29 PgC yr-1) occurred in boreal forests (Extended Data Fig.5). While this −0.36 result confirms that wildfires are a key source of forest disturbance in boreal forests relative to other ecosystems50, other types of forest area disturbances within this biome were ~ 1.7 times larger than those caused by stand-replacing fires.” (Page 15 lines 284–292). [Reviewer #3 Comment 11] Line 284-287: \"... indicating that Landsat forest loss data cannot capture small-scale fire-burned forest loss patches and/or that disturbances such as drought and low-intensity understory wildfires may not always result in stand replacement but cause biomass carbon losses, consistent with the studies of Fan et al. (2023)14 and Rogers et al. (2015)50.\" I agree that most wildfires are not stand-replacing. But I am curious that if 30m data cannot capture such fire-burned forest loss, how can 25km data capture the fire-burned forest loss? OK, I see you are using MODIS fire data (Ref # 42, 250m?), then I suggest add one sentence in main text saying that MODIS has a better way to capture fire events. (And it is multiple year data with global coverage, better than Landsat single year data) [Response]: Thank you very much for your comments. Our attribution results indicate that the gross biomass carbon losses caused by stand-replacing fires (calculated as the forest area loss due to stand-replacing fires multiplied by the biomass carbon density before fire) were significantly lower than when considering all fire events. This suggests that the Landsat forest loss data may fail to capture small-scale fire-burned forest loss patches and/or that disturbances such as drought and low-intensity understory wildfires may not always result in stand replacement but still contribute to biomass carbon losses. To emphasize that this is a hint and to make the manuscript more rigorous, we have revised this statement. In Page 16 lines 307–312, we have revised the original sentence as: “L-VOD-derived gross biomass carbon losses in partly burned areas of the boreal regions were 2.7 times higher than those estimated solely from stand-replacing fires (−0.30−0.27 PgC yr-1), implying that Landsat forest loss data may −0.34 not capture all small-scale fire-burned forest loss patches and/or that disturbances such as drought and low-intensity understory wildfires may not always result in stand replacement but cause biomass carbon losses” Additionally, following your suggestion, we have included a description of the advantages of MODIS burned area data: “This widely used dataset provides global coverage, making it ideal for monitoring fire activity across diverse ecosystems, while its high temporal resolution enables timely detection over large areas13,14.” (Page 28 lines 710–712) [Reviewer #3 Comment 12] L295-297: \"The downward trend from 2016 to 2022 is reflected by a net loss of biomass carbon over the northern temperate biome (−0.26−0.33−0.17 PgC yr-1) and boreal Eurasia (−0.03−0.04−0.02 PgC yr-1).\" Here it is obvious that temperate biome has bigger biomass carbon loss than boreal, right? But in Line 254 you mentioned -2.96 PgC for total loss, in Line 142 you mentioned -1.5 PgC. That makes me thinking boreal might lost ~1.46 PgC, almost equal to temperate region. Maybe I read the data wrong, but please double check. [Response]: Thanks for the comment. Following your suggestion, we have rechecked each figure and found no errors. First, from a net change perspective, your understanding is entirely correct: the net loss in the northern ecosystem from 2016 to 2022 is primarily reflected in the northern temperate biome (−0.26−0.17 PgC yr-1), exceeding the net loss in boreal Eurasia (−0.03−0.02 PgC −0.33 −0.04 yr-1). Second, it is important to note that this net change results from the offset between gross losses and gross gains. Although the boreal region’s annual gross losses (−1.29−1.22 PgC yr-1) are nearly −1.47 comparable to those in the temperate region (−1.50−1.26 PgC yr-1) during 2016-2022, the annual −1.77 gross gains in the boreal region amount to +1.32+1.46 PgC yr-1, surpassing those of the temperate +1.26 region at +1.24+1.44 PgC yr-1. +1.09 [Reviewer #3 Comment 12] Line 300-302: \"Notably, 2022, the year most affected 300 by drought, showed the sharpest decline in biomass carbon stocks over the 2010-2022 period, with a net change of −0.67−0.69−0.41 PgC.\" ~2% biomass loss in the year 2022 is still a huge number at global scale. [Response]: We agree with the reviewer that the net decrease of −0.67−0.41 PgC in 2022 −0.69 represents a substantial figure, highlighting the vulnerability of vegetation carbon stocks to emerging climate-induced disturbances. We have rechecked and confirmed this calculation, noting that this value represents 0.69% of the total biomass carbon stock in the northern ecosystem (97.27113.83 PgC in 2016). 84.35 [Reviewer #4 Comment 1] I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts. [Response]: We appreciate Reviewer 4’s time and effort in reviewing the manuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "My review is included as PDF file.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "Dear authors, I am satisfied with your revision. You have added substantial explanations and new info to the revision. All my questions are answered. Great work. Recommendation: Accept", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We sincerely appreciate your valuable comments and suggestions, which have helped us refine the manuscript. We have fully addressed each of your requests, and our point-by-point responses are provided below. The requests are reiterated in italics, with our responses and revisions highlighted in blue and red text, respectively. [Reviewer #1 Comment 1] I would suggest critically revisiting the formulations around ‘carbon sinks’, ‘stock changes’ and ‘emissions’. To avoid that readers confuse changes in biomass stocks with atmospheric carbon sinks or sources, I would like to suggest the following edits: [Response]: Thank you for your helpful suggestions. Besides addressing them in the responses below, we have also reviewed the entire manuscript to ensure clarity in terminology. - Abstract: I urge the authors to use the word “biomass carbon stock” only together with “live” since that is what was observed. Soil carbon is not part of the story, and neither is coarse woody debris (CWD) or harvested timber. For a layperson your finding might imply that the lost biomass stocks are now in the atmosphere and we have reached a NH source, which is not the case. [Response]: Adopted! We fully agree that adding live makes the terminology more precise. - Line 51: I suggest rephrasing “the sustainability of carbon sinks in this key region” to “the sustainability of the current large terrestrial carbon sink in this key region (1.4- 2.0 PgC/yr Friedlingstein et al., 2024)” to reinforce the fact that we are still drawing down substantial amounts of atmospheric carbon in the NH biosphere. [Response]: Thanks. We have made the revisions as per your suggestion and cited the paper you mentioned to enhance the background information. - Line 102: “a net mean sink”: I suggest to not refer to biomass changes as sources or sinks, since they have only an indirect relation to uptake or loss from the atmosphere, and this might lead to confusion. Same in line 282. [Response]: Thank you for pointing this out. Your suggestions have been carefully considered and adopted. For Line 102, we have revised \"boreal regions were a net mean sink of +0.15+0.15 PgC yr-1, while +0.16 northern temperate regions were almost at a state of neutrality with a net flux of +0.03+0.04 PgC +0.02 yr-1 \" to \"boreal regions had a net increase of +0.15+0.15 PgC yr-1 in live biomass carbon, while +0.16 northern temperate regions remained nearly neutral with a net change of +0.03+0.04 PgC yr-1\" +0.02 For Line 282, we have revised \" in reducing the carbon sink in boreal biomes \" to \"in reducing the live biomass carbon stocks in boreal biomes\" - Line 235 “…partly due to limitations in their prognostic fire modules in realistically reflecting wildfire occurrence and emissions.” I think the leap to emissions might have been too fast here, since these were not assessed in this manuscript. I would suggest keeping it at “loss of live biomass”, which is what you retrieved from L-VOD. The hypothesis that their emissions are all wrong is likely to be true, but there might be canceling errors in emission factors, burning completeness, and even in delayed breakdown of CWD. Also here, substantial salvage logging post-fire might contribute to the observed live biomass change, without direct emissions to the atmosphere. [Response]: We appreciate your patient guidance. We agree with your point regarding the compensating errors in DGVMs. We have now revised the sentence to: \"partly due to limitations in their prognostic fire modules in realistically reflecting wildfire occurrence and intensity as well as the associated loss of live biomass\" - Line 371: The final paragraph, or one early in the paper, should mention the current size of the NH land sink (from GCB for example) to give some context to the statements on reversing or disappearing sinks. [Response]: Thanks. Following your suggestion, we have added the GCB estimate of the current NH land sink in the introduction, as: \"Northern ecosystems (north of 30°N) hold about 41% of the world’s forest area1, and are critical for mitigating global warming2-3, contributing an estimated 1.4–2.0 PgC annually to the global terrestrial carbon sink4.\" [Reviewer #1 Comment 2] As mentioned, the RF-SHAP approach is sufficiently explained now. The additional information, including the attempt to perform Multiple Linear Regression, is useful and appreciated. However, the figures of the analysis (3 and 4) still raise a few questions: - In Figure 3h, the Boreal non-fire influence seems very large on the y-axis, and it has a very steep slope. Does that suggest it is a very sensitive parameter (process) to include in modeling of forecasting the Δbiomass? Radiation exceeds it in the SHAP contribution %, but it’s relation with Δbiomass seems weaker in magnitude and less interpretable (U-shape)? It would help the reader if you interpret the outcomes in panels b-e and g-j of this figure, to what we need to understand or focus on as a community. That would truly leverage the power of the RF-SHAP analysis in my opinion. [Response]: We agree that further exploring the implications of partial dependencies plots (PDPs) (i.e., Figs. 3b-e and g-j) is valuable. At the same time, we would like to emphasize that PDP illustrates the model’s response to a given variable under specific conditions, assuming all other factors remain unchanged. The steep slope and large y-axis range observed in the PDP of boreal non-fire influence indicate that the model’s predicted Δbiomass changes more sharply when this variable varies, keeping others fixed. However, this represents only a localized effect and does not necessarily mean it is the most influential variable overall. In contrast, SHAP rankings (Fig. 3a,f) quantify each variable’s total contribution across all predictions, including direct effects and interactions. Radiation ranks highest in SHAP importance, meaning it has the greatest overall influence on Δbiomass predictions. However, its PDP exhibits a narrower y-axis range, suggesting that its effect remains relatively consistent across its range rather than fluctuating sharply at specific values. This suggests that while Δbiomass in boreal biomes may exhibit strong localized sensitivity to non-fire drivers, the impact of radiation is more consistent and widespread. Therefore, both analyses should be considered together. In consideration of your next comment, we have adopted SHAP-based PDPs (model’s response to a given variable capturing both their direct effects and interactions with other variables) instead of the smoother traditional PDPs (model’s response to a given assuming all other factors remain unchanged) in the revised manuscript. In addition, in order to clearly express the purpose of PDP use, we have added: \"To elucidate the nonlinear relationships between Δbiomass in the two ecosystems and their four key explanatory variables, SHAP-based partial dependence plots (PDPs) were generated, showing how SHAP values change with each variable while accounting for both its direct effects and interactions with other predictors68. By combining SHAP rankings with SHAP-based PDPs, we not only identify variable importance but also gain insights into how the top variables shape Δbiomass patterns.\" (Page 25 lines 714-719) We have also added the following clarification to the caption of Fig.3: \"The partial dependence plots illustrate how individual variables contribute to the predicted Δbiomass, capturing both their direct effects and interactions with other variables (Methods).\" - Fig 3 has changed from the original, now showing nearly no spread in the relations because the plot shows 1000 lines (mostly overlapping) instead of the dots in the previous version. What has changed, and why is this decision taken? The previous figure to me fit better with the idea that relations were not perfectly known, and that both explanatory power and uncertainty remained. Now they show that your optimal solution is very rigid and all your efforts lead to the same set of relations. But that does not make them more true or credible, it just shows that your setup does not contribute to the spread in the outcome. Would you agree with this more nuanced interpretation? If so, please consider reverting back to the previous design with dots. [Response]: Thank you for your valuable comments. We initially switched from SHAP-based PDPs (previous version) to traditional PDPs (current version) to make the y-axis representation more intuitive, as the former displays SHAP values while the latter directly displays the modeled Δbiomass. However, we acknowledge that this change reduced the variability previously conveyed by SHAP-based PDP dots, which inherently captured interactions with other variables. The smoother appearance of traditional PDPs results from averaging predictions across all data points, thereby reducing variability. Given that SHAP-based PDPs better reflect variable influence by incorporating feature interactions, we agree with your suggestion to revert to the previous dot-based design. To further aid readers in interpreting SHAP PDP values, we have also added the following explanation: \"The Y-axis represents the SHAP value for the corresponding predictor (X-axis), where positive values contribute to higher predicted Δbiomass values, and negative values to lower ones.\" - Would it help the interpretation of the RF-SHAP analysis if you split the dataset into gains and losses before training? It is likely that these two are not associated with the same drivers, nor responses. This is visible from the panels with steep changes in the slope, where a parameter only contributes for either a positive or negative anomaly. And contributions to negatives (losses) more abundant than positives (gains) it seems. Did you try this, would it help? I would understand if this is too large of a methodological challenge to do, and in that case your considerations and some reflections of them in the text would suffice. [Response]: Thank you for this thoughtful suggestion. We agree that biomass gains and losses likely reflect different ecological processes and drivers. For example, gains are often associated with secondary forest growth or post-disturbance recovery, whereas losses are more commonly linked to deforestation, degradation, or other disturbances. Moreover, a recent study has highlighted that understanding biomass gains is an emerging research focus (Feng et al. 2024), as most work emphasizes losses. We acknowledge that the traditional PDPs previously used may have given the misleading impression that variables contribute more strongly to biomass losses, due to averaging effects across the dataset. The updated SHAP-based PDPs better preserve the asymmetry between gains and losses, and this pattern is no longer observed. Following your suggestion, we have added the following description in the revised manuscript: \"Here, we aim to understand the underlying processes driving net changes in live biomass carbon density (Δbiomass) in the northern region during 2016–2022, considering both recovery and disturbance factors. While this study focuses on net changes, a separate analysis of biomass carbon losses56 and gains57 may offer additional insights and merits more detailed investigation in future work.\" (Page 25, lines 693-696) - I am likewise surprised by the shown relation with soil moisture: it seems less nonlinear than I would expect based on theory (a threshold, with an S-shaped curve is more typical), and both positive and negative soil moisture changes are associated with biomass loss (but less so for increased soil moisture). Can you comment on this? I don’t expect a full revision of the analysis, but would like to read your views. [Response]: We appreciate your feedback. The relationship between soil moisture (SM) and biomass in the previous PDP may have been influenced by the limitations of traditional PDPs, which fix the effects of other variables in an averaged form and may not fully capture nonlinear threshold behaviors. To address this, we have now adopted SHAP-based PDPs, which more accurately reflect the actual contributions of variables by considering both direct effects and interactions. As shown in the revised Fig. 3(above), the updated curve exhibits a more pronounced non-linear pattern and shows signs of saturation at high SM levels. Moreover, the SHAP PDPs revealed that decreases in SM were generally associated with lower predicted biomass, while increases in SM led to higher predictions, aligning better with theoretical expectations. In the revised manuscript, we have updated the following descriptions and references to emphasize the importance of SM. \"Additionally, our interpretable machine learning analysis identified change in radiation and soil moisture as the two most important climatic contributors influencing the spatial variability of live biomass carbon changes in the northern temperate biomes. Knowing that increased incoming radiation can further exacerbate drought formation, the high ranking of the latter underscores the important role of soil moisture stress in vegetation growth and regulating carbon fluxes in the region55-56. \" (Page 14, lines 356-360) - The caption of Figure 4d states “soil moisture trends”, and the labels on the legend should say the same, as now it suggests there is negative soil moisture (red shaded). [Response]: Thank you for pointing this out. We have corrected the label issue accordingly. - Extended data Figure 1: For (a), the log-log plot feels not so useful, while the second one seems leveraged by three large points. What if you leave these points out, does the relation hold? For (a), could you normalize the change to the total country stock, or forested area and then plot on a linear scale? [Response]: Thank you for your valuable comments. We agree that using a linear scale for Fig.1a would be more intuitive. We would also like to note that this figure presents the total live biomass carbon stock across countries in the study area, rather than temporal changes. This type of national-level comparison has been used in previous studies to assess the spatial consistency between remote sensing estimates and FAO data (Xu et al., 2021; Spawn et al., 2020). We tried applying the linear scale as suggested (figure below), but large cross-country differences made the plot difficult to interpret. For these reasons, we have chosen to retain the original scaling in Fig. 1a. To address your comment, we have revised the text as follows to better clarify our intent. \"Scatter plot of national total carbon stocks in live vegetation biomass estimated by L-VOD and reported by FAO for 2015, as done in Xu et al. (2021)11.\" Regarding Fig. 1b, the three countries indeed have a strong influence on the results, as they account for approximately 70% of the study area. However, even if the correlation is lower, it remains statistically significant even when the three countries are excluded (see Fig.1b below). Moreover, the tendency of L-VOD to underestimate UNFCCC-reported biomass changes in many countries may be explained by national forest inventories (NFI) primarily focusing on mature forests, while overlooking biomass carbon losses from degradation and deforestation in surrounding areas, as discussed in detail in Text S3: \"This underestimation can be partly explained by the fact that NFI data represent long-term biomass carbon change, since several years (typically >5 years) are required to revisit all plots within a national network. In addition, NFI primarily captures local biomass carbon changes in intact, mature, or old-growth forests, and does not account for carbon losses due to deforestation or degradation in the surrounding areas. Another reason could be that the NFI only samples within managed forest areas. We used intact forest boundaries as a proxy for managed forest, which may differ from what is actually applied by NFI in different countries (Grassi et al., 2023). \" [Reviewer #1 Comment 3] I would like to ask the authors to reflect on the following questions in the discussion: - In Fig 3j I notice a very steep slope for VPD, suggesting that biomass gains were only possible in regions that did not become drier and hotter. If that is the correct interpretation, it could perhaps be highlighted in the discussion? [Response]: Thank you for your valuable suggestion. We have now highlighted the role of VPD in this region in the discussion section, as suggested. \"…, notably the role of VPD in this region27, where increasing VPD was associated with decreasing biomass.\" (Page 15, lines 384-385) - Logging is mentioned and implicitly incorporated in the analysis as ‘non-fire forest loss’. However, the discussion misses a section on the role of logging or timber harvest on the observed decrease in yearly net L-VOD biomass change. What would be the relevance of timber harvest in light of your results? Moreover, and this relates closely to comment #1, how to avoid that readers think that all losses observed in the data are reductions of the atmospheric sinks? A concrete suggestion to answer this is to assess the reported harvest quantities by countries within the study area, for instance to the FAO or UNFCCC. Do these reported quantities correspond to the trend observed in the L-VOD data or to the observed ‘non-fire related forest loss’ data by Hansen et al. (2013)? [Response]: Thank you for this important comment. We agree that logging (timber harvest) may contribute to the observed biomass carbon decrease captured by L-VOD. While our analysis implicitly includes timber harvest under 'non-fire forest loss', we acknowledge the need to better contextualize its role. Following your suggestion, we compiled timber harvest volume data for countries in the study area based on FAO reports, including both industrial roundwood and fuelwood. Although timber harvest volumes fluctuated between 2016 and 2022, the overall increase was limited and did not exhibit a clear trend, whereas L-VOD estimates showed a decreasing trend in live biomass carbon over the same period. Moreover, the average annual carbon extracted via timber harvest was estimated at 0.27 PgC for industrial roundwood and 0.04 PgC for fuelwood, accounting for less than half of the L-VOD-derived biomass carbon losses attributed to non-fire forest loss (−0.76−0.68 PgC yr-1). This suggests that logging may not be the dominant factor of non-fire forest loss, but this loss is caused by other types of disturbances, consistent with studies in North American forests (Ma et al., 2012; Wang et al.2021; Kurz et al., 2008). In response to your comment, we have added: \"The average annual carbon removed through timber harvest reported by FAO is less than half of the losses from non-fire stand-replacing disturbances (Fig. S8), suggesting that harvesting may not be the main factor of these losses but rather other types of disturbances18,27, warranting further investigation.\" (Page 12, lines 289-292) Since there is currently a lack of harvest boundary data to explicitly isolate biomass loss caused by timber harvesting, we also added: \"Future studies could apply the approach of this study to isolate live biomass carbon loss from timber harvest, provided spatial data defining harvest disturbances or reliable proxies (e.g., forest management boundaries) become available59. \" (lines 371-373) In addition, to avoid the potential misunderstanding that all observed biomass losses correspond to reductions in atmospheric carbon sinks, we have added the following clarification: \"Moreover, live biomass carbon losses from non-fire disturbances may not be immediately emitted into the atmosphere in the year they occur but could be delayed (e.g., through the production of harvested wood products) or transferred to other carbon pools (e.g., from live biomass to dead wood and litter)35. In addition, forest fires can have low combustion completeness13,50-51, e.g., in Siberian larch, Siberian fir, or temperate broadleaved forests. Therefore, future studies should further quantify the fate of carbon to better contextualize net carbon stock changes within the framework of the global carbon budget and terrestrial carbon emissions and removals4,11. \"(Page 15, lines 373-379) - The RF-SHAP analysis nicely outlines the effect of various explanatory variables on net biomass changes, but the total variance explained by all explanatory variables combined is only ~50-60%. Could the authors briefly speculate about processes that this study did not consider and that could be underlying the unexplained variance? Reading the authors’ thoughts on this may be relevant to direct further studies on this subject. I would like to read at least a hypothesis on which processes remain unexplained with the chosen set of explanatory variables. [Response]: Thank you for your insightful feedback. We agree that the unexplained variance could be related to additional climatic and biotic factors. For example, stand-replacing disturbances caused by insect outbreaks and other non-fire, non-harvest factors are found to be the dominant drivers of biomass carbon loss in North American forests (Wang et al., 2021; Kurz et al., 2008), whereas in the tropics, fire radiative power and canopy structural indices are considered key drivers of biomass carbon gain (Yang et al., 2022; Feng et al., 2024). Additionally, processes operating at finer spatial scales might not be well captured by the current dataset. In response, we have incorporated the following perspective into our discussion: \"However, it is noteworthy that although multiple variables were included to explain spatial variations in live biomass carbon changes, 49% and 39% of the variance remain unexplained in temperate and boreal biomes, respectively. This suggests that additional climatic and biotic factors, such as canopy structure57, fire radiative power58, and insect disturbances27, merit exploration in future research, together with enhanced L-VOD spatial resolution to better capture small-scale drivers.\" (Page 15, lines 360-365) - In your rebuttal you mention that “This discrepancy is partly due to the fact that DGVMs do not simulate accurately actual wildfire occurrences and miss extreme fires.” Can you substantiate somewhere that extreme fires are indeed the ones that lead to largest carbon loss, perhaps with a reference? For the modeling community this would signify a stronger focus on the tails of the fire distribution instead of on the mean/median. If you agree, consider making that recommendation. [Response]: We agree with your comments. To address this, we have added a relevant recommendation and supporting references in the revised text below. \"our results showed that the DGVMs were not able to reproduce the effects of wildfires on the observed live biomass carbon evolution in boreal regions, possibly due to insufficient consideration of fire disturbance regimes in the models, particularly inaccurate representation of extreme fire occurrences60-61\"(Page 15, lines 380-383) - I appreciate the extended Table 1 and comparison to other estimates in the new text starting on line 324. When I suggested this I had in mind a comparison also to individual extreme events. On line 347 you mention the European summer drought and say your 0.25-0.28 Pg/yr (was the sign on the top number correct?) is “consistent” with refs 22 and 27. Can you specify more closely what this means? Again, the Δbiomass is not the same as the atmospheric sink anomaly, so what is it consistent with? Estimates of reduced NPP plus salvage logging? Or is “consistent” simply saying that your product also sees severe drought impacts qualitatively? [Response]: Thank you for your valuable suggestions and guidance. We fully agree with and understand your point. Unfortunately, no other up-to-date live biomass datasets are currently available to support a direct comparison, particularly for extreme events such as the 2022 European summer drought. Therefore, the use of “consistent” here refers to qualitative consistency, indicating that our L-VOD-derived biomass losses reflect the broader impacts of the event. To ensure clarity, we have revised the wording accordingly. We have also double-checked the reported values and did not find any errors. \" Of this reduction, 37% (−0.25−0.25 PgC) occurred in Europe, −0.28 which is qualitatively consistent with the widespread impacts of drought in the region observed during the summer of 2022, including reduced net carbon uptake by forests21,26.\"(lines 352-354)"}]

Decision letter