Variation in CO₂ and CH₄ fluxes among land cover types in heterogeneous Arctic tundra in northeastern Siberia

Arctic tundra is facing unprecedented warming, resulting in shifts in the vegetation, thaw regimes, and potentially in the ecosystem–atmosphere exchange of carbon (C). However, the estimates of regional carbon dioxide (CO₂) and methane (CH₄) budgets are highly uncertain. We measured CO₂ and CH₄ fluxes, vegetation composition and leaf area index (LAI), thaw depth, and soil wetness in Tiksi (71^∘ N, 128^∘ E), a heterogeneous site located within the prostrate dwarf-shrub tundra zone in northeastern Siberia. Using the closed chamber method, we determined the net ecosystem exchange (NEE) of CO₂, ecosystem respiration in the dark (ER), ecosystem gross photosynthesis (Pg), and CH₄ flux during the growing season. We applied a previously developed high-spatial-resolution land cover map over an area of 35.8 km² for spatial extrapolation. Among the land cover types varying from barren to dwarf-shrub tundra and tundra wetlands, the NEE and Pg at the photosynthetically active photon flux density of 800 µmol m⁻² h⁻¹ (NEE₈₀₀ and Pg₈₀₀) were greatest in the graminoid-dominated habitats, i.e., streamside meadow and fens, with NEE₈₀₀ and Pg₈₀₀ of up to −21 (uptake) and 28 mmol m⁻² h⁻¹, respectively. Vascular LAI was a robust predictor of both NEE₈₀₀ and Pg₈₀₀ and, on a landscape scale, the fens were disproportionately important for the summertime CO₂ sequestration. Dry tundra, including the dwarf-shrub and lichen tundra, had smaller CO₂ exchange rates. The fens were the largest source of CH₄, while the dry mineral soil tundra consumed atmospheric CH₄, which on a landscape scale amounted to −9 % of the total CH₄ balance during the growing season. The largest seasonal mean CH₄ consumption rate of 0.02 mmol m⁻² h⁻¹ occurred in sand- and stone-covered barren areas. The high consumption rate agrees with the estimate based on the eddy covariance measurements at the same site. We acknowledge the uncertainty involved in spatial extrapolations due to a small number of replicates per land cover type. This study highlights the need to distinguish different land cover types including the dry tundra habitats to account for their different CO₂ and CH₄ flux patterns, especially the consumption of atmospheric CH₄, when estimating tundra C exchange on a larger spatial scale.

1 Introduction

It is uncertain whether the Arctic tundra is a sink or a source of atmospheric carbon (C). The current estimates suggest a sink of 13–110 Tg C yr⁻¹, but their uncertainty range crosses the zero balance (McGuire et al., 2012; Virkkala et al., 2021). Improving these estimates is vital, because the Arctic tundra covers a vast area of 7.6 million km² (Walker, 2000) that is experiencing substantial warming (IPCC, 2013; Chen et al., 2021). Warming can alter C exchange, either amplifying or mitigating climate change through ecosystem–atmosphere interactions. Some local-scale studies suggest that the Arctic tundra is shifting from a small sink to a source of C (Webb et al., 2016; Euskirchen et al., 2017). It is likely that the climate change response of the ecosystem carbon dioxide (CO₂) sink strength and methane (CH₄) emissions, whether an increase or a decrease, depends on site-specific changes in thawing, wetness, temperature, and vegetation (McGuire et al., 2018). Dynamics of C exchange need to be quantified across the arctic habitats to improve the upscaling of arctic CO₂ and CH₄ balances and to monitor how ecosystems respond to environmental changes.

The uncertainty in the Arctic C balance estimates arises from the sparse and uneven observation network, which provides poor support for model-based spatial extrapolation (cf. McGuire et al., 2018; Virkkala et al., 2021; Kuhn et al., 2021). On a local scale, landscape heterogeneity and the related difficulty of mapping the spatial distribution of habitats and their C fluxes add to this uncertainty (McGuire et al., 2012; Treat et al., 2018; Saunois et al., 2020). Furthermore, year-to-year variations in seasonal features, particularly the timing of spring, summer temperatures, and snow depth have been found to cause substantial variation in the annual CO₂ and CH₄ balances (Aurela et al., 2004; Humphreys and Lafleur, 2011; Zhang et al., 2019). Fine-scale spatial heterogeneity in soil water saturation, thaw depth, vegetation characteristics, and soil organic content is typical of the tundra landscape (e.g., Virtanen and Ek, 2014; Mikola et al., 2018; Lara et al., 2020). These factors control CO₂ and CH₄ exchange, and on an annual scale, tundra wetlands typically act as net CO₂ sinks while upland tundra areas have a close-to-neutral CO₂ balance (e.g., Marushchak et al., 2013; Virkkala et al., 2021). While tundra wetlands are substantial sources of CH₄, dry tundra acts as a small sink or small source of atmospheric CH₄ (Bartlett and Harriss, 1993; Kuhn et al., 2021).

