# Statistical upscaling of ecosystem CO 2 fluxes across the terrestrial tundra and boreal domain: Regional patterns and uncertainties

Anna-maria Virkkala, Juha Aalto, Brendan Rogers, Torbern Tagesson, Claire

Treat, Susan Natali, Jennifer Watts, Stefano Potter, Aleksi Lehtonen,
Marguerite Mauritz, et al.

► To cite this version:

Anna-maria Virkkala, Juha Aalto, Brendan Rogers, Torbern Tagesson, Claire Treat, et al. Statistical upscaling of ecosystem CO 2 fluxes across the terrestrial tundra and boreal domain: Regional patterns and uncertainties. Global Change Biology, 2021, 27 (17), pp.4040-4059. 10.1111/gcb.15659. hal-03260396

HAL Id: hal-03260396

https://hal.science/hal-03260396

Submitted on 16 Jun 2021

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Statistical upscaling of ecosystem CO₂ fluxes across the terrestrial tundra and boreal domain: Regional patterns and uncertainties

Anna-Maria Virkkala<fnref n="1" />,<fnref n="2" /> ID | Juha Aalto<fnref n="1" />,<fnref n="3" /> ID | Brendan M. Rogers<fnref n="2" /> ID | Torbern Tagesson<fnref n="4" />,<fnref n="5" /> ID |

Claire C. Treat<fnref n="6" /> ID | Susan M. Natali<fnref n="2" /> ID | Jennifer D. Watts<fnref n="2" /> | Stefano Potter<fnref n="2" /> |

Aleksi Lehtonen<fnref n="7" /> ID | Marguerite Mauritz<fnref n="8" /> ID | Edward A. G. Schuur<fnref n="9" /> ID |

John Kochendorfer<fnref n="10" /> ID | Donatella Zona<fnref n="11" />,<fnref n="12" /> ID | Walter Oechel<fnref n="11" />,<fnref n="13" /> ID | Hideki Kobayashi<fnref n="14" /> ID |

Elyn Humphreys<fnref n="15" /> ID | Mathias Goeckede<fnref n="16" /> ID | Hiroki Iwata<fnref n="17" /> ID | Peter M. Lafleur<fnref n="18" /> ID |

Eugenie S. Euskirchen¹⁹ ID | Stef Bokhorst²⁰ ID | Maija Marushchak²¹,²² ID |

Pertti J. Martikainen²² ID | Bo Elberling²³ ID | Carolina Voigt²²,²⁴ ID | Christina Biasi²² ID |

Oliver Sonnentag²⁴ ID | Frans-Jan W. Parmentier<fnref n="4" />,²⁵ ID | Masahito Ueyama²⁶ ID |

Gerardo Celis²⁷ ID | Vincent L. St. Louis²⁸ ID | Craig A. Emmerton²⁸ ID | Matthias Peichl²⁹ ID |

Jinshu Chi²⁹ ID | Järvi Järveoja²⁹ ID | Mats B. Nilsson²⁹ ID | Steven F. Oberbauer³⁰ ID |

Margaret S. Torn³¹ ID | Sang-Jong Park³² ID | Han Dolman³³ ID | Ivan Mammarella³⁴ ID |

Namyi Chae³⁵ ID | Rafael Poyatos³⁶,³⁷ ID | Efrén López-Blanco³⁸,³⁹ ID |

Torben Røjle Christensen³⁹ ID | Min Jung Kwon⁴⁰,⁴¹ ID | Torsten Sachs⁴² ID |

David Holl⁴³ ID | Miska Luoto<fnref n="1" /> ID

<footnote n="1">¹Department of Geosciences and Geography, Faculty of Science, University of Helsinki, Helsinki, Finland</footnote>

<footnote n="2">²Woodwell Climate Research Center, Falmouth, MA, USA</footnote>

<footnote n="3">³Weather and Climate Change Impact Research, Finnish Meteorological Institute, Helsinki, Finland</footnote>

<footnote n="4">⁴Department of Physical Geography and Ecosystem Science, Lund University, Lund, Sweden</footnote>

<footnote n="5">⁵Department of Geosciences and Natural Resource Management, Copenhagen University, Copenhagen, Denmark</footnote>

<footnote n="6">⁶Alfred Wegener Institute Helmholtz Center for Polar and Marine Research, Potsdam, Germany</footnote>

<footnote n="7">⁷Natural Resources Institute Finland, Helsinki, Finland</footnote>

<footnote n="8">⁸University of Texas at El Paso, El Paso, TX, USA</footnote>

<footnote n="9">⁹Center for Ecosystem Science and Society, Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, USA</footnote>

<footnote n="10">¹⁰Atmospheric Turbulence and Diffusion Division of NOAA's Air Resources Laboratory, Oak Ridge, TN, USA</footnote>

<footnote n="11">¹¹San Diego State University, San Diego, CA, USA</footnote>

<footnote n="12">¹²University of Sheffield, Sheffield, UK</footnote>

<footnote n="13">¹³University of Exeter, Exeter, UK</footnote>

<footnote n="14">¹⁴Research Institute for Global Change, Japan Agency for Marine-Earth Science and Technology, Yokoama, Japan</footnote>

<footnote n="15">¹⁵Carleton University, Ottawa, ON, Canada</footnote>

<footnote n="16">¹⁶Dept. Biogeochemical Signals, Max Planck Institute for Biogeochemistry, Jena, Germany</footnote>

<footnote n="17">¹⁷Department of Environmental Science, Shinshu University, Matsumoto, Japan</footnote>

<footnote n="18">¹⁸School of the Environment, Trent University, Peterborough, ON, Canada</footnote>
¹⁹Institute of Arctic Biology, University of Alaska Fairbanks, Fairbanks, AK, USA

²⁰Vrije Universiteit Amsterdam, Amsterdam, the Netherlands

²¹Department of Biological and Environmental Science, University of Jyväskylä, Jyväskylä, Finland

²²Department of Environmental and Biological Sciences, University of Eastern Finland, Kuopio, Finland

²³Center for Permafrost, Department of Geoscience and Natural Resource Management, University of Copenhagen, Copenhagen, Denmark

²⁴Département de géographie, Université de Montréal, Montréal, QC, Canada

²⁵Centre for Biogeochemistry in the Anthropocene, Department of Geosciences, University of Oslo, Oslo, Norway

²⁶Graduate School of Life and Environmental Sciences, Osaka Prefecture University, Sakai, Japan

²⁷Agronomy Department, University of Florida, Gainesville, FL, USA

²⁸Department of Biological Sciences, University of Alberta, Edmonton, AB, Canada

²⁹Department of Forest Ecology and Management, Swedish University of Agricultural Sciences, Umeå, Sweden

³⁰Department of Biological Sciences, Florida International University, Miami, FL, USA

³¹Berkeley Lab and UC Berkeley, Berkeley, CA, USA

³²Division of Atmospheric Sciences, Korea Polar Research Institute, Incheon, Republic of Korea

³³Department of Earth Sciences, Free University Amsterdam, Amsterdam, the Netherlands

³⁴Institute for Atmospheric and Earth System Research/Physics, Faculty of Science, University of Helsinki, Helsinki, Finland

³⁵Institute of Life Science and Natural Resources, Korea University, Seoul, Republic of Korea

³⁶CREAF, Catalonia, Spain

³⁷Universitat Autònoma de Barcelona, Barcelona, Spain

³⁸Department of Environment and Minerals, Greenland Institute of Natural Resources, Nuuk, Greenland

³⁹Department of Bioscience, Arctic Research Center, Aarhus University, Roskilde, Denmark

⁴⁰Laboratoire des Sciences du Climat et de l'Environnement, Gif-sur-Yvette, France

⁴¹Division of Life Sciences, Korea Polar Research Institute, Incheon, Republic of Korea

⁴²GFZ German Research Centre for Geosciences, Potsdam, Germany

⁴³Center for Earth System Research and Sustainability (CEN), University of Hamburg, Hamburg, Germany

### Correspondence

Anna-Maria Virkkala, Woodwell Climate Research Center, Falmouth, MA, USA.
Email: avirkkala@woodwellclimate.org

### Funding information

Humboldt Fellowship for Experienced Researchers; European Commission, Grant/Award Number: H2020-BG-09-2016 and 727890; Helsingin Yliopisto; Jenny ja Antti Wihurin Rahasto; Väisälä fund; Natural Sciences and Engineering Research Council of Canada; Alfred Kordelinin Säätiö; Arctic Challenge for Sustainability, Grant/Award Number: JPMXD1420318865; Suomen Kulttuurirahasto; Netherlands Earth System Science Centre; Korean government, Grant/Award Number: KOPRI-PN21011, NRF-2021M1A5A1065679, NRF-2021M1A5A1065425 and NRF-2018R1D1A1B07047778; Gordon and Betty Moore Foundation, Grant/Award Number: 8414; Skogssällskapet, Grant/Award Number: 2018-485-Steg 2 2017; Office of Biological and Environmental Research; Natural Sciences and Engineering Research Council; Svenska Forskningsrådet Formas, Grant/Award Number: 2016-01289 and 942-2015-49; Vetenskapsrådet, Grant/Award Number: 2017-05268; Norges Forskningsråd, Grant/Award Number: 274711; NASA,

### Abstract

The regional variability in tundra and boreal carbon dioxide ($CO_2$) fluxes can be high, complicating efforts to quantify sink-source patterns across the entire region. Statistical models are increasingly used to predict (i.e., upscale) $CO_2$ fluxes across large spatial domains, but the reliability of different modeling techniques, each with different specifications and assumptions, has not been assessed in detail. Here, we compile eddy covariance and chamber measurements of annual and growing season $CO_2$ fluxes of gross primary productivity (GPP), ecosystem respiration (ER), and net ecosystem exchange (NEE) during 1990–2015 from 148 terrestrial high-latitude (i.e., tundra and boreal) sites to analyze the spatial patterns and drivers of $CO_2$ fluxes and test the accuracy and uncertainty of different statistical models. $CO_2$ fluxes were up-scaled at relatively high spatial resolution (1 km$^2$) across the high-latitude region using five commonly used statistical models and their ensemble, that is, the median of all five models, using climatic, vegetation, and soil predictors. We found the performance of machine learning and ensemble predictions to outperform traditional regression methods. We also found the predictive performance of NEE-focused models to be low, relative to models predicting GPP and ER. Our data compilation and ensemble predictions showed that $CO_2$ sink strength was larger in the boreal biome (observed and predicted average annual NEE –46 and –29 $g C m^{-2} yr^{-1}$, respectively) compared to tundra (average annual NEE +10 and –2 $g C m^{-2} yr^{-1}$). This pattern was associated with large spatial variability, reflecting local heterogeneity in soil organic carbon stocks, climate, and vegetation productivity. The terrestrial ecosystem $CO_2$ budget, estimated
Grant/Award Number: NNH17ZDA001N, NNX15AT81A and NNX17AE13G; Danmarks Grundforskningsfond, Grant/Award Number: CENPERM DNRF100; Nordic Center of Excellence; Arctic Data Center; EU FP7-ENV, Grant/Award Number: 282700; Natural Sciences and Engineering Research Council Discovery Grants; Suomen Akatemia, Grant/Award Number: 286950, 312912, 314630, 317054, 325680, 332196, 337549, 33761 and 337552; Greenland Research Council, Grant/Award Number: 80.35; Canada Research Chairs; US Geological Survey; Nordenskiöld-samfundet; Swedish National Space Board; NSF Research, Synthesis, and Knowledge Transfer in a Changing Arctic: Science Support for the Study of Environmental Arctic Change, Grant/Award Number: 1331083; NSF PLR Arctic System Science Research Networking Activities (Permafrost Carbon Network: Synthesizing Flux Observations for Benchmarking Model Projections of Permafrost Carbon Exchange), Grant/Award Number: 1931333; NSF grant, Grant/Award Number: 1203583, 1204263, 1702797, 1702798, DEB-1636476, PLR1504381, PLR1836898, AON 856864, 1304271, 0632264 and 1107892; EU 6th Framework Programme, Grant/Award Number: 036993; Danish National Research Foundation, Grant/Award Number: DNRF100; Research Council of Norway; Swedish Research Council, Grant/Award Number: contract #2018-03966; the national research infrastructures SITES and ICOS, funded by VR and partner institutes

using the annual NEE ensemble prediction, suggests the high-latitude region was on average an annual CO₂ sink during 1990–2015, although uncertainty remains high.

### KEYWORDS

Arctic, CO₂ balance, empirical, greenhouse gas, land, permafrost, remote sensing

# 1 | INTRODUCTION

The terrestrial ecosystem carbon dioxide (CO₂) balance is one of the largest uncertainties in the global carbon budget (Friedlingstein et al., 2020), with high latitudes (i.e., tundra and boreal biomes) representing one of the least-constrained budgets (López-Blanco et al., 2019; Schuur et al., 2015; Zscheischler et al., 2017). Moreover, due to polar amplification and large carbon stocks, the high latitudes have the potential for substantial positive feedbacks to climate warming (Abbott et al., 2016; Gasser et al., 2018; Schuur et al., 2008; Turetsky et al., 2020). Currently, in the absence of major disturbances (e.g., fire), boreal forests are generally CO₂ sinks (Bradshaw & Warkentin, 2015; Pan et al., 2011), while regional estimates of tundra vary from sinks (McGuire et al., 2009, 2012, 2016) to sources (Belshe et al., 2013). Both the growing and non-growing seasons are important for these annual budget estimates. A recent synthesis found that non-growing season soil CO₂ emissions from the northern permafrost region are larger than previously estimated (Natali et al., 2019). However, CO₂ uptake by plants over the growing season can be substantial and is often the dominant component of the annual CO₂ budget (Alekseychik et al., 2017; Kolari et al., 2009; Lafleur et al., 2012). The current state

of the annual terrestrial high-latitude CO₂ budget (net sink or source) remains highly uncertain. A key research priority is to develop robust data-driven quantitative frameworks to constrain regional boreal and tundra CO₂ budgets at annual and seasonal time scales.

Estimating high-latitude CO₂ fluxes across large areas and over long timescales is challenging due to their high spatiotemporal variability (Ai et al., 2018; Wilkman et al., 2018) that is controlled by a range of environmental variables (Camps-Valls et al., 2015; Lund et al., 2010). The ecosystem CO₂ balance (i.e., net ecosystem CO₂ exchange; NEE) is the relatively small difference between the two large CO₂ fluxes of photosynthesis (gross primary production; GPP) and ecosystem respiration (ER; comprising autotrophic and heterotrophic respiration). Although NEE can be measured with the eddy covariance (EC) and chamber techniques (Baldocchi et al., 1988; Lundegårdh, 1927), GPP and ER are estimated indirectly using environmental light and temperature measurements for EC sites (Lasslop et al., 2010; Reichstein et al., 2005) and dark chamber measurements for chamber sites (Shaver et al., 2007). Field studies have shown that GPP, ER, and NEE depend on climatic conditions (e.g., temperature, precipitation, and radiation) (López-Blanco et al., 2017; Nobrega & Grogan, 2008; Zhang et al., 2018), vegetation (Cahoon et al., 2012; Fox et al., 2008; Järveoja et al., 2018),
<figure><img src="image_2.png" /></figure>

