Wildfires offset the increasing but spatially heterogeneous Arctic–boreal CO₂ uptake

Estimating the terrestrial net ecosystem CO₂ exchange (NEE) of the Arctic–Boreal Zone (ABZ) poses a major challenge^1,2,3,4 due to complex biogeophysical dynamics^4,5 and a limited network of field measurements^6,7. As a result, models show a wide range of CO₂ budgets, from substantial net atmospheric sinks (−1,800 Tg C yr⁻¹) to moderate atmospheric sources (600 Tg C yr⁻¹)^1,4,8,9,10, a concerning discrepancy as northern permafrost soils hold nearly half of the global soil organic carbon stocks¹¹. The release of this soil carbon to the atmosphere as CO₂ could considerably exacerbate climate change¹². There is thus an urgent need to improve CO₂ budget estimates across the ABZ.

The rapid climate change of the ABZ makes this discrepancy even more critical¹³. Increasing air and soil temperatures in both summer and non-summer seasons are causing changes in the CO₂ budget that remain poorly understood⁹. Furthermore, it is not known how the widespread but spatially heterogeneous increase in vegetation productivity and greening¹⁴ impacts the annual CO₂ balance, although links to enhanced CO₂ sinks during the spring–summer period have been found¹⁵. Some of the enhanced uptake might be offset by CO₂ losses associated with increasing non-summer season respiration, vegetation dieback (‘browning’) and the escalating frequency and intensity of disturbances such as abrupt permafrost thaw (for example, thermokarst), drought and fires, further complicating the understanding of ABZ carbon dynamics and climate feedbacks^16,17,18.

Current evidence on recent ABZ CO₂ budget trends and their main drivers is limited to few in situ data-driven synthesis and modelling studies without a regional perspective on where and why CO₂ budgets are changing^1,8,9,10. These studies have focused primarily on ecosystem CO₂ fluxes (that is, not incorporating fire emissions), coarse annual or seasonal CO₂ fluxes (that is, overlooking the intra-annual dynamics), and spatial patterns in CO₂ fluxes with data from only one to two decades. Most importantly, earlier studies have not extended into the 2020s, a period when warming has further accelerated and fires have burned large areas¹⁹. We thus lack a comprehensive understanding of the regional and seasonal patterns in recent ecosystem CO₂ fluxes (including fire emissions), their multidecadal trends and the links to changing environmental conditions across the ABZ.

Here we address this knowledge gap using an ABZ CO₂ flux dataset that includes monthly terrestrial photosynthesis (gross primary productivity (GPP)), ecosystem respiration (R_eco) and NEE data from 200 terrestrial eddy covariance and flux chamber sites (4,897 site-months). This dataset is at least four times larger than ones used in earlier upscaling efforts and covers a longer period, with data extending to 2020. The same dataset was previously used to analyse in situ CO₂ flux trends in permafrost versus non-permafrost regions, with the conclusion that the annual net uptake is increasing in the non-permafrost region but not in the permafrost region²⁰. Here we extend that study from the site level to the full ABZ region by combining flux observations with meteorological, remote sensing and soil data, together with random forest models, to estimate CO₂ budgets across the ABZ. We do this upscaling over two periods, 2001–2020 (1-km resolution) and 1990–2016 (8-km resolution); the results in the main text are based on the 1-km models unless stated otherwise. We then assess regional and seasonal patterns and trends in ABZ ecosystem CO₂ fluxes and their environmental drivers. We also integrate annual fire emissions from 2002 to 2020²¹ to provide near-complete terrestrial CO₂ budget estimates (referred to as NEE + fire).

Using machine learning models that had a high predictive performance (up to two times higher cross-validated R² than earlier efforts^8,9), we found that from 2001–2020 circumpolar tundra was on average CO₂ neutral without accounting for fire emissions (in situ NEE, −4 ± 44 g C m⁻² yr⁻¹; upscaled NEE, 7 ± 3 g C m⁻² yr⁻¹; upscaled budget, 45 ± 53 Tg C yr⁻¹ (mean ± standard deviation); Table 1; a negative sign indicates net uptake, and a positive sign indicates net emissions). In contrast, the boreal was a strong sink (in situ NEE, −42 ± 82 g C m⁻² yr⁻¹; upscaled NEE, −43 ± 7 g C m⁻² yr⁻¹; upscaled budget, −593 ± 101 Tg C yr⁻¹). Including fire emissions (on average 237 Tg C yr⁻¹ (ref. ²¹)—that is, 2% of R_eco and 43% of the net CO₂ uptake budgets in the ABZ) changed the budget to −383 ± 101 Tg C yr⁻¹ in the boreal and to 64 ± 53 Tg C yr⁻¹ in the tundra. With fire emissions included, the permafrost region became CO₂ neutral (NEE, −249 ± 123 Tg C yr⁻¹, NEE + fire, −24 ± 123 Tg C yr⁻¹).

