Figure 2 . Circumpolar percentage coverage of the five adapted BAWLD terrestrial land cover types (Boreal Forests, Non-permafrost Wetlands, Permafrost Bogs, Dry Tundra, and Tundra Wetlands) used for ecosystem-based upscaling of GHG flux budgets in this study. Note that these maps show the distributions across the full BAWLD domain as presented by Olefeldt et al (2021), not the more limited extent of the RECCAP2 permafrost BAWLD domain used in this study.
The land cover mean GHG flux (Fjx) were obtained for each of the five terrestrial land cover classes after homogenising and analysing three comprehensive GHG flux datasets: Virkkala et al. (2022) for CO2 fluxes; Kuhn et al. (2021a) for CH4 fluxes; and Voigt et al. (2020) for N2O fluxes. Additional data was extracted from literature for Boreal Forest N2O fluxes (Schiller and Hastie, 1996; Simpson et al., 1997; Kim and Tanaka, 2003; Morishita et al. 2007; Matson et al. 2009; Ullah et al. 2009; Köster et al. 2018a), since the N2O flux dataset from Voigt et al. (2020) does not cover Boreal Forest ecosystems. These datasets comprise roughly 1000in-situ growing-season and annual observations (including multiple observations from some sites) of terrestrial fluxes obtained from more than 200 sites using chamber (for CH4, N2O, and CO2), diffusion (for CH4 and CO2), and eddy covariance (for CO2 and CH4) methods. The growing season length was defined as June to August (90 days) for the tundra and permafrost bogs sites, and May to September (150 days) for the Boreal Forests and Non-Permafrost Wetlands. The CO2 dataset comprises year-round measurements of net ecosystem exchange (NEE), which we used to calculate growing season and annual NEE. Average fluxes were calculated based on 93 sites and 403 observations for growing season NEE and 54 sites and 222 observations for annual NEE. The CH4 and N2O datasets provide growing-season measurements based on 98 sites and 458 observations of CH4 exchange and 47 sites and 91 observations of N2O exchange. For sites with incomplete growing season measurements, we multiplied average daily fluxes to the length of the growing season. Annual CH4 fluxes were estimated assuming that growing season emissions accounted for 64% of annual emissions (Treat et al. 2018), except for boreal forests were we assumed growing season emissions accounted for 100% of annual emissions as the sites averaged net CH4 growing season uptake and available data for winter season fractions only covers CH4-emitting ecosystems (Treat et al. 2018). Our Boreal Forest annual estimate should therefore be considered conservative. Annual N2O fluxes were estimated assuming that growing season emissions accounted for 50% of annual emissions as reported in Voigt et al. (2020). For all three GHGs, only sites with no record of large-scale upland hillslope abrupt thaw disturbance in the metadata were included in the flux estimates to avoid double-counting emissions from upland hillslope abrupt thaw (see methodology for disturbances). However, although scarce, we included other disturbed sites in our CO2 estimates to account for ecosystem CO2 losses following disturbances and their different successional stages (e.g., 4 sites reporting thermokarst; Virkkala et al. 2022). Sites from the above-mentioned GHG flux datasets were classified into one of the five terrestrial land cover classes using the metadata provided in each of the datasets. More details on how ecosystem flux upscaling was performed can be found in the supporting information.
While the focus of this study is the period 2000-2020, we include allin-situ measurements obtained between 1991 and 2020 in order to overcome the limited amount of flux measurements in some of the ecosystems and therefore ensure adequate spatial representation of ecosystem fluxes. A separate analysis of decadal CO2fluxes from 1991 to 2020 revealed no differences, suggesting that the extension of time series to 1991 does not impact our findings (Table A2).

2.4 GHG fluxes from inland waters

Similarly to the method used to calculate GHG emissions from terrestrial land cover types, GHG fluxes from inland waters were calculated by upscaling mean GHG fluxes from lakes and rivers (see below) using the estimated surface area of these aquatic classes from the BAWLD classification (Olefeldt et al. 2021), adjusted to the study region (see supplementary Table A1 for estimated aerial extent of inland waters).

GHG budgets for rivers