Mineral soil tundra barrens, however, have been found to have high consumption rates of atmospheric CH₄, which is due to the high-affinity methane oxidizing bacteria (Emmerton et al., 2014; Jørgensen et al., 2014; D'Imperio et al., 2017; Oh et al., 2020). These bacteria can utilize atmospheric CH₄ as energy source at low atmospheric concentrations, opposite to the low-affinity methane oxidizers that require higher CH₄ concentrations and occur in wetlands (e.g., Oh et al., 2020). A modeling exercise that introduced the high-affinity methanotrophy for mineral-rich soils resulted in a doubling of the circumpolar soil CH₄ sink above 50^∘ N compared to previous estimates (Oh et al., 2020). Thus, distinguishing dry and wet tundra with their moisture and vegetation characteristics is crucial when mapping C exchange within the tundra biome. Treat et al. (2018) tested spatial resolution requirements for such mapping on a landscape level and found that a 20 m pixel size captured the spatial variation in a reasonable manner, while a coarser resolution resulted in underestimation of both the landscape-scale CO₂ uptake and CH₄ emissions. In addition, understanding the spatial heterogeneity of ecosystem C exchange substantially improves analyses of eddy covariance (EC) measurements that, while in principle representing spatially integrated fluxes, may provide biased gas flux balances in a highly heterogeneous source and sink environment, as the spatial integration of EC involves non-uniform weighting of the surface elements that contribute to the measured flux (Tuovinen et al., 2019).

The aim of this study was to assess the spatial patterns and magnitudes of CO₂ and CH₄ fluxes within heterogeneous prostrate dwarf-shrub tundra in Tiksi, located in northeastern Russia. Growing season fluxes of CO₂ (ecosystem net exchange, photosynthesis, and respiration) and CH₄ were determined using the chamber method to answer the questions: (i) what is the magnitude of these fluxes in different land cover types and (ii) how do they depend on vegetation characteristics and soil wetness? In addition, we extrapolated the plot-level measurements in space and compared them with the ecosystem-level data measured with the EC technique.

2 Materials and methods

2.1 Study site

The study site is located near the Tiksi Observatory in Sakha (Yakutia) (see Uttal et al., 2016), northeastern Russia (71.5943^∘ N, 128.8878^∘ E), 500 m inland from the Laptev Sea coast and, on average, 7 m above sea level (Fig. 1a). The area belongs to the middle-Arctic prostrate dwarf-shrub tundra subzone (Walker, 2000) and has continuous permafrost. At the end of the growing season, the maximum thaw depth is ca. 40 cm (Mikola et al., 2018). The climate in Tiksi is defined by cold winters and cool summers. The long-term mean annual temperature and mean annual precipitation were −12.7 ^∘C and 232 mm, respectively, during the climate normal period 1981–2010. The growing season lasts about 3 months, the soils typically freeze in the end of September, and the permanent snow falls in October and thaws in June (AARI, 2018).

Figure 1 (a) Location of the study area in Tiksi, Yakutia, Russia. (b) Land cover map with the chamber flux measurement points (dark dots) and the EC mast (cross), and photos of the land cover types: (c) lichen tundra with barren ground and patches of vegetation, (d) dwarf-shrub tundra, (e) bog, (f) wet and dry fen, (g) graminoid tundra, and (h) meadow by the stream. See Tuovinen et al. (2019) for the EC footprint climatology.

Bedrock is alkaline, resulting in high plant species richness. Vegetation consists of mosses, lichens, grasses, sedges, prostrate dwarf-shrubs such as willows (Salix spp.), dwarf birch (Betula nana) and Diapensia lapponica, and forb species (Table 1). The average height of dwarf-shrub species is 4–6 cm, and the leaf area index (LAI) of vascular plants reaches up to 1 m² m⁻² in the fen and meadow habitats with graminoid vegetation (Juutinen et al., 2017). The land cover at the site has been classified a priori and mapped based on a combination of field inventories and high-spatial-resolution satellite images (Mikola et al., 2018). The a priori land cover types (LCTs) consist of wet fen, dry fen, graminoid tundra, bog, meadow at the stream bank, dwarf-shrub tundra, and lichen tundra that consists of barren ground with rocks, sand, and patches of vegetation (Table 1, Fig. 1c–h; for a closer view see Appendix Fig. A1). The depth of organic layer is negligible in lichen tundra and a few centimeters in dwarf-shrub tundra, meadow, and graminoid tundra. In bog, dry fen, and wet fen, the organic layer depth is at least the maximum depth of the active layer, ca. 30–40 cm. Soil organic content can reach ca. 40 % in tundra wetlands (Mikola et al., 2018). A section of the wet and dry fen within the EC footprint area is disturbed by vehicle tracks that create open water surfaces, and there is also an area of eroded bare-peat surface on a dry fen.

Table 1 Soil and vegetation characteristics of the land cover types (LCT) and their proportions in the EC impact area (90 % of the cumulative footprint).

¹ Combines the bare and lichen tundra LCTs defined in Juutinen et al. (2017), Mikola et al. (2018), and Tuovinen et al. (2019). ² Proportion within the 90 % coverage of the mean EC footprint area during the growing season of 2014 (Tuovinen et al., 2019).

2.2 CO₂ and CH₄ flux measurements

Fluxes of CO₂ and CH₄ were measured using static chambers equipped with a fan and set on pre-installed collars of 50 cm × 50 cm in area. The measurement points (collars) were set to cover the heterogeneity in land cover, and in each study year, there were 1–4 measurement points per LCT (Table 2). Most of the data were collected during a study campaign during 15 July–16 August 2014 (12 collars). The growing season had started earlier due to a warm period, and the daily mean air temperature had stayed over 5 ^∘C since 5 July (Fig. 2 and Tuovinen et al., 2019). The net ecosystem exchange of CO₂ (NEE) and ecosystem respiration of CO₂ in the dark (ER) were measured using transparent and opaque chambers (transparent chamber covered with a hood), respectively, allowing the partitioning of ecosystem gross photosynthesis (Pg) and ER. Fluxes of CH₄ were determined from closures of both transparent and opaque chambers, but because there was no difference between them when performed consecutively, the data from opaque chamber measurements were used for flux calculations. In addition, CH₄ fluxes were measured during shorter campaigns in 2012, 2013, 2016, and 2019 (Table 2). These data also included the plots disturbed by vehicle tracks and an eroded bare-peat surface, which were measured in 2019.

Figure 2 (a) Air temperature accumulation with the threshold T_air and surface T_soil of 0 ^∘C, (b) seasonal dynamics of NDVI in the study area (16 d MODIS data), (c) weekly means of soil temperature at a depth of 5 cm in wet fen and (d) in lichen tundra, (e) water-table level relative to the ground surface in wet fen, and (f) LCT-specific means of thaw depth in the measurement collars in 2014. Rocks in the ground prevented detecting the thaw depth of lichen tundra.

In 2012 and 2013, four air samples were taken from the chambers using syringes. The samples were stored in glass vials prior to the analysis. First, a vial was flushed with the sample and then filled to overpressure. The samples were analyzed for CH₄ concentration using a TSVET 500 M gas chromatograph (Chromatek, Russia) with a flame ionization detector at the laboratory of the Voeikov Main Geophysical Observatory within a month of sampling. Each measurement was accompanied by calibration using standard gas mixtures with known CH₄ concentrations (the NOAA2004 scale). The vials were tested prior to the field sampling using a standard gas: after 2 weeks, the vials were still overpressurized and the sample CH₄ concentrations were within ±3 ppb of the initial standard gas concentration. Since July 2014, CH₄ and CO₂ concentrations inside the chambers were recorded every second during closures of about 5 min using a gas analyzer (DLT-100, Los Gatos Research, Inc., San Jose, CA, USA) (see Appendix Fig. A2 for examples).

Gas fluxes between the ecosystem and the atmosphere were calculated from the linear concentration change in the chamber head space over time, accounting for temperature, volume, and atmospheric pressure. The concentration change during each chamber closure was evaluated visually to determine the closure start time and to remove cases showing nonlinearity due to leaks, ebullition, or saturation. The first data points were generally neglected when determining the slope of concentration change over time. The linearity was screened also on the basis of the coefficient of determination of the fit (R²). An R² greater than 0.9 was required, except for nearzero fluxes. There were a few ebullition cases at the vehicle track measurement points that had only sparse or no vegetation cover, and those measurements were included in the final data. When determining the NEE fluxes measured using the transparent chamber, data were screened for variation in photosynthetically active photon flux density (PPFD), measured during the chamber closure, and the flux measurement was rejected if the PPFD variation exceeded 100 µmol m⁻² s⁻¹ during the closure.

The fluxes of CO₂ and CH₄ were also measured by the micrometeorological EC method, which provides continuous data of the atmosphere–biosphere fluxes averaged on an ecosystem scale. The EC system consisted of a three-dimensional sonic anemometer (USA-1, METEK GmbH, Elmshorn, Germany), a closed-path CH₄ analyzer (RMT-200, Los Gatos Research, Inc., San Jose, CA, USA), and a closed-path CO₂ $/$ H₂O analyzer (LI-7000, LI-COR, Inc., Lincoln, NE, USA). The fluxes were calculated as 30 min averages and processed using standard methods (Aubinet et al., 2012). The EC measurement system and the post-processing procedures have been presented in more detail by Tuovinen et al. (2019).

Supporting meteorological measurements including air temperature (T_air) (HMP, Vaisala), soil temperature (T_soil) (IKES, Nokeval), PPFD (PQS1, Kipp & Zonen), and water table level relative to the ground surface (WT) (8438.66.2646, Trafag) were collected by a Vaisala QML datalogger as 30 min averages. We also present meteorological data for the period 2011–2019 to relate the conditions during the measurement campaigns in 2012, 2013, 2014, 2016, and 2019, to longer-term variations.

2.3 Vegetation and topographic wetness index

On a site level, vegetation and soil characteristics were inventoried in plots assigned into a systematic grid outside the area covered by the gas flux measurement points in 2014 (see Juutinen et al., 2017; Mikola et al., 2018). The projection cover (%) of plant species and species groups, and the mean canopy height of each species group were recorded. Eight species groups were included in the inventory: Sphagnum mosses, feather mosses, brown mosses, dwarf shrubs, Betula nana, Salix species, forbs, and graminoids. A subset of the plots was harvested, and vascular plant leaves were scanned to determine the one-sided LAI for empirical relationships between LAI and %-cover and canopy height, which were used to estimate the LAI in the collars (see Juutinen et al., 2017). In the collars, the projection cover and canopy height of each species group were recorded weekly during the gas flux measurement campaign during 15 July–16 August 2014. Because there were no observational vegetation data for the other years than 2014, the green chromatic coordinate (GCC) calculated from digital photographs was used as a proxy for the amount of green above-ground vascular plants (e.g., Richardson, 2019). The GCC was calculated from the digital numbers of red (R), green (G), and blue (B) color channels as the proportion of green in the RGB images, GCC = G $/$ (R + G + B), of the vegetation inside the collars. The photographs were taken at the time of measurements. We determined an empirical relationship between LAI and GCC by using a data set of harvested plots with digital photographs and measured LAI data (n=91). For the LAI estimation, we used a linear relationship (R²=0.46, p<0.001) between LAI and GCC determined using the entire data set (see Fig. A3 for the data and equation).

Table 2 Measurement periods, measured fluxes (CH₄, ER, NEE), and number of measurement points in each land cover type (LCT) across the study years.

^* Measured only on 30 May 2016.

To quantify the potential soil wetness at each measurement point, we calculated the mean topographic wetness index (TWI) based on a 2 m spatial resolution digital elevation model (Mikola et al., 2018). To characterize differences between growing seasons as manifested by vegetation greenness, the MODIS Normalized Difference Vegetation Index (NDVI) with 16 d temporal and 500 m spatial resolution was calculated for a circular area with a 300 m radius from the flux tower using Google Earth Engine (Gorelick et al., 2017). NDVI was derived for 2011–2019 to place the measurement years in the context of year-to-year variation in weather and plant growth.

2.4 Data analyses

When examining the role of the LCTs in CO₂ and CH₄ exchange, we applied the land cover classification presented by Mikola et al. (2018). The data collected in 15 July–16 August 2014 were used for examining gas exchange in relation to the variation in LAI, GCC, WT, and TWI among the collars. The light-normalized Pg and NEE at PPFD = 800 µmol m⁻² s⁻¹ (Pg₈₀₀ and NEE₈₀₀, respectively), were estimated by fitting a hyperbolic response function of CO₂ vs PPFD utilizing the ER and NEE flux data:

$$\begin{array}{}\text{(1)}& \mathrm{NEE}=\mathrm{ER}-{\mathrm{Pg}}_{\max }\times \mathrm{PPFD}/(\mathit{\beta }+\mathrm{PPFD}),\end{array}$$

where Pg_max is the asymptotic maximum of photosynthetic CO₂ uptake rate, and β is the half-saturation PPFD. Fluxes of CH₄ are expressed as temporal averages for each collar. We used a sign convention where a positive flux means net release to the atmosphere and a negative flux denotes net uptake by the ecosystem. Fluxes of CH₄ measured over all study years, 2012–2019, were averaged for each LCT.

Regression analyses were used to test the relationships between gas flux estimates and vascular LAI, GCC, WT, and TWI. All CH₄ flux data from the years 2012–2014, 2016, and 2019 were used to quantify the mean growing season CH₄ flux for each LCT and examine the relationship between CH₄ and GCC and TWI. To find the main factors and gradients in the plant community, gas flux, and environmental variables data measured in the flux collars in 2014, we performed a detrended correspondence analysis (DCA) of the species group data with a post-hoc fit of environmental variables, including gas fluxes, WT, LAI, GCC, elevation, and thaw depth as supplementary variables. The DCA was performed on logarithmically transformed, centered species data (species or species groups) using Canoco 5 (Šmilauer and Leps, 2014).

We compared the LCT-specific flux estimates obtained from the chamber measurements with the estimates based on EC measurements during the same period (15 July–16 August 2014). Partitioning of the EC-based CO₂ fluxes to Pg and ER and the estimates of Pg₈₀₀ and NEE₈₀₀ were calculated similarly to those derived from the chamber data (Eq. 1). The EC flux data were classified into five wind sectors (30–125, 125–185, 185–239, 239–310, 310–360^∘) based on the mean EC flux footprint, modeled for the growing season of 2014 by Tuovinen et al. (2019). The sectors distinguished areas dominated by different LCTs, especially tundra heaths and wetlands, and similarly those with a large and small vascular LAI. For each sector, the footprint-weighted areal proportions of LCTs and mean vascular LAI were derived from the high-spatial-resolution LCT and LAI maps (Mikola et al., 2018). For this comparison, sector averages of Pg₈₀₀, ER, NEE₈₀₀, and CH₄ flux were calculated from the chamber data by weighting the LCT-specific flux estimates with the above-mentioned LCT proportions in each sector. Because there were no chamber measurement points within graminoid tundra, we applied wet fen (for CO₂) and dry fen (for CH₄) flux estimates for the graminoid tundra based on the observed similarities in LAI and soil wetness, respectively. Overall, graminoid tundra can be considered part of the fen continuum in terms of soil characteristics (notably high organic content) and CH₄ exchange (Mikola et al., 2018; Tuovinen et al., 2019).

Finally, to synthesize the CO₂ and CH₄ exchange variability across the tundra, we upscaled the LCT-specific NEE_800, Pg₈₀₀, ER, and CH₄ flux (2014 data) estimate averages to the 35.8 km² area surrounding our study site, for which a LCT map was produced by Mikola et al. (2018).

3 Results

3.1 Environmental conditions

In 2014, when we collected most of the flux data, temperature sum accumulation (with a 0 ^∘C T_air threshold) took place at a near-average rate during the thaw period (the period when soil surface temperature was continuously above 0 ^∘C), but the spring and mid-growing season were warmer than on average (Fig. 2a). The average air temperature was 15 ^∘C during the gas flux measurements. Accordingly, the MODIS NDVI showed an early start of greening (Fig. 2b), and vegetation development had already started at the beginning of the measurement period. In 2011–2019, which includes all the CH₄ measurement years, the thaw period lasted for 74–124 d, creating a temperature sum range of 642–1003 ^∘C days (Fig. 2a). Surface soils thawed between 28 May and 9 July and froze again between 21 September and 1 October. Among the observation years, the years 2012 and 2019 had notably longer and warmer thaw periods than the other years. The driest habitat, lichen tundra, with least snow accumulation, thawed 10–15 d earlier than the other habitats and had an about 3 ^∘C higher soil temperature at the depth of 5 cm than wet fen (Fig. 2c–d). Water table level, measured at a wet fen location, showed only subtle interannual variation (Fig. 2e). In 2014, the active layer depth, measured close to the collars during the flux measurement period, was deepest in the end of August, reaching ca. 40 cm in wet fen but remaining <30 cm in the dwarf-shrub tundra (Fig. 2f). Lichen tundra had rocks underneath the loose surface layer, which made it impossible to measure the actual thaw depth.

3.2 Exchange of CO₂ and CH₄

Among different LCTs, the estimates of Pg₈₀₀ varied from about 5 mmol m⁻² h⁻¹ in lichen tundra to 22 and 27 mmol m⁻² h⁻¹ in wet fen and meadow, respectively (Table 3a). Pg₈₀₀ was strongly and positively correlated with the vascular plant LAI and the greenness index GCC (Fig. 3). There was also a positive correlation between Pg₈₀₀ and both WT and TWI, possibly because the highest LAI occurred at the wet fen and meadow plots. However, the TWI values for the two meadow plots located on an elevated bank of the stream were disproportionately high in relation to the WT at these plots, probably because of insufficient spatial accuracy or an artifact of the digital elevation model. Ecosystem respiration was highest in the two meadow plots, on average 18 mmol m⁻² h⁻¹. The relationship between ER and LAI was weaker than between Pg₈₀₀ and LAI (Fig. 3). NEE₈₀₀ varied from about zero in the lichen tundra plots to a net CO₂ uptake of 16 mmol m⁻² h⁻¹ in the meadow and wet fen plots (Table 3a). NEE₈₀₀ was more tightly linked to Pg₈₀₀ than to ER and it was correlated with LAI, GCC, WT, and TWI (Fig. 3).

Figure 3 Variation in the estimates of Pg₈₀₀, ER, and NEE₈₀₀ (Eq. 1) and collar means of CH₄ fluxes in relation to variation in the collar means of LAI, GCC, WT, and TWI during 6 July–16 August 2014. Error bars denote the standard error of estimate (n=15 or 16). Fitted regression lines and adjusted coefficients of determination ($R_{\mathrm{adj}.}^{\mathrm{2}}$) are included for the significant linear relationships. The two meadow plots were not included in the TWI regressions.

Table 3 (a) Means, medians, and standard deviations (sd) of LCT specific CO₂ and CH₄ fluxes in 2014 calculated from collar specific estimates (CO₂) or seasonal means (CH₄). There were 15 or 16 data points per each collar and 1–3 collars per LCT (see Table 2). (b) Proportions of LCTs in landscape totals of area, and Pg₈₀₀, NEE₈₀₀, and CH₄ fluxes based on the LCT means (part a). (c) LCT specific means, medians, and standard deviations of CH₄ fluxes based on multiyear data (n is number of the observations).

¹ Graminoid tundra contribution estimated using values for wet fen (CO₂) and dry fen (CH₄), ² Area-weighted mean. ³ Water not shown.

There was substantial consumption of atmospheric CH₄ in the barren tundra, where the mean of all measured fluxes was −0.018 mmol m⁻² h⁻¹, and in the vegetated lichen tundra with a mean of −0.005 mmol m⁻² h⁻¹ (Table 3c, Figs. 4 and 5). Minor consumption occurred in the bog, meadow, and dwarf-shrub tundra plots (mean fluxes from −0.0002 to −0.001 mmol m⁻² h⁻¹), while efflux to the atmosphere was observed in the dry fen (mean 0.04 mmol m⁻² h⁻¹) and wet fen plots (mean 0.17 mmol m⁻² h⁻¹). Fluxes were also high in the eroded bare-peat plot within the dry fen habitat and the vehicle-track plots in wet fen (Table 3c).

Variation among the plot means (Fig. 3 for 2014) and LCT means (Fig. 5 for all years) of CH₄ flux was related to WT, and CH₄ emissions occurred when TWI was greater than 4. The two meadow plots that showed net consumption of CH₄ had an unrealistically high TWI relative to their WT (see above and Figs. 3 and 5). Variation in CH₄ fluxes was incoherently related to the variation in LAI and GCC because of the high emission cases in plots with little vegetation, including the wettest wet fen, vehicle-track, and bare-peat plots (Fig. 5).

Figure 4 Instantaneous CH₄ fluxes in each LCT. The data are a composite of all study years. Barren surfaces are indicated among the lichen tundra data. The eroded bare-peat and vehicle-track plots (crosses) are plotted as part of the dry fen and wet fen data, respectively. Note that the panel groups have different y-axis scales.

Figure 5 LCT mean (± SE) CH₄ fluxes in relation to the corresponding mean (± SE) TWI (excluding the meadow) and mean GCC. Data from years 2012–2019; see Table 3c for the number of measurements.

The DCA ordination of species groups with a post-hoc fit of environmental variables (elevation, WT, thaw depth, LAI, GCC, and CO₂ and CH₄ exchange) showed that species distributed along a moisture gradient. Axis 1 explained 49 % of the variation in the species data and distinguished the wet and dry LCTs (Fig. 6). Graminoids and brown mosses occurred in the wet end of the gradient, while evergreen dwarf-shrubs, Betula nana, and lichens occurred in the dry end of it. The barren plot (the other lichen tundra plot) with its negligible vegetation differed most from the other plots. Axis 2 explained additional 14 % of the variation in the species data (Fig. 6). The supplementary variables WT, vascular plant LAI, thaw depth, GCC, Pg₈₀₀, NEE₈₀₀, and CH₄ flux correlated positively with Axis 1 having post-hoc correlations (r) of 0.6–0.9, as derived from the DCA-weighted correlation matrix, while plot's elevation and ER had positive correlations with Axis 2 (r=0.8 and 0.4, respectively).

Figure 6 DCA ordination diagram based on species (species groups) data from the measurement collars in 2014. The explained variation in the species data is indicated for the axes. The scores are indicated for species groups (crosses), sample plots (small circles), and post-hoc fits of the supplementary variables (blue arrows). Land cover types of the sample plots are indicated in gray and the plots assigned to each LCT are circled (dashed ellipses).

In both the southern (125–185^∘) and southwestern (185–239^∘) wind sectors, vegetation mainly consisted of graminoids, as the LCTs dry fen, wet fen, graminoid tundra, and meadow comprised 80 % of the total EC footprint-weighted area (Fig. 7a). The northern sector (310–360^∘) was characterized by lichen tundra and bare ground that accounted for 68 % of the footprint-weighted LCT areas, while all the other LCTs covered less than 18 % in total. The other wind direction sectors had more even LCT distributions. The differences between the sectors were similar in the EC-based and spatially weighted chamber-based averages of CO₂ exchange (Fig. 7b–d). Both Pg₈₀₀ and NEE₈₀₀ were largest in the southern and south-western sectors and clearly smallest in the barren–lichen tundra-dominated sector in the north. The chamber-based estimates of CO₂ exchange were, however, lower: on average, Pg₈₀₀ was 57 %, ER 93 %, and NEE₈₀₀ 44 % of the mean EC-based fluxes among the wind direction sectors.

The southern and southwestern wind sectors with abundant dry and wet fens and graminoid tundra had clearly the largest CH₄ fluxes (Fig. 7f). The estimate based on chamber measurements was 30 % and 50 % larger than the mean EC-based flux in the east sector (dominated by dry fen and bog) and south sector (dominated by dry fen and wet fen), respectively. In contrast, the chamber-based estimate was smaller than the EC flux for the other sectors, which were dominated by graminoid tundra, lichen tundra, and barren ground. Both the EC- and chamber-based measurements showed consumption of atmospheric CH₄ in the northernmost sector, of which barren ground and lichen tundra covered 50 % and 20 %, respectively. The mean EC flux was three times the chamber-based estimate.

Figure 7 Footprint-weighted mean contribution of each LCT to the EC measurements divided into wind direction sectors (a), and comparison of EC and chamber-based sector means of CO₂ exchange (Pg₈₀₀, ER, and NEE₈₀₀) (b–d) vascular plant LAI (e), and CH₄ fluxes (f). The chamber-based data are weighted by the LCT proportions shown in panel a. All data were measured in 2014. Map of LAI (Tuovinen et al., 2019) and the LAI measured in the collars were used to estimate the EC – and chamber – related sector means, respectively, in panel (e).

Within the extended study area of 35.8 km², the LCT-weighted mean NEE₈₀₀ was −4.6 mmol m⁻² h⁻¹ (uptake relative to the atmosphere). The corresponding mean Pg₈₀₀ was 11 mmol m⁻² h⁻¹, and the mean CH₄ flux was 0.05 mmol m⁻² h⁻¹ (Table 3a). Relative to their spatial cover (28 % in total), wet and dry fens were disproportionally important for the landscape-level Pg₈₀₀, NEE₈₀₀, and CH₄ emissions, because the fens contributed 47 % of total Pg₈₀₀ and 74 % of NEE₈₀₀, and were the largest source of CH₄ (Table 3b). Consumption of CH₄ by barren and lichen tundra, dwarf-shrub tundra, and meadow tundra soils contributed −9 % of the CH₄ balance, and the barren ground dominated the sink. It should be noted that these data represent the growing season conditions when both the CH₄ emissions and consumption of atmospheric CH₄ were at their highest during the year due to high temperatures, thawed soils and active vegetation.

4 Discussion

The studied tundra site in Tiksi in northeastern Siberia has heterogeneous land cover, which is reflected as equally heterogeneous CO₂ and CH₄ exchange. We found that the LAI of vascular plants was a robust predictor of Pg₈₀₀ and NEE₈₀₀ across the LCTs. On the one hand, the tundra wetlands had a disproportionate role in the landscape-level CO₂ uptake capacity. The fens also dominated the landscape's CH₄ emissions. On the other hand, our results highlight the substantial CH₄ consumption of atmospheric CH₄ within the dry tundra areas, particularly in barrens. The CH₄ consumption by dry tundra contributed −9 % of the total CH₄ balance estimated for this landscape from the data collected during the growing season. This finding is in agreement with other studies and suggests distinguishing non-vegetated dry tundra habitats when upscaling plot-scale CH₄ fluxes (Table 4). In Tiksi, the barren was characterized by sand and rocks underlain by schists (Fig. A1). The consumption of CH₄ was smaller if the sand and stones were partly covered with vegetation (Figs. 5 and A1).

Table 4 Summary of reported CH₄ fluxes in mineral soil dry tundra.

¹ Mean estimated from a figure. ² Minimum and maximum estimated from a figure. ³ Three day measurement. ⁴ Estimated from EC measurements with a statistical model.

The land-cover-categorical approach serves to distinguish the basic features of spatial variation in CO₂ and CH₄ fluxes, and the extreme ends of the moisture and vegetation gradients from barren to wet fen are clearly distinguishable, also in terms of CO₂ and CH₄ exchange (Fig. 6). Overall, microrelief, moisture gradient, vegetation types, and ecosystem functions are connected. For instance, barren areas are wind swept and thus have minimal snow accumulation, while in wet depressions snow accumulation further increases soil moisture (Fig. 6; Callaghan et al., 2011). The spatial extrapolation of fluxes, however, is here sensitive to a small number of chamber measurement points as there is large within-LCT variation in fluxes and LAI, as observed in the wet fen and meadow data. Moreover, the LCTs share common features and form a continuum as shown by the DCA ordination (Fig. 6). Mikola et al. (2018) used a larger soil and vegetation data set from Tiksi and also found that the neighboring LCTs overlapped in terms of soil and plant attributes. Despite the limited number of observations, our conclusions drawn from the chamber data are corroborated by the temporally matching section of EC data, which show high similarity to the chamber data (Fig. 7). Furthermore, the statistical analysis of EC data by Tuovinen et al. (2019) showed that it is possible to find significant differences between different LCT categories representing high and low CH₄ emitters and CH₄ sinks. However, for spatial modeling of ecosystem functions, maps of key variables, such as LAI and WT, that drive CO₂ and CH₄ exchange would be preferable to categorical LCT classification (Räsänen et al., 2021).

The spatial pattern of the growing season Pg₈₀₀ and NEE₈₀₀ was strongly related to the corresponding pattern of the LAI of vascular plants (Fig. 3). Hence, the abundance of graminoid (Cyperaceae and Poaceae) vegetation was associated with a large NEE₈₀₀, which varied from near zero in lichen tundra up to −25 mmol m⁻² h⁻¹ in wet fen. Ecosystem respiration had a smaller role than Pg in determining NEE, but we note that our data cover only a section of the growing season with warmer temperatures and half- to full-grown vegetation. The importance of ER is likely to be different when considering the full annual balance (e.g., Hashemi et al., 2021). While our data represent only the growing season, a similar relationship has also been found between the annual NEE and LAI at a tundra site with a mixture of wet and dry tundra in northeastern Europe (Marushchak et al., 2013), in a multi-site EC study in Alaskan tundra (McFadden et al., 2003), in Canadian low arctic tundra wetlands (Lafleur et al., 2012), and across tundra sites (Street et al., 2007; Shaver et al., 2007).

The magnitude of Pg₈₀₀ and NEE₈₀₀ in the fen and meadow plots of this study were similar to the maximum Pg and NEE found in a tundra wetland in Seida in northeastern Europe (Marushchak et al., 2013), at low tundra wetland sites in eastern Canada (Lafleur et al., 2012), and at a wetland-dominated but more continental site (with an equally long growing season) in northeastern Siberia (van der Molen et al., 2007). The vegetation and Pg₈₀₀ of lichen tundra and dwarf-shrub tundra in our study resembled those observed within the polygon rim habitat of the polygon tundra in the Lena River delta, while those of meadow, dry fen, and wet fen resembled the wet polygon center habitats (Eckhardt et al., 2019). In our study, the spatial variation of ecosystem respiration resulted from the variation in vascular plant LAI, soil organic content, and water saturation: the highest ER occurred in the mineral soil meadow plots with a high LAI, suggesting substantial autotrophic respiration and likely deep rooting and large root biomass contributing to ecosystem respiration (Fig. 3).

Our chamber-based estimate of the average CH₄ flux within the 35.8 km² upscaling area was 0.05 mmol m⁻² h⁻¹, which is close to 0.04 mmol m⁻² h⁻¹ obtained by Tuovinen et al. (2019), who combined EC data with footprint modeling to statistically determine LCT group-specific CH₄ fluxes. Within this upscaling area, we estimate that 28 % of the area emitted CH₄, while the other habitats either consumed atmospheric CH₄ (barren and lichen tundra, dwarf-shrub tundra, meadow) or were close to neutral (bog) relative to the atmosphere (Fig. 4, Table 3a–b). The relationship between the vascular plant LAI and CH₄ flux was confused by the occurrence of large CH₄ fluxes in plots with little or no vegetation. Those fluxes were observed at the wettest fen plot and the bare-peat and vehicle track plots (Figs. 4–5). A high LAI, a high WT, and a high CH₄ emission systematically co-occurred in wet fen (Fig. 6). In addition, in the bare-peat and vehicle-track plots, erosion or anthropogenic disturbance may have created CH₄ flux hotspots due to permafrost scars, water saturation, and recently thawed organic matter (e.g., Bubier et al., 1995; McCalley et al., 2014; Wickland et al., 2020). These are small-scale landscape features, while on a larger scale our data encourage applying indices of wetness and vegetation as a means of CH₄ flux upscaling in a tundra environment.

The recognition of CH₄ consuming tundra habitats is important for accurately estimating the net CH₄ balance of tundra. The substantial uptake of atmospheric CH₄ by lichen tundra (here a mixture of bare ground and sparse vegetation) in Tiksi was inferred by Tuovinen et al. (2019) based on a source allocation analysis of EC data: the average flux of the consuming area was estimated at −0.03 mmol m⁻² h⁻¹, which corresponded to −22 % of the total upscaled CH₄ balance. In this study, the average seasonal CH₄ flux was −0.02 mmol m⁻² h⁻¹ in the barren tundra and an order of magnitude lower in meadow and dwarf-shrub tundra. This difference between the estimates likely originates from the LCT-weighting and the small sample of the chamber-based data and, in general, demonstrates the inherent sensitivity involved in upscaling of fluxes of opposite direction.

High consumption of atmospheric CH₄ in barrens is associated with the high affinity methanotrophs (Emmerton et al., 2014; Jørgensen et al., 2014; D'Imperio et al., 2017; St Pierre et al., 2019). Our summary of the CH₄ fluxes in mineral-rich dry tundra (Table 4) shows that the consumption rates in Tiksi are higher than those observed elsewhere. This may be due to a local feature associated with the parent material of the ground. Similar rates, however, have been recorded at other dry tundra sites with little or no vegetation. For instance, on Disko Island, Greenland, which consists of similar land cover types to Tiksi, CH₄ uptake by bare ground was 0.005–0.01 mmol m⁻² h⁻¹ during the growing season, while a mean uptake of 0.003–0.004 mmol m⁻² h⁻¹ was observed in dry tundra heath (D'Imperio et al., 2017). These consumption rates associated with tundra barrens and high-affinity methanotrophs can be even higher than those measured on north-boreal forest soils (e.g., Lohila et al., 2016).

5 Conclusions

Our results provide new observations of carbon exchange for the prostrate dwarf-shrub tundra sub-zone, which covers a substantial area of the Arctic. These data augment the knowledge on the functional diversity, namely the distribution of different land cover types and their emission factors, across the vast Arctic tundra and will lend support to bottom-up and top-down extrapolations across the Arctic. Graminoid vegetation that favored the wet and moist habitats, such as wet fens, was characterized by large CO₂ uptake and CH₄ emissions. In addition, our data support the observation of notable consumption of atmospheric CH₄ in barren tundra that has substantial coverage across the Arctic. The heterogeneity of landscape and the related large spatial variability of CO₂ and CH₄ fluxes observed in this study encourage to monitor the Arctic sites for changes in habitat type distribution. Such changes can include the forming of meadows and wet fens and appearance of new vegetation communities, such as erect shrubs, that benefit from warming-induced changes in thaw depth and soil wetness. The spatial extrapolation based on a small number of measurement points involves inherent uncertainty but still allowed us to identify key relationships between CO₂ and CH₄ fluxes and vegetation and moisture features, which can be utilized in more robust upscaling studies that make use of EC measurements.

Appendix A

Figure A1 Examples of the barren (a, c) and lichen tundra (b, d) plots with close views (c, d). Vegetation consists of lichens Flavocetraria sp., Thamnolia sp., Alectoria sp., dwarf-shrubs Dryas octopetala, Vaccinium vitis-idaea, Cassiope tetragona, and graminoids and forbs such as Carex spp. and Polygonum viviparum.

Figure A2 Examples of gas concentration variations during chamber closures measured using the gas analyzer (DLT-100, Los Gatos Research, Inc., San Jose, CA, USA). The examples represent lichen tundra, barren, and wet fen.

Figure A3 Relationship between GCC and vascular plant LAI in the harvested plots. LCTs are indicated with symbols. In the LCT-specific regressions (not shown), the coefficient of determination ($R_{\mathrm{adj}.}^{\mathrm{2}}$) was lowest for dry fen (0.06) and highest for wet fen (0.54). Regression slopes varied from 8.3 for dry fen to 17.8 for the combined graminoid tundra and meadow LCT.