<figure><img src="image_1.png" /></figure>

HAL
open science

# Wildfires offset the increasing but spatially heterogeneous Arctic-boreal CO2 uptake

Anna-Maria Virkkala, Brendan M Rogers, Jennifer D Watts, Kyle A Arndt, Stefano Potter, Isabel Wargowsky, Edward a G Schuur, Craig R See, Marguerite Mauritz, Julia Boike, et al.

► To cite this version:

Anna-Maria Virkkala, Brendan M Rogers, Jennifer D Watts, Kyle A Arndt, Stefano Potter, et al. Wildfires offset the increasing but spatially heterogeneous Arctic-boreal CO2 uptake. Nature Climate Change, 2025, 15 (2), pp.188 - 195. 10.1038/s41558-024-02234-5. hal-04950986

HAL Id: hal-04950986

https://hal.science/hal-04950986v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

Received: 24 January 2024

Accepted: 3 December 2024

Published online: 21 January 2025

<figure><img src="image_4.png" /></figure>

Check for updates

A list of authors and their affiliations appears at the end of the paper

The Arctic-Boreal Zone is rapidly warming, impacting its large soil carbon stocks. Here we use a new compilation of terrestrial ecosystem CO₂ fluxes, geospatial datasets and random forest models to show that although the Arctic-Boreal Zone was overall an increasing terrestrial CO₂ sink from 2001 to 2020 (mean ± standard deviation in net ecosystem exchange, -548 ± 140 Tg C yr⁻¹; trend, -14 Tg C yr⁻¹; *P* < 0.001), more than 30% of the region was a net CO₂ source. Tundra regions may have already started to function on average as CO₂ sources, demonstrating a shift in carbon dynamics. When fire emissions are factored in, the increasing Arctic-Boreal Zone sink is no longer statistically significant (budget, -319 ± 140 Tg C yr⁻¹; trend, -9 Tg C yr⁻¹), and the permafrost region becomes CO₂ neutral (budget, -24 ± 123 Tg C yr⁻¹; trend, -3 Tg C yr⁻¹), underscoring the importance of fire in this region.

Estimating the terrestrial net ecosystem CO₂ exchange (NEE) of the Arctic-Boreal Zone (ABZ) poses a major challenge¹⁻⁴ due to complex biogeophysical dynamics⁴,⁵ and a limited network of field measurements⁶,⁷. 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⁻¹)¹,⁴,⁸⁻¹⁰, 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¹⁶⁻¹⁸.

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¹,⁸⁻¹⁰. 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<sub>eco</sub>) 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

<figure><img src="image_5.png" /></figure>

e-mail: avirkkala@woodwellclimate.org
<table><caption>Table 1 | Average Arctic-boreal CO₂ fluxes</caption><thead><tr><th rowspan="2">Class<br>Flux and unit</th><th colspan="3">In situ average</th><th colspan="3">Upscaled per-area average</th><th colspan="3">Average regional budget</th><th rowspan="2">Proportion of summer (June–August) net uptake budget of non-summer (September–May) net emissions (%)</th><th rowspan="2">Average regional budget with fire<br>NEE + fire (TgCyr<sup>-1</sup>)</th><th rowspan="2">Area (×10<sup>6</sup> km<sup>2</sup>)</th></tr><tr><th>NEE (gCm<sup>-2</sup>yr<sup>-1</sup>)</th><th>GPP (gCm<sup>-2</sup>yr<sup>-1</sup>)</th><th>R<sub>eco</sub> (gCm<sup>-2</sup>yr<sup>-1</sup>)</th><th>NEE (gCm<sup>-2</sup>yr<sup>-1</sup>)</th><th>GPP (gCm<sup>-2</sup>yr<sup>-1</sup>)</th><th>R<sub>eco</sub> (gCm<sup>-2</sup>yr<sup>-1</sup>)</th><th>NEE (TgCyr<sup>-1</sup>)</th><th>GPP (TgCyr<sup>-1</sup>)</th><th>R<sub>eco</sub> (TgCyr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>ABZ</td><td>-32 (±76)</td><td>618 (±396)</td><td>588 (±385)</td><td>-26 (±5)</td><td>482 (±20)</td><td>460 (±15)</td><td>-548 (±140)</td><td>9,970 (±144)</td><td>9,525 (±90)</td><td>1.4</td><td>-319</td><td>20.7</td></tr><tr><td>Tundra</td><td>-4 (±44)</td><td>302 (±125)</td><td>312 (±133)</td><td>7 (±3)</td><td>300 (±14)</td><td>306 (±12)</td><td>45 (±53)</td><td>2,049 (±49)</td><td>2,090 (±33)</td><td>0.9</td><td>64</td><td>6.8</td></tr><tr><td>Boreal</td><td>-42 (±82)</td><td>705 (±402)</td><td>664 (±398)</td><td>-43 (±7)</td><td>572 (±24)</td><td>537 (±17)</td><td>-593 (±101)</td><td>7,920 (±106)</td><td>7,435 (±74)</td><td>1.6</td><td>-383</td><td>13.9</td></tr><tr><td>Permafrost region</td><td>-21 (±62)</td><td>458 (±197)</td><td>445 (±171)</td><td>-15 (±5)</td><td>416 (±20)</td><td>405 (±16)</td><td>-249 (±123)</td><td>6,918 (±109)</td><td>6,719 (±69)</td><td>1.2</td><td>-24</td><td>16.6</td></tr></tbody></table>

Mean GPP, R<sub>eco</sub> and NEE fluxes and budgets over 2001–2020 and NEE + fire budgets over 2002–2020 are shown. The uncertainties are standard deviations across sites or pixels (for the mean fluxes) or across 20 bootstrapped budget estimates. Positive numbers for NEE indicate net CO₂ loss to the atmosphere, and negative numbers indicate net CO₂ uptake by the ecosystem. Mismatches in the site-level versus upscaled CO₂ fluxes are probably related to sites being biased to certain regions and years, while upscaled summaries should provide more representative regional estimates but are influenced by model performance.

<figure><img src="image_6.png" /><figcaption>**Fig. 1 |** Spatial variability in Arctic-boreal CO₂ fluxes. 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</figcaption></figure>

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⁻¹.

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).