Atmospheric riverine GHG fluxes were calculated in different ways for each GHG, depending on available source data, and when possible scaled across the region using riverine area from the permafrost region (0.12 x 106 km2), reported in BAWLD.
Estimates of river and stream CO2 flux were calculated from gridded monthly flux data estimated by Liu et al. (2021; https://doi.org/10.5061/dryad.d7wm37pz9; Dryad) from river water dissolved CO2 pressure and gas transfer velocity. We combined the monthly fluxes from the start of May to the end of October, assuming that this corresponds to the ice-free season (when water-to-air gas transfer can occur). This time extent (184 days) is nine days longer than the duration Liu et al. (2022) cite as the mean ice-free period for Arctic lakes (175 days). This data is delivered as unprojected global grids with a 0.0083 degree resolution (which is ca 1*0.2 km pixels in the high Arctic). The global grids were clipped to the extent of BAWLD and then reprojected to an equal area grid at 100*100 m resolution. Calculations from this data yields a total stream and river area of 0.069*106 km2, and a total flux of 94 Tg CO2-C yr-1. Assuming the mean river flux (1,370 g C m-2 yr-1) can be scaled also to smaller streams and rivers, we applied the area of streams and rivers in BAWLD (0.12*106 km2). Because spatially explicit estimates of uncertainty are not available, we report a coefficient of variation proportional to the global uncertainty reported by Liu et al. (2022). Riverine CH4 emissions were determined using the mean CH4 diffusive flux reported in the MethDB (Stanley et al. 2016). Stanley et al. (2016) found that diffusive CH4 emissions did not statistically differ across latitudes and scaled global river CH4emissions using one mean value. Given the limited number of reported CH4 fluxes for rivers in the Arctic (e.g. Zolkos et al. 2020), we used the same approach as Stanley et al. (2016) and applied a global mean diffusive flux of 135 mg CH4m-2 d-1 to the river area. Because there are few studies that measure CH4 emissions upon ice-out, we applied for CH4 a conservative estimate that 17% of annual fluxes occur during the ice-free period (Denfield et al. 2018; consistent with the approach by Liu et al. 2022). Ebullition was not included for river CH4 emission estimates due to few available measurements in the literature for this region (Stanley et al. 2016). Estimates of river N2O flux were derived from gridded annual N2O flux estimated by a mechanistic mass balance model developed globally for inland waters by Maavara et al. (2019). These data was reprojected from an original 0.5 degree unprojected grid to an equal area grid at 1 km resolution and clipped to the BAWLD extent. As the original lake and river surface area was not known, no correction of inland water surface area was made. Uncertainties for river GHG budgets were determined using the standard error and coefficient of variance reported by Liu et al. (2022), Stanley et al. (2016) and Maavara et al. (2019), respectively, for CO2, CH4, and N2O.

GHG budgets for lakes

CH4 fluxes (diffusion and ebullition) were extracted from the BAWLD-CH4 aquatic ecosystem dataset and classified based on classes (yedoma lakes, peatland ponds, and glacial/post-glacial organic poor lakes and ponds) and sizes, from large (> 10 km2) to midsize (0.1 to 10 km2) to small lakes (< 0.1 km2) (Kuhn et al. 2021a; total area = 1.255 106 km2; Table A3). Notably, no minimum size for lakes was considered in the BAWLD dataset, as the dataset gives an estimate of the overall area covered by lakes in each size-class (Olefeldt et al. 2021). Conceptually, any area which is likely to be inundated >50% of the growing season period (long term average) is considered part of the lake land cover classes. Ice-free days were determined based on averages of reported ice-free days for each lake type and this information was used to determine ice-free season fluxes (supplementary Table A1). In addition to ice-free emissions, spring ice-out emissions (i.e. winter contribution) were considered to be 23% of the annual total (Wik et al. 2016).
Estimated lake CO2 fluxes were compiled from multiple available sources based on a literature search made in May 2022 (Humborg et al. 2010; Rocher-Ros et al. 2017; Karlsson et al. 2013; Sepulveda-Jauregui et al. 2015; Pelletier et al. 2014; Rasilo et al. 2014; Korteliane et al. 2006) and are summarised in Table A4). The studies report lake CO2 fluxes as mean flux values for various binned lake surface areas. We took these averages and grouped them by the lake size classes included in BAWLD (<0.1, 0.1-10, >10 km2). We found no statistical differences in fluxes between the size groups and thus used one mean lake CO2 flux to scale across the year and the region (315 ± 196 mg C m-2 d-1). We applied the same number of ice-free days used to scale lake CH4emissions (ice-free days reported in the literature for each lake class).
To estimate lake fluxes of N2O, gridded global data of annual flux from Lauerwald et al. (2019) were used. This estimate is based on the nitrous oxide (N2O) emission model developed by Maavara et al. (2019) and the HydroLAKES database and was reprojected from an original 0.5 degree unprojected grid to an equal area grid at 1 km resolution and clipped to the BAWLD extent. As the original lake and river surface area was not known, no correction of inland water surface area was made. Uncertainties for lake N2O were determined using the coefficient of variance reported for regions north of 50 deg latitude in Lauerwald et al. (2019).