and soil properties (e.g., soil nutrients and moisture) (Arens et al., 2008; Dagg & Lafleur, 2011; Lund et al., 2009). However, our understanding of the influence of these drivers on GPP and ER, and particularly on NEE, across the entire high-latitude region remains limited (see e.g., Belshe et al., 2013; Lund et al., 2010).

Knowledge of the contemporary high-latitude terrestrial CO₂ budget is further limited by an increasing, but still relatively sparse, flux measurement network (Alton, 2020; Chu et al., 2017; Virkkala et al., 2018). The majority of flux sites are concentrated within a few intensively studied regions, particularly Alaska and Fennoscandia (Metcalfe et al., 2018; Pastorello et al., 2020; Virkkala et al., 2019), with just a few sites in other large regions such as Siberia and northern Canada. Consequently, issues related to the temporal, geographical and environmental representativeness of the measurements need to be considered to accurately estimate high-latitude carbon budgets and their uncertainties. Previous studies have used a variety of synthesis approaches (Belshe et al., 2013; McGuire et al., 2012), and statistical (Natali et al., 2019), process-based (López-Blanco et al., 2019; McGuire et al., 2018; Rawlins et al., 2015; Wania et al., 2009) and atmospheric inversion models (McGuire et al., 2012), yielding highly different CO₂ budgets. Most of these modeling studies have been conducted at coarse spatial resolutions (25–100 km; Natali et al., 2019; Rawlins et al., 2015; López-Blanco et al., 2019) that do not fully capture the heterogeneity in high-latitude environments despite their importance for the regional CO₂ budgets (Raynolds et al., 2019; Treat et al., 2018). New efforts synthesizing the current distribution of flux data and developing models at high spatial resolution are required to improve our understanding on the spatial patterns and magnitudes of CO₂ fluxes.

Models that rely on the statistical relationships between CO₂ flux and predictor variables have been increasingly employed to constrain global and high-latitude CO₂ budgets (e.g., Jung et al., 2020; Natali et al., 2019; Warner et al., 2019). These statistical models are useful for predicting fluxes across larger areas (i.e., upscaling) because they directly draw upon relationships between fluxes and environmental variables, can account for environmental variability across space and time at high resolutions, and are able to handle biases in the geographic representation of the data (Jung et al., 2020; Natali et al., 2019; Warner et al., 2019). A broad range of statistical models and data sources are available for upscaling, but not all of these have been fully utilized. For example, many past studies have upscaled high-latitude fluxes using a single model (Natali et al., 2019; Peltola et al., 2019; Ueyama, Ichii, et al., 2013), but how different models compare with each other is not well known (with exception of Jung et al., 2017 and Tramontana et al., 2016). Further, most of these studies have primarily used machine learning models due to their ability to capture non-linear relationships and interactions in data (Elith et al., 2008). However, traditional regression methods can be a powerful tool in upscaling high-latitude ground conditions due to their ability to extrapolate beyond the range of data used for training, and due to their generalizability and ease of interpretation (Aalto et al., 2018). Finally, many of the recent upscaling studies have relied on EC flux measurements only, neglecting chamber

measurements despite their importance as additional data sources (with exception of Natali et al., 2019). Chambers are useful especially in remote, sparsely measured treeless tundra where they can capture the entire ecosystem CO₂ balance and directly measure NEE and ER (Sørensen et al., 2019). Thus, a compilation of both EC and chamber flux measurements and the comparison of available modeling techniques is clearly required to ensure accurate CO₂ flux estimates from existing data and models.

Here, we synthesize annual and growing season CO₂ fluxes from EC and chamber measurements across the high-latitude terrestrial tundra and boreal region. We then use this new database to upscale annual average ecosystem CO₂ fluxes at relatively high spatial resolution (1 km²) across the high-latitude domain using several statistical models. We compare our new database of *in situ* CO₂ fluxes to past tundra syntheses (Belshe et al., 2013; McGuire et al., 2012), provide a detailed assessment of model performance, analyze the spatial patterns and drivers of CO₂ fluxes, and discuss the resulting CO₂ budget estimates and recommendations for future work. We focus on understanding the spatial variability in average CO₂ fluxes instead of a temporal analysis of CO₂ flux change; however, our modeling framework also considers the interannual variability in fluxes.

## 2 | MATERIAL AND METHODS

### 2.1 | Data collection

#### 2.1.1 | Collection of CO₂ flux data

Our study area was defined by the high-latitude tundra and boreal biomes (>45°N) based on global ecoregions (20.6 × 10⁶ km²; Figure 1; Dinerstein et al., 2017). We first conducted a literature survey to identify existing EC and chamber-based terrestrial CO₂ flux observations of GPP, ER, and NEE over annual and growing season periods across the domain. Potential sites were identified from previous studies (Ichii et al., 2017; Marushchak et al., 2013; McCallum et al., 2013; Watts et al., 2014) and prior synthesis efforts (Belshe et al., 2013; McGuire et al., 2012; Virkkala et al., 2018). We augmented the resulting site list using a Web of Science search with key words ("tundra" or "boreal" or "arctic") and ("CO₂ flux" or "CO₂ exchange" or "CO₂ budget"). Additionally, a community call was solicited through a CO₂ flux synthesis workshop (Parmentier et al., 2019), whereby investigators contributed their most current unpublished data. Additional EC data were downloaded from FLUXNET2015 (Pastorello et al., 2020). The compiled dataset represents all natural terrestrial vegetation types (categorized by needle- or broadleaf forest, shrubland, grassland, wetland, and sparse vegetation) present in the high-latitude region.

We included studies and sites with NEE, GPP, and ER estimates over a full growing season or calendar year (i.e., cumulative flux). Growing season flux measurements are provided by EC and chambers. Non-growing season flux measurements include a variety of methods in addition to EC and chambers (e.g., a gas diffusion method by Björkman et al., 2010, soda lime by Welker et al., 2004, or an
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Measured median annual (a-c) and growing season (d-f) fluxes of GPP (gross primary production), ER (ecosystem respiration), and NEE (net ecosystem exchange) in the study domain (>45°N). The color of the point defines the median flux of the site (i.e., a sampling location), and the size of the point the number of observations that was measured (i.e., number of years). The background map represents the high-latitude region (dark gray = boreal biome, light gray = tundra biome). In all panels, sites that had only eddy covariance measurements are shown with black outline color around the point, and chamber measurements are without outline. One site had both eddy covariance and chamber measurements, but this is shown with black outline color. Positive numbers for NEE indicate net ecosystem CO₂ loss to the atmosphere (i.e., CO₂ source) and negative numbers indicate net ecosystem CO₂ gain (i.e., CO₂ sink)</figcaption></figure>

empirical model by Vogel et al., 2009). Growing season length and measurement period were defined in multiple ways at individual sites. To allow inter-site comparison, we filtered out measurements that did not represent the entire growing season and standardized the remaining measurements (see Supplementary Text Section 1.1 and a similar approach in Belshe et al., 2013). From this filtered dataset, we calculated average growing season daily flux rates based on the reported measurement length and standardized the fluxes based on a common growing season length. The final list of sites having

representative annual or growing season measurements is provided in Table S1, sites that were excluded from our analysis are in Table S2.

The resulting dataset included 148 sites with CO₂ fluxes from 1990 to 2015 from variable measurement periods (Figure 1). We compiled 1448 cumulative annual and growing season flux values (when chamber measurements were aggregated per site; Figure 1); 82% of the aggregated observations are from EC and 18% are from chambers. Annual and growing season NEE were the most widely reported fluxes in the dataset (Figure 1). Unlike McGuire et al. (2012)
and Belshe et al. (2013) we also included data from the boreal biome, additional tundra sites, and wetlands (not synthesized in Belshe et al., 2013; Figure S1). Similar to McGuire et al. (2012) and Belshe et al. (2013), our database primarily represents undisturbed environments. However, it also includes measurements from ca. 10 sites that have experienced high natural, anthropogenic or anthropogenically induced disturbances, such as permafrost thaw (Bäckstrand et al., 2010; Cassidy et al., 2016; Trucco et al., 2012), fires (Iwata et al., 2011; Ueyama et al., 2019), insect outbreaks (Heliasz et al., 2011; López-Blanco et al., 2017; Lund et al., 2017), or extensive harvesting practices (Coursole et al., 2012; Machimura et al., 2005). Throughout the text, positive numbers for NEE indicate net CO₂ loss to the atmosphere (i.e., CO₂ source) and negative numbers indicate net CO₂ gain (i.e., CO₂ sink). GPP and ER are always given as positive numbers.

### 2.1.2 | Gridded predictors and reference flux data

We acquired 10 eco-physiologically relevant predictors at 1-km² resolution (0.0083°) representing climate, vegetation, topographic, and soil conditions: growing degree days (GDD3; °C), freezing degree days (FDD; °C), water balance (WAB; mm), maximum growing season normalized difference vegetation index (NDVI), topographic wetness index (TWI), potential incoming direct annual solar radiation (RAD; MJ cm⁻² yr⁻¹), soil organic carbon stocks in the upper 2 meters (SOC; tons per ha), topsoil (0–5 cm) pH, topsoil clay content (CLAY; %), and land cover (LC; classes were mixed or broadleaved forest, needle-leaved forest, grassland and shrubland, wetland, sparse vegetation; see Supplementary Text Section 1.2 and Figure S2 for more information about the predictors). These predictors characterize previously identified key relationships between CO₂ fluxes and summer and winter temperatures, radiation, precipitation, local hydrology and soil conditions, soil carbon stocks, and vegetation properties (i.e., see Beer et al., 2010; Belshe et al., 2013; Lund et al., 2010; Natali et al., 2019; Ueyama, Iwata, et al., 2013). NDVI further reflects disturbances as it can show spectral browning signals related to drought, harvesting, or fires (Myers-Smith et al., 2020; Figure S3; Supplementary Text Section 2.5). We recognize that GPP and ER partitioning and gap filling rely on supporting environmental data (e.g., temperature and radiation), and consequently these fluxes already include some information about variables that we also used as predictors in our statistical models. We used annual (1990–2015) data for GDD3, FDD, WAB, and NDVI; the remaining predictors were considered to be static. All predictor datasets were masked to only include high-latitude tundra and boreal biomes (Dinerstein et al., 2017), and to exclude permanent water bodies, urban areas, and croplands based on a land cover dataset developed by ESA (2017).

We compared our annual ecosystem NEE predictions and budgets (see Section 2.2.1) with FLUXCOM, a global product derived from FLUXNET EC towers and machine learning at 0.5° resolution (Baldocchi et al., 2001; Jung et al., 2017; Tramontana

et al., 2016) and an ensemble of global Earth system models from the Coupled Model Intercomparison Project Phase 5 (CMIP5) at 1.92 × 1.5° resolution (Taylor et al., 2012) (Supplementary Text Section 1.3).

## 2.2 | Data analysis

### 2.2.1 | Statistical modeling

Our main response variables were annual and growing season cumulative GPP, ER, and NEE, but we also modeled daily average GPP, ER, and NEE during the growing season. Annual and growing season CO₂ fluxes were linked to the environmental predictors using a range of different statistical modeling methods (Figure S4). We used five statistical models; two were extensions of linear regression models, and three were based on machine-learning. All of these models have been widely used in empirical CO₂ flux upscaling studies (Bond-Lamberty & Thomson, 2010; Hursh et al., 2016; Tramontana et al., 2016; Ueyama, Ichii, et al., 2013). Specifically, we examined generalized linear models (GLMs); generalized additive models (GAMs); generalized boosted regression trees (GBMs); random forest (RF models); and support vector machines (SVMs).

We used several model approaches because individual models have inherent strengths and weaknesses (Supplementary Text Section 2). For example, machine learning methods might suffer from overfitting, whereas regression methods might result in unrealistic values when extrapolated outside the model data range. Further, individual models may detect different patterns in the data, and the best performing models are not always the same for different response variables (Segurado & Araújo, 2004). We also produced an ensemble prediction by calculating a median prediction over the five predictions from the individual modeling methods (see also Tramontana et al., 2016). We used the median instead of the mean to avoid extreme predicted values inflating the ensemble prediction. In this procedure, the uncertainty of the ensemble is expected to be lower than the uncertainty of a single model (Aalto et al., 2018). Consequently, we produced six model predictions for each of our response variables.

To determine the main drivers of the spatial patterns of response variables, the relative contribution of predictors in the models was assessed using a prediction re-shuffling approach (Niittynen & Luoto, 2018). We first fit the model and developed predictions using the original data, and then repeated this procedure with the values for one predictor randomly permuted. The contribution of a variable was calculated as a correlation between these two predictions (i.e., original model and the model with a shuffled predictor) subtracted from one:

$$ \text{Relative contribution} = 1 - \text{correlation} \\ \text{(Prediction}_{\text{original data}}, \text{Prediction}_{\text{Randomly permuted data}}) $$

Values close to 1 indicate that the two predictions were different, indicating high variable importance of the predictor variable.
Each predictor was randomly permuted 100 times for each flux with each of the modelling methods, and an ensemble contribu- tion was derived by taking a mean of the values. To visualize a predictor's effect on a response variable after controlling for the effects of other predictors, partial dependence plots were derived from the random forest model. For both variable importance and partial dependence plot analyses, we used daily average growing season fluxes because the growing season length estimates that were used to calculate growing season fluxes are not independent from GDD3. We found that the daily average fluxes correlated strongly with the growing season fluxes (Pearson's correlation 0.93–0.94), so they can be assumed to reflect the same relation- ships with the predictors.

To extrapolate across the study domain, we fit the models using
the entire dataset to produce annual flux predictions and their
ensembles that were subsequently averaged to 1990–2015 mean
values. Because the ensemble predictions were among the most
accurate and least uncertain predictions across all response vari-
ables, and because their use is generally recommended in predic-
tive efforts (Araújo & New, 2007), our final flux maps and budgets
were based on the flux ensemble. In addition to annual and growing
season budgets, we also calculated a non-growing season budget
(see Table S4). We had different numbers of observations and sites
available for each flux and model, and consequently observed and
predicted ER and GPP fluxes and budgets do not sum up to NEE.