## CO₂ budgets across the ABZ

Using machine learning models that had a high predictive performance (up to two times higher cross-validated R² than earlier efforts⁸,⁹), 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. 21)—that is, 2% of R<sub>eco</sub> 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
<figure><img src="image_7.png" /></figure>

**Fig. 2 | Trends in CO₂ budgets.** 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*<sub>eco</sub> trend maps identifying how trends in GPP and *R*<sub>eco</sub> 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.

(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.

## Model performance and comparison

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²³,²⁴), 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*<sub>eco</sub> 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).

## Temporal trends in upscaled ABZ CO₂ budgets

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*<sub>eco</sub>, probably associated with disturbed ecosystems (for example, forest fires and harvesting) with high *R*<sub>eco</sub> 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
<figure><img src="image_8.png" /><figcaption>**Fig. 3** | **Seasonal shifts in CO₂ flux dynamics.** **a–c**, Average upscaled monthly NEE (a), GPP (b) and *R*<sub>eco</sub> (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.</figcaption></figure>

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.

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*<sub>eco</sub> 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¹⁶,²⁸, 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²⁹,³⁰. 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*<sub>eco</sub>) 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).

## Drivers of ABZ CO₂ fluxes

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*<sub>eco</sub> (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.

## Continental and regional patterns in CO₂ budgets and their trends

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¹,³². 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
<figure><img src="image_9.png" /><figcaption>**Fig. 4 | Regional variability in CO₂ budget trends.** 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</figcaption></figure>

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.

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).

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¹⁶,²⁸—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.

## Discussion