2.5 Disturbances - fires and abrupt thaw

Monthly GHG fire emissions were extracted for the study region from the Global Fire Emission Database version 4s (GFED; van der Werf et al. 2017). The GFED4s spans from 1997-2016 and estimates of burned areas are based on remote sensing data at a spatial resolution of 0.25 degrees (van der Werf et al. 2017). GHG emissions in the GFED4s are derived from the multiplication of burned area and fuel consumption per unit burned area, the latter being the product of modelled fuel loads per unit area and combustion completeness. For our purpose, we extracted mean annual GHG emissions from burned areas for the period 2000-2016 and assumed similar rates for the period 2016-2020.
Localised, but widespread, disturbances associated with abrupt thaw are thought to contribute significantly to GHG emissions from permafrost (Abbott and Jones, 2015, Yang et al. 2018, Walker et al. 2019, Turetsky et al. 2020; Holloway et al. 2020, Marushchak et al. 2021, Runge et al. 2022). Abrupt thaw includes thawing processes that affect permafrost soils in periods of days to several years (Grosse et al. 2011), and is typically associated with thermokarst and thermoerosion processes that lead to the formation of hillslope erosional features (thaw slumps, thermo-erosion gullies and active layer detachments), thermokarst lakes, and thermokarst wetlands (i.e., collapse scar bogs and fens). We report abrupt thaw areas and derived annual CO2 and CH4 emissions using the inventory-based abrupt thaw model by Turetsky et al. (2020), in which atmospheric emissions are estimated for three generalised types of abrupt thaw terrains: mineral-rich lowlands, upland hillslopes, and organic-rich wetlands. In the abrupt thaw model, abrupt thaw areas are based on synthesised field observations and remote sensing measurements. GHG emissions from abrupt thaw were synthesised for each ecosystem state within each abrupt thaw type from the literature (ca. 20 published papers). The abrupt thaw model was initialised for a historical assessment period (1900-2000) to provide the model with a spin up and prevent the regional carbon fluxes starting at zero at the beginning of the dynamic measurement period. Thaw rates were generally in equilibrium with succession and recovery of surface permafrost during this initialization period. Changes in the area of each successional state were tracked over time by multiplying initial starting areas by transition rates. Estimates of abrupt thaw GHG emissions following the historical assessment period were done by increasing rates of abrupt thaw through time. This increase in thaw rate was prescribed to follow the average output of ‘permafrost-enabled’ land surface models, all of which were forced by atmospheric climate anomalies from the Community Climate System Model 4 (CCSM4) Earth system model under an RCP8.5 projection. For our purpose, we ran the abrupt thaw model for the period 2000-2020 and extracted cumulative CO2 and CH4 emissions from active and stabilised abrupt thaw features, and derived annual fluxes for each abrupt thaw terrain for the time period 2000-2020. We used the reported uncertainty ranges of ± 40% on the upland hillslope areas, ± 30% on the mineral-rich lowland areas, and ± 35% on the organic-rich wetland areas as in Turetsky et al. (2020). Additional details on the inventory model can be found in Turetsky et al. (2020). Since GHG datasets that we used for ecosystem upscaling partly account for abrupt thaw and to prevent double counting GHG fluxes, CO2 and CH4 fluxes from abrupt thaw were added as a sub-flux (not added to the total) of terrestrial and inland water land cover fluxes and their contribution to the total GHG budget is discussed. Due to the lack of in situ observations of abrupt thaw impacts on N2O fluxes in the used datasets, no N2O budget is presented for abrupt thaw.

2.6 Lateral fluxes and geological emissions