Although the entire ABZ ___domain was a terrestrial CO₂ sink across all years during 2001–2020 with an average NEE of −548 ± 140 Tg C yr⁻¹, our upscaling of NEE revealed a large areal fraction of annual ecosystem CO₂ sources across the ___domain (34% of the total region; Fig. 1). These upscaled CO₂ source regions were located particularly in the northern parts of the ABZ characterized by low temperatures and normalized difference vegetation index (NDVI) values (Extended Data Fig. 1). For the permafrost ___domain, the fraction of annual CO₂ sources was even higher (41% of the region). This large fraction is also seen in our in situ CO₂ flux database, with 29% of sites being CO₂ sources (NEE between 0 and 142 g C m⁻² yr⁻¹). These CO₂ source sites were mostly in Alaska (44%) but also in northern Europe (25%), Canada (19%) and Siberia (13%). One key factor driving CO₂ sources is the long and persistent non-summer-season (September–May) emissions in the tundra that, on average, exceed the short summer (June–August) net CO₂ uptake (Table 1). In the boreal, longer summers with strong uptake still dominate over non-summer emissions.

a,b, Maps showing the mean annual terrestrial NEE (a) and its trends (b) based on site-level data, our upscaling, the atmospheric inversion ensemble and the CMIP6 process model ensemble. The in situ trends in b are based on sites that have more than seven years of data. Supplementary Fig. 2 shows the uncertainty in upscaled NEE and the significance of the trends. While the average upscaled NEE values go up to 116 g C m⁻² yr⁻¹, most of the values are below 60 g C m⁻² yr⁻¹.

Full size image

We observed moderate correlation of our upscaled NEE results with an ensemble (that is, mean) of atmospheric inversions²² across space (Pearson’s correlation, 0.5; P < 0.001), but the correlation between the temporal trends was weaker (Pearson’s correlation, 0.2; P < 0.001) (Fig. 1). However, the ensemble net uptake budgets from the inversions, as well as from a global machine-learning-based upscaling product (FLUXCOM-X-BASE^23,24), were 1.5 to 3 times larger than our upscaled budgets (Supplementary Section 5). Moreover, the global Coupled Model Intercomparison Project Phase 6 (CMIP6) process model ensemble²⁵ had barely any annual CO₂ sources across the ABZ, indicating that the process models may not accurately simulate CO₂ source situations (Fig. 1), especially given the prevalence of site-level sources. The cross-validated predictive performances of our random forest models for GPP, R_eco and NEE showed high correlations between observed and predicted fluxes (R² varied from 0.57 to 0.73, and the root mean square error (RMSE) varied from 19.4 to 37.3 g C per m² per month; Extended Data Fig. 2), but upscaling uncertainties remain. For example, areas with the most extensive strong sink or source estimates rarely had in situ data and were thus largely extrapolated (for example, sources in central Siberia or sinks in southern Siberia; Supplementary Fig. 1). These areas also had the highest uncertainties in our analysis (approximately twice as large as in the more densely measured areas; Supplementary Fig. 2).

The ABZ has been an increasing terrestrial CO₂ sink on the basis of NEE alone from 2001 to 2020 (temporal trend, −14 Tg C yr⁻¹; P < 0.001) (Fig. 2). However, the increasing sink strength was no longer statistically significant when fire emissions were added to NEE (average NEE + fire budget trend, −9 Tg C yr⁻¹ over 2002–2020). In the permafrost region, the NEE + fire trend was only −3.3 Tg C yr⁻¹. Nevertheless, on the basis of our NEE upscaling, 23% of the region showed a statistically significant (P < 0.05) increase in net CO₂ uptake from 2001 to 2020 (Fig. 1), with increasing net sink pixels occurring across all key regions. Most of the increasing net sink activity was driven by an increase in GPP, especially in Siberia (Fig. 2). Some of the trends were also related to a declining R_eco, probably associated with disturbed ecosystems (for example, forest fires and harvesting) with high R_eco during the first post-disturbance years now recovering²⁶. However, evidence for the increasing overall net uptake trend from the in situ data is limited due to the low number of long-term sites (more than seven years of year-round measurements; nine sites), of which only one site showed a statistically significant trend (increasing uptake at a boreal forest site in Finland). Some of the relationships in our model are probably thus influenced by spatial differences across the sites rather than temporal and truly causal patterns, creating some uncertainty in upscaled trends²⁷. However, the model reproduces temporal changes at individual sites well (Supplementary Figs. 3 and 4), and our upscaled trends are similar to a recent in situ time-series analysis²⁰ and somewhat similar to those estimated from the inversion ensemble (Fig. 1), providing confidence in our trend results.