Our results show that the ABZ was on average an increasing terrestrial CO₂ sink (GPP is increasing more than *R*<sub>eco</sub> + 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³⁻³⁶, 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<sub>eco</sub> instead³⁹⁻⁴³. 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⁴⁶⁻⁴⁸ 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⁵⁰,⁵¹, rapid recovery of ecosystems after fire⁵², and high cover of larch forests that can rapidly take up CO₂ (Supplementary Table 1)⁷,⁵³. 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.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-024-02234-5.

## References

1. Watts, J. D. et al. Carbon uptake in Eurasian boreal forests dominates the high-latitude net ecosystem carbon budget. *Glob. Change Biol.* **29**, 1870–1889 (2023).

2. Fan, L. et al. Siberian carbon sink reduced by forest disturbances. *Nat. Geosci.* **16**, 56–62 (2022).

3. McGuire, A. D. et al. Dependence of the evolution of carbon dynamics in the northern permafrost region on the trajectory of climate change. *Proc. Natl Acad. Sci. USA* **115**, 3882–3887 (2018).

4. McGuire, A. D. et al. Variability in the sensitivity among model simulations of permafrost and carbon dynamics in the permafrost region between 1960 and 2009. *Glob. Biogeochem. Cycles* **30**, 1015–1037 (2016).

5. Treat, C. C. et al. Permafrost carbon: progress on understanding stocks and fluxes across northern terrestrial ecosystems. *J. Geophys. Res. Biogeosci.* **129**, e2023JG007638 (2024).

6. Pallandt, M. M. et al. Representativeness assessment of the pan-Arctic eddy covariance site network and optimized future enhancements. *Biogeosciences* **19**, 559–583 (2021).

7. Virkkala, A.-M. et al. The ABCflux database: Arctic-boreal CO₂ flux observations and ancillary information aggregated to monthly time steps across terrestrial ecosystems. *Earth Syst. Sci. Data* **14**, 179–208 (2021).

8. Virkkala, A.-M. et al. Statistical upscaling of ecosystem CO₂ fluxes across the terrestrial tundra and boreal domain: regional patterns and uncertainties. *Glob. Change Biol.* https://doi.org/10.1111/gcb.15659 (2021).

9. Natali, S. M. et al. Large loss of CO₂ in winter observed across the northern permafrost region. *Nat. Clim. Change* **9**, 852–857 (2019).

10. Hugelius, G. et al. Permafrost region greenhouse gas budgets suggest a weak CO₂ sink and CH₄ and N₂O sources, but magnitudes differ between top-down and bottom-up methods. *Glob. Biogeochem. Cycles* **38**, e2023GB007969 (2024).

11. Hugelius, G. et al. Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps. *Biogeosciences* **11**, 6573–6593 (2014).

12. Schuur, E. A. G. et al. Permafrost and climate change: carbon cycle feedbacks from the warming Arctic. *Annu. Rev. Environ. Resour.* **47**, 343–371 (2022).

13. Rantanen, M. et al. The Arctic has warmed nearly four times faster than the globe since 1979. *Commun. Earth Environ.* **3**, 168 (2022).

14. Berner, L. T. et al. Summer warming explains widespread but not uniform greening in the Arctic tundra biome. *Nat. Commun.* **11**, 4621 (2020).

15. Lund, M. et al. Variability in exchange of CO₂ across 12 northern peatland and tundra sites. *Glob. Change Biol.* **16**, 2436–2448 (2010).

16. Euskirchen, E. S., Bret-Harte, M. S., Shaver, G. R., Edgar, C. W. & Romanovsky, V. E. Long-term release of carbon dioxide from Arctic tundra ecosystems in Alaska. *Ecosystems* **20**, 960–974 (2017).

17. Varner, R. K. et al. Permafrost thaw driven changes in hydrology and vegetation cover increase trace gas emissions and climate forcing in Stordalen Mire from 1970 to 2014. *Phil. Trans. R. Soc. A* **380**, 20210022 (2022).

18. Arndt, K. A., Hashemi, J., Natali, S. M., Schiferl, L. D. & Virkkala, A.-M. Recent advances and challenges in monitoring and modeling non-growing season carbon dioxide fluxes from the Arctic Boreal Zone. *Curr. Clim. Change Rep.* **9**, 27–40 (2023).
19. Watts, J. Regional hotspots of change in northern high latitudes informed by observations from space. Preprint at ESS Open Archive https://doi.org/10.22541/au.170497370.03373595/v1 (2024).

20. See, C. R. et al. Decadal increases in carbon uptake offset by respiratory losses across northern permafrost ecosystems. *Nat. Clim. Change* **14**, 853–862 (2024).

21. van Wees, D. et al. Global biomass burning fuel consumption and emissions at 500 m spatial resolution based on the Global Fire Emissions Database (GFED). *Geosci. Model Dev.* **15**, 8411–8437 (2022).

22. Friedlingstein, P. et al. Global carbon budget 2022. *Earth Syst. Sci. Data* **14**, 4811–4900 (2022).

23. Jung, M. et al. Scaling carbon fluxes from eddy covariance sites to globe: synthesis and evaluation of the FLUXCOM approach. *Biogeosciences* **17**, 1343–1365 (2020).

24. Nelson, J. A. et al. X-BASE: the first terrestrial carbon and water flux products from an extended data-driven scaling framework, FLUXCOM-X. *Biogeosciences* **21**, 5079–5115 (2024).

25. Eyring, V. et al. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.* **9**, 1937–1958 (2016).

26. Ueyama, M. et al. Carbon dioxide balance in early-successional forests after forest fires in interior Alaska. *Agric. For. Meteorol.* **275**, 196–207 (2019).

27. Damgaard, C. A critique of the space-for-time substitution practice in community ecology. *Trends Ecol. Evol.* **34**, 416–421 (2019).

28. Celis, G. et al. Tundra is a consistent source of CO₂ at a site with progressive permafrost thaw during 6 years of chamber and eddy covariance measurements. *J. Geophys. Res. Biogeosci.* **122**, 1471–1485 (2017).

29. Forkel, M. et al. Enhanced seasonal CO₂ exchange caused by amplified plant productivity in northern ecosystems. *Science* **351**, 696–699 (2016).

30. Graven, H. D. et al. Enhanced seasonal exchange of CO₂ by northern ecosystems since 1960. *Science* **341**, 1085–1089 (2013).

31. McGuire, A. D. et al. An assessment of the carbon balance of Arctic tundra: comparisons among observations, process models, and atmospheric inversions. *Biogeosciences* **9**, 3185–3204 (2012).

32. Lin, X. et al. Siberian and temperate ecosystems shape northern hemisphere atmospheric CO₂ seasonal amplification. *Proc. Natl Acad. Sci. USA* **117**, 21079–21087 (2020).

33. Jorgenson, M. T., Racine, C. H., Walters, J. C. & Osterkamp, T. E. Permafrost degradation and ecological changes associated with a warming climate in central Alaska. *Climatic Change* **48**, 551–579 (2001).

34. Euskirchen, E. S., Edgar, C. W., Turetsky, M. R., Waldrop, M. P. & Harden, J. W. Differential response of carbon fluxes to climate in three peatland ecosystems that vary in the presence and stability of permafrost. *J. Geophys. Res. Biogeosci.* **119**, 1576–1595 (2014).

35. Harris, L. I. et al. Permafrost thaw causes large carbon loss in boreal peatlands while changes to peat quality are limited. *Glob. Change Biol.* **29**, 5720–5735 (2023).

36. Euskirchen, E. S. et al. Persistent net release of carbon dioxide and methane from an Alaskan lowland boreal peatland complex. *Glob. Change Biol.* **30**, e17139 (2024).

37. Hayes, D. J. et al. Is the northern high-latitude land-based CO₂ sink weakening? *Glob. Biogeochem. Cycles* https://doi.org/10.1029/2010gb003813 (2011).

38. Commane, R. et al. Carbon dioxide sources from Alaska driven by increasing early winter respiration from Arctic tundra. *Proc. Natl Acad. Sci. USA* **114**, 5361–5366 (2017).

39. Zona, D. et al. Pan-Arctic soil moisture control on tundra carbon sequestration and plant productivity. *Glob. Change Biol.* **29**, 1267–1281 (2023).

40. Zona, D. et al. Earlier snowmelt may lead to late season declines in plant productivity and carbon sequestration in Arctic tundra ecosystems. *Sci. Rep.* **12**, 3986 (2022).

41. Byrne, B. et al. Multi-year observations reveal a larger than expected autumn respiration signal across northeast Eurasia. *Biogeosciences* **19**, 4779–4799 (2022).

42. Pulliainen, J. et al. Early snowmelt significantly enhances boreal springtime carbon uptake. *Proc. Natl Acad. Sci. USA* **114**, 11081–11086 (2017).

43. Parker, T. C., Unger, S. L., Moody, M. L., Tang, J. & Fetcher, N. Intraspecific variation in phenology offers resilience to climate change for *Eriophorum vaginatum*. *Arct. Sci.* https://doi.org/10.1139/as-2020-0039 (2021).

44. Brovkin, V. et al. Comparative carbon cycle dynamics of the present and last interglacial. *Quat. Sci. Rev.* **137**, 15–32 (2016).

45. Oechel, W. C. et al. Recent change of Arctic tundra ecosystems from a net carbon dioxide sink to a source. *Nature* **361**, 520–523 (1993).

46. Lawrence, D. M., Koven, C. D., Swenson, S. C., Riley, W. J. & Slater, A. G. Permafrost thaw and resulting soil moisture changes regulate projected high-latitude CO₂ and CH₄ emissions. *Environ. Res. Lett.* **10**, 094011 (2015).

47. Parker, T. C., Subke, J.-A. & Wookey, P. A. Rapid carbon turnover beneath shrub and tree vegetation is associated with low soil carbon stocks at a subarctic treeline. *Glob. Change Biol.* **21**, 2070–2081 (2015).

48. Voigt, C. et al. Warming of subarctic tundra increases emissions of all three important greenhouse gases—carbon dioxide, methane, and nitrous oxide. *Glob. Change Biol.* **23**, 3121–3138 (2017).

49. Piao, S. et al. Characteristics, drivers and feedbacks of global greening. *Nat. Rev. Earth Environ.* **1**, 14–27 (2019).

50. Frost, G. V. & Epstein, H. E. Tall shrub and tree expansion in Siberian tundra ecotones since the 1960s. *Glob. Change Biol.* **20**, 1264–1277 (2014).

51. Kharuk, V. I., Ranson, K. J., Im, S. T. & Petrov, I. A. Climate-induced larch growth response within the central Siberian permafrost zone. *Environ. Res. Lett.* **10**, 125009 (2015).

52. Schulze, E.-D. et al. Factors promoting larch dominance in central Siberia: fire versus growth performance and implications for carbon dynamics at the boundary of evergreen and deciduous conifers. *Biogeosciences* **9**, 1405–1421 (2012).

53. Hiyama, T. et al. Lessons learned from more than a decade of greenhouse gas flux measurements at boreal forests in eastern Siberia and interior Alaska. *Polar Sci.* **27**, 100607 (2021).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025
Anna-Maria Virkkala Ⓓ1, Brendan M. Rogers Ⓓ1, Jennifer D. Watts<fnref n="1" />, Kyle A. Arndt Ⓓ1, Stefano Potter Ⓓ1, Isabel Wargowsky Ⓓ1, Edward A. G. Schuur Ⓓ2,3, Craig R. See², Marguerite Mauritz Ⓓ4, Julia Boike<fnref n="5" />,⁶, M. Syndonia Bret-Harte Ⓓ7, Eleanor J. Burke Ⓓ8, Arden Burrell<fnref n="1" />, Namyi Chae Ⓓ9, Abhishek Chatterjee Ⓓ10, Frederic Chevallier Ⓓ11, Torben R. Christensen Ⓓ12,13, Roisin Commane Ⓓ14, Han Dolman Ⓓ15,16, Colin W. Edgar Ⓓ7, Bo Elberling Ⓓ17, Craig A. Emmerton Ⓓ18, Eugenie S. Euskirchen Ⓓ7, Liang Feng¹⁹, Mathias Göckede Ⓓ20, Achim Grelle Ⓓ21, Manuel Helbig Ⓓ22,23,24, David Holl Ⓓ25, Järvi Järveoja Ⓓ26, Sergey V. Karsanaev Ⓓ27, Hideki Kobayashi Ⓓ28, Lars Kutzbach Ⓓ25, Junjie Liu Ⓓ10, Ingrid T. Luijkx Ⓓ29, Efrén López-Blanco Ⓓ12,30, Kyle Lunneberg³¹, Ivan Mammarella Ⓓ32, Maija E. Marushchak Ⓓ33, Mikhail Mastepanov Ⓓ12,13, Yojiro Matsuura³⁴, Trofim C. Maximov Ⓓ27, Lutz Merbold Ⓓ35, Gesa Meyer Ⓓ23,36, Mats B. Nilsson Ⓓ26, Yosuke Niwa Ⓓ37,38, Walter Oechel Ⓓ31, Paul I. Palmer Ⓓ19, Sang-Jong Park Ⓓ39, Frans-Jan W. Parmentier Ⓓ40, Matthias Peichl Ⓓ26, Wouter Peters Ⓓ29,41, Roman Petrov Ⓓ27, William Quinton Ⓓ42, Christian Rödenbeck Ⓓ20, Torsten Sachs Ⓓ24,43, Christopher Schulze Ⓓ23,44, Oliver Sonnentag Ⓓ23, Vincent L. St. Louis Ⓓ18, Eeva-Stiina Tuittila Ⓓ45, Masahito Ueyama Ⓓ46, Andrej Varlagin Ⓓ47, Donatella Zona Ⓓ31 & Susan M. Natali<fnref n="1" />

<footnote n="1">¹Woodwell Climate Research Center, Falmouth, MA, USA. ²Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, AZ, USA.</footnote>

<footnote n="3">³Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, USA. ⁴Biological Sciences, University of Texas at El Paso, El Paso, TX, USA.</footnote>

<footnote n="5">⁵Permafrost Research Section, Alfred Wegener Institute Helmholtz Center for Polar and Marine Research, Potsdam, Germany. ⁶Department of Geography, Humboldt University, Berlin, Germany. ⁷Institute of Arctic Biology, University of Alaska Fairbanks, Fairbanks, AK, USA. ⁸Met Office Hadley Centre, Exeter, UK. ⁹Institute of Life Science and Natural Resources, Korea University, Seoul, South Korea. ¹⁰NASA Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA. ¹¹Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France. ¹²Department of Ecoscience, Arctic Research Center, Aarhus University, Roskilde, Denmark. ¹³Water, Energy and Environmental Engineering Research Unit, University of Oulu, Oulu, Finland. ¹⁴Department of Earth and Environmental Sciences, Columbia University, New York, NY, USA.</footnote>

<footnote n="15">¹⁵Royal Netherlands Institute for Sea Research, Texel, the Netherlands. ¹⁶Vrije Universiteit, Amsterdam, the Netherlands. ¹⁷Department of Geosciences and Natural Resource Management, University of Copenhagen, Copenhagen, Denmark. ¹⁸Department of Biological Sciences, University of Alberta, Edmonton, Alberta, Canada. ¹⁹National Centre for Earth Observation, University of Edinburgh, Edinburgh, UK. ²⁰Max Planck Institute for Biogeochemistry, Jena, Germany. ²¹Linnaeus University, Växjö, Sweden. ²²Department of Physics and Atmospheric Science, Dalhousie University, Halifax, Canada.</footnote>

<footnote n="23">²³Département de Géographie, Université de Montréal, Montréal, Québec, Canada. ²⁴GFZ German Research Centre for Geosciences, Potsdam, Germany.</footnote>

<footnote n="25">²⁵Institute of Soil Science, Center for Earth System Research and Sustainability (CEN), Universität Hamburg, Hamburg, Germany. ²⁶Department of Forest Ecology and Management, Swedish University of Agricultural Sciences, Umeå, Sweden. ²⁷Institute for Biological Problems of Cryolithozone of the Siberian Branch of the RAS, Division of Federal Research Centre “The Yakut Scientific Centre of the Siberian Branch of the Russian Academy of Sciences”, Yakutsk, Russia. ²⁸Research Institute for Global Change, Japan Agency for Marine-Earth Science and Technology, Yokosuka, Japan. ²⁹Department of Meteorology and Air Quality, Wageningen University, Wageningen, the Netherlands. ³⁰Department of Environment and Minerals, Greenland Institute of Natural Resources, Nuuk, Greenland. ³¹Department Biology, San Diego State University, San Diego, CA, USA. ³²Institute for Atmospheric and Earth System Research / Physics, Faculty of Science, University of Helsinki, Helsinki, Finland. ³³Department of Environmental and Biological Sciences, University of Eastern Finland, Kuopio, Finland. ³⁴Forestry and Forest Products Research Institute, Tsukuba, Ibaraki, Japan. ³⁵Integrative Agroecology, Agroecology and Environment, Agroscope, Zurich, Switzerland. ³⁶Climate Research Division, Environment and Climate Change Canada, Victoria, British Columbia, Canada. ³⁷Earth System Division, National Institute for Environmental Studies, Tsukuba, Japan. ³⁸Department of Climate and Geochemistry Research, Meteorological Research Institute, Tsukuba, Japan. ³⁹Division of Ocean and Atmosphere Sciences, Korea Polar Research Institute, Incheon, Republic of Korea. ⁴⁰Centre for Biogeochemistry in the Anthropocene, Department of Geosciences, University of Oslo, Oslo, Norway. ⁴¹Centre for Isotope Research, Energy and Sustainability Research Institute Groningen, Groningen University, Groningen, the Netherlands. ⁴²Cold Regions Research Centre, Wilfrid Laurier University, Waterloo, Ontario, Canada. ⁴³Institute of Geoecology, Technische Universität Braunschweig, Braunschweig, Germany. ⁴⁴Department of Renewable Resources, University of Alberta, Edmonton, Alberta, Canada. ⁴⁵School of Forest Sciences, University of Eastern Finland, Joensuu, Finland. ⁴⁶Graduate School of Agriculture, Osaka Metropolitan University, Sakai, Japan. ⁴⁷A.N. Severtsov Institute of Ecology and Evolution, Russian Academy of Sciences, Moscow, Russia. ✉e-mail: avirkkala@woodwellclimate.org</footnote>
# Methods

## 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*<sub>eco</sub> 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⁵⁴,⁵⁵, 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*<sub>eco</sub> 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*<sub>eco</sub> models (Supplementary Table 4). These differences in data distribution resulted in slight mismatches in the upscale NEE versus GPP – *R*<sub>eco</sub> 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. 60), 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 upscale 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*<sub>eco</sub> 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<sub>eco</sub>; 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²⁶,⁶⁶. 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. 68), 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<sub>eco</sub> 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⁷⁰,⁷¹ 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 perma-
frost region⁷³ and vegetation types⁷. Supplementary Table I 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²⁵, atmos-
pheric inversions used in the Global Carbon Project's Global Carbon
Budget 2022²²,⁷⁴ 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⁵,⁷⁵ (Sup-
plementary Table 6). The final ensemble outputs used here represent
terrestrial NEE (GPP – R<sub>eco</sub>) 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<sup>-1</sup>, respectively,
in the northern permafrost region (that is, greater than the NEE + fire
budget calculated in this study for the same domain)⁷⁶.

Data availability

The in situ data used here can be accessed via ORNL DAAC (https://doi.
org/10.3334/ORNLDAAC/1934)⁷⁷, and the geospatial data are available
via the links and references provided in Supplementary Tables 3 and 6.
The 1-km upscaled rasters of NEE, GPP and R<sub>eco</sub> are available via ORNL
DAAC (https://doi.org/10.3334/ORNLDAAC/2377)⁷⁸.

Code availability

The analysis code is available via Zenodo at https://doi.org/10.5281/
zenodo.13691584 (ref. 61).

References

54. Pallandt, M. M. et al. Representativeness assessment of the pan-Arctic eddy covariance site network and optimized future enhancements. *Biogeosciences* **19**, 559–583 (2022).

55. Pallandt, M. M. T. A., et al. High-latitude eddy covariance temporal network design and optimization. *J. Geophys. Res. G* **129**, e2024JG008406 (2024).

56. Mu, C. et al. Ecosystem CO₂ exchange and its economic implications in northern permafrost regions in the 21st century. *Global Biogeochem. Cycles* **37**, e2023GB007750 (2023).

57. Giglio, L., Boschetti, L., Roy, D. P., Humber, M. L. & Justice, C. O. The Collection 6 MODIS burned area mapping algorithm and product. *Remote Sens. Environ.* **217**, 72–85 (2018).

58. Giglio, L., Schroeder, W. & Justice, C. O. The Collection 6 MODIS active fire detection algorithm and fire products. *Remote Sens. Environ.* **178**, 31–41 (2016).

59. Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. *Science* **342**, 850–853 (2013).

60. Kuhn, M. Building predictive models in R using the caret package. *J. Stat. Softw.* **28**, 1–26 (2008).

61. Virkkala, A.-M. et al. Code for ‘An increasing Arctic-boreal CO₂ sink offset by wildfires and source regions’ (version 1). *Zenodo* https://doi.org/10.5281/zenodo.13691585 (2024).

62. Greenwell, B. M. pdp: an R package for constructing partial dependence plots. *R J.* **9**, 421–436 (2017).

63. Greenwell, B. M. & Boehmke, B. C. Variable importance plots—an introduction to the vip package. *R J.* **12**, 343–366 (2020).

64. Breiman, L. Random forests. *Mach. Learn.* **45**, 5–32 (2001).

65. Tramontana, G. et al. Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms. *Biogeosciences* **13**, 4291–4313 (2016).

66. Amiro, B. D., Ian MacPherson, J., Desjardins, R. L., Chen, J. M. & Liu, J. Post-fire carbon dioxide fluxes in the western Canadian boreal forest: evidence from towers, aircraft and remote sensing. *Agric. For. Meteorol.* **115**, 91–107 (2003).

67. Pinheiro, J., Bates, D. & R Core Team. nlme: Linear and nonlinear mixed effects models. R package version 3.1-149 (2020).

68. Dinerstein, E. et al. An ecoregion-based approach to protecting half the terrestrial realm. *Bioscience* **67**, 534–545 (2017).

69. Hijmans, R. J. terra: Spatial data analysis. R package version https://doi.org/10.32614/cran/package.terra (2020).

70. Sen, P. K. Estimates of the regression coefficient based on Kendall's tau. *J. Am. Stat. Assoc.* **63**, 1379–1389 (1968).

71. Bronaugh, D., Schoeneberg, A. & Zeman, L. zyp: Zhang + Yue-Pilon trends package. R package version 0.11-1 (2023).

72. Zhang, X., Vincent, L. A., Hogg, W. D. & Niitsoo, A. Temperature and precipitation trends in Canada during the 20th century. *Atmos. Ocean* **38**, 395–429 (2000).

73. Heginbottom, J., Brown, J., Ferrians, O. & Melnikov, E. S. *Circum-Arctic Map of Permafrost and Ground-Ice Conditions, Version 2* (NSIDC, 2002); https://doi.org/10.7265/SKBG-KF16

74. Luijkx, I. T. et al. Global CO₂ Gridded Flux Fields from 9 Atmospheric Inversions in GCB2022 (ICOS Data Portal, 2023); https://doi.org/10.18160/7AH8-K1X4

75. Liu, Z. et al. Respiratory loss during late-growing season determines the net carbon dioxide sink in northern permafrost regions. *Nat. Commun.* **13**, 5626 (2022).

76. Ramage, J. L. et al. The net GHG balance and budget of the permafrost region (2000–2020) from ecosystem flux upscaling. *Glob. Biogeochem. Cycles* **38**, e2023GB007953 (2023).

77. Virkkala, A.-M. et al. The ABCflux Database: Arctic-Boreal CO₂ Flux and Site Environmental Data, 1989–2020 (ORNL DAAC, 2021); https://doi.org/10.3334/ORNLDAAC/1934

78. Virkkala, A.-M. et al. Machine Learning-Based Arctic-Boreal Terrestrial Ecosystem CO₂ Fluxes, 2001–2020 (ORNL DAAC, 2024); https://doi.org/10.3334/ORNLDAAC/2377

Acknowledgements

This work was supported by funding from the Gordon and Betty
Moore Foundation (grant number 8414; A.-M.V., B.M.R., J.D.W., S.P.,
S.M.N., K.A.A. and I.W.) and funding catalysed by the TED Audacious
Project (Permafrost Pathways; A.-M.V., B.M.R., J.D.W., S.P., S.M.N.,
K.A.A. and I.W.). We additionally acknowledge funding from the NASA
Arctic-Boreal Vulnerability Experiment and Carbon Cycle Science
programmes (grant number NNX17AE13G; B.M.R.); NSF PLR Arctic
System Science Research Networking Activities (grant number
1931333; E.A.G.S. and M. Mauritz); Minderoo Foundation (E.A.G.S.);
NSF (grant numbers DEB LTREB 1354370, 2011257, DEB-0425328,
DEB-0724514 and DEB-0830997; E.S.E. and C.W.E.); the US Geological
Survey Climate R&D Program (E.S.E., M.S.B.-H. and C.W.E.); NSF Arctic
Observatory Network (grant number 1936752; M.S.B.-H.; and grant
numbers 1503912 and 1107892; E.S.E., M.S.B.-H. and C.W.E.); Office
of Polar Programs of the NSF (award numbers 2149988 and 1932900;
D.Z.); KAKENHI (grant number 19H05668; M.U.); the Danish National
Research Foundation (CENPERM DNRF100; B.E.); EU HORIZON
GreenFeedBack, grant agreement number 101056921 (I.M., M.G.,
E.L.-B., T.R.C. and M. Mastepanov) and LiweFor number 101079192
(I.M.); ICOS-FI via University of Helsinki funding (I.M.); Danish Arctic
Climate support through Greenland Ecosystem Monitoring and ICOS grants (E.L.-B., T.R.C. and M. Mastepanov); Natural Sciences and Engineering Research Council (NSERC; V.L.S.L. and C.A.E.); the Deutsche Forschungsgemeinschaft (German Research Foundation) under Germany's Excellence Strategy—EXC 2037 'CLICCS—Climate, Climatic Change, and Society'—project number 390683824 (L.K. and D.H.); Met Office Hadley Centre Climate Programme funded by DSIT and European Union's Horizon 2020 research and innovation programme under Grant Agreement Number 101003536 (ESM2025—Earth System Models for the Future; E.J.B.); NASA Grant/Cooperative Agreement No. NNX17AD69A (A.C.); the Research Council of Norway (BioGov, project no. 323945; F.-J.W.P.); European Research Council (ERC synergy project Q-Arctic, grant agreement number 951288; M.G.); the Copernicus Atmosphere Monitoring Service, implemented by ECMWF on behalf of the European Commission (grant number CAMS2 55; F.C.); the Environment Research and Technology Development Fund of the Environmental Restoration and Conservation Agency of Japan (grant number JPMEERF21S20810; Y.N.); ArCSII (grant no. JPMXD142031886; M.U. and H.K.); the Swedish Research Council (VR) and consortium partners to ICOS Sweden (grant numbers 2015–06020 and 2019–00205; M.P. and M.B.N.) and SITES (grant number 2017–00635; M.P. and M.B.N.); VR grant numbers 2019-04676 and 2018-03966 (M.P. and M.B.N.); ArcticNet and NSERC (W.Q.); NOAA Cooperative Agreement NA16SEC4810008 (W.O.); Research Council of Finland (project numbers 341349, 330840, 349503 ICOS-FIRI and ACCC—Atmosphere and Climate Competence Center; I.M., M.E.M. and E.-S.T.); NRF-2021M1A5A1065425 (KOPRI-PN24011; S.-J.P.); NRF-2021M1A5A1065679 (N.C.); NRF-2021R1I1A1A01053870 (N.C.); the Dutch Research Council (NWO) (project number VI.Vidi.213.143; I.T.L.); the UK National Centre for Earth Observation funded by the National Environment Research Council (grant number NE/R016518/1) and the UK Space Agency (L.F. and P.I.P.); Ministry of Economic Development and Commerce of the Russian Federation (registration number 123030300031-6; A.V., R.P., T.C.M. and S.V.K.); NSF OPP 1936752 (M.S.B.-H.); NSERC Discovery Grant (M.H.); Helmholtz Association in the framework of MOSES (Modular Observation Solutions for Earth Systems: J.B.); and the Natural Environment Research Council through the National Centre for Earth Observation (grant number NE/R000115/1; L.F.). Part of the research was carried out at the Jet Propulsion Laboratory,

California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004; A.C. and J.L.). Part of the inverse analyses were performed on the supercomputer systems at the National Institute for Environmental Studies and Meteorological Research Institutes (NEC SX-Aurora TSUBASA and FUJITSU PRIMERGY CX2550M5; Y.N.) and at the HPC cluster Aether at the University of Bremen, financed by the Deutsche Forschungsgemeinschaft within the scope of the Excellence Initiative (Germany; I.T.L.). We further thank the Arctic Data Center Working Group: Reconciling historical and contemporary trends in terrestrial carbon exchange of the northern permafrost-zone, and we thank R. Treharne, T. Smith and Y. Yang for help with the analysis.

## Author contributions

A.-M.V. together with B.M.R., J.D.W. and S.M.N. designed the work. A.-M.V. led the paper, wrote the draft and, together with S.P. and K.A.A., analysed the data. B.M.R., K.A.A., S.M.N. and J.D.W. substantively revised the draft. All other authors acquired and interpreted the data and edited the drafts.

## Competing interests

The authors declare no competing interests.

## Additional information

**Extended data** is available for this paper at
https://doi.org/10.1038/s41558-024-02234-5.

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s41558-024-02234-5.

**Correspondence and requests for materials** should be addressed to Anna-Maria Virkkala.

**Peer review information** *Nature Climate Change* thanks Shaorun Lin, Shawn Pedron and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at
www.nature.com/reprints.
a

<figure><img src="image_10.png" /></figure>

Clusters

b

<figure><img src="image_11.png" /></figure>

Extended Data Fig. 1 | Clustering of the environmental conditions controlling CO2 fluxes across the ABZ. Results showing the clusters on a map (a), the percentage of vegetation types within each cluster (b), and the variability in key

<figure><img src="image_12.png" /></figure>

environmental conditions across the clusters based on a random spatial sample of 10,000 pixels per cluster (c). A description of the clusters can be found in Supplementary Methods Section 4.2.
a

<figure><img src="image_13.png" /></figure>

b

<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 2** | Predictive performance of the NEE model. The monthly NEE variability of observed fluxes and predicted fluxes, based on the same train and test dataset (no cross validation; no CV) and a leave-one-site-out CV (CV) approach (a). The model training dataset comprised a total of 4765 samples. Monthly sample sizes were distributed as follows: 288 in January; 289 in February; 329 in March; 377 in April; 449 in May; 497 in June; 527 in July; 522 in August; 475 in September; 386 in October; 320 in November; and 306 in December. The lower and upper hinges in the boxes correspond to the first and third quartiles (the 25th and 75th percentiles) and the line between the median. The whiskers extend from the quartiles to the furthest data points within 1.5 \* IQR, where IQR is the interquartile range; data beyond the end of the whisker are visualized with dots. Subplots b and c

c

<figure><img src="image_15.png" /></figure>

show the correlation between observed and predicted fluxes, colored by the density of observations (b) or the deviance from average site-level monthly flux (c). Subfigure c indicates that the model struggles the most when observations from individual sites have a large deviance from the mean. Coefficient of determination (R²) describes the strength of the linear relationship between the observed and predicted fluxes. Mean bias error (MBE) characterizes the average bias between prediction and observation, with negative values indicating the model to underestimate NEE (that is, overestimate net CO₂ sinks or underestimate net CO₂ sources). Mean absolute error (MAE) describes the absolute bias between prediction and observation, with larger values describing larger errors. For a similar figure for GPP and R<sub>eco</sub>, see Supplementary Figs. 18 and 19.
Variable

Month
Solar radiation
Land surface temperature (daytime)

NDVI

Soil temperature

Snow cover

Vapor pressure deficit

Vegetation type

Percent tree cover

Soil organic carbon stock

Soil pH

Atmospheric CO2

Compound topographic index

Percent non-tree vegetation

Percent non-vegetated

Permafrost probability

Volumetric soil water

<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 3** | Variable importance and partial dependence plots for the most important predictors of the 1-km NEE model. The variable importance boxes show the scores from 100 simulations, with the lower and upper hinges corresponding to the first and third quartiles (the 25th and 75th percentiles) and the line between the median. The whiskers extend from the quartiles to the furthest data points within 1.5 * IQR, where IQR is the interquartile range; data beyond the end of the whisker are visualized with dots. The values on the y axis of each partial dependence plot can be interpreted as followed: yhat is conditional

on other predictors in the model and their relationships with the predictor in the plot in question. Therefore, yhat values should not be directly compared with observed or predicted values, rather the patterns in yhat should be explored more generally. The x-axis represents the actual predictor values and can be used to infer, for example, conditions that lead to changes in yhat (tipping points), and the strength and direction of the relationship. For a similar figure for GPP and R<sub>eco</sub>, see Supplementary Figs. 20 and 21.
Fire emissions Tg C yr⁻¹

Alaskan boreal

Trend = -0.15, p=0.363

Average = 12

<figure><img src="image_17.png" /></figure>

Alaskan tundra

Trend = -0.13, p=0.363

Average = 4

<figure><img src="image_18.png" /></figure>

Canadian boreal

Trend = -0.52, p=0.762

Average = 56

<figure><img src="image_19.png" /></figure>

Canadian tundra

Trend = -0.02, p=0.780

Average = 4

<figure><img src="image_20.png" /></figure>

European boreal

Trend = -0.02, p=0.363

Average = 1

<figure><img src="image_21.png" /></figure>

European tundra

Trend = 0, p=0.649

Average = 0

<figure><img src="image_22.png" /></figure>

Siberian boreal

Trend = 9.7, p=0.023

Average = 148

<figure><img src="image_23.png" /></figure>

Siberian tundra

Trend = 0.06, p=0.834

Average = 12

<figure><img src="image_24.png" /></figure>

Year
---
500-m fire emissions from GFED (2002-2020)

**Extended Data Fig. 4** | Annual fire emission budgets across the key domains. The magnitude of the trend was computed using the Theil-Sen approach, and the p-value determined using the Mann-Kendall test.
<figure><img src="image_25.png" /><figcaption>Extended Data Fig. 5 | Trends in upscaled shoulder season mean NEE fluxes across the key domains from 2001 to 2020. The magnitude of the trend was computed using the Theil-Sen approach, and the p-value determined using the Mann-Kendall test.</figcaption></figure>
<figure><img src="image_26.png" /></figure>

**Extended Data Fig. 6** | Model performance across sites in disturbed vs. non-disturbed conditions. Model fit (that is, no cross validation) and cross-validated predictive performance estimates for each site across sites with recent (<20 years since disturbance) fire disturbance (n = 10), other disturbance (for example, permafrost thaw, harvest; n = 6), or no disturbance or information about disturbance (n = 159). The lower and upper hinges in the boxes correspond to the first and third quartiles (the 25th and 75th percentiles) and the line between the median. The whiskers extend from the quartiles to the furthest data points within

<figure><img src="image_27.png" /></figure>

1.5 *IQR*, where IQR is the interquartile range; data beyond the end of the whisker are visualized with dots. Coefficient of determination (R²) describes the strength of the linear relationship between the observed and predicted fluxes. Mean bias error (MBE) characterizes the average bias between prediction and observation, with negative values indicating the model to underestimate NEE (that is, overestimate net CO₂ sinks or underestimate net CO₂ sources). Mean absolute error (MAE) describes the absolute bias between prediction and observation, with larger values describing larger errors.

Recent disturbance