Lateral C and N fluxes from riverine transport and coastal erosion (i.e. DOC and DON losses from the permafrost region to the ocean) are taken from Terhaar et al. (2021), representative for all land north of 60° N. They estimated riverine lateral fluxes for the six largest Arctic rivers (Mackenzie, Yukon, Kolyma, Lena, Ob, Yenisei) from the Arctic Great River Observatory (ArcticGRO) dataset and extrapolated to the entire Arctic catchment. Emissions from coastal erosion were calculated by multiplying spatially resolved estimates of coastal erosion rates by estimates of C content in coastal soils provided in Lantuit et al. (2012).
Estimates of geological emissions of CH4 (from subsurface fossil hydrocarbon reservoirs) are taken from an upscaled circumpolar permafrost region estimate for gas seeps along permafrost boundaries and lake beds made by Walter Anthony et al. (2012). We note that there is some risk of double counting such fluxes, especially in sites where eddy covariance flux towers may have unknowingly been placed close to seeps of geological CH4 emissions. No separate estimates of geological emission for CO2 or N2O are available for the permafrost region. For CO2, the full global geological emissions are estimated to 0.16 Pg CO2-C yr-1 (Mörner and Etiope 2002).

3 Results and Discussion

3.1 Net GHG exchange from terrestrial land cover types