2.2.2 | Model fit, predictive performance and uncertainty

To evaluate model fit, we predicted fluxes over the entire model
training data. To assess the predictive performance of the models,
we used a leave-one-site-out cross validation scheme in which each
site was iteratively left out from the dataset, and the remaining
data were used to predict fluxes for the excluded site (Bodesheim
et al., 2018). For both model fit and predictive performance, we
calculated bias as an average of the absolute error between predic-
tion and actual observations, Pearson correlation (r) to determine
the strength of the linear relationship between the observed and
predicted fluxes, and root mean squared error (RMSE) to estimate
the model error. We use the terms "observed" and "predicted" to
distinguish between field measurements and model predictions but
acknowledge that some of these observed values represent indirect
estimates of fluxes (e.g., GPP).

We evaluated the prediction uncertainty of all flux models and
the budget uncertainty of annual and growing season NEE mod-
els using a repeated random resampling procedure (Aalto et al.,
2018). Prediction uncertainty was calculated to characterize the
spatial variability in flux predictions across the high-latitude region,
whereas budget uncertainty quantified the range of potential NEE
budget values. We used bootstrapping (fractional resampling with
replacement based on LC classes) to subset the model training data
into 200 different datasets, all of which had the same number of

observations as the original flux data itself. These 200 datasets were
then used to produce 200 individual predictions with all five sta-
tistical models and their ensemble for each flux and for each year
from 1990 to 2015 to assess prediction uncertainty which was
summarized using the prediction interval (PI; 95th percentile – 5th
percentile). Uncertainty for annual and growing season NEE budgets
was estimated by calculating the range of budgets from the 50 first
ensemble predictions out of the 200 predictions for each year from
1990 to 2015, due to computational constraints. For more details,
see Supplementary Text Section 2.4 and Figure S5.

3 | RESULTS

3.1 | Observed flux variation

Flux measurements showed considerable variation in magnitudes and signs (CO₂ sink vs source) across the high-latitude environments (Figure 1 and Table 1). Observed annual NEE (no upscaling) was on average a small source of CO₂ in the most northern parts of the study domain (tundra: +10 g C m⁻² yr⁻¹, 42 sites; northern permafrost region: +6 g C m⁻² yr⁻¹, 63 sites) and in drier environments (tundra upland: +16 g C m⁻² yr⁻¹, 36 sites), whereas the boreal biome (−46 g C m⁻² yr⁻¹, 41 sites), and in particular boreal uplands (−47 C m⁻² yr⁻¹, 34 sites), and non-permafrost regions (−90 g C m⁻² yr⁻¹, 20 sites) were net ecosystem CO₂ sinks. All environmental categories were, on average, net CO₂ sinks during the growing season, with the average NEE ranging from −37 to −115 g C m⁻² period⁻¹ (Table 1). Tundra upland and non-permafrost regions had the lowest average growing season sink strength. The non-permafrost region sink was greatly reduced by one disturbed site that had large source values up to +600 g C m⁻² period⁻¹ (Petrone et al., 2014), but this was not apparent in the annual averages because the same site did not report annual fluxes. Although the environmental conditions at the sites were fairly representative of the entire high-latitude region (Figure S6), colder environments with low NDVI and GDD3 as well as high FDD were less well represented (e.g., large areas of Siberia; Figure 1). Some chamber sites were located in conditions that would have otherwise remained undersampled (Figure S6). These included sites with relatively high soil organic carbon stocks in Hudson Bay Lowland and northwestern Canada, and wet climates in Greenland and northern Fennoscandia.

3.2 | Predictive performance of the models

The model fit and predictive performance analyses indicated that
the GBM, RF and SVM (machine learning) methods outperformed
the GLM and GAM (regression model) approaches across most of
the response variables (in particular with NEE, but also with GPP
and ER; model fit of annual machine learning models: *r* = 0.69–0.99
vs. regression models: *r* = 0.6–0.92; predictive performance of an-
nual machine learning methods: *r* = 0.2–0.73 vs. regression models:
<table><caption>TABLE 1 Summary statistics of observed and predicted (using the average ensemble prediction) annual and growing season GPP (gross primary productivity), ER (ecosystem respiration), and NEE (net ecosystem exchange) fluxes (g C m⁻² yr⁻¹ for annual and g C m⁻² period⁻¹ for growing season fluxes) in different environments across the high-latitude region over 1990–2015. The time-series of the sites were averaged prior to calculating the observed mean flux (i.e., one flux value from one site was used when the regional averages were calculated). Positive numbers for NEE indicate net CO₂ loss to the atmosphere (i.e., CO₂ source) and negative numbers indicate net CO₂ gain (i.e., CO₂ sink). Note that ER and GPP do not sum up to NEE as different numbers of observations and sites were available for each flux and model. Moreover, some plant uptake occurs outside of our defined growing season, and consequently growing season GPP and annual GPP do not equal to each other. The average fluxes were calculated based on the extent of the high-latitude tundra and boreal biomes (Dinerstein et al., 2017), permafrost zones (Brown et al., 2002), and land cover (i.e., wetlands, and everything else is upland; ESA, 2017). The confidence intervals for the observed fluxes and the uncertainty ranges for the predicted fluxes can be found in the Table S6</caption><thead><tr><th>Category</th><th>Annual GPP</th><th>Annual ER</th><th>Annual NEE</th><th>Growing season GPP</th><th>Growing season ER</th><th>Growing season NEE</th></tr></thead><tbody><tr><td><strong>Observed mean flux</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>High-latitude</td><td>482</td><td>456</td><td>-17</td><td>317</td><td>262</td><td>-63</td></tr><tr><td>Boreal</td><td>624</td><td>605</td><td>-46</td><td>420</td><td>347</td><td>-87</td></tr><tr><td>Tundra</td><td>250</td><td>259</td><td>10</td><td>232</td><td>192</td><td>-44</td></tr><tr><td>Boreal upland</td><td>676</td><td>647</td><td>-47</td><td>432</td><td>350</td><td>-84</td></tr><tr><td>Boreal wetland</td><td>406</td><td>381</td><td>-38</td><td>347</td><td>330</td><td>-102</td></tr><tr><td>Tundra upland</td><td>250</td><td>259</td><td>16</td><td>232</td><td>192</td><td>-37</td></tr><tr><td>Tundra wetland</td><td></td><td></td><td>-24</td><td></td><td></td><td>-115</td></tr><tr><td>No permafrost</td><td>831</td><td>773</td><td>-90</td><td>405</td><td>370</td><td>-37</td></tr><tr><td>Permafrost</td><td>342</td><td>350</td><td>6</td><td>302</td><td>241</td><td>-67</td></tr><tr><td><strong>Predicted mean flux</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>High-latitude</td><td>554</td><td>508</td><td>-20</td><td>343</td><td>283</td><td>-50</td></tr><tr><td>Boreal</td><td>638</td><td>594</td><td>-29</td><td>396</td><td>327</td><td>-52</td></tr><tr><td>Tundra</td><td>378</td><td>326</td><td>-2</td><td>230</td><td>192</td><td>-46</td></tr><tr><td>Boreal upland</td><td>653</td><td>604</td><td>-30</td><td>399</td><td>328</td><td>-51</td></tr><tr><td>Boreal wetland</td><td>437</td><td>458</td><td>-18</td><td>358</td><td>303</td><td>-64</td></tr><tr><td>Tundra upland</td><td>378</td><td>326</td><td>-1</td><td>229</td><td>191</td><td>-45</td></tr><tr><td>Tundra wetland</td><td>367</td><td>347</td><td>-29</td><td>281</td><td>242</td><td>-71</td></tr><tr><td>No permafrost</td><td>805</td><td>736</td><td>-56</td><td>447</td><td>375</td><td>-53</td></tr><tr><td>Permafrost</td><td>489</td><td>448</td><td>-11</td><td>315</td><td>259</td><td>-49</td></tr></tbody></table>

r = 0.12–0.72; Figure 2g-i). We found that the machine learning-based methods were less uncertain (Figure S7) and predicted values within the range of the observed fluxes as opposed to regression models. However, the machine learning method that performed best and had the least uncertainties varied depending on the flux response variable.

Ensemble predictions were among the best performing models (model fit of annual and growing season ensemble models: *r* = 0.68–0.94; predictive performance of annual and growing season ensemble models: *r* = 0.21–0.73; Figure 2 and Figure S8). However, similar to the individual models, model fit and predictive performance was lower for annual and growing season NEE compared to GPP and ER (model fit for GPP and ER: *r* = 0.89–0.94 vs. NEE: *r* = 0.68–0.77; predictive performance for GPP and ER: *r* = 0.53–0.71 vs. NEE: *r* = 0.21–0.27; Figure 2 and Figure S8). Annual models for ER and NEE exhibited a better fit and predictive performance than the growing season models (based on *r*), whereas the opposite was true for GPP (Figure 2 and Figure S8). The growing season GPP model fit and predictive performance was higher than that of the ER models,

but annual GPP and ER models performed equally well. Model fit and predictive performance were similar in models trained with and without chambers (Table S3). In most predictive performance analyses, the lowest and highest observed fluxes were over- and underestimated, respectively, indicating overall poor predictive performance at the extremes (Figures S9 and S10).

Average predicted and observed fluxes were of similar magnitude (Table 1). However, there was a tendency for the average predicted values to have slightly larger GPP and ER values (e.g., observed and predicted annual GPP in the tundra: 250 g C m⁻² yr⁻¹ and 378 g C m⁻² yr⁻¹, respectively) and stronger net CO₂ sink values than what was observed (e.g., observed and predicted annual NEE in the tundra: +10 g C m⁻² yr⁻¹ and -2 g C m⁻² yr⁻¹, respectively). Our cross-comparison of annual and growing season flux ensemble predictions showed there was a mismatch between annual and growing season component fluxes in approximately 2% of the pixels (growing season GPP/ER larger than annual GPP/ER) and that unrealistic flux values (negative GPP or ER) were found in less than 0.01% of the pixels in the ensemble predictions.
<figure><img src="image_4.png" /></figure>

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