a,b, Terrestrial CO₂ budgets for 1-km (blue; 2001–2020) and 8-km (grey; 1990–2016) NEE as well as 1-km NEE + fire emissions (red; 2002–2020) across the ABZ (a) and the permafrost region (b). c, An overlay analysis of NEE, GPP and R_eco trend maps identifying how trends in GPP and R_eco relate to trends in NEE over 2001–2020 (includes significant and non-significant trends). d, Pixels burned during 2002–2020. The central values (that is, annual budgets; solid lines) in a and b are derived from the outputs of the final model using the complete training dataset. The standard deviations (shaded areas) are calculated from the outputs of 20 different models, each trained on a unique bootstrapped sample of the original training data. The magnitude of each trend was computed using the Theil–Sen approach, and the P value determined using the Mann–Kendall test.

Full size image

Parts of the ABZ also show increasing annual net CO₂ emissions over time (Fig. 1). Such trends have been observed at six long-term sites (2 to 17 g C m⁻² yr⁻¹, P > 0.05) and in 2% of the upscaled region (P < 0.05) from 2001 to 2020. Most of the increasing net emission trends were driven by an increase in R_eco instead of a decline in GPP (Fig. 2). Regions experiencing increased net CO₂ emissions in upscaling were found especially in (1) northern Europe and Canada (dominated by evergreen needleleaf forests with mild and moderately wet climates), (2) parts of central Alaska and northern Siberia (sparse boreal ecosystems and graminoid tundra with permafrost and high soil carbon stocks), and (3) Hudson Bay and Siberian lowlands (wetlands with some permafrost and high soil organic carbon stocks). Some sites in central Alaska have increasing net emissions of CO₂ due to permafrost thaw^16,28, but it is unclear whether similar changes are occurring in other regions with increasing net CO₂ emissions.

We calculated an overall 25% increase in the seasonal amplitude of CO₂ fluxes from the upscaled NEE time series from 2001 to 2020 across the ABZ, on par with earlier atmospheric and modelling studies^29,30. Both increasing summer uptake and non-summer season emissions—the key dynamics driving increasing annual sinks and sources—were evident in the tundra and boreal biomes (Fig. 3). However, over the 2001–2020 period, the increasing uptake (GPP) during summer months dominated the increasing net emissions (R_eco) during non-summer months across most of the ___domain. On average across both biomes, net uptake increased the most during July (an average upscaled increase of −5 g C per m² per month in the boreal and −3 g C per m² per month in the tundra in 2011–2020 compared with 2001–2010), and increasing net emissions occurred throughout the entire non-summer season, with no clear peaks (0.1–0.9 g C per m² per month). Although increases in the early growing season (May–June) uptake were evident, trends in the late growing season (September) were absent or minimal (Fig. 3).

a–c, Average upscaled monthly NEE (a), GPP (b) and R_eco (c) in boreal and tundra biomes during the past two decades. Negative NEE values represent net uptake, and positive values indicate net release. The central mean values in the figures are derived from the outputs of the final model using the complete training dataset. The standard deviations (error bars) are calculated from the outputs of 20 different models, each trained on a unique bootstrapped sample of the original training data. Error bars are shown only for the 2011–2020 period but are similar for the 2001–2010 period. Note that NEE was 1.4 g C m⁻² month⁻¹ lower in September in 2011–2020 than in 2001–2010 in the boreal biome, but this is not shown in the figure.

Full size image

There are several environmental conditions driving CO₂ budgets across the ABZ. Our variable importance analysis showed that CO₂ fluxes, and thus the overall increasing sink strength, are explained by dynamic variables of air or surface temperatures, solar radiation and NDVI, and partially also by soil temperature, snow cover and the vapour pressure deficit (Extended Data Fig. 3). Other less important dynamic variables were vegetation cover and atmospheric CO₂ concentration. Volumetric soil water content was not important in our models, probably due to the large uncertainties and coarse spatial resolution in the gridded product, although in situ studies have shown drier soils to be linked to larger net CO₂ emissions and wetter soils to enhanced plant growth due to the lack of water limitation³¹. Static variables (primarily vegetation type, soil carbon stock and soil pH) were also important in explaining spatial differences.

The most important dynamic variables had a positive overall effect on net uptake, GPP and R_eco (Extended Data Fig. 3); however, these relationships are more nuanced in reality. In fact, the recent permafrost in situ trend analysis of CO₂ fluxes using the same database suggests that the CO₂ flux response to warmer temperatures ranges from positive to negative, depending on the availability of water and nutrients at the site²⁰. Consequently, strong warming or greening trends did not always translate into increasing net CO₂ sinks in our upscaling (Supplementary Fig. 5). For example, while 49% of the region experienced greening (June–August average NDVI; on the basis of MODIS NDVI, P < 0.05), only 12% of those greening pixels showed an annual increasing net CO₂ uptake trend, and 29% showed an increasing June–August net uptake.