Terrestrial ecosystems represented a decadal-scale sink for CO2, and source for CH4 and N2O (Table 1, Fig. 3). The mean annual CO2 flux was a net sink, but could not be distinguished from CO2 neutral when the 95% confidence interval was considered (-339.6 (-835.5, 156.3) Tg CO2-C y-1). The broad uncertainty interval can be attributed both to the large natural variability in CO2 fluxes across sites and to the heterogeneity of ecosystem types included in each of the land cover classes defined in the BAWLD classification. Boreal Forests and Non-permafrost Wetlands were CO2sinks (-270.3 and -69.4 Tg CO2-C y-1, respectively) while Tundra Wetlands and Permafrost Bogs were close to neutral (-2.7 and -0.05 Tg CO2-C y-1, respectively). Dry Tundra was the only ecosystem type classified as an annual ecosystem CO2 source (2.9 Tg CO2-C y-1), but the very broad uncertainty range (-147.6, 153.5 Tg CO2-C y-1) indicates low confidence in the sign of this flux. Terrestrial ecosystems were overall a net sink of CO2 during the growing season (-1611 (-2148, -1074) Tg CO2-C gs-1), with the strongest sink in the boreal forest (-1034 (-1305, -763) Tg CO2-C gs-1) (Table 2).
Annual terrestrial CO2 flux budgets have been reported for high-latitudes in recent papers using different upscaling approaches. While closely related due to overlap in flux data, a higher NEE uptake is reported by both Virkkala et al. (2021) and Watts et al. (2023) (-419 (95% CI of -559 to -189) Tg CO2-C y-1 and -601 (standard error of ± 1138) Tg CO2-C y-1, respectively). However the estimated NEE uptakes for the permafrost region solely are weaker, with an uptake of -181 (-305, 32) Tg CO2-C y-1 and -230 (± 22) Tg CO2-C y-1, respectively). The difference between the later NEE uptakes and our results relates to the subset of data included in the analyses (exclusively eddy covariance tower fluxes in Watts et al. (2023)), the different years covered in the analyses (Virkkala et al. 2021: 1990-2015, Watts et al. 2023: 2003-2015), the different spatial extents, and the upscaling approach applied (Arctic Terrestrial Carbon Flux Model (TCFM-Arctic) in Watts et al. (2023), and statistical upscaling in Virkkala et al. (2021)). Both of these studies as well as the previous RECCAP synthesis (1990-2006, McGuire et al. 2012) report the tundra as a weak CO2 sink (-13 (-81, 62); -16 (±84–270); and -16 (-42, 10) Tg CO2-C y-1, respectively) although they also show that annual tundra budgets cannot be distinguished from CO2 neutral when taking into account the uncertainty range. Dry Tundra CO2 budget was also identified as a source of 10 (-27, 47) Tg CO2-C y-1 in McGuire et al. (2012).
Our estimated annual net CH4 source of 25.6 (14.7, 36.4) Tg CH4-C y-1 from terrestrial ecosystems (Table 1) was largely driven by emissions from Non-permafrost Wetlands (20.6 (14.3, 26.9) Tg CH4-C y-1). As in Treat et al. (2018), Non-permafrost Wetlands emitted more than Tundra Wetlands. Annual CH4flux estimates for Tundra Wetlands (3.3 (2.7, 3.9) Tg CH4 y-1) and Dry Tundra (2.1 (-0.4, 4.5 Tg CH4-C y-1) were in the lower range from the previous estimates provided in McGuire et al. (2012), in which the tundra was estimated to release 11 (0, 22) Tg CH4-C y-1 (between 1990 and 2006). Our growing season CH4 budget was a source of 16 (8.6, 23.3) Tg CH4-C gs-1 (Table 2) with Non-permafrost Wetlands contributing 83%. All terrestrial ecosystems except Boreal Forests were net CH4 emitters. Boreal Forests were a net sink of CH4 (-1.1 (-2.3, 0.0) Tg CH4-C gs-1). Our CH4annual budget was lower than the ones estimated for the northern high latitude wetlands (>45°N) at 31, 32, and 35 Tg CH4-C y-1 (depending on wetland distribution maps) by Peltola et al. (2019) and 38 Tg CH4-C y-1 by Watts et al. (2023). However, our CH4 growing season budget estimate was higher than the budget based on 93 observations presented in Treat et al. (2018) except for the Tundra Wetlands where they remain within the same range. Despite their large spatial coverage, Dry Tundra was a small source of CH4 during the growing season (1.4 (-0.3, 2.9) Tg CH4-C gs-1), although the low end of the CI suggests that it could remain a sink. More measurements from these drier ecosystems are needed.
Our N2O annual budget estimate of 0.55 (-0.03, 1.1) Tg N2O-N y-1 (Table 1) suggests that terrestrial ecosystems were a N2O source, although the uncertainty range around N2O fluxes extends from a small sink to a larger source. These high uncertainties partly relate to the limited number of observations of N2O fluxes (47 sites and 91 observations), which only includes growing-season observations. Our estimated annual N2O budget is within the range of the one previously reported by Voigt et al. (2020)(0.14-1.27 Tg N2O-N y-1 median-mean-based estimate). In our study, Dry Tundra was the largest N2O source (0.23 (0.04, 0.42) Tg N2O-N y-1). Boreal Forests were the second largest N2O source (0.14 (-0.01, 0.30) Tg N2O-N y-1) due to their large area, although their fluxes per unit area were small (Table A5, 52.43 ug N2O m-2d-1). Although they occupy a small portion of the landscape (5%), Permafrost Bogs were the largest N2O emitters per unit area (Table A5, 645.14 ug N2O m-2 d-1) and their contribution to the regional budget was 18%. The estimate for Permafrost Bogs includes emissions from barren peat surfaces, where vascular plants are absent - surfaces previously identified as N2O hot spots in the Arctic due to ideal conditions for N2O production (Repo et al. 209; Marushchak et al., 2011; Gil et al. 2017). A challenge remains regarding the mapping of Permafrost Bogs and barren ground and integration within land cover classifications. Therefore, we did not differentiate between vegetated and non-vegetated Permafrost Bog areas when upscaling. N2O emissions from Tundra Wetlands were negligible (0.01 (0.00, 0.02) Tg N2O-N y-1), which can be explained by the lack of nitrate supply as an N2O precursor in reduced conditions and reduction of N2O to N2 during denitrification when the water table is high (Butterbach-Bahl et al. 2011; Voigt et al. 2017). Recent observations not included in the N2O review dataset (Voigt et al 2020) show that wetlands may also function as net N2O sinks in the Arctic (Schulze et al. 2023).
Table 1. Greenhouse gas (GHGs - CO2, CH4, and N2O) budget for the permafrost region based on ecosystem upscaling. Negative GHG emissions represent an uptake while positive emissions represent a release. GHG emissions from terrestrial ecosystems are reported as mean fluxes with 2.5 and 97.5% confidence intervals (CI). GHG emissions from inland waters and fires are reported with 5 and 95% CI. GHG emissions from abrupt thaw are reported with +-40% uncertainty range. *these fluxes are estimated using the abrupt thaw model from Turetsky et al (2020) and are considered as additive to the total for these categories (to avoid double counting of fluxes). **includes CO2, CH4 and lateral fluxes.