FIGURE 2 Observed and predicted annual fluxes of GPP (gross primary production; a and d), ER (ecosystem respiration; b and e), and NEE (net ecosystem exchange; c and f) based on model fit (a–c) and predictive performance (d–e). Model fit was evaluated by predicting fluxes over the entire model training data, while predictive performance was assessed using a leave-one-site-out cross validation scheme in which each site was iteratively left out from the dataset, and the remaining data were used to predict fluxes for the excluded site. Model fit and predictive performance statistics (r = Pearson's correlation between observed and predicted fluxes, g; Bias = mean absolute bias, h; RMSE = root mean square error, i across annual fluxes and five modeling methods; GLM = generalized linear model; GAM = generalized additive model; GBM = generalized boosted regression tree; RF = random forest; SVM = support vector machine) and their median ensemble (ENS) are shown in subfigures g–i. The black line indicates a 1:1 relationship

3.3 | Predicted flux variation

Predicted fluxes showed high spatial variability across the region
with a general trend towards decreasing fluxes and sink strength
with increasing latitude for GPP, ER, and NEE (Figure 3 and Figure
S11). The variability was related to differences in climate (GDD3
and FDD), solar radiation (RAD) and vegetation greenness (NDVI),
which had the strongest influence on most of the fluxes (Figure 4).
Moreover, SOC, CLAY, and LC were important variables for annual
NEE; CLAY and SOC both had a positive yet saturating relationship
(Figure S12). The relationship between LC and NEE suggested that

the annual and growing season net sink strength was largest in wet-
lands and smallest in sparse vegetation (Figures S12 and S13). Some
variables had a very low variable importance for most of the fluxes
(e.g., TWI, soil pH).

Our predictions revealed regional hot spots in annual and grow-
ing season NEE, GPP, and ER. Strong annual and growing season CO2
sinks, having low ER and high GPP, were found in forested regions
with high GDD3, NDVI, RAD, and low FDD across Fennoscandia and
European Russia, southern Canada, and southern Siberia (Figure 3
and Figure S11). Annual CO2 sources were identified within northern
and central Siberia, Greenland, northern and central Alaska, as well
<figure><img src="image_6.png" /></figure>

FIGURE 3 Average predictions of annual CO₂ fluxes at 1-km² resolution over 1990–2015. Annual predictions (a–c), associated uncertainties (d–f) and mean fluxes and uncertainties along latitudes (g–i) of GPP (gross primary production), ER (ecosystem respiration), and NEE (net ecosystem exchange) of the statistical model ensembles over 1990–2015. The uncertainty (prediction interval, PI; 90% uncertainty range) is quantified as the variability of predictions over a random subset of pixels (n = 10,000) interpolated across the study domain based on a repeated (n = 200) bootstrap sampling procedure. It demonstrates how robust the relationships in the models are and how differences in model training data influence the predictions. The gray lines in a–f represent the borders of northern countries and points in g–i site locations

as northern Canada. These regions were located mainly in the tundra, characterized by high FDD, and low GDD3 and NDVI. Growing season CO₂ sources were located in southeastern Siberia, northern Siberia and some parts of southern and northern Canada. Largest uncertainties in flux predictions were found in areas with relatively strong CO₂ sinks in the boreal biome, such as in Fennoscandia and eastern Canada, but also in the tundra (e.g., Canadian Arctic Archipelago; Figure 3 and Figure S11). The largest differences across our annual NEE, and CMIP5 and FLUXCOM predictions were found in European Russia, Fennoscandia, and southeastern Canada (Figure 5a–d).

## 3.4 | Terrestrial ecosystem NEE budget for the high-latitude region

Our ensemble predictions showed that the high-latitude tundra and boreal region was on average an annual terrestrial ecosystem CO₂ sink over the 26-year (1990–2015) study period (Table 2). The annual NEE budget (based on upscaled NEE data) averaged –419 Tg C yr⁻¹ (90% uncertainty range: –559 to –189 Tg C yr⁻¹; range of budgets across the study period: –449 to –366 Tg C yr⁻¹). When estimating annual NEE according to the separately modeled annual GPP

(11,344 Tg C yr⁻¹) and ER (10,397 Tg C yr⁻¹) budgets, we obtain an NEE budget of –948 Tg C yr⁻¹. The average high-latitude growing season NEE budget over the period of 1990–2015 was –1018 Tg C yr⁻¹ (–1332 to –455 Tg C yr⁻¹, 90% uncertainty range), which was supported by the difference between the average growing season ER (5800 Tg C yr⁻¹) and GPP (7016 Tg C yr⁻¹) budgets. For the regional budgets, see Table 2.

The average annual NEE budgets over the study period from CMIP5 and FLUXCOM were –488 and –1056 Tg C yr⁻¹, respectively (Table S5). In the boreal biome, average annual GPP in our study was 8850 compared to 8561 Tg C yr⁻¹ in FLUXCOM. In the tundra biome, the average annual GPP in this study was twice as high as in FLUXCOM (2495 and 1267 Tg C yr⁻¹, respectively). Differences were larger for annual ER. Our annual ER budget for the boreal and tundra biomes was 8241 and 2156 Tg C yr⁻¹, respectively, but the same budgets were only 6363 and 1200 Tg C yr⁻¹ in FLUXCOM. For the regional NEE budgets estimated with CMIP5 and FLUXCOM, see Table S5.

## 4 | DISCUSSION

This study provides a conceptual and methodological framework to bridge the gap between local, regional, and high-latitude scales in
<figure><img src="image_7.png" /></figure>

FIGURE 4 Variable importance for annual and growing season fluxes of GPP (gross primary production), ER (ecosystem respiration), and NEE (net ecosystem exchange). Explanatory variables are GDD3 (growing degree days), FDD (freezing degree days), WAB (water balance), NDVI (normalized difference vegetation index), TWI (topographic wetness index), RAD (potential incoming direct annual radiation), SOC (soil organic carbon stocks up to 2 m), pH (topsoil pH), CLAY (topsoil clay content), and LC (land cover). Variable importance was calculated by assessing how a randomly permuted predictor influences the predictions across all five statistical models. Values close to 0 and 1 indicate low and high importance of the predictor variable, respectively. The box corresponds to the 25th and 75th percentiles. The lines denote the 1.5 IQR of the lower and higher quartile, where IQR is the inter-quartile range, or distance between the first and third quartiles

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

statistical flux upscaling. Our framework is unique in that it (a) com-
piles a new dataset of growing season and annual fluxes using EC and
chamber data and investigates the drivers of these fluxes; (b) quanti-
fies the performance of different statistical models; and (c) provides
the first spatially continuous high-latitude maps of CO₂ fluxes and
their uncertainties at high spatial resolution, capturing the inherent
spatial heterogeneity in predictors and fluxes and minimizing biases
in upscaling compared to coarser scale models (Figure 5e). The bet-
ter geographical and environmental coverage of the flux measure-
ments compared to past efforts improves our understanding of the
spatial patterns and regional budgets of terrestrial ecosystem CO₂
fluxes, however, uncertainties in our direct model estimates of NEE
remained rather high.

## 4.1 | Drivers and spatial patterns of GPP, ER, and NEE

Our results suggest that climatic, vegetation, and soil variables
were all important predictors for terrestrial ecosystem CO₂ fluxes.
However, almost all CO₂ fluxes were strongly driven by the broad

climatic gradients and spatiotemporal variability in solar radiation,
growing and non-growing season climatic conditions, water balance,
and the resulting vegetation greenness patterns, supporting the find-
ings of previous syntheses (Belshe et al., 2013; Lund et al., 2010;
Natali et al., 2019). Even though these climatic variables are not in-
dependent of our GPP and ER estimates (see Section 4.2), confidence
in these results can be drawn from the underlying mechanistic rela-
tionships between the climate drivers and fluxes. For example, GPP
across large scales is dependent on growing season temperatures,
length of season, and radiation, which regulate and provide resources
for plant growth (López-Blanco et al., 2017; Lund et al., 2010), and
ER is largely driven by enzymatic processes, which are tightly linked
with temperatures (Davidson et al., 2006) as well as plant growth (La
Puma et al., 2007). In general, we found that warmer, moderately
wet, and greener conditions (i.e., environments of higher biomass as
indicated by NDVI) increased the magnitude of annual GPP and ER.
However, our results also indicate that the overall net sink strength
increases with larger greenness, warmer and shorter winters, and
wetter climate. These results suggest that GPP and ER respond rather
similarly to changes in climate and vegetation conditions across the
high-latitude region, although GPP might increase even more due to
<figure><img src="image_9.png" /><figcaption>FIGURE 5 Complementing annual NEE predictions averaged over 1990–2015. Mean annual NEE derived by subtracting annual ER (ecosystem respiration) from GPP (gross primary production) in this study (a), from a global upscaling product FLUXCOM (b), and from a process model ensemble CMIP5 (Coupled Model Intercomparison Project Phase 5) (c), and the standard deviation of these and the independently modeled annual NEE in this study (visualized in Figure 3c) (d). A regional-scale example of the spatial variation of annual NEE in our prediction in northern Alaska, with black outlines depicting the size of the pixel in one of the highest resolution (smallest pixel size) models in the CMIP5 ensemble (1.92 × 1.5°) (e)</figcaption></figure>

<table><caption>TABLE 2 Annual and growing season average GPP, ER, and NEE budgets (Tg C yr⁻¹) over 1990–2015 across the environments and the spatial extent of each environmental category when permanent water bodies, urban areas, and croplands were masked away. The NEE budgets are based on upscaled NEE data and include an uncertainty range derived by bootstrapping. Note that ER and GPP do not sum up to NEE as different numbers of observations and sites were available for each flux and model. For the non-growing season CO₂ budgets estimated based on annual and growing season budgets, see Table S4</caption><thead><tr><th>Category</th><th>Annual GPP</th><th>Annual ER</th><th>Annual NEE</th><th>Growing season GPP</th><th>Growing season ER</th><th>Growing season NEE</th><th>Area ×10<sup>6</sup> km<sup>2</sup></th></tr></thead><tbody><tr><td>High-latitude</td><td>11,344</td><td>10,397</td><td>-419 (-559 to -189)</td><td>7016</td><td>5800</td><td>-1018 (-1332 to -455)</td><td>20.6</td></tr><tr><td>Boreal</td><td>8850</td><td>8241</td><td>-406 (-499 to -239)</td><td>5496</td><td>4531</td><td>-715 (-1037 to -224)</td><td>13.9</td></tr><tr><td>Tundra</td><td>2495</td><td>2156</td><td>-13 (-81 to 62)</td><td>1520</td><td>1269</td><td>-303 (-338 to -224)</td><td>6.7</td></tr><tr><td>Boreal upland</td><td>8437</td><td>7808</td><td>-389 (-475 to -226)</td><td>5158</td><td>4245</td><td>-655 (-973 to -196)</td><td>12.9</td></tr><tr><td>Boreal wetland</td><td>412</td><td>433</td><td>-17 (-28 to -10)</td><td>338</td><td>287</td><td>-60 (-70 to -29)</td><td>0.9</td></tr><tr><td>Tundra upland</td><td>2451</td><td>2115</td><td>-9 (-78 to 64)</td><td>1486</td><td>1240</td><td>-294 (-330 to -218)</td><td>6.6</td></tr><tr><td>Tundra wetland</td><td>44</td><td>41</td><td>-4 (-3 to -1)</td><td>34</td><td>29</td><td>-8 (-9 to -6)</td><td>0.1</td></tr><tr><td>No permafrost</td><td>3407</td><td>3116</td><td>-238 (-288 to -185)</td><td>1895</td><td>1587</td><td>-223 (-353 to -45)</td><td>4.2</td></tr><tr><td>Permafrost</td><td>7924</td><td>7269</td><td>-181 (-305 to 32)</td><td>5114</td><td>4207</td><td>-793 (-1000 to -414)</td><td>16.3</td></tr></tbody></table>

increases in vegetation greenness (Berner et al., 2020) and changing climate (Lund et al., 2010). However, differences in these relationships might occur in different regions and land cover types (Baldocchi et al., 2018; Belshe et al., 2013; Lafleur et al., 2012).

In addition to the climate and greenness variables operating mostly at large scales, other more local-scale variables such as soil organic carbon stock and land cover helped explain CO₂ fluxes. Soil organic carbon stock was the most important predictor for annual NEE, and it had a positive relationship with it, demonstrating that

areas with high carbon stocks might lose more CO₂ to the atmosphere. However, this result was not supported by the annual ER models, which would represent the main process behind this positive relationship (i.e., larger carbon stocks have more potential for increased CO₂ emissions, particularly in dry conditions (Voigt et al., 2019)). The lack of this relationship might be due to annual ER models not covering the full range of conditions represented by the annual NEE models, or spurious causal relationships being identified by the relatively poorly performing NEE models. The importance of land
<figure><img src="image_10.png" /></figure>

cover was expected as it summarizes many key processes related to carbon cycling (e.g., the carbon uptake capacity, temperature sensitivity, as well as quantity and quality of carbon inputs into the soil; Sørensen et al., 2019) and other environmental characteristics (e.g., soil moisture is likely higher in wetlands than in sparse vegetation).

Our ensemble prediction suggested that most of the southern portion of the high-latitude terrestrial region was an annual net ecosystem CO₂ sink while the central and northern regions were neutral or small net CO₂ sources. Observed and predicted spatial patterns in fluxes were similar to those described by most previous studies. For example, our compiled field observations and predictions are consistent with the majority of Alaskan tundra being an annual ecosystem CO₂ source on average, similar to the average observed fluxes in McGuire et al. (2012) or the prediction in Ueyama, Ichii, et al. (2013). The strongest annual ecosystem CO₂ sinks in our study were located in southern European Russia, Fennoscandia, and southern Canada, as also observed in the FLUXCOM product (Jung et al., 2017; Tramontana et al., 2016).

For some regions, our ensemble prediction differed from the predictions of previous studies. The distribution of annual net CO₂ sources across the tundra biome was larger in our prediction compared to FLUXCOM, particularly in Siberia and Canada. This was likely explained by our models including some tundra sites from Canada, Greenland, European Russia, and Siberia, which were not covered by the FLUXCOM model training data. Some of the sites in these regions were annual net CO₂ sources on some years (Emmerton et al., 2016; Karelin et al., 2013). A similar disagreement was found between an Asia-focused statistical upscaling analysis by Ichii et al. (2017) which suggested stronger sink strength across large parts of Siberia, likely due to a limited number of northern eddy covariance sites used to train their models. The largest regional differences between our predictions, CMIP5, and FLUXCOM occurred in central Siberia, Fennoscandia, European Russia, and eastern Canada and the Canadian Arctic Archipelago, and these differences were primarily driven by the fact that CMIP5 showed these regions to be sources whereas they were sinks in FLUXCOM and our analysis (Figure 5). These regional differences demonstrate that these particular areas should be studied further to understand the sink-source patterns more accurately in the future.

Our uncertainty estimation suggests that CO₂ flux predictions should be interpreted carefully in areas that lack sampling locations or have large variability in fluxes that cannot be captured by the predictor variables. Such areas are particularly concentrated in European Russia, eastern Canada, and the Canadian Arctic Archipelago. As the accuracy of the prediction can usually be improved with increases in the quantity and quality of data, new measurements in these regions or better predictors would likely improve the performance of high-latitude CO₂ flux models.

## 4.2 | Key sources of uncertainty in our modeling approach

No single best model could be identified across the five modeling methods. However, the three machine learning methods

outperformed the two regression models, particularly for NEE, as demonstrated by the improved model performance, lower uncertainty, and the lack of unrealistically high or low flux values in predictions. The better performance of the machine learning methods was likely related to their flexibility and capability to find complex structures in the flux data (Elith et al., 2008). Our results demonstrate that several machine learning methods should be tested to produce the most accurate high-latitude flux predictions and that ensemble methods provide robust predictions (Araújo & New, 2007). Our results also indicate that an ensemble prediction based on machine learning methods alone would likely lead to higher model accuracy and transferability (see also Tramontana et al., 2016).

Our models performed well when predicting to the same data that the models were trained with, but the models had challenges when tested against independent validation data. The predictive performance of our ensemble predictions was comparable to (annual GPP and ER) or less than (growing season GPP, ER, NEE, and annual NEE) that of in other global and regional upscaling studies (Ichii et al., 2017; Natali et al., 2019; Peltola et al., 2019; Tramontana et al., 2016; Ueyama, Ichii, et al., 2013). However, comparisons of cross-validation results are hampered by different cross-validation techniques used in studies, with some of the studies including observations from the same site both in the model training and validation data, therefore providing overly optimistic accuracy estimates based on non-independent data. Moreover, these other studies primarily focused on a smaller area and/or shorter time period (with the exception of Tramontana et al., 2016), and used a different set of predictors, further complicating this comparison. In these other studies, the correlation (r) between observed and predicted fluxes (derived with cross validation), measured mostly throughout the year as daily-to-monthly fluxes, was roughly 0.65–0.7 for NEE and 0.7–0.8 for GPP and ER. There are several reasons for why some of our models performed more poorly than these previous studies, which we explain below.

The lower quantity of measurements and weaker comparability of fluxes derived with EC and chamber techniques and with variable measurement lengths might explain the lower predictive performance in our study compared to the other upscaling studies. As we used aggregated fluxes over the growing season and annual time scales, the sample size in our models was smaller than in other efforts which all used daily or monthly fluxes (a few hundred observations versus thousands of observations). A larger sample size usually increases the predictive performance of the models, particularly when these measurements cover variable environmental conditions that can be captured by the predictors. For example, FLUXCOM models (Jung et al., 2017, 2020; Tramontana et al., 2016) might have had a higher predictive performance than our models because they use a global FLUXNET database (Pastorello et al., 2020), which covers broad environmental gradients. However, FLUXNET data originates mostly from lower latitudes (e.g., only five sites from the Arctic and 34 from the boreal out of 224 global sites in total used in Tramontana et al., 2016). This could explain the larger net sink strength in FLUXCOM compared to our predictions. The higher
predictive performance of FLUXCOM compared to our prediction
might also be explained by the fact that FLUXNET is based on a
single flux measurement technique (EC) with standardized filtering,
gap-filling, and partitioning procedures. We included chambers to
our analysis as they covered conditions that were not covered by
the EC network even though we acknowledge that using both cham-
ber and EC measurements, and different partitioning methods for
EC, increased the number of different flux measurement techniques
and study designs, and may have made the comparison of fluxes
across sites more uncertain (Fox et al., 2008; Tramontana et al.,
2016). However, we observed no significant differences in fluxes
estimated with the two approaches indicative of these mismatches
(Figure S6d), and the performance of models did not change when
chambers were excluded from model training data. These results
suggest that the relatively low performance of some models is re-
lated to the high variability in both EC and chamber-derived CO2
flux estimates that is not captured by our predictors. Further, it
demonstrates that including chamber measurements, despite op-
erating at different spatial and temporal resolutions than the EC
technique, did not decrease the model performance. It is also possi-
ble that the lower predictive performance of growing season models
compared to annual models was related to the variable growing sea-
son measurement periods used across the studies. We accepted this
variability because our goal was to use as many published fluxes as
possible to improve the geographical and environmental coverage
of sites.

The accuracy of our ensemble predictions varied depending on
the flux, with the predictive performance being lowest for NEE mod-
els (r = 0.21–0.27). The predictive performance of our GPP and ER
models was higher (r = 0.53–0.73) and is comparable to past efforts
(Ichii et al., 2017; Natali et al., 2019; Tramontana et al., 2016; Ueyama,
Ichii, et al., 2013) because these fluxes represent the ecophysiologi-
cal and biogeochemical processes describing CO2 uptake and loss, re-
spectively. GPP and ER also already included some information about
temperature and radiation variables that we used as predictors in our
statistical models, which may introduce some circularity and artificially
inflate the model performance. Our NEE models over- and underesti-
mated low and high (i.e., large negative and positive) values, respec-
tively, by approximately 100–200 g C m⁻² yr⁻¹, which has also been
demonstrated with NEE and other fluxes in previous upscaling studies
(Ichii et al., 2017; Tramontana et al., 2016; Warner et al., 2019). These
extreme values were often from disturbed sites experiencing for ex-
ample, permafrost thaw or extreme forest management practices, or
represented an observation that was notably different from the site
mean. Based on the cross validation results of the individually modeled
annual NEE, a substantial fraction (53%) of annual source observations
were predicted to be sinks (similar to the pattern observed in Ichii et al.,
2017; Figure 3b), but some sink observations (24%) were also predicted
as sources. We also discovered that the observed average annual NEE
was often larger (more positive) than the individually predicted aver-
age NEE, which was either a result of the model not being able to pre-
dict sources accurately, or of the distribution of flux sites being biased
towards environments with larger CO2 source observations than the

entire region on average (see the large number of sites with source ob-
servations originating primarily only from Alaska in Figure 1). These re-
sults demonstrate that the predictors included in our analyses did not
fully represent the spatial gradients and dynamic temporal variability
in environmental conditions that influence carbon cycle processes, and
particularly the high and low NEE conditions. Further research should
explore improvements offered by other current and potential future
predictors related to the disturbance and permafrost conditions, snow
cover duration and snow depth, soil moisture and nutrient availabil-
ity, and phenology, root properties, and microbial communities (Illeris
et al., 2003; Järveoja et al., 2018; Nobrega & Grogan, 2007).

Even though the geographical and environmental coverage of the flux sites was improved in our study compared to previous efforts, our models included only ca. 10 sites from heavily disturbed conditions (see Section 2.1.1). Consequently, our sites did not cover the full range of disturbance and post-disturbance recovery conditions and the associated impacts on CO2 fluxes. For example, rapidly thawing permafrost and burned landscapes remained largely under-sampled across Siberia. These disturbances have a substantial impact on carbon cycling in high-latitude ecosystems (Abbott et al., 2016; Walker et al., 2019), including direct emissions from the disturbance (not estimated with our models) and typically increased net CO2 emissions for several years to decades after the disturbance (Courolle et al., 2012; Kittler et al., 2017; Lund et al., 2017; Turetsky et al., 2020) which should ideally be captured by our models. The lack of flux data representing disturbed and post-disturbance recovery conditions likely leads to underestimations in net ecosystem CO2 emissions, and is generally thought as one of the key limitations in statistical upscaling efforts (Jung et al., 2020; Zscheischler et al., 2017).

## 4.3 | Terrestrial ecosystem CO₂ budget and its uncertainty

Although our models may be biased towards sinks, our results sug-
gest that high-latitude terrestrial ecosystems were on average an an-
nual net CO2 sink during 1990–2015. The uncertainty of this budget
was high, as demonstrated by the low predictive performance of the
annual NEE model, and the fact that budgets derived from differ-
ent predictions (individual NEE predictions and ER-GPP predictions)
differed by ca. 500 Tg C yr⁻¹ – the latter most likely being linked to
the different numbers of observations and sites available for each
flux and model (Figure 1). Nevertheless, the annual NEE budget was
of similar magnitude to the one estimated by CMIP5 models and
larger (less negative) than the one estimated by FLUXCOM (Table
S5). The boreal biome was responsible for most of this sink strength
(–406 Tg C yr⁻¹, from –499 to –239 Tg C yr⁻¹; 13.9 × 10⁶ km²). In
contrast, the tundra biome was on average a small sink (–13 Tg C
yr⁻¹, from –81 to +62 Tg C yr⁻¹; 6.7 × 10⁶ km²) or a small source
(+10 g C m⁻² yr⁻¹), based on our predictions and observations. This
suggests that the tundra biome was on average close to CO2 neu-
tral even though the large soil organic carbon stocks of this region
would indicate larger historical CO2 sink strength (Hugelius et al.,
2014). Our tundra budget is within the range (though on average more positive, indicating stronger source) of the one comprising process and inversion models, and field-based estimates by McGuire et al. (2012) (–103 Tg C yr⁻¹, from –297 to +89 Tg C yr⁻¹). However, it differs from the source budget (+462 Tg C yr⁻¹, from +94 to +840 Tg C yr⁻¹; 10.5 × 10⁶ km²; wetlands not included) estimated by Belshe et al. (2013). The divergence of average annual NEE across our and Belshe et al. (2013) study is likely explained by our inclusion of fluxes from wetlands, which were on average annual net ecosystem CO₂ sinks (Table 1). The discrepancy between our and the McGuire et al. (2012) study can be explained by a 50% increase in new annual tundra source observations in our dataset (see e.g., Celis et al., 2017; Euskirchen et al., 2014), which were not included in the McGuire et al. (2012) analysis. Further, there are some differences in the study domain boundaries (e.g., the tundra domain in Belshe et al., 2013 was larger than in this study) which might explain some of the discrepancies between these studies, although the general patterns of these boundaries were rather similar (see e.g., Figure 1 in McGuire et al., 2012 vs. our tundra domain in Figure 1).

Our findings suggest that both the boreal and tundra biomes were relatively strong CO₂ sinks during the growing season. Our growing season CO₂ budgets estimated for the same seasons as in previous studies (see Supplementary Text Section 2.3), derived both by predicting NEE as well as subtracting GPP from ER suggest that the growing season net uptake is stronger than or similar to the estimates in Belshe et al. (2013) and Natali et al. (2019). The growing season NEE budget calculated for 100 days in the tundra was –296 Tg C yr⁻¹ in this study, compared to –137 ± 80 Tg C yr⁻¹ in Belshe et al. (2013). The growing season NEE budget estimated for 153 days in the northern permafrost region in this study was –1122 Tg C yr⁻¹, whereas the process model estimates varied between –687 and –1647 Tg C yr⁻¹ in Natali et al. (2019). Further, the observed daily average growing season NEE in tundra demonstrated a stronger sink strength than the average growing season NEE reported in McGuire et al. (2012) and Belshe et al. (2013) (–0.6, –0.3, and –0.2 g C m⁻² d⁻¹, respectively). Even though we acknowledge that some plant uptake and CO₂ emissions occur outside of our defined growing season (i.e., our growing season estimates did not capture the spring and autumn seasons), our results demonstrate that growing season CO₂ uptake might be larger than previously thought.

## 4.4 | Summary and next steps in high-latitude CO₂ flux upscaling

Overall, our findings suggest that statistical predictions aimed at describing high-latitude CO₂ flux patterns provide new insights into the understanding of broad GPP and ER patterns but have uncertainty with NEE. Furthermore, this study demonstrates that machine learning models are a robust and accurate empirical approach to predicting high-latitude terrestrial CO₂ fluxes, and that no individual machine learning model outperformed the others. This therefore supports the use of ensemble predictions to reduce uncertainties associated with a

single method and to produce more robust predictions. Nevertheless, the building of better models with an improved flux measurement network remains the highest research priority. Our results suggest that the next steps for future high-latitude upscaling efforts are to (a) measure fluxes over the entire year in as many sites as possible, (b) establish new sites in data-poor regions and regions where CO₂ predictions were most uncertain, such as in European Russia, Siberia, eastern Canada, and Canadian Arctic Archipelago, and specifically in disturbed and high-Arctic conditions, (c) develop better geospatial predictors (e.g., describing soil moisture and nutrients or permafrost thaw) to explain fluxes, (d) conduct detailed sensitivity tests of the importance of the flux measurement method, data distribution, and different predictor datasets influencing the budgets, and (e) build models at a finer temporal resolution than annual and growing season, to capture rapidly changing transition periods and bypass issues associated with temporal aggregation and varying definitions of seasons. High-latitude specific models are needed to more accurately monitor current emissions and improve understanding of the role of high-latitude regions in the global carbon cycle, as large changes in carbon cycling are likely in the near future.

## ACKNOWLEDGEMENTS

AMV was supported by Nordenskiöld-samfundet, The Finnish Cultural Foundation, Alfred Kordelin Foundation, Väisälä fund, and Jenny and Antti Wihuri Foundation. AMV and ML were also funded by the Academy of Finland (grant 286950). JA acknowledges the funding by Academy of Finland (grants 33761, 337552), while AL acknowledges strategic research funding by the Academy of Finland for SOMPA project (grants 312912 and 325680). TT was funded by the Swedish National Space Board (SNSB Dnr 95/16). BR was supported by the NASA Carbon Cycle Science and Arctic-Boreal Vulnerability Experiment programs (ABOVE grant NNX17AE13G), SMN by NASA ABoVE (grant NNX15AT81A) and JDW by NNX15AT81A and NASA NIP grant NNH17ZDA001N. AMV, BR, SN, and JDW were also funded by the Gordon and Betty Moore foundation (grant #8414). EAGS acknowledges NSF Research, Synthesis, and Knowledge Transfer in a Changing Arctic: Science Support for the Study of Environmental Arctic Change (grant #1331083) and NSF PLR Arctic System Science Research Networking Activities (Permafrost Carbon Network: Synthesizing Flux Observations for Benchmarking Model Projections of Permafrost Carbon Exchange; grant #1931333). JK acknowledges NSF grant 1203583, DZ NSF 1204263 and 1702797 and WO NSF 1204263, and 1702798. WO and DZ further acknowledge NOAA NA16SEC4810008, NASA NNX-15AT74A and NNX16AF94A, EU Horizon 2020 727890, and UK NERC NE/ P002552/1. HK, MU, and HI were funded by Arctic Challenge for Sustainability II grant JPMXD1420318865, and EH and PL by Natural Sciences and Engineering Research Council. MG acknowledges European Commission (INTAROS project, H2020-BG-09-2016, project 727890) and ESE NSF grants DEB-1636476, AON 856864, 1304271, 0632264, and 1107892, and the US Geological Survey. MM was funded by Academy of Finland (grant 317054) and MM and PJM were funded by the EU 6th Framework Programme project
CARBO-North (grant 036993). CB and CV were funded by the EU FP7-ENV project PAGE21 (grant 282700) and CB, CV, and PJM by the Nordic Center of Excellence project DEFROST. CB was further funded by the Academy of Finland (grant 314630), and CV by Academy of Finland (grant 332196). BE acknowledges Danish National Research Foundation (CENPERM DNRF100) and FJWP Research Council of Norway (Winterproof, grant 274711) and Swedish Research Council (WinterGap, project 2017-05268). JC was funded by FORMAS (grant 942-2015-49). VLSL and CE were funded by the Natural Sciences and Engineering Research Council and MP by FORMAS #2016-01289. JJ was funded by the Swedish Forest Society Foundation (2018-485-Steg 2 2017) and SFO by NSF grants PLR1504381 and PLR1836898. MST acknowledges Office of Biological and Environmental Research, DOE Office of Science; SJP the Korean government (NRF-2021M1A5A1 065425, KOPRI-PN21011); and NC the Korean government (MSIP) (NRF-2018R1D1A1B07047778 and NRF-2021M1A5A1065679). HS was funded by Netherlands Earth System Science Centre (NESSC), and IM by Academy of Finland Flagship funding (project 337549) and ICOS-Finland by University of Helsinki funding. RP was funded by Humboldt Fellowship for Experienced Researchers, MBN by Swedish Research Council, contract #2018-03966 and the national research infrastructures SITES and ICOS, funded by VR and partner institutes, and ELB by Greenland Research Council grant number 80.35, financed by the "Danish Program for Arctic Research". OS was supported through the Canada Research Chairs and Natural Sciences and Engineering Research Council Discovery Grants programs. The authors would also like to acknowledge Liangzhi Chen for his help with the literature review. Funding for the CO₂ flux synthesis workshop was provided by the Arctic Data Center.

## AUTHOR CONTRIBUTIONS

AMV and ML designed the study. AMV extracted the flux data from the literature and the data from the community call was designed and gathered by MM, TS et al. AMV, JA, and SP prepared the gridded datasets. ML, JA, and AMV developed the modeling framework. TT, CT, BR, JDW, and SMN commented on the analysis and AMV, with the help of JA and ML, conducted the analysis. Other authors contributed data and all authors were involved in the writing.

## DATA AVAILABILITY STATEMENT

Data are archived and freely available at Zenodo. The synthesis dataset is available at http://doi.org/10.5281/zenodo.4519583. Annual and averaged flux predictions are available at http://doi.org/10.5281/zenodo.4521852. The codes to run the statistical models and predictions together with the uncertainty estimation can be found in an R Markdown file as a supplement (Virkkalaetal_CO2flux_upscaling.pdf).

## ORCID

Anna-Maria Virkkala ID https://orcid.org/0000-0003-4877-2918

Juha Aalto ID https://orcid.org/0000-0001-6819-4911

Brendan M. Rogers ID https://orcid.org/0000-0001-6711-8466

Torbern Tagesson ID https://orcid.org/0000-0003-3011-1775

Claire C. Treat ID https://orcid.org/0000-0002-1225-8178

Susan M. Natali ID https://orcid.org/0000-0002-3010-2994

Aleksi Lehtonen ID https://orcid.org/0000-0003-1388-0388

Marguerite Mauritz ID https://orcid.org/0000-0001-8733-9119

Edward A. G. Schuur ID https://orcid.org/0000-0002-1096-2436

Mathias Goeckede ID https://orcid.org/0000-0003-2833-8401

Hiroki Iwata ID https://orcid.org/0000-0002-8962-8982

Peter M. Lafleur ID https://orcid.org/0000-0003-0347-9128

Stef Bokhorst ID https://orcid.org/0000-0003-0184-1162

Maija Marushchak ID https://orcid.org/0000-0002-2308-5049

Pertti J. Martikainen ID https://orcid.org/0000-0003-0415-8449

Bo Elberling ID https://orcid.org/0000-0002-6023-885X

Carolina Voigt ID https://orcid.org/0000-0001-8589-1428

Christina Biasi ID https://orcid.org/0000-0002-7413-3354

Oliver Sonnentag ID https://orcid.org/0000-0001-9333-9721

Frans-Jan W. Parmentier ID https://orcid.org/0000-0003-2952-7706

Masahito Ueyama ID https://orcid.org/0000-0002-4000-4888

Vincent L. St. Louis ID https://orcid.org/0000-0001-5405-1522

Craig A. Emmerton ID https://orcid.org/0000-0001-9511-9191

Matthias Peichl ID https://orcid.org/0000-0002-9940-5846

Jinshu Chi ID https://orcid.org/0000-0001-5688-8895

Järvi Järveoja ID https://orcid.org/0000-0001-6317-660X

Mats B. Nilsson ID https://orcid.org/0000-0003-3765-6399

Sang-Jong Park ID https://orcid.org/0000-0002-6944-6962

Ivan Mammarella ID https://orcid.org/0000-0002-8516-3356

Rafael Poyatos ID https://orcid.org/0000-0003-0521-2523

Efrén López-Blanco ID https://orcid.org/0000-0002-3796-8408

Torben Røjle Christensen ID https://orcid.org/0000-0002-4917-148X

Min Jung Kwon ID https://orcid.org/0000-0002-7330-2320

Torsten Sachs ID https://orcid.org/0000-0002-9959-4771

David Holl ID https://orcid.org/0000-0002-9269-7030

Miska Luoto ID https://orcid.org/0000-0001-6203-5143

## REFERENCES

Aalto, J., Karjalainen, O., Hjort, J., & Luoto, M. (2018). Statistical forecasting of current and future circum-arctic ground temperatures and active layer thickness. *Geophysical Research Letters*, 45, 4889-4898.
https://doi.org/10.1029/2018GL078007

Abbott, B. W., Jones, J. B., Schuur, E. A. G., Chapin III, F. S., Bowden, W. B., Bret-Harte, M. S., Epstein, H. E., Flannigan, M. D., Harms, T. K., Hollingsworth, T. N., Mack, M. C., McGuire, A. D., Natali, S. M., Rocha, A. V., Tank, S. E., Turetsky, M. R., Vonk, J. E., Wickland, K. P., Aiken, G. R., ... Zimov, S. (2016). Biomass offsets little or none of permafrost carbon release from soils, streams, and wildfire: An expert assessment. *Environmental Research Letters*, 11(3), 1-13.
https://doi.org/10.1088/1748-9326/11/3/034014

Ai, J., Jia, G., Epstein, H. E., Wang, H., Zhang, A., & Hu, Y. (2018). MODIS-based estimates of global terrestrial ecosystem respiration. *Journal of Geophysical Research: Biogeosciences*, 123(2), 326-352. https://doi.org/10.1002/2017JG004107

Alekseychik, P., Mammarella, I., Karpov, D., Dengel, S., Terentieva, I., Sabrekov, A., & Lapshina, E. (2017). Net ecosystem exchange and energy fluxes in a West Siberian bog. *Atmospheric Chemistry and Physics Discussions*, 17, 9333-9345. https://doi.org/10.5194/acp-17-9333-2017

Alton, P. B. (2020). Representativeness of global climate and vegetation by carbon-monitoring networks; implications for estimates of gross and net primary productivity at biome and global levels.
Agricultural and Forest Meteorology, 290(April), 108017. https://doi.org/10.1016/j.agrformet.2020.108017

Araújo, M. B., & New, M. (2007). Ensemble forecasting of species distributions. *Trends in Ecology and Evolution*, 22(1), 42–47. https://doi.org/10.1016/j.tree.2006.09.010

Arens, S. J. T., Sullivan, P. F., & Welker, J. M. (2008). Nonlinear responses to nitrogen and strong interactions with nitrogen and phosphorus additions drastically alter the structure and function of a high arctic ecosystem. *Journal of Geophysical Research-Biogeosciences*, 113(G3), 10. https://doi.org/10.1029/2007jg000508

Bäckstrand, K., Crill, P. M., Jackowicz-Korczyński, M., Mastepanov, M., Christensen, T. R., & Bastviken, D. (2010). Annual carbon gas budget for a subarctic peatland. *Biogeosciences*, 7(1), 95–108. https://doi.org/10.5194/bg-7-95-2010

Baldocchi, D., Chu, H., & Reichstein, M. (2018). Inter-annual variability of net and gross ecosystem carbon fluxes: A review. *Agricultural and Forest Meteorology*, 249, 520–533. https://doi.org/10.1016/j.agrformet.2017.05.015

Baldocchi, D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C. H., Davis, K., Evans, R., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X., Malhi, Y., Meyers, T., Munger, W., Oechel, W., ... Wofsy, S. (2001). FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. *Bulletin of the American Meteorological Society*, 82(11), 2415–2434.

Baldocchi, D. D., Hicks, B. B., & Meyers, T. P. (1988). Measuring biosphere-atmosphere exchanges of biologically related gases with micrometeorological methods. *Ecology*, 69(5), 1331–1340. https://doi.org/10.2307/1941631

Beer, C., Reichstein, M., Tomelleri, E., Ciais, P., Jung, M., Carvalhais, N., Rodenbeck, C., Arain, M. A., Baldocchi, D., Bonan, G. B., Bondeau, A., Cescatti, A., Lasslop, G., Lindroth, A., Lomas, M., Luyssaert, S., Margolis, H., Olesen, K. W., Roupsard, O., ... Papale, D. (2010). Terrestrial gross carbon dioxide uptake: Global distribution and covariation with climate. *Science*, 329(5993), 834–838. https://doi.org/10.1126/science.1184984

Belshe, E. F., Schuur, E. A. G., & Bolker, B. M. (2013). Tundra ecosystems observed to be CO₂ sources due to differential amplification of the carbon cycle. *Ecology Letters*, 16(10), 1307–1315. https://doi.org/10.1111/ele.12164

Berner, L. T., Massey, R., Jantz, P., Forbes, B. C., Macias-Fauria, M., Myers-Smith, I., Kumpula, T., Gauthier, G., Andreu-Hayles, L., Gaglioti, B. V., Burns, P., Zetterberg, P., D'Arrigo, R., & Goetz, S. J. (2020). Summer warming explains widespread but not uniform greening in the Arctic tundra biome. *Nature Communications*, 11(1), 1–12. https://doi.org/10.1038/s41467-020-18479-5

Björkman, M. P., Morgner, E., Björk, R. G., Cooper, E. J., Elberling, B., & Klemedtsson, L. (2010). A comparison of annual and seasonal carbon dioxide effluxes between sub-Arctic Sweden and High-Arctic Svalbard. *Polar Research*, 29(1), 75–84. https://doi.org/10.1111/j.1751-8369.2010.00150.x

Bodesheim, P., Jung, M., Gans, F., Mahecha, M. D., & Reichstein, M. (2018). Upscaled diurnal cycles of land-Atmosphere fluxes: A new global half-hourly data product. *Earth System Science Data*, 10(3), 1327–1365. https://doi.org/10.5194/essd-10-1327-2018

Bond-Lamberty, B., & Thomson, A. (2010). Temperature-associated increases in the global soil respiration record. *Nature*, 464, 579–582. https://doi.org/10.1038/nature08930

Bradshaw, C. J. A., & Warkentin, I. G. (2015). Global estimates of boreal forest carbon stocks and flux. *Global and Planetary Change*, 128, 24–30. https://doi.org/10.1016/j.gloplacha.2015.02.004

Brown, J., Ferrians, O., Heginbottom, J. A., & Melnikov, E. (2002). *Circum-arctic map of permafrost and ground-ice conditions, version 2*. NSIDC National Snow and Ice Data Center.

Cahoon, S. M. P., Sullivan, P. F., Shaver, G. R., Welker, J. M., & Post, E. (2012). Interactions among shrub cover and the soil microclimate

may determine future Arctic carbon budgets. *Ecology Letters*, 15, 1415–1422. https://doi.org/10.1111/j.1461-0248.2012.01865.x

Camps-Valls, G., Jung, M., Ichii, K., Papale, D., Tramontana, G., Bodesheim, P., & Reichstein, M. (2015). Ranking drivers of global carbon and energy fluxes over land. *International Geoscience and Remote Sensing Symposium*, 4416–4419. https://doi.org/10.1109/IGARSS.2015.7326806

Cassidy, A. E., Christen, A., & Henry, G. H. R. (2016). The effect of a permafrost disturbance on growing-season carbon-dioxide fluxes in a high Arctic tundra ecosystem. *Biogeosciences*, 13(8), 2291–2303. https://doi.org/10.5194/bg-13-2291-2016

Celis, G., Mauritz, M., Bracho, R., Salmon, V. G., Webb, E. E., Hutchings, J., & Schuur, E. A. G. (2017). Tundra is a consistent source of CO₂ at a site with progressive permafrost thaw during 6 years of chamber and eddy covariance measurements. *Journal of Geophysical Research: Biogeosciences*, 122(6), 1471–1485. https://doi.org/10.1002/2016JG003671

Chu, H., Baldocchi, D. D., John, R., Wolf, S., & Reichstein, M. (2017). Fluxes all of the time? A primer on the temporal representativeness of FLUXNET. *Journal of Geophysical Research: Biogeosciences*, 122(2), 289–307. https://doi.org/10.1002/2016JG003576

Coursolle, C., Giasson, M.-A., Margolis, H. A., & Bernier, P. Y. (2012). Moving towards carbon neutrality: CO₂ exchange of a black spruce forest ecosystem during the first 10 years of recovery after harvest. *Canadian Journal of Forest Research*, 42(11), 1908–1918. https://doi.org/10.1139/x2012-133

Dagg, J., & Lafleur, P. (2011). Vegetation community, foliar nitrogen, and temperature effects on tundra CO₂ exchange across a soil moisture gradient. *Arctic Antarctic and Alpine Research*, 43(2), 189–197. https://doi.org/10.1657/1938-4246-43.2.189

Davidson, E. A. J. E. A., Janssens, I. A. I. A., Marks, D., Murdock, M. D., Ahl, R. S., Woods, S. W., & Loffler, J. (2006). Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. *Nature*, 440(7081), 165–173. https://doi.org/10.1038/nature04514

Dinerstein, E., Olson, D., Joshi, A., Vynne, C., Burgess, N. D., Wikramanayake, E., Hahn, N., Palminteri, S., Hedao, P., Noss, R., Hansen, M., Locke, H., Ellis, E. C., Jones, B., Barber, C. V., Hayes, R., Kormos, C., Martin, V., Crist, E., ... Saleem, M. (2017). An ecoregion-based approach to protecting half the terrestrial realm. *BioScience*, 67(6), 534–545. https://doi.org/10.1093/biosci/bix014

Elith, J., Leathwick, J. R., & Hastie, T. (2008). A working guide to boosted regression trees. *Journal of Animal Ecology*, 77(4), 802–813. https://doi.org/10.1111/j.1365-2656.2008.01390.x

Emmerton, C. A., St. Louis, V. L., Humphreys, E. R., Gamon, J. A., Barker, J. D., & Pastorello, G. Z. (2016). Net ecosystem exchange of CO₂ with rapidly changing high Arctic landscapes. *Global Change Biology*, 22(3), 1185–1200. https://doi.org/10.1111/gcb.13064

ESA (2017). Land cover CCI product user guide version 2. Tech. Rep. Retrieved from http://maps.elie.ucl.ac.be/CCI/viewer/download/ESACCI-LC-Ph2-PUGv2_2.0.pdf

Euskirchen, E. S., Edgar, C. W., Turetsky, M. R., Waldrop, M. P., & Harden, J. W. (2014). Differential response of carbon fluxes to climate in three peatland ecosystems that vary in the presence and stability of permafrost. *Journal of Geophysical Research G: Biogeosciences*, 119(8), 1576–1595. https://doi.org/10.1002/2014JG002683

Fox, A. M., Huntley, B., Lloyd, C. R., Williams, M., & Baxter, R. (2008). Net ecosystem exchange over heterogeneous Arctic tundra: Scaling between chamber and eddy covariance measurements. *Global Biogeochemical Cycles*, 22(2), 1–15. https://doi.org/10.1029/2007G003027

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Hauck, J., Olsen, A., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Quéré, C., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S., Aragão, L. E. O. C., Arneth, A., Arora, V., Bates, N. R., ... Zaehle, S. (2020). Global carbon budget 2020. *Earth System Science Data*, 12, 3269–3340.
<figure><img src="image_11.png" /></figure>

Gasser, T., Kechiar, M., Ciais, P., Burke, E. J., Kleinen, T., Zhu, D., Huang, Y., Ekici, A., & Obersteiner, M. (2018). Path-dependent reductions in CO₂ emission budgets caused by permafrost carbon release. *Nature Geoscience*, *11*(11), 830–835. https://doi.org/10.1038/s41561-018-0227-0

Heliasz, M., Johansson, T., Lindroth, A., Mölder, M., Mastepanov, M., Friborg, T., Callaghan, T. V., & Christensen, T. R. (2011). Quantification of C uptake in subarctic birch forest after setback by an extreme insect outbreak. *Geophysical Research Letters*, *38*(1), 1–5. https://doi.org/10.1029/2010GL044733

Hugelius, G., Strauss, J., Zubrzycki, S., Harden, J. W., Schuur, E. A. G., Ping, C.-L., Schirrmeister, L., Grosse, G., Michaelson, G. J., Koven, C. D., O'Donnell, J. A., Elberling, B., Mishra, U., Camill, P., Yu, Z., Palmtag, J., & Kuhry, P. (2014). Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps. *Biogeosciences*, *11*(23), 6573–6593. https://doi.org/10.5194/bg-11-6573-2014

Hursh, A., Ballantyne, A., Cooper, L., Maneta, M., Kimball, J., & Watts, J. (2016). The sensitivity of soil respiration to soil temperature, moisture, and carbon supply at the global scale. *Global Change Biology*, *23*(5), 2090–2103. https://doi.org/10.1111/gcb.13489

Ichii, K., Ueyama, M., Kondo, M., Saigusa, N., Kim, J., Alberto, M. C., & Zhao, F. (2017). New data-driven estimation of terrestrial CO₂ fluxes in Asia using a standardized database of eddy covariance measurements, remote sensing data, and support vector regression. *Journal of Geophysical Research: Biogeosciences*, *122*(4), 767–795. https://doi.org/10.1002/2016JG003640

Illeris, L., Michelsen, A., & Jonasson, S. (2003). Soil plus root respiration and microbial biomass following water, nitrogen, and phosphorus application at a high arctic semi desert. *Biogeochemistry*, *65*(1), 15–29. https://doi.org/10.1023/a:1026034523499

Iwata, H., Ueyama, M., Harazono, Y., Tsuyuzaki, S., Kondo, M., & Uchida, M. (2011). Quick recovery of carbon dioxide exchanges in a burned black spruce forest in Interior Alaska. *Scientific Online Letters on the Atmosphere*, *7*(1), 105–108. https://doi.org/10.2151/sola.2011-027

Järveoja, J., Nilsson, M. B., Gažovič, M., Crill, P. M., & Peichl, M. (2018). Partitioning of the net CO₂ exchange using an automated chamber system reveals plant phenology as key control of production and respiration fluxes in a boreal peatland. *Global Change Biology*, *24*(8), 3436–3451. https://doi.org/10.1111/gcb.14292

Jung, M., Reichstein, M., Schwalm, C. R., Huntingford, C., Sitch, S., Ahlström, A., Arneth, A., Camps-Valls, G., Ciais, P., Friedlingstein, P., Gans, F., Ichii, K., Jain, A. K., Kato, E., Papale, D., Poulter, B., Raduly, B., Rödenbeck, C., Tramontana, G., ... Zeng, N. (2017). Compensatory water effects link yearly global land CO₂ sink changes to temperature. *Nature*, *541*(7638), 516–520. https://doi.org/10.1038/nature20780

Jung, M., Schwalm, C., Migliavacca, M., Walther, S., Camps-Valls, G., Koirala, S., & Walker, A. (2020). Scaling carbon fluxes from eddy covariance sites to globe: Synthesis and evaluation of the FLUXCOM approach. *Biogeosciences*, *17*, 1343–1365. https://doi.org/10.5194/bg-2019-368

Karelin, D. V., Zamolodchikov, D. G., Zukert, N. V., Chestnykh, O. V., Pochikalov, A. V., & Krayev, G. N. (2013). Interannual changes in PAR and soil moisture during the warm season may be more important for directing of annual carbon balance in tundra than temperature fluctuations. *Zhurnal Obshchei Biologii*, *74*(1), 3–22. https://doi.org/10.1134/S2079086413050058

Kittler, F., Heimann, M., Kolle, O., Zimov, N., Zimov, S., & Göckede, M. (2017). Long-term drainage reduces CO₂ uptake and CH₄ emissions in a Siberian permafrost ecosystem. *Global Biogeochemical Cycles*, *31*(12), 1704–1717. https://doi.org/10.1002/2017GB005774

Kolari, P., Kulmala, L., Pumpanen, J., Launiainen, S., Ilvesniemi, H., Hari, P., & Nikinmaa, E. (2009). CO₂ exchange and component CO₂ fluxes of a boreal Scots pine forest. *Boreal Environment Research*, *14*, 761–783.

La Puma, I. P., Philippi, T. E., & Oberbauer, S. F. (2007). Relating NDVI to ecosystem CO₂ exchange patterns in response to season length and soil warming manipulations in arctic Alaska. *Remote Sensing of Environment*, *109*(2), 225–236. https://doi.org/10.1016/j.rse.2007.01.001

Lafleur, P. M., Humphreys, E. R., St. Louis, V. L., Myklebust, M. C., Papakyriakou, T., Poissant, L., Barker, J. D., Pilote, M., & Swystun, K. A. (2012). Variation in peak growing season net ecosystem production across the Canadian arctic. *Environmental Science and Technology*, *46*(15), 7971–7977. https://doi.org/10.1021/es300500m

Lasslop, G., Reichstein, M., Papale, D., Richardson, A. D., Arneth, A., Barr, A., Stoy, P., & Wohlfahrt, G. (2010). Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: Critical issues and global evaluation. *Global Change Biology*, *16*(1), 187–208. https://doi.org/10.1111/j.1365-2486.2009.02041.x

López-Blanco, E., Exbrayat, J.-F., Lund, M., Christensen, T. R., Tamstorf, M. P., Slevin, D., Hugelius, G., Bloom, A. A., & Williams, M. (2019). Evaluation of terrestrial pan-Arctic carbon cycling using a data-assimilation system. *Earth System Dynamics*, *10*(2), 233–255. https://doi.org/10.5194/esd-10-233-2019

López-Blanco, E., Lund, M., Williams, M., Tamstorf, M. P., Westergaard-Nielsen, A., Exbrayat, J. F., & Christensen, T. R. (2017). Exchange of CO₂ in Arctic tundra: Impacts of meteorological variations and biological disturbance. *Biogeosciences*, *14*(19), 4467–4483. https://doi.org/10.5194/bg-14-4467-2017

Lund, M., Christensen, T. R., Mastepanov, M., Lindroth, A., & Ström, L. (2009). Effects of N and P fertilization on the greenhouse gas exchange in two northern peatlands with contrasting N deposition rates. *Biogeosciences*, *6*(10), 2135–2144. https://doi.org/10.5194/bg-6-2135-2009

Lund, M., Lafleur, P. M., Roulet, N. T., Lindroth, A., Christensen, T. R., Aurela, M., Chojnicki, B. H., Flanagan, L. B., Humphreys, E. R., Laurila, T., Oechel, W. C., Olejnik, J., Rinne, J., Schubert, P., & Nilsson, M. B. (2010). Variability in exchange of CO₂ across 12 northern peatland and tundra sites. *Global Change Biology*, *16*(9), 2436–2448. https://doi.org/10.1111/j.1365-2486.2009.02104.x

Lund, M., Raudrupp, K., Westergaard-Nielsen, A., López-Blanco, E., Nymand, J., & Aastrup, P. (2017). Larval outbreaks in West Greenland: Instant and subsequent effects on tundra ecosystem productivity and CO₂ exchange. *Ambio*, *46*, 26–38. https://doi.org/10.1007/s13280-016-0863-9

Lundegårdh, H. (1927). Carbon dioxide evolution of soil and crop growth. *Soil Science*, *23*(6), 417–453.

Machimura, T., Kobayashi, Y., Iwahana, G., Hirano, T., Lopez, L., Fukuda, M., & Fedorov, A. N. (2005). Change of carbon dioxide budget during three years after deforestation in eastern Siberian larch forest. *Journal of Agricultural Meteorology*, *60*(5), 653–656. https://doi.org/10.2480/agmet.653

Marushchak, M. E., Kiepe, I., Biasi, C., Elsakov, V., Friborg, T., Johansson, T., Soegaard, H., Virtanen, T., & Martikainen, P. J. (2013). Carbon dioxide balance of subarctic tundra from plot to regional scales. *Biogeosciences*, *10*(1), 437–452. https://doi.org/10.5194/bg-10-437-2013

McCallum, I., Franklin, O., Moltchanova, E., Merbold, L., Schmullius, C., Shvidenko, A., Schepaschenko, D., & Fritz, S. (2013). Improved light and temperature responses for light-use-efficiency-based GPP models. *Biogeosciences*, *10*(10), 6577–6590. https://doi.org/10.5194/bg-10-6577-2013

McGuire, A. D., Anderson, L. G., Christensen, T. R., Dallimore, S., Guo, L., Hayes, D. J., Heimann, M., Lorenson, T. D., Macdonald, R. W., & Roulet, N. (2009). Sensitivity of the carbon cycle in the Arctic to climate change. *Ecological Monographs*, *79*(4), 523–555.

McGuire, A. D., Christensen, T. R., Hayes, D., Heroult, A., Euskirchen, E., Kimball, J. S., Koven, C., Lafleur, P., Miller, P. A., Oechel, W., Peylin,
<figure><img src="image_12.png" /></figure>

P., Williams, M., & Yi, Y. (2012). An assessment of the carbon balance of Arctic tundra: comparisons among observations, process models, and atmospheric inversions. *Biogeosciences*, 9(8), 3185-3204. https://doi.org/10.5194/bg-9-3185-2012

McGuire, A. D., Koven, C., Lawrence, D. M., Clein, J. S., Xia, J., Beer, C., Burke, E., Chen, G., Chen, X., Delire, C., Jafarov, E., MacDougall, A. H., Marchenko, S., Nicolsky, D., Peng, S., Rinke, A., Saito, K., Zhang, W., Alkama, R., ... Zhuang, Q. (2016). Variability in the sensitivity among model simulations of permafrost and carbon dynamics in the permafrost region between 1960 and 2009. *Global Biogeochemical Cycles*, 30, 1015–1037. https://doi.org/10.1002/2016GB005405

McGuire, A. D., Lawrence, D. M., Koven, C., Clein, J. S., Burke, E., Chen, G., Jafarov, E., MacDougall, A. H., Marchenko, S., Nicolsky, D., Peng, S., Rinke, A., Ciais, P., Gouttevin, I., Hayes, D. J., Ji, D., Krinner, G., Moore, J. C., Romanovsky, V., ... Zhuang, Q. (2018). Dependence of the evolution of carbon dynamics in the northern permafrost region on the trajectory of climate change. *Proceedings of the National Academy of Sciences of the United States of America*, 115(15), 3882–3887. https://doi.org/10.1073/pnas.1719903115

Metcalfe, D. B., Hermans, T. D. G., Ahlstrand, J., Becker, M., Berggren, M., Björk, R. G., Björkman, M. P., Blok, D., Chaudhary, N., Chisholm, C., Classen, A. T., Hasselquist, N. J., Jonsson, M., Kristensen, J. A., Kumordzi, B. B., Lee, H., Mayor, J. R., Prevéy, J., Pantazatou, K., ... Abdi, A. M. (2018). Patchy field sampling biases understanding of climate change impacts across the Arctic. *Nature Ecology and Evolution*, 2(9), 1443–1448. https://doi.org/10.1038/s41559-018-0612-5

Myers-Smith, I. H., Kerby, J. T., Phoenix, G. K., Bjerke, J. W., Epstein, H. E., Assmann, J. J., John, C., Andreu-Hayles, L., Angers-Blondin, S., Beck, P. S. A., Berner, L. T., Bhatt, U. S., Bjorkman, A. D., Blok, D., Bryn, A., Christiansen, C. T., Cornelissen, J. H. C., Cunliffe, A. M., Elmendorf, S. C., ... Wipf, S. (2020). Complexity revealed in the greening of the Arctic. *Nature Climate Change*, 10, 106–117. https://doi.org/10.1038/s41558-019-0688-1

Natali, S. M., Watts, J. D., Rogers, B. M., Potter, S., Ludwig, S. M., Selbmann, A.-K., Sullivan, P. F., Abbott, B. W., Arndt, K. A., Birch, L., Björkman, M. P., Bloom, A. A., Celis, G., Christensen, T. R., Christiansen, C. T., Commane, R., Cooper, E. J., Crill, P., Czimczik, C., ... Zona, D. (2019). Large loss of CO₂ in winter observed across the northern permafrost region. *Nature Climate Change*, 9, 852–857. https://doi.org/10.1038/s41558-019-0592-8

Niittynen, P., & Luoto, M. (2018). The importance of snow in species distribution models of arctic vegetation. *Ecography*, 41(6), 1024–1037. https://doi.org/10.1111/ecog.03348

Nobrega, S., & Grogan, P. (2007). Deeper snow enhances winter respiration from both plant-associated and bulk soil carbon pools in birch hummock tundra. *Ecosystems*, 10(3), 419–431. https://doi.org/10.1007/s10021-007-9033-z

Nobrega, S., & Grogan, P. (2008). Landscape and ecosystem-level controls on net carbon dioxide exchange along a natural moisture gradient in Canadian low arctic tundra. *Ecosystems*, 11(3), 377–396. https://doi.org/10.1007/s10021-008-9128-1

Pan, Y., Birdsey, R. A., Fang, J., Houghton, R., Kauppi, P. E., Kurz, W. A., & Hayes, S. (2011). A large and persistent carbon sink in the world's forests. *Science*, 333, 1239–1243. https://doi.org/10.1002/clen.200720011

Parmentier, F.-J.-W., Sonnentag, O., Mauritz, M., Virkkala, A.-M., & Schuur, E. A. G. (2019). Is the northern permafrost zone a source or a sink for carbon? *EOS*, 100.

Pastorello, G., Trotta, C., Canfora, E., Chu, H., Christianson, D., Cheah, Y.-W., Poindexter, C., Chen, J., Elbashandy, A., Humphrey, M., Isaac, P., Polidori, D., Reichstein, M., Ribeca, A., van Ingen, C., Vuichard, N., Zhang, L., Amiro, B., Ammann, C., ... Papale, D. (2020). The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data. *Scientific Data*, 7(1), 225. https://doi.org/10.1038/s41597-020-0534-3

Peltola, O., Vesala, T., Gao, Y., Räty, O., Alekseychik, P., Aurela, M., Chojnicki, B., Desai, A. R., Dolman, A. J., Euskirchen, E. S., Friborg, T., Göckede, M., Helbig, M., Humphreys, E., Jackson, R. B., Jocher, G., Joos, F., Klatt, J., Knox, S. H., ... Aalto, T. (2019). Monthly gridded data product of northern wetland methane emissions based on upscaling eddy covariance observations. *Earth System Science Data*, 11, 1263–1289. https://doi.org/10.5194/essd-11-1263-2019

Petrone, R. M., Chasmer, L., Hopkinson, C., Silins, U., Landhäuser, S. M., Kljun, N., & Devito, K. J. (2014). Effects of harvesting and drought on CO₂ and H₂O fluxes in an aspen-dominated western boreal plain forest: Early chronosequence recovery. *Canadian Journal of Forest Research*, 45(1), 87–100. https://doi.org/10.1139/cjfr-2014-0253

Rawlins, M. A., McGuire, A. D., Kimball, J. S., Dass, P., Lawrence, D., Burke, E., & Sueyoshi, T. (2015). Assessment of model estimates of land-atmosphere CO₂ exchange across Northern Eurasia. *Biogeosciences*, 12(14), 4385–4405. https://doi.org/10.5194/bg-12-4385-2015

Raynolds, M. K., Walker, D. A., Balser, A., Bay, C., Campbell, M., Cherosov, M. M., Daniëls, F. J. A., Eidesen, P. B., Ermokhina, K. A., Frost, G. V., Jedrzejek, B., Jorgenson, M. T., Kennedy, B. E., Kholod, S. S., Lavrinenko, I. A., Lavrinenko, O. V., Magnússon, B., Matveyeva, N. V., Metúsalemsson, S., ... Troeva, E. (2019). A raster version of the Circumpolar Arctic Vegetation Map (CAVM). *Remote Sensing of Environment*, 232, 111297. https://doi.org/10.1016/j.rse.2019.111297

Reichstein, M., Falge, E., Baldocchi, D., Papale, D., Aubinet, M., Berbigier, P., Bernhofer, C., Buchmann, N., Gilmanov, T., Granier, A., Grunwald, T., Havrankova, K., Ilvesniemi, H., Janous, D., Knohl, A., Laurila, T., Lohila, A., Loustau, D., Matteucci, G., ... Valentini, R. (2005). On the separation of net ecosystem exchange into assimilation and ecosystem respiration: Review and improved algorithm. *Global Change Biology*, 11(9), 1424–1439. https://doi.org/10.1111/j.1365-2486.2005.001002.x

Schuur, E. A. G., Bockheim, J., Canadell, J. G., Euskirchen, E., Field, C. B., Goryachkin, S. V., Hagemann, S., Kuhry, P., Lafleur, P. M., Lee, H., Mazhitova, G., Nelson, F. E., Rinke, A., Romanovsky, V. E., Shiklomanov, N., Tarnocai, C., Venevsky, S., Vogel, J. G., & Zimov, S. A. (2008). Vulnerability of permafrost carbon to climate change: Implications for the global carbon cycle. *BioScience*, 58(8), 701. https://doi.org/10.1641/B580807

Schuur, E. A. G., McGuire, A. D., Schädel, C., Grosse, G., Harden, J. W., Hayes, D. J., Hugelius, G., Koven, C. D., Kuhry, P., Lawrence, D. M., Natali, S. M., Olefeldt, D., Romanovsky, V. E., Schaefer, K., Turetsky, M. R., Treat, C. C., & Vonk, J. E. (2015). Climate change and the permafrost carbon feedback. *Nature*, 520(7546), 171–179. https://doi.org/10.1038/nature14338

Segurado, P., & Araújo, M. B. (2004). An evaluation of methods for modelling species distributions. *Journal of Biogeography*, 31(10), 1555–1568.

Shaver, G. R., Street, L. E., Rastetter, E. B., Van Wijk, M. T., & Williams, M. (2007). Functional convergence in regulation of net CO₂ flux in heterogeneous tundra landscapes in Alaska and Sweden. *Journal of Ecology*, 95(4), 802–817. https://doi.org/10.1111/j.1365-2745.2007.01259.x

Sørensen, M. V., Graae, B. J., Classen, A., Enquist, B. J., & Strimbeck, R. (2019). Drivers of C cycling in three arctic-alpine plant communities. *Arctic, Antarctic, and Alpine Research*, 51(1), 128–147. https://doi.org/10.1080/15230430.2019.1592649

Taylor, K. E., Stouffer, R. J., & Meehl, G. A. (2012). An overview of CMIP5 and the experiment design. *Bulletin of the American Meteorological Society*, 93(4), 485–498. https://doi.org/10.1175/BAMS-D-11-00094.1

Tramontana, G., Jung, M., Schwalm, C. R., Ichii, K., Camps-Valls, G., Ráduly, B., & Papale, D. (2016). Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms. *Biogeosciences*, 13(14), 4291–4313. https://doi.org/10.5194/bg-13-4291-2016
<figure><img src="image_13.png" /></figure>

Treat, C. C., Marushchak, M. E., Voigt, C., Zhang, Y. U., Tan, Z., Zhuang, Q., Virtanen, T. A., Räsänen, A., Biasi, C., Hugelius, G., Kaverin, D., Miller, P. A., Stendel, M., Romanovsky, V., Rivkin, F., Martikainen, P. J., & Shurpali, N. J. (2018). Tundra landscape heterogeneity, not interannual variability, controls the decadal regional carbon balance in the Western Russian Arctic. *Global Change Biology*, 24(11), 5188-5204. https://doi.org/10.1111/gcb.14421

Trucco, C., Schuur, E. A. G. G., Natali, S. M., Belshe, E. F., Bracho, R., & Vogel, J. (2012). Seven-year trends of CO₂ exchange in a tundra ecosystem affected by long-term permafrost thaw. *Journal of Geophysical Research: Biogeosciences*, 117(2), 1-12. https://doi.org/10.1029/2011JG001907

Turetsky, M. R., Abbott, B. W., Jones, M. C., Anthony, K. W., Olefeldt, D., Schuur, E. A. G., Grosse, G., Kuhry, P., Hugelius, G., Koven, C., Lawrence, D. M., Gibson, C., Sannel, A. B. K., & McGuire, A. D. (2020). Carbon release through abrupt permafrost thaw. *Nature Geoscience*, 13(2), 138-143. https://doi.org/10.1038/s41561-019-0526-0

Ueyama, M., Ichii, K., Iwata, H., Euskirchen, E. S., Zona, D., Rocha, A. V., Harazono, Y., Iwama, C., Nakai, T., & Oechel, W. C. (2013). Upscaling terrestrial carbon dioxide fluxes in Alaska with satellite remote sensing and support vector regression. *Journal of Geophysical Research: Biogeosciences*, 118(3), 1266-1281. https://doi.org/10.1002/jgrg.20095

Ueyama, M., Iwata, H., Harazono, Y., Euskirchen, E. S., Oechel, W. C., & Zona, D. (2013). Growing season and spatial variations of carbon fluxes of Arctic and boreal ecosystems in Alaska (USA). *Ecological Applications*, 23(8), 1798-1816. https://doi.org/10.1890/11-0875.1

Ueyama, M., Iwata, H., Nagano, H., Tahara, N., Iwama, C., & Harazono, Y. (2019). Carbon dioxide balance in early-successional forests after forest fires in interior Alaska. *Agricultural and Forest Meteorology*, 275, 196-207. https://doi.org/10.1016/j.agrformet.2019.05.020

Virkkala, A.-M., Abdi, A. M., Luoto, M., & Metcalfe, D. (2019). Identifying multidisciplinary research gaps across Arctic terrestrial gradients. *Environmental Research Letters*, 14, 124061. https://doi.org/10.1088/1748-9326/ab4291

Virkkala, A.-M., Virtanen, T., Lehtonen, A., Rinne, J., & Luoto, M. (2018). The current state of CO₂ flux chamber studies in the Arctic tundra: A review. *Progress in Physical Geography*, 42(2), 162-184. https://doi.org/10.1177/0309133317745784

Vogel, J., Schuur, E. A. G. G., Trucco, C., & Lee, H. (2009). Response of CO₂ exchange in a tussock tundra ecosystem to permafrost thaw and thermokarst development. *Journal of Geophysical Research-Biogeosciences*, 114(G4), 1-14. https://doi.org/10.1029/2008JG000901

Voigt, C., Marushchak, M. E., Mastepanov, M., Lamprecht, R. E., Christensen, T. R., Dorodnikov, M., Jackowicz-Korczyński, M., Lindgren, A., Lohila, A., Nykänen, H., Oinonen, M., Oksanen, T., Palonen, V., Treat, C. C., Martikainen, P. J., & Biasi, C. (2019). Ecosystem carbon response of an Arctic peatland to simulated permafrost thaw. *Global Change Biology*, 25(5), 1746-1764. https://doi.org/10.1111/gcb.14574

Walker, X. J., Baltzer, J. L., Cumming, S. G., Day, N. J., Ebert, C., Goetz, S., Johnstone, J. F., Potter, S., Rogers, B. M., Schuur, E. A. G., Turetsky, M. R., & Mack, M. C. (2019). Increasing wildfires threaten historic carbon sink of boreal forest soils. *Nature*, 572(7770), 520-523. https://doi.org/10.1038/s41586-019-1474-y

Wania, R., Ross, I., & Prentice, I. C. (2009). Integrating peatlands and permafrost into a dynamic global vegetation model: 2. Evaluation and sensitivity of vegetation and carbon cycle processes. *Global Biogeochemical Cycles*, 23, GB3015. https://doi.org/10.1029/2008G003413

Warner, D. L., Bond-Lamberty, B., Jian, J., Stell, E., & Vargas, R. (2019). Spatial predictions and associated uncertainty of annual soil respiration at the global scale. *Global Biogeochemical Cycles*, 33(12), 1733-1745. https://doi.org/10.1029/2019GB006264

Watts, J. D., Kimball, J. S., Parmentier, F. J. W., Sachs, T., Rinne, J., Zona, D., & Aurela, M. (2014). A satellite data driven biophysical modeling approach for estimating northern peatland and tundra CO₂ and CH₄ fluxes. *Biogeosciences*, 11(7), 1961-1980. https://doi.org/10.5194/bg-11-1961-2014

Welker, J. M., Fahnestock, J. T., Henry, G. H. R., O'Dea, K. W., & Chimner, R. A. (2004). CO₂ exchange in three Canadian High Arctic ecosystems: Response to long-term experimental warming. *Global Change Biology*, 10(12), 1981-1995. https://doi.org/10.1111/j.1365-2486.2004.00857.x

Wilkman, E., Zona, D., Tang, Y., Gioli, B., Lipson, D. A., & Oechel, W. (2018). Temperature response of respiration across the heterogeneous landscape of the Alaskan arctic tundra. *Journal of Geophysical Research: Biogeosciences Research*, 123, 2287-2302. https://doi.org/10.1029/2017JG004227

Zhang, W., Jansson, P. E., Schurgers, G., Hollesen, J., Lund, M., Abermann, J., & Elberling, B. (2018). Process-oriented modeling of a high arctic tundra ecosystem: Long-term carbon budget and ecosystem responses to interannual variations of climate. *Journal of Geophysical Research: Biogeosciences*, 123(4), 1178-1196. https://doi.org/10.1002/2017JG003956

Zscheischler, J., Mahecha, M. D., Avitabile, V., Calle, L., Carcalhais, N., Ciais, P., & Reichstein, M. (2017). Reviews and syntheses: An empirical spatiotemporal description of the global surface - Atmosphere carbon fluxes: Opportunities and data limitations. *Biogeosciences*, 14, 3685-3703. https://doi.org/10.5194/bg-14-3685-2017

## SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section.

**How to cite this article:** Virkkala A-M, Aalto J, Rogers BM, et al. Statistical upscaling of ecosystem CO₂ fluxes across the terrestrial tundra and boreal domain: Regional patterns and uncertainties. *Glob Change Biol.* 2021;00:1-20. https://doi.org/10.1111/gcb.15659
PRIMARY RESEARCH ARTICLE

# Statistical upscaling of ecosystem CO₂ fluxes across the terrestrial tundra and boreal domain: Regional patterns and uncertainties

Anna-Maria Virkkala<fnref n="1" />,<fnref n="2" /> ID | Juha Aalto<fnref n="1" />,<fnref n="3" /> ID | Brendan M. Rogers<fnref n="2" /> ID |

Torbern Tagesson<fnref n="4" />,<fnref n="5" /> ID | Claire C. Treat<fnref n="6" /> ID | Susan M. Natali<fnref n="2" /> ID | Jennifer D. Watts<fnref n="2" /> ID |

Stefano Potter<fnref n="2" /> ID | Aleksi Lehtonen<fnref n="7" /> ID | Marguerite Mauritz<fnref n="8" /> ID | Edward A. G. Schuur<fnref n="9" /> ID |

John Kochendorfer<fnref n="10" /> ID | Donatella Zona<fnref n="11" />,<fnref n="12" /> ID | Walter Oechel<fnref n="11" />,<fnref n="13" /> ID | Hideki Kobayashi<fnref n="14" /> ID |

Elyn Humphreys<fnref n="15" /> ID | Mathias Goeckede<fnref n="16" /> ID | Hiroki Iwata<fnref n="17" /> ID | Peter M. Lafleur<fnref n="18" /> ID |

Eugenie S. Euskirchen¹⁹ ID | Stef Bokhorst²⁰ ID | Maija Marushchak²¹,²² ID | Pertti

J. Martikainen²² ID | Bo Elberling²³ ID | Carolina Voigt²²,²⁴ ID | Christina Biasi²² ID |

Oliver Sonnentag²⁴ ID | Frans-Jan W. Parmentier<fnref n="4" />,²⁵ ID | Masahito Ueyama²⁶ ID |

Gerardo Celis²⁷ ID | Vincent L. St. Louis²⁸ ID | Craig A. Emmerton²⁸ ID | Matthias Peichl²⁹ ID |

Jinshu Chi²⁹ ID | Järvi Järveoja²⁹ ID | Mats B. Nilsson²⁹ ID | Steven F. Oberbauer³⁰ ID |

Margaret S. Torn³¹ ID | Sang-Jong Park³² ID | Han Dolman³³ ID | Ivan Mammarella³⁴ ID |

Namyi Chae³⁵ ID | Rafael Poyatos³⁶,³⁷ ID | Efrén López-Blanco³⁸,³⁹ ID | Torben

Røjle Christensen³⁹ ID | Min Jung Kwon⁴⁰ ID | Torsten Sachs⁴¹ ID | David Holl⁴² ID |

Miska Luoto<fnref n="1" /> ID

<footnote n="1">¹Department of Geosciences and Geography, Faculty of Science, University of Helsinki, Helsinki, Finland</footnote>

<footnote n="2">²Woodwell Climate Research Center, Falmouth, MA, USA</footnote>

<footnote n="3">³Weather and Climate Change Impact Research, Finnish Meteorological Institute, Helsinki, Finland</footnote>

<footnote n="4">⁴Department of Physical Geography and Ecosystem Science, Lund University, Lund, Sweden</footnote>

<footnote n="5">⁵Department of Geosciences and Natural Resource Management, Copenhagen University, Copenhagen, Denmark</footnote>

<footnote n="6">⁶Alfred Wegener Institute Helmholtz Center for Polar and Marine Research, Potsdam, Germany</footnote>

<footnote n="7">⁷Natural Resources Institute Finland, Helsinki, Finland</footnote>

<footnote n="8">⁸University of Texas at El Paso, El Paso, TX, USA</footnote>

<footnote n="9">⁹Center for Ecosystem Science and Society, Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, USA</footnote>

<footnote n="10">¹⁰Atmospheric Turbulence and Diffusion Division of NOAA's Air Resources Laboratory, Oak Ridge, TN, USA</footnote>

<footnote n="11">¹¹San Diego State University, San Diego, CA, USA</footnote>

<footnote n="12">¹²University of Sheffield, Sheffield, UK</footnote>

<footnote n="13">¹³University of Exeter, Exeter, UK</footnote>

<footnote n="14">¹⁴Research Institute for Global Change, Japan Agency for Marine-Earth Science and Technology, Yokoama, Japan</footnote>

<footnote n="15">¹⁵Carleton University, Ottawa, ON, Canada</footnote>

<footnote n="16">¹⁶Dept. Biogeochemical Signals, Max Planck Institute for Biogeochemistry, Jena, Germany</footnote>

<footnote n="17">¹⁷Department of Environmental Science, Shinshu University, Matsumoto, Japan</footnote>

<footnote n="18">¹⁸School of the Environment, Trent University, Peterborough, ON, Canada</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>

<footnote>© 2021 The Authors. *Global Change Biology* published by John Wiley & Sons Ltd.</footnote>