Our upscaling showed clear continental patterns in NEE budgets and trends (Fig. 4), with the boreal biome primarily driving the budget and trend differences between the continents^1,32. The increasing net uptake trend was more pronounced in Eurasia (−11 Tg C yr⁻¹, P < 0.001) than in North America (−3 Tg C yr⁻¹, P < 0.05), which corresponds with the smaller area and weaker warming, declining snow cover and greening trends in North America (Supplementary Figs. 6–10). We found statistically significant declining summer soil moisture trends in the Siberian boreal (Supplementary Fig. 9), but this did not translate into stronger net emissions. When fire emissions were added, continental differences were less pronounced due to the much larger and more rapidly increasing CO₂ emissions from Siberian fires (on average 160 compared with 76 Tg C yr⁻¹ in North America; Extended Data Fig. 4). Fire emissions even reversed some NEE trends: the strong increasing sink in Siberia became a source when fire emissions were included (trend, +0.7 Tg C yr⁻¹; P > 0.05; Fig. 4). However, Siberian ecosystems have the largest uncertainty for both the upscaled fluxes and inversion-based estimates due to the lack of in situ observations, making it challenging to accurately determine the magnitude of continental differences (Fig. 4 and Supplementary Fig. 11).

a,b, Terrestrial CO₂ budgets for NEE and NEE + fire in key regions across the boreal (a) and tundra (b). Terrestrial CO₂ budgets are shown for 1-km (blue; 2001–2020) and 8-km (grey; 1990–2016) NEE as well as 1-km NEE + fire emissions (red; 2002–2020). The inset in the Alaskan boreal plot in a shows the time series with a narrower y axis than that in the main figure to better depict interannual variability. The central lines in the figures are derived from the outputs of the final model using the complete training dataset. The standard deviations (shaded areas) are calculated from the outputs of 20 different models, each trained on a unique bootstrapped sample of the original training data. The magnitude of the trend was computed using the Theil–Sen approach, and the P value determined using the Mann–Kendall test.

Full size image

Alaska is an important contributor to the weaker North American CO₂ sink. On the basis of our analysis, Alaska as a whole was consistently CO₂ neutral or a source over 2002–2020 (NEE + fire emissions), in both the boreal (budget, +5 Tg C yr⁻¹) and tundra (budget, +7 Tg C yr⁻¹). Alaska is therefore different from the other ABZ regions, where boreal regions remain on average CO₂ sinks. Alaska has a relatively high density of observations, making this result more certain than those for other regions. Potential reasons for the Alaskan CO₂ source include Alaska having the most rapidly warming autumns and declining autumn snow cover, which also have high interannual variability (Supplementary Figs. 8 and 10). Furthermore, field measurements suggest that many of the observed changes in Alaskan ecosystems can be attributed to permafrost thaw^16,28—a phenomenon that has accelerated notably in response to Alaska’s pronounced warming trend since the 1950s³³. However, we were unable to incorporate permafrost thaw into our models, as high-resolution geospatial data from 1990 to 2020 were not available. The question of whether analogous trends will manifest in other regions across the northern permafrost region remains an important research priority.

Our results show that the ABZ was on average an increasing terrestrial CO₂ sink (GPP is increasing more than R_eco + fire), indicating that the region still creates an important negative feedback to global warming. However, our study also suggests some positive feedbacks to climate change that have been more regional and of shorter duration in recent decades. We show that the presence of annual sources was large, as indicated by several site-level and regional studies^34,35,36, and even larger with fire emissions included³⁷. There were also extreme years when fire emissions exceeded annual net CO₂ uptake (for example, 2003 in the Siberian boreal, 2012 in the Canadian boreal and several years in the permafrost region; Fig. 2). We also observed increasing shoulder-season net emission trends, particularly in Alaska³⁸ (Extended Data Fig. 5). Moreover, while the summer net uptake increase still dominates over non-summer CO₂ emissions, net CO₂ uptake is increasing only in the early and peak growing season (May–August in the boreal and June–August in the tundra) and not in the late growing season (September), because GPP does not increase later in the season due to plant physiological limitations, and drier and warmer conditions cause enhanced R_eco instead^{39,40,41,42,43}. A better understanding of how soil moisture and hydrology have been and will be changing, and the impact of these changes on CO₂ fluxes, is critical for more accurate ABZ CO₂ budgets.

Our findings reveal a noteworthy shift in carbon dynamics in the tundra. While the tundra region has been a carbon sink for millennia⁴⁴, our results suggest that many tundra regions may now have started to function as CO₂ sources. This transition from an ecosystem CO₂ sink to a CO₂ source may have begun prior to 1990⁴⁵, but the precise timing of this transformation remains uncertain. The main drivers of this pattern may be related to warming-induced permafrost thaw, the drying of soils or vegetation shifts^46,47,48 but remain unresolved. Tundra regions are also progressing towards conditions where average annual soil (0–7 cm) temperatures are above freezing, resulting in more soil organic material being susceptible to decomposition (Supplementary Fig. 6). Overall, the primary reasons behind the annual CO₂ emissions from tundra ecosystems are the limited duration of the high net CO₂ uptake period and the substantial non-summer-season net emissions. However, we observed lower in situ and upscaled October–April NEE fluxes and budgets than did Natali et al.⁹ throughout the entire period (Supplementary Fig. 12).

Our results demonstrate the need to further study Siberian CO₂ flux trends. Our upscaling indicated that some of the strongest net sources and sinks and the strongest increasing sink trends occur in the Siberian boreal. Increasing sink trends in the Siberian tundra were also the strongest across tundra regions. The Siberian sink trend might be explained by strong greening trends⁴⁹, earlier growing season starts and increasing carbon uptake due to declining spring snow cover (Supplementary Fig. 8), increases in tree growth and distribution^50,51, rapid recovery of ecosystems after fire⁵², and high cover of larch forests that can rapidly take up CO₂ (Supplementary Table 1)^7,53. However, the large inversion model spread, the sparse measurement network and our upscaling uncertainties indicate that it remains challenging to determine the magnitude of the Siberian CO₂ balance². This is a major problem given that Siberia stores more than half of the permafrost region’s carbon stocks and is now warming more rapidly than other ABZ regions.

In summary, our study reveals distinct spatial and temporal patterns in CO₂ budgets across the ABZ and underscores the importance of three decades’ worth of data. Relatively robust spatial patterns can be seen, such as the Alaskan CO₂ sources and southern Eurasian boreal sinks, while the temporal trends remain more uncertain. While CO₂ fluxes can be relatively well modelled using machine learning and existing gridded datasets, gaps persist, such as the incomplete characterization of fire, thermokarst and harvesting disturbances and their links to ecosystem CO₂ fluxes (Extended Data Fig. 6), the lack of accurate predictors describing soil moisture, and the need to quantify landscape heterogeneity and carbon dynamics at even higher spatial (metres) and temporal resolutions (days). Sustaining long-term sites is crucial to accurately tracking trends in ABZ CO₂ balance, while establishing new year-round sites in data-poor areas such as Siberia and the Canadian Arctic is vital to filling knowledge gaps and enhancing our understanding of carbon dynamics.

In situ data overview

We used a recently compiled dataset of in situ Arctic–boreal terrestrial ecosystem CO₂ flux measurements (ABCflux, led by Virkkala et al.⁷) within the ABZ (Supplementary Sections 1 and 2 and Supplementary Table 2). The synthesized data were cumulative flux densities of gross NEE, GPP and R_eco aggregated at the monthly timescale (3,968 to 4,897 site-months depending on the flux). In addition to eddy covariance data, we included fluxes estimated with chamber methods to increase data coverage, especially during the growing season. The majority of the data for the 1-km models were based on eddy covariance (55% of sites and 88% of months). Each site had from 1 to 213 months of measurements in our database; the average number of months per site was 25. Long-term sites with more than seven years of year-round data included boreal forest sites (FI-Hyy, FI-Sod, CA-Oas, CA-Obs, CA-Gro, US-Uaf and SE-Deg), a wetland site (FI-Kaa) and a tundra site (US-EML).

In total, 14% of the sites in ABCflux had experienced some level of natural or anthropogenic disturbances. This proportion might be less than the overall proportion of disturbances across the entire ABZ. For example, 7% of the ABZ was burned during the 2002–2020 period²¹, post-fire succession can take many decades, and the areas experiencing thermokarst and harvesting are also extensive. The flux site distribution might thus be biased towards non-disturbed or only moderately disturbed sites^54,55, leading to potential underestimations of the effects of disturbance on CO₂ emissions. We had 21 sites that reported fire disturbance. Only four of those were longer-term sites (operating for more than three years) with recent fire history (less than ten years since burn). All four of these sites in young recovering ecosystems were measured year-round and originally had an evergreen (black spruce) forest cover, which underwent a shift to a more deciduous shrub- and tree-dominated cover after a stand-replacing fire. These include (1) CA-NS7, with four years of data starting four years after the fire; (2) CA-SF3, with six years of data starting three years after the fire; (3) US-Rpf, with a six-year time series starting four years after the fire; and (4) US-CR-Fire, with a four-year time series starting the next year after the fire. Of all the sites, three reported thermokarst, but gradual permafrost thaw was present in many more sites. At least five forest sites had been harvested.

The ABCflux dataset is more comprehensive than the ones used in earlier upscaling studies, as it represents monthly fluxes from the entire year if available, while Virkkala et al. focused on coarse seasonal or annual fluxes⁸, Natali et al. on monthly winter fluxes⁹ and Mu et al. on a more limited temporal period (2002–2017)⁵⁶. Furthermore, we included more data from recent years (805 monthly observations from 2015–2020 compared with 32 fluxes in Virkkala et al.⁸ and 95 fluxes in Natali et al.⁹), and the sample size in our models was 4 to 25 times larger than that in the earlier upscaling efforts.

Geospatial data

We used data from geospatial products as predictor variables to upscale fluxes. Our models had the following predictors: month, incident solar radiation, vapour pressure deficit, atmospheric CO₂ concentration, vegetation type, snow cover (the fraction covered by snow), soil temperature (0–7 cm), soil moisture (0–7 cm), NDVI (MODIS- or AVHRR-based), land surface temperature (or air temperature; MODIS- or ERA5 Land-based), compound topographic index (that is, topographic wetness index), continuous vegetation fields describing per cent non-tree vegetation and non-vegetated fraction and per cent tree cover (MODIS- or AVHRR-based), soil pH (0–5 cm), soil organic carbon stock in 0–2 metres, and permafrost probability (Supplementary Table 3 and Supplementary Section 3). In our analysis, NDVI was the primary predictor describing disturbances, with declines in NDVI being related to disturbances⁸ (Supplementary Figs. 13 and 15). The data were in daily, weekly, monthly, annual and static formats (that is, no temporal changes, such as in the compound topographic index). If data were of higher temporal resolution than monthly, they were aggregated to monthly time steps. Gaps in MODIS and AVHRR NDVI time series were filled to produce a continuous time series. The data were reprojected to North Pole Lambert Azimuthal Equal Area Projection at 1-km and 8-km spatial resolution and extracted at the flux sites. See Supplementary Section 3 for further descriptions and data sources.

We used a 500-m fire dataset based on the Global Fire Emissions Database (GFED) series of products²¹ to calculate fire emissions and burned fractions. The product uses the MODIS-based MCD64A1 product⁵⁷ for burned area at 500-m resolution but also adds a partial burned fraction to pixels that experienced fires according to MODIS active fire acquisitions⁵⁸ and forest loss according to Hansen et al.⁵⁹. The model calculates fire emissions at 500 m on the basis of the GFED4 model but incorporates an updated field database of fuel load and consumption for model calibration, including over 800 sites from boreal forests. The higher resolution of the 500-m model than that of earlier coarser models improved the detection of small-scale fires and understanding of landscape heterogeneity, and reduced the scale mismatch in comparing field measurements to model grid cells.

Machine learning modelling

We used random forest models to upscale GPP, R_eco and NEE to the ABZ from 1990 to 2020, the period with in situ flux measurements. Two sets of predictive models were developed: (1) models using primarily predictors with a spatial resolution ≤1 km from 2001 to 2020 (that is, the MODIS era) at 1-km spatial resolution (referred to as 1-km models) and (2) models using coarser-resolution predictors from the AVHRR GIMMS era (1990–2016) from 1990 to 2016 at 8-km spatial resolution (referred to as 8-km models) (Supplementary Table 3). Each model included all available monthly fluxes from the entire year (that is, no separate models for individual months or seasons were developed), as this approach resulted in the best predictive performance. All models included 17 predictors, but the sample sizes varied depending on the amount of data available for each flux and time period; NEE models had the highest amount of model training data and more sites than GPP and R_eco models (Supplementary Table 4). These differences in data distribution resulted in slight mismatches in the upscaled NEE versus GPP − R_eco estimates.

Model parameter tuning was performed separately for each response variable on the basis of leave-one-site-out cross validation to achieve the lowest RMSE values. The models were run using the caret package in R version 4.2 (ref. ⁶⁰), and the analysis code is available in Virkkala et al.⁶¹. The tuned parameter (the number of variables to randomly sample as candidates at each split) varied from 2 to 17 in the final models. To interpret the models, we used partial dependence plots (that is, response graphs) estimated with the pdp package⁶² and calculated importance scores for the predictors from each of the models using the vip package⁶³. Variable importance scores were estimated by randomly permuting the values of the predictor in the training data and exploring how this influenced model performance on the basis of RMSE values, with the idea that random permutation would decrease model performance⁶⁴. We used 100 simulations to calculate 100 importance scores. We assessed the area of extrapolation and clustering of environmental conditions to further understand the patterns in upscaled fluxes (Supplementary Section 4).

We assessed the predictive performance of the final models using (1) R², (2) the RMSE, (3) the mean absolute error and (4) the mean bias error between predicted and observed values using the cross-validation data. Larger RMSE and mean absolute error values indicate larger errors, and positive mean bias error values indicate overestimation. The predictive performance of our models was good or high, with R² ranging from 0.57 to 0.73 and RMSE from 19.4 to 37.3 g C per m² per month. However, our performance metrics also indicate that strong sinks and sources as well as high GPP and R_eco were underestimated—a common issue in any kind of modelling⁶⁵. As reflected by the small and positive mean bias error values, the NEE models had a small tendency to overestimate fluxes (that is, predict too-small net sinks or too-high net emissions). Furthermore, the distribution of NEE residuals was slightly skewed towards negative residuals (that is, NEE was overestimated; Supplementary Fig. 14). Overestimation was particularly evident with the model struggling with strong sink observations (observations of around −180 to −80 g C per m² per month were predicted to be −80 to −20 g C per m² per month in cross validation; Extended Data Fig. 1). However, the model also underestimated strong net CO₂ sources (with a deviation of up to 80 g C per m² per month). Potential reasons for these biases are probably related to (1) our model not being able to identify landscape heterogeneity, with nearby sites showing large differences in CO₂ fluxes (for example, a forest and a wetland site), and (2) our model not capturing interannual variability at individual sites, both of which are probably attributed to the coarse, uncertain and missing predictors characterizing such conditions (for example, soil moisture and disturbances) (Supplementary Figs. 1–3). However, compared with earlier ABZ upscaling efforts, our cross-validated performance metrics indicate better performance. For example, the R² of our models ranged from 0.57 to 0.73, whereas Natali et al.⁹ had an R² of 0.49 for winter NEE, and Virkkala et al.⁷ had an R² of 0.07 for annual NEE and an R² of 0.5 for annual GPP and annual R_eco; note, though, that the cross validation in Natali et al.⁹ was not based on a leave-one-site-out approach.

We evaluated the uncertainty of predictions by creating 20 bootstrapped model training datasets (with replacement and with the same sample size as in the original model training data) and using those to develop 20 individual models and predictions. For these bootstrapped datasets and models, we did not include the categorical month and land cover datasets as predictors due to bootstrapping resulting in situations where a factor level was entirely missing from the model training data (for example, for a barren class that had little data), which prevented us from predicting fluxes across the entire ___domain. Using the 20 predictions, we calculated the standard deviation to represent prediction uncertainty. Similar to the predictive performance metrics, the uncertainty analysis also points towards the highest uncertainties in areas with strong sinks, such as in northern Europe and southwestern Russia. However, when the uncertainty estimates were presented relative to the average flux, uncertainties were the highest in tundra regions and parts of northern boreal Canada, which generally have low in situ flux data coverage.

Model performance in burned ecosystems

In addition to direct fire emissions from combustion (that is, burning) derived with GFED, fires have a profound impact on carbon budgets by modulating post-fire ecosystem CO₂ fluxes^26,66. Across all the burned sites in our model training data, the in situ flux data and remotely sensed NDVI time series showed a clear pattern of July NDVI values, GPP and net carbon uptake steadily increasing after the fire (Supplementary Fig. 13). This post-fire recovery signal was captured by our upscaling, as our upscaled GPP and net uptake dropped after a fire and then returned to higher levels (Supplementary Fig. 15). However, while our random forest model fitted the time series of the longer-term sites with recent fire history relatively well, the predictions based on cross validation (that is, model training data excluding each site) were variable (Supplementary Fig. 4), indicating that our model might have struggled in extrapolating post-fire ecosystem CO₂ fluxes in other areas. The model performance at sites experiencing recent fire or other disturbances was also lower than at sites without disturbance or disturbance information, as the model had a lower R² and a tendency to underestimate NEE values (that is, underestimate net CO₂ emissions or overestimate net CO₂ uptake) (Extended Data Fig. 6). The model bias based on cross validation was up to 75 g C per m² per month at the burned sites (Extended Data Fig. 6 and Supplementary Fig. 14). Such biases were high but not exceptional at other disturbed or undisturbed sites. When looking at the models trained with or without disturbed sites, there were no major differences in predictive performance estimates (Supplementary Table 5).

We further evaluated the NEE model bias as a function of years since fire (0 to 18 years) at burned sites. There was no evidence of a nonlinear relationship during this period (Supplementary Fig. 16); thus, we developed linear mixed-effects models for each month with model bias from cross validation as a response variable, years since fire based on site-level information as a predictor and site ID as a random effect (both slope and intercept)⁶⁷. The population-level slopes and intercepts were variable across months, with May, June and July being the only months with statistically significant slope and intercept parameters. During the May–July months, there was a positive bias (that is, net CO₂ sinks were overestimated) during the first 7–10 years after the burn and a negative bias (that is, net CO₂ sinks were underestimated) during the last years starting 8–11 years after the burn. The number of years with a positive bias varied depending on the month (ten for May, seven for June and nine for July), which were modelled separately.

Due to this relationship between fire history and model bias, we applied a sensitivity test to assess whether correcting for the bias would change the annual NEE budget in 2020. The bias correction was done for May–July only, and for 2020 alone to ensure the maximum years-since-fire time series across the ABZ (18 years; the GFED 500-m product begins in 2002). We summed the original upscaled monthly NEE value by the monthly bias predicted separately with the mixed model from 0 to 18 years since fire. The bias correction was applied to the pixels that had burned in 2002–2020 on the basis of the GFED product. To account for the burns not covering the entire pixel, the bias was scaled by the fraction of the pixel that was burned. The bias correction resulted in a minor increase in the annual net uptake budget across the ABZ for 2020, with the original budget changing from −809 Tg C yr⁻¹ to −822 Tg C yr⁻¹ with the bias correction. The correction thus increased net uptake, probably due to the largest fire years over 2002–2020 occurring in the early 2000s (Extended Data Fig. 4), during which the model was underestimating net CO₂ sinks.

Spatial upscaling of fluxes

We upscaled fluxes across the Arctic–boreal terrestrial area ≥49 °N (ref. ⁶⁸), which comprises 20.7 × 10⁶ km² of land (excluding glaciers and ice sheets; Fig. 1) with lake areas removed. The models were applied at a monthly time step from 2001 to 2020 for the 1-km models and from 1990 to 2016 for the 8-km models. In total, we produced 1,692 upscaled flux data layers. The 8-km upscaled layers were further multiplied by the terrestrial surface fraction within each 8-km pixel on the basis of the 1-km ESA CCI + CAVM land cover dataset to remove fluxes from water bodies. A comparison of the 1-km and 8-km average annual NEE maps for 2001–2016 showed that NEE was similar across the two pixel resolutions (Supplementary Fig. 17). We also compared upscaled NEE maps and budgets from two approaches: one based on modelling NEE directly and one based on deriving it indirectly from the upscaled GPP and R_eco maps. NEE from these two approaches yielded similar results, providing confidence in our results (Table 1 and Supplementary Fig. 17). Overall, our upscaling results revealed a latitudinal pattern of average CO₂ fluxes, with stronger sinks in the south and weaker sinks or sources in the north (Fig. 1). However, the correlation between latitude and average NEE was moderate (Pearson’s correlation for in situ NEE, 0.26; P = 0.053; for upscaled NEE, 0.55; P < 0.001), suggesting that the latitudinal climate and radiation gradients were not the sole drivers of spatial CO₂ flux patterns.

We analysed the upscaled flux data layers as well as fire emission and environmental predictor rasters for temporal trends using the terra package⁶⁹ and the non-parametric Mann–Kendall test from the zyp package^70,71 with pre-whitening (the Zhang method⁷²) to remove autocorrelation. We report the significance of Kendall’s correlation coefficient (the strength of the time series) and the Theil–Sen slope to describe trends over time. Finally, we calculated zonal statistics of average annual, seasonal and monthly fluxes and trends across key regions (Siberia was defined as all land east of the Ural mountains, including a small portion of Mongolia; the rest of Eurasia, including Greenland, was grouped within the European classes), biomes (tundra and boreal)⁶⁸, the permafrost region⁷³ and vegetation types⁷. Supplementary Table 1 shows the budgets for different vegetation types and regions.

Comparison to process models and atmospheric inversions

We compared our estimates with the CMIP6 process models²⁵, atmospheric inversions used in the Global Carbon Project’s Global Carbon Budget 2022^22,74 and a global upscaling product, FLUXCOM-X-BASE²⁴ (Supplementary Section 5). We included a subset of CMIP6 process models (13 in total) that had soil thermal processes at several depths to assure that they had some information about the freeze–thaw patterns in the permafrost region. We included inversions with data from the whole 2001–2020 period (that is, we included five inversions and excluded four). Fire CO₂ emissions²¹ were subtracted from the inversions. CMIP6 process model outputs were available only for the 2001–2014 period. Model outputs are presented here as ensembles (that is, means of the individual models). There is some heterogeneity between individual inversions and CMIP6 models within the ensembles, but overall the ensemble results can be considered robust^5,75 (Supplementary Table 6). The final ensemble outputs used here represent terrestrial NEE (GPP − R_eco) in a similar way across the models except for inversions that also include vertical CO₂ fluxes from water bodies. None of the models include lateral transport of carbon. For reference, lateral transport of carbon and vertical lake and river CO₂ emissions were recently summarized to be 93, 66 and 164 Tg C yr⁻¹, respectively, in the northern permafrost region (that is, greater than the NEE + fire budget calculated in this study for the same ___domain)⁷⁶.