Growing-season variability and environmental controls of CH₄ and N₂O fluxes across high-latitude terrestrial ecosystems

Rapid warming and accelerating ecological change across high-latitude (>66° N) terrestrial environments are reshaping greenhouse gas dynamics. However, spatiotemporal patterns and environmental drivers of CH₄ and N₂O fluxes in high-latitude ecosystems remain poorly understood.

In this study we describe the seasonal and spatial variation in CH₄ and N₂O fluxes across key subarctic and boreal environmental gradients. A secondary aim is to examine the environmental drivers underlying these flux patterns. We measured fluxes during the snow-free season of 2024 and 2025 in northern Finland, starting early June (spring snowmelt) and ending late September (first snowfall). We conducted 1272 chamber measurements at 144 locations spanning eight high-latitude land cover classes: bare tundra, dwarf-shrub tundra, deciduous forest, evergreen forest, fen, bog, permafrost bog and tundra wetland.

Bare tundra, dwarf-shrub tundra, deciduous forest and evergreen forest sites acted as net CH₄ sinks (mean flux of −0.16 to −0.05 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, depending on land cover class) and as both N₂O sources and sinks (−0.69 to 1.69 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). Fen, bog, permafrost bog and tundra wetland sites were net CH₄ sources (mean flux of 0.93 to 27.81 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) and net N₂O sinks (−0.89 to −6.53 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). Highest mean CH₄ fluxes occurred in fens (23.4 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) and permafrost bogs (27.8 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) which often coincided with substantial N₂O sinks (mean flux of −6.53 in fens and −5.40 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ in permafrost bogs). Seasonal trends were present for CH₄ but were land cover-dependent, with both sinks and sources typically peaking during July–August. Land cover class, soil temperature, and soil pH had the strongest effects on CH₄ fluxes while variation in N₂O was best explained by soil pH and land cover class. However, a significant portion of the variability in N₂O remained unexplained.

These results demonstrate that strong spatial heterogeneity and pronounced seasonal dynamics during the snow-free period generate substantial variability in CH₄ and N₂O fluxes across high-latitude ecosystems. Our findings suggest that monthly measurements of insitu gas fluxes and soil conditions in combination with accurate land cover mapping are essential for modelling and upscaling of GHG fluxes. This is especially important in high-latitude regions, where rapid warming is altering both ecosystem distributions and soil conditions with significant consequences for current and future GHG budgets.

1 Introduction

High-latitude (>66° N) ecosystems exhibit highly dynamic greenhouse gas (GHG) fluxes shaped by the interactions of soil, microbial, and vegetation processes (Whalen, 2005; Riutta et al., 2007; Olefeldt et al., 2013; Martins et al., 2017). Climate change can alter all these factors by reshaping abiotic conditions such as precipitation, temperature, and snow cover. Changing climate also drives biological shifts, including changes in microbial communities (Deslippe et al., 2012) and broader ecosystem transitions such as shrubification, paludification, and permafrost thaw (Myers-Smith et al., 2011; Biskaborn et al., 2019). Arctic regions store a major fraction of global soil carbon reserves (Batjes, 2016; Hugelius et al., 2020) and are warming faster than the global average (Rantanen et al., 2022), underscoring their importance in the global carbon cycle. However, the interplay of complex GHG processes and the heterogeneous landscape makes both measurements and modelling particularly challenging.

Soil carbon storage and GHG dynamics vary widely across high-latitude ecosystems (Hugelius et al., 2024). Peatlands (wetlands with peat soils), for example, function as major long-term carbon sinks due to ongoing peat accumulation, yet they are also the largest natural source of atmospheric methane (CH₄) because of their persistently anoxic soils (Saunois et al., 2025; Kübert et al., 2026). Most high-latitude wetlands can be further categorized into bogs, permafrost bogs and fens based on their hydrology and permafrost extent (Olefeldt et al., 2021). We have additionally included a tundra wetland category as both their ecology and hydrology significantly differ from the other categories (see Sect. 2.3 for detailed land cover classifications). In addition to wetlands, we also studied dry tundra and forest soils which were further categorized into bare tundra, shrub tundra, evergreen forests and deciduous forests (see Sect. 2.3). In this study, dry tundra soils specifically refer to non-wetland upland mineral soils with no tree cover. Although dry tundra and forest ecosystems store less carbon per unit area than peatlands, their vast extent, up to 80 % of the boreal-Arctic area (Raynolds et al., 2019; Olefeldt et al., 2021), means they exert a major influence on high-latitude carbon cycling through CH₄ uptake, biomass and soil carbon storage, and gross-primary production (GPP) (Virkkala et al., 2021; Voigt et al., 2023; Li et al., 2025; Jiang et al., 2025). In addition to CH₄ dynamics, nitrous oxide (N₂O) sources and sinks have been reported from high-latitude ecosystems, but their magnitudes, spatial patterns, and environmental drivers remain poorly understood (Chapuis-Lardy et al., 2007; Marushchak et al., 2011; Voigt et al., 2020). CH₄ and N₂O are the second and third most influential greenhouse gases after CO₂. Over a 100-year time horizon, CH₄ has a global warming potential (GWP100) that is 27–29 times greater than that of CO₂, whereas N₂O has a GWP100 of 273 (IPCC, 2023).

High latitude peatlands have a disproportionately strong effect on carbon dynamics due to their intense CH₄ emissions and large soil carbon stocks and because of that, they are overrepresented in many flux databases (Kuhn et al., 2021; Virkkala et al., 2026). Yet, high-latitude landscapes are heterogenous with a mosaic composition of wetlands, forests and tundra of which tundra soils are especially underrepresented in literature. Furthermore, it is typical to focus on peak growing season fluxes (Olefeldt et al., 2013, Kuhn et al., 2021) as it makes measurements and modelling more straightforward. However, this ignores the strong seasonal dynamics of the environmental drivers in high-latitude ecosystems which control GHG fluxes (Rautakoski et al., 2024; Kübert et al., 2026). Multi-month measurements across ecosystems are therefore required to resolve the seasonal progression of GHG fluxes and their underlying controls, which can shift substantially across environmental gradients and between spring, summer peak productivity, and autumn transition.

This study aims to characterize the seasonal and spatial variability of CH₄ and N₂O fluxes across key subarctic and boreal land cover classes, focusing on representative forest, wetland, and tundra ecosystems. A secondary aim is to examine the key environmental drivers and their role in shaping these flux patterns. To capture ecosystem-level heterogeneity, we classify these biomes into eight land cover classes following the BAWLD dataset framework (Olefeldt et al., 2021; Kuhn et al., 2021): deciduous and evergreen forests, bare tundra, dwarf-shrub tundra, fens, bogs, permafrost bogs and tundra wetlands (see Sect. 2.3 for details). Specifically, we address three research questions: (i) How do CH₄ and N₂O fluxes differ among the land cover classes? (ii) How do these fluxes develop seasonally during the snow-free period (June–September)? (iii) How do various environmental variables explain the variance in GHG fluxes? The analysis is based on 1273 in situ chamber measurements collected at 144 locations spanning eight high-latitude land cover classes at two research sites in Northern Finland (Kilpisjärvi and Pallas). The measurements cover the full duration of growing seasons of years 2024 and 2025, from spring snowmelt to late autumn.

2 Materials and methods

2.1 Location and area descriptions

The two study locations are in the Pallas-Yllästunturi national park in the Kittilä and Muonio municipalities in western Finnish Lapland, and the Kilpisjärvi region in the Enontekiö municipality, in northwestern Finnish Lapland. The Pallas study sites (67.9° N, 24.0° E, 250–570 m a.s.l.) are located at the edge of the northern-boreal and subarctic climatic zones with a long-term mean annual temperature of −1.0 °C and annual precipitation of 644 mm (Pirinen et al., 2012). It is characterized by a mosaic of northern boreal ecosystems consisting of mixed and coniferous forests (Picea abies, Pinus sylvestris, Betula pubescens), open fen wetlands and open tundra on the fells. The evergreen forest sites are only found in the Pallas area (Fig. 1, map 2) as the Kilpisjärvi sites are located above the coniferous forest limit. The Kilpisjärvi (Gilbbesjávri in Northern Sámi language) study sites are located in the subarctic climate zone (69.0° N, 20.8° E, 474–737 m a.s.l.) with a long-term mean annual temperature of −1.9 °C and annual precipitation of 487 mm (Pirinen et al., 2012). The Kilpisjärvi region lies at the ecotone between the boreal and tundra biomes and is characterized by dwarf-shrub dominated tundra highlands, Carex spp. dominated open fens and short stature mountain birch (Betula pubescens var pumila) forests at lower altitudes. The permafrost bog, tundra wetland, and deciduous forest sites occur only in the Kilpisjärvi research area (Fig. 1a and b). The Pallas area lies too far south and at too low an elevation to support permafrost bogs and extensive mountain birch forests, and comparable tundra wetlands are absent from the vicinity of the Pallas research area. The Pallas area is underlain by acidic Precambrian bedrock, resulting in low soil pH, whereas Kilpisjärvi is characterized by a geological mosaic of Paleozoic sedimentary rocks and acidic Precambrian bedrock (Simonen, 1990), yielding intermediate soil pH conditions. This geological contrast provides a robust framework for assessing the influence of pH on CH₄ and N₂O fluxes. Snowmelt and associated flooding and nutrient pulses typically occur in June, with rapidly increasing daylight hours progressing into a ∼2-month period of continuous daylight from late-May to late-July. Afterwards, a senescing period with occasional sub-zero temperatures occurs already during September. All measurement areas are relatively undisturbed with no ditching, recent logging or other major human influences.

Figure 1 Maps of the two research sites. Subarctic Kilpisjärvi site (1a and 1b) and boreal Pallas site (2). Colored dots indicate the measurement plots (n=144) and their eight different land cover classifications.

Figure 2 Representative in situ photographs of the eight land-cover classes included in the study, illustrating the range of environments sampled across the research areas (Photos: Johannes Niemi, Saku Mäki-Ikola).

2.2 Field campaign

The measurement campaign was conducted in two parts (years 2024 and 2025). The first part took place between 9 June and 1 October (2024) and the second part between 3 June and 21 September (2025). A total of 144 measurement plots were repeatedly sampled once a month (June–July–August–September) each year. Additional repeated measurements were taken from roughly a third of the research plots. In total, we conducted 1273 chamber measurements (n=1273 for each GHG) while continuously measuring microclimate and green vegetation cover at each measurement plot. Total chamber measurements were distributed between land cover classes as follows: shrub tundra (n=266, 20.9 %), fen (n=248, 19.5 %), evergreen forest (n=186, 14.6 %), deciduous forest (n=180, 14.2 %), bare tundra (n=106, 8.3 %), tundra wetland (n=106, 8.3 %), bog (n=100, 7.9 %) and permafrost bog (n=80, 6.3 %).

2.3 Land cover classification and site descriptions

The measurement plots were equally distributed between the Pallas (boreal) and Kilpisjärvi (subarctic) regions with three main ecosystem types of dry tundra, forests and wetlands, which were further divided between eight different key land covers: bare tundra, shrub tundra, evergreen forest, deciduous forest, tundra wetland, bog, fen and permafrost bog (Fig. 2). These land cover categorizations are based on the BAWLD dataset (Olefeldt et al., 2021), with some minor changes to more accurately describe our research sites. Specifically, the differences to the BAWLD categorization: our tundra wetland sites are not permafrost affected, and we separated dry upland soils into bare tundra and shrub tundra as we expected them to have different GHG flux profiles due to the lack of vegetation cover and shallower organic layer thickness in bare tundra plots.

The following descriptions are specifically for our research sites and based on our in situ observations. Shrub tundra plots (n=28) are located above the tree line and typically have a fairly thin dry organic soil layer (mean organic layer thickness, OLT=5.0 cm). Vegetation is generally dominated by Betula nana, Empetrum nigrum and Vaccinium myrtillus, although vegetation can occasionally be more diverse in wind protected valleys and creeks. Bare tundra plots (n=11) lack vascular plants and typically result from repeating freeze-melt cycles (cryoturbation) and persistent winds at higher altitudes. These plots may occasionally support small amounts of lichen and mosses and have a thin OLT (mean OLT=3.2 cm).

Evergreen forest plots (n=21) are located at the northern edge of the boreal zone with a tree cover dominated by Picea abies or Pinus sylvestris but can also have some deciduous trees mixed in, i.e. Betula pubescens or Populus tremula. Field layer is dominated by dwarf-shrubs (e.g. Vaccinium myrtillus, Vaccinium vitis-ideae) and various graminoids (e.g. Avenella spp. and Poa spp.). Evergreen forests had the highest OLT of the dry land cover classes (mean OLT=9.82 cm). The deciduous forest plots (n=20) are primarily located at the edges of the forest limit, and the tree cover consists solely of short growing Betula pubescens. var. pumila (mountain birch). The field layer is mixed with various dwarf-shrubs, graminoids (Avenella spp.) and forbs (e.g. Cornus suecica, Viola palustris). Most deciduous forests occur on well-drained glaciofluvial eskers and therefore have shallow organic layers (mean OLT=2.7 cm).

Fen plots (n=31) are minerotrophic, wet and treeless wetlands with peat depths up to ∼4.4 m. Field layer is dominated by Carex spp. and Sphagnum spp. Bog plots (n=11) are all of the dry and “hummocky” (permafrost free) wetland forms, i.e. if a wetland as a whole would be considered a fen, it might still have elevated dry Sphagnum hummocks with bog-type shrub vegetation which have considerably different greenhouse gas profiles compared to the wetter fen-type areas. Permafrost bog plots (n=10) are similar to the previously mentioned bogs but with permafrost cores that are in various stages of degradation. The permafrost bog research sites include two different peatlands representing varying degrees of permafrost degradation: one with largely intact but degrading palsas, the other with only remnants of old permafrost left. Permafrost extent was confirmed using soil cores to distinguish permafrost bogs from bogs without permafrost. Lastly, tundra wetland plots (n=12) are Carex spp., Trichoforum spp. and moss dominated wetlands (note: generally, no Sphagnum spp.). The tundra wetland sites are situated on hillsides, with small creeks running through them. They have a shallow organic soil layer compared to the other wetland types (mean organic layer thickness, OLT=11.2 cm). The upper soil layer consists of a mixture of fluvial minerogenic and organic deposits, including mud, gyttja, litter, peat, and roots. They have variable seasonal hydrology, i.e. water table was typically above the soil surface during spring snowmelt and autumn rains but below the surface during summer. Despite the low summer water table, soil moisture remained high throughout the measurement campaigns.

2.4 Environmental variables

The field layer vegetation cover of each plot was estimated by image analysis with Canopeo (Patrignani and Ochsner, 2015). Canopeo measures the relative number of green pixels from a photo and converts them into percentages, i.e. it acts as a measurement for green vegetation cover (“vegetation cover”). A photo was taken of each measurement plot every time chamber measurements were conducted. The pictures covered roughly a 1 m² square around the measurement plot and thus included vegetation from outside the chamber collar.

Microclimatic variables were measured with TOMST TMS-4 sensors (Wild et al., 2019), which are located 0.5–1 m from the chamber collar and measure continuously throughout the year. Loggers record microclimatic variables in 15 min intervals: soil temperature at a depth of −6 cm and soil moisture at depth of −14 cm. Each logger operated with a temperature resolution of 0.0625 °C and an accuracy of ±0.5 °C. For the soil temperature, the data quality control and processing procedures followed the approach outlined in Aalto et al. (2022). First the temperature time series were visually inspected to remove readings that clearly reflected sensor malfunction, such as abrupt shifts in soil moisture or temperature indicative of the logger being displaced. After the first screening of the data, the temperature sensors were then calibrated based on their individual response differences when stored under stable conditions (i.e. data logged while in storage indoors). Finally, any abnormal spikes in the data were filtered out, for example, a data point was flagged as erroneous if soil temperature change between consecutive time steps exceeded 5 °C or change in air temperature exceeded 10 °C.

For the soil moisture data, the quality control and preprocessing of the data followed the same outline as in Kemppinen et al. (2023). Similarly to the temperature time series, the raw soil moisture data was screened manually for clearly erroneous data where values were removed if the logger had been displaced, or the raw soil moisture value was <200, as it is outside the range of the sensor and usually indicates of poorly placed logger. The TOMST logger records raw soil moisture data based on electromagnetic pulses, producing time-domain transmission measurements (see Wild et al., 2019, for details). These were subsequently converted to volumetric water content (VWC, %) using the generic calibration function for unknown soils by Kopecký et al. (2021). The function for unknown soils was chosen, as it provided more realistic ranges of VWC values compared to soil-specific calibrations, however it does limit the maximum possible VWC to 60 %. This method has already been validated in Finland by Kemppinen et al. (2023), supporting its applicability to this research. After quality control and preprocessing, microclimate data was coupled with flux measurements by selecting the nearest 15 min logger reading corresponding to the flux measurement time. In cases where the closest time point contained erroneous or missing data, the next closest valid reading was used.

Both water table depth and soil moisture have been identified as key drivers for GHG fluxes (Moore and Dalva, 1993; Olefeldt et al., 2013; Rinne et al., 2018; Voigt et al., 2020; Virkkala et al., 2024). We chose to use soil moisture over water table depth for three reasons: (a) determining water table depth in mineral soils typically requires invasive methods that disturb soil structure and alter hydrological conditions (b) water table in dry tundra and forest sites is generally very deep, thus it has minimal impact on biogeochemistry occurring at the surface or at intermediate depths (c) we wanted model parameters that can be measured at all sites (i.e. both mineral soils and wetlands) for robust modelling across ecosystem gradients.

Soil samples for pH data were collected once during the measurement campaign from the organic layer underneath the vegetation (outside the chamber collars). The laboratory analysis was conducted by the Finnish Natural Resources Institute (LUKE) according to the ISO 10390 soil Quality-Determination of pH standards, in a suspension of soil in water (pH in H₂O).

2.5 Chamber measurements and flux calculations

The gas measurements were conducted as closed opaque chamber measurements with a duration of 3 min, taken with two portable gas analyzers (LI-COR-LI7810 and LI-COR-LI7820) that measure gas concentrations in 1 s intervals. Gas fluxes for CH₄ and N₂O were both measured simultaneously from the same chamber.

The chamber had a built-in fan and a wireless Bluetooth temperature sensor (Ruuvi tag pro) for measuring chamber temperatures. The volume of the chamber was ∼7.85 L (cylindrical, 20 cm×25 cm) and was placed on top of a steel collar. As the collar height changed between measuring plots (depending on how soft and deep the soil was), the height was always measured separately for each gas measurement to get accurate chamber volumes. For wetlands, the collars were installed into the soil roughly 7–14 d before the measurement campaign started to avoid disturbances in the peat that could cause ebullition or other artefacts in the flux. For mineral soils, the collar was removed between measurements to minimize changes in microclimatic conditions such as sheltering from wind. If unusually high fluxes were measured, they were repeated multiple times in situ to make sure the emissions were not caused by any influence of the measurer (i.e. ebullition due to the pressure on surrounding soil). If the repeat measurement showed similar fluxes, then the initial measurement was used (assuming linear change of gas concentration).

Fluxes were calculated as shown in Eq. (1), based on a linear fit between starting gas concentration at the start of the measurement and at the end of the measurement, utilizing the ideal gas law

$$\begin{array}{}\text{(1)}& \text{flux}={\displaystyle \frac{(ct-c\mathrm{0})\cdot V\cdot M\cdot p}{t\cdot R\cdot (T+\mathrm{273.15})\cdot A}}\end{array}$$

where $(ct-c\mathrm{0})/t$ is the slope (expressed as ppm s⁻¹), V is chamber volume, M is molar mass, p is atmospheric pressure, R is the ideal gas constant, T is chamber temperature and A is chamber area. Output of the formula without conversions is $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. Fluxes were converted from s⁻¹ to h⁻¹ (multiply by seconds in an hour; 3600 s) and from micrograms to milligrams for CH₄ (by dividing by 1000). Final units used for fluxes are mg ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ for CH₄ and $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ for N₂O.

Each chamber measurement was visually graphed to make sure the gas concentrations were changing linearly and to remove any possible measurement errors, such as anomalies at the start and end of the measurements that come from pressure changes or from handling the chamber (Pavelka et al., 2018). The same flux calculation method has been used in previous studies by Rissanen et al. (2023), Karjalainen et al. (2025) and Suominen et al. (2026).

Python script for flux calculations is also available: https://github.com/janivaltteri/djp-mcds/blob/main/djpmcds/mcds/fluxcalc.py (last access: 2 November 2026).

We followed Nickerson (2016) to calculate a minimum detectable flux (MDF) for chamber measurements which is an extended formula from Christiansen et al. (2015a) specifically designed for continuous high-resolution laser-based measurements. The calculation is based on the LI-COR sampling frequency (1 Hz) and measurement precision (0.6 ppb for CH₄ and 0.4 ppb for N₂O as reported by LI-COR), chamber volume, closure time and ideal gas law. For CH₄, 1.18 % of measurements were under the MDF (i.e. not confidently distinguishable from zero) while 10.6 % were below the MDF for N₂O.

Clear measurement errors (i.e. flat zero measurements that are sometimes caused by LI-COR malfunction) were removed, but all other measurements which had linear gas concentration change were kept in the data. Following quality control, 15 N₂O measurements (0.12 % of all N₂O measurements) and 9 CH₄ measurements (0.01 % of all CH₄ measurements) were excluded from the analyses. Non-linear fluxes are typically caused by the increasing concentration gradient at high gas concentrations, which slows down the gas diffusion rates from the soil to the chamber causing a log-curve (Pihlatie et al., 2013). In such cases, only the linear part of the measurement was used for the calculation to avoid underestimating fluxes. Chamber measurements that showed exponential concentration growth were typically repeated in situ until a linear trend was observed, as this kind of flux behavior is usually caused by influence of the measurer (Pavelka et al., 2018).

2.6 Statistical analysis and modelling

We used the lme4 package (Bates et al., 2015) in R (R: The R Project for Statistical Computing, 2026) to fit linear mixed-effects models examining the relationships between GHG fluxes and environmental variables. We used a joint model for both wetlands and dry land cover classes to investigate useful parameters for modelling across broad environmental gradients. The final model for CH₄ consisted of four fixed effects, three two-way interactions terms and a random effect. As the main effects, we included soil temperature (continuous), soil moisture (continuous), pH (continuous) and land cover (categorical, 8 levels). Interaction terms were land cover paired with each of the other fixed variables. The land cover interaction term was used to give the model flexibility for the direction of the GHG flux response, e.g. increasing temperature can have both positive and a negative effect on GHG flux depending on the land cover class. Measurement point identification was used as a random effect (random intercept) to account for the repeated measurements design (e.g. pseudoreplication). Green vegetation cover (continuous) and measurement year (categorical, 2 levels) were also tested as fixed effects but were dropped due to worse model performance.

The final model for N₂Oconsisted of five fixed effects, two two-way interactions and a random effect. As main effects we included soil temperature (continuous), pH (continuous), Green vegetation cover (continuous), year (categorical, 2 levels) and land cover (categorical, 8 levels). Interaction terms were land cover paired with pH and vegetation cover. Measurement point identification was used as a random effect (random intercept). Soil moisture, soil moisture-land cover interaction and soil temperature-land cover interaction were also tested as fixed effects but were dropped due to worse model performance.

Organic layer thickness (OLT) was considered as a potential explanatory variable but was not included in the final models as it was strongly correlated with land-cover class and soil moisture. Because these variables already captured much of the variation associated with OLT, its inclusion would have increased collinearity and model complexity without providing substantial additional explanatory power.

Model validation was done through Akaike information criterion (AIC) and model parameter effects sizes were described with partial eta-squared (η²p). The GHG flux data was skewed and included both negative and positive values and as such was normalized with a YeoJohnson-transformation (Yeo and Johnson, 2000) for CH₄ and with arcsinh-transformation for N₂O. Optimal data transformation options were investigated with the bestNormalize R-package (Peterson, 2021). Statistical differences between months and between land cover classes were tested with Aligned Rank Transformed Anova, (ART-Anova, Wobbrock et al., 2011) using the ARTool R-package (Kay et al., 2025). This allows testing for differences between groups (months or land cover class) while accounting for repeated measurements and uneven group sizes. Post-Hoc tests were done following the ART-C contrast testing procedure (Elkin et al., 2021) with Tukey's HSD using the art-contrasts R-package (Elkin et al., 2021). Differences between groups were considered statistically significant at p<0.05.

3 Results

3.1 CH₄ between land cover classes

All dry land cover classes (shrub and bare tundra, deciduous and evergreen forest) were small net sinks (mean −0.12 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) for CH₄ and conversely all wetland classes were net sources (mean 17.1 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) of CH₄ (Table A1), fluxes across all measurements varied from −0.55 to 1008 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (Table A1). Highest positive mean CH₄ fluxes were measured from permafrost bogs (mean 27.8 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; range −0.24 to 602), followed by fens (23.4 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; 0.13 to 1008), bogs (10.1 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −0.11 to 668) and tundra wetlands (0.93 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −0.04 to 7.59). Largest negative mean CH₄ fluxes (sinks) were measured from shrub tundra (−0.16 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −0.55 to 0.001) followed by bare tundra (−0.13 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −0.52 to 0.01), deciduous forests (−0.12 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −0.50 to 3.02) and lastly evergreen forests (−0.045 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −0.40 to 0.57). Although both deciduous and evergreen forest were net sinks, they also included some small positive fluxes which all came from locations in ecotones adjacent to wetlands or streams.Pairwise comparisons showed that 21 out of 28 land cover classes differed significantly in CH₄ fluxes (Table A7).

3.2 CH₄ temporal trends

Across all land cover types, both negative and positive fluxes were small in early June, increased to their highest values during July–August, and declined again in September. However, the magnitude of monthly differences was dependent on land cover type. Monthly mean fluxes can be found in Table A2 and statistical differences between months in Table A5. Wetland sites remained net sources throughout the measurement period, with positive fluxes persisting until the end of the growing season. (Fig. 3a, Table A2). No significant effect of measurement year on CH₄ fluxes was detected. Within-plot variability across months was often high, but extreme values typically originated from the same few measurement locations. For example, one permafrost bog plot, across multiple months, produced repeated CH₄ outliers of 50, 116, 366, and 602 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, resulting in a within-plot range of over 500 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$.

Figure 3 Growing season variability of (a) CH₄ and (b) N₂O (both log₁₀-transformed) with significant monthly differences in the eight land cover classes. Positive values are emissions and negative values are sinks. Boxplots have 25th and 75th quantiles, median (horizontal line), 1.5 IQR from the lower and higher quantiles (whiskers) and individual data points (black dots). In addition, significant differences between land cover classes for (c) CH₄, and (b) N₂O (${}^{\ast \ast \ast }$ p<0.001, ${}^{\ast \ast }$ p<0.01, ^∗ p<0.05).

3.3 N₂O between land cover classes

N₂O showed high variance within land cover classes, with fluxes ranging from −141 to 59.5 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ across all sites. Notably, all wetland types were on average, net sinks of N₂O, although all land covers exhibited both positive and negative fluxes (Fig. 3b). Positive mean N₂O fluxes were observed in deciduous forests (mean 1.70 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; range −15.6 to 39.2) and shrub tundra (0.44 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −36.9 to 23.0). The strongest negative mean fluxes (net sinks) occurred in fens (−6.53 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −142 to 5.13) and permafrost bogs (−5.40 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; −92.8 to 16.6). N₂O sinks were also measured in bogs (−3.51 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), tundra wetlands (−0.89 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), bare tundra (−0.69 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), and evergreen forests (−0.52 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). Pairwise comparisons showed that 15 out of 28 land cover classes differed significantly in N₂O fluxes (Table A8).

Mean CH₄ and N₂O fluxes were significantly inversely correlated (Spearman $R=-\mathrm{0.70}$, p<0.001), with positive CH₄ fluxes generally associated with strong negative N₂O fluxes (sinks) (Fig. A2).

3.4 N₂O Temporal trends

Unlike CH₄, most land cover classes showed no significant trends between months for N₂O (Table A6). However, there was a significant difference in N₂O fluxes between years, where year 2025 had smaller positive N₂O fluxes and larger negative fluxes (Fig. 7c for model results and Fig. A1 for observed yearly fluxes). Furthermore, within plot variability was often very high and some plots would measure both negative and positive fluxes in monthly repeated measurements.

3.5 Environmental variables

Soil temperature during the conducted chamber measurements varied between −0.02 and 22.18 °C (mean 10.23 °C). September was the coldest month (soil temperature) with the exceptions of evergreen forests, deciduous forests and permafrost bogs, where the coldest month was June. Mean soil temperatures peaked in July or August depending on the land cover (Fig. 4b), with the highest single measurements generally occurring during August. Soil moisture varied between 4.7 % and 60.7 % across all plots. Mean soil moisture was highest in fens (56.5 %), followed by tundra wetlands (55.4 %), bogs (40.8 %), permafrost bogs (34.2 %), evergreen forests (32.5 %), deciduous forests (31.9 %), shrub tundra (26.9 %) and bare tundra (25.2 %). Generally, soil moisture was highest in June due to snow melt waters and then steadily dropped towards August. Depending on the land cover, August or September was generally the driest month (Fig. 4a). Green vegetation cover ranged from 0.01 % to 85.6 %, with mean values highest in deciduous forest (37.5 %) and lowest in bare tundra (3.2 %); other land covers were between 12 %–25 %. Green vegetation cover peaked in July or August depending on the land cover while the lowest values were recorded in June and September (Fig. 4c). Soil pH was acidic across all plots except for a few locations with dolomite runoff. Mean pH was highest in tundra wetlands (5.53), likely reflecting the influence of dolomite runoff at these sites, followed by bare tundra (4.99), deciduous forests (5.07), fens (4.78), shrub tundra (4.49), bogs (4.22), evergreen forests (4.34) and permafrost bogs (4.12) (Fig. 4d).

Figure 4 Observed variation in soil and vegetation characteristics. Boxplots of (a) soil moisture, (b) soil temperature and (c) green vegetation cover between different land covers and months and (d) pH between land covers. Boxplots have 25th and 75th quantiles, median (horizontal line), mean (white dot), 1.5 IQR from the lower and higher quantiles (whiskers).

3.6 Multivariate models

3.6.1  CH₄

The mixed-effects model for CH₄ showed high explanatory power (R2 cond.=0.96, R2 marg. 0.81). Moisture, temperature and pH all influenced CH₄ fluxes, however the direction of the response was dependent on the land cover (Figs. 5 and 6). Highest effect variables based on partial eta-squared (η²p) were soil temperature (interaction) and pH (interaction) (Fig. 5, final model). For land cover main effect, the CH₄ fluxes were highest in fens, followed by bogs, permafrost bogs, tundra wetlands, evergreen forests, bare tundra, shrub tundra and deciduous forests (Table A3). All dry land cover classes modelled as net sinks (negative fluxes).

Figure 5 Model performance indicators and effect sizes for the CH₄ linear mixed-effects models. The baseline model, which included only main effects, showed lower explanatory power compared to the interaction model that incorporated interactions with land cover class. The final model was obtained after removing the least influential predictors based on AIC to achieve a more parsimonious model structure. Partial eta-squared (Eta2 part.) describes the effect size of each model parameter.

Figure 6 Land cover specific modelled (linear mixed-effect model) responses of CH₄ to (a) soil temperature, (b) soil moisture, (c) pH. All CH₄ model variables included interaction terms, so no simple main effects are shown. Dashed lines represent wetlands. The figure includes only predicted (normalized) values for GHG fluxes.

3.6.2  N₂O

The mixed-effects model for N₂O showed modest explanatory power (R2 cond.=0.38, R2 marg.=0.30, Fig. 7). Modelled N₂O fluxes were highest in deciduous forests, followed by shrub tundra, bare tundra, evergreen forests, tundra wetlands, bogs, fens and permafrost bogs (Fig. 8a). Based on η²p, land cover and its interaction with pH were the strongest predictors while both soil moisture and temperature had noticeably less effect (Fig. 7). Unlike CH₄, years 2024 and 2025 also showed a noticeable difference for N₂O with 2025 generally having smaller positive fluxes and larger sinks (Fig. 8, Table A7). In addition, vegetation cover had a small but significant effect on N₂O fluxes depending on the land cover type (Figs. 7 and 8). The point identification random effect also had a weaker effect compared to the CH₄ model.

Figure 7 Model performance indicators and effect sizes for the N₂O linear mixed-effects models. The baseline model, which included only main effects, showed lower explanatory power compared to the interaction model that incorporated interactions with land cover class. The final model was obtained after removing the least influential predictors based on AIC to achieve a more parsimonious model structure. Partial eta-squared (Eta2 part.) describes the effect size of each model parameter.

Figure 8 Modelled (linear mixed-effects model) responses of N₂O to (a) land cover, (b) soil temperature, (c) year, (d) pH-land cover interaction and (e) green vegetation cover-land cover interaction. Note differences in the scale of y axis between subfigures (use η²p for comparison of effect size between parameters). Dashed lines represent wetlands. The figure includes only predicted (normalized) values for GHG fluxes.

4 Discussion

Our findings on the seasonal and spatial patterns of CH₄ and N₂O fluxes provide several insights into how environmental drivers shape gas exchange across high-latitude subarctic and boreal landscapes. Both GHGs showed strong variation among land-cover classes (research question i), with 21 of 28 pairwise comparisons being significant for CH₄ (Table A7) and 15 of 28 for N₂O fluxes (Table A8). The seasonal development of CH₄ fluxes (research question ii) was strongly dependent on land-cover class, with significant differences between months observed mainly in fen, deciduous forest, and shrub tundra sites (Table A5). In contrast, N₂O showed no clear seasonal trends (Table A6). Regarding environmental controls (research question iii), the strongest predictors of CH₄ fluxes were soil temperature and pH, with a weaker effect of soil moisture. For N₂O, pH was the strongest predictor, with weaker effects of year and green vegetation cover. For both GHGs, the direction and strength of these relationships often varied among land-cover classes. This suggests that the influence of environmental drivers is strongly context dependent and that the same change in environmental conditions can lead to contrasting flux responses across ecosystems. Consequently, broad-scale models based on uniform driver-response relationships may fail to capture important spatial variation in GHG exchange across subarctic and boreal landscapes.

4.1 GHG flux magnitudes

Dry land cover classes acted as net CH₄ sinks (−0.16 to −0.05 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, depending on land cover class) and as both N₂O sources and sinks (−0.69 to 1.69 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). Wetlands were net CH₄ sources (0.93 to 27.81 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) and net N₂O sinks (−0.89 to −6.53 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), with fens and permafrost bogs having the highest CH₄ fluxes and largest N₂O sinks (see full descriptions of GHG magnitudes in Tables A1 and A2).

Overall, the wetland CH₄ emissions were quite high but within maximum ranges reported in published literature (Turetsky et al., 2014; Tupek et al., 2015; Bao et al., 2021; Kuhn et al., 2021; Virkkala et al., 2026). Similar to our measurements, the ABCFlux CH₄ database (Virkkala et al., 2026) also had peak fluxes in August. Our forest and tundra CH₄ sinks were larger (more negative) than many previously reported values (Emmerton et al., 2014; Voigt et al., 2023; Virkkala et al., 2026), where some studies have even reported positive fluxes from upland tundra soils (Kuhn et al., 2021). Also, our mean net CH₄ sinks peaked in August, whereas in the ABCFlux database CH₄ sinks peaked in July (Virkkala et al., 2026). Our N₂O sources were generally at the lower end of the reported ranges, while the N₂O sinks were at the higher end (more negative) (Repo et al., 2009; Lohila et al., 2010; Marushchak et al., 2011; Tupek et al., 2015; Voigt et al., 2020; Virkkala et al., 2024; Triches et al., 2025).

4.2 Environmental drivers of CH₄

We observed that CH₄ emissions increased with higher soil moisture across most of the study sites. This pattern is consistent with the underlying microbial processes: methanogenesis occurs primarily in anaerobic soils (Le Mer and Roger, 2001; Lai, 2009), whereas methanotrophy dominates under aerobic conditions (Whalen, 2005; Lai, 2009). Soil oxygen content follows soil water content (Mainiero and Kazda, 2005) and as such soil moisture should have a direct causal effect on CH₄ emissions. Our results show that in dry land cover classes, soil moisture was a good indicator for CH₄. The few positive CH₄ measurements in dry land cover classes (Fig. 3a) were always from ecotones close to wetlands or streams or because of spring snow melt waters saturating the soils for extended periods of time. Conversely, CH₄ sinks in dry soils also grew stronger (more negative) with decreasing soil moisture. Notably, all mean and median fluxes for CH₄ in dry land cover classes were negative. Christiansen et al. (2015b), Voigt et al. (2023) and Virkkala et al. (2024) reported similar findings, with soil moisture being one of the major drivers for CH₄ sinks at high-latitude dry upland soils and Lohila et al. (2016) found that a prolonged wet autumn period changed a net CH₄ sink forest into a net source for several months. Although the CH₄ sinks are relatively small when compared to wetland emissions, dry tundra and boreal forest soils cover extensive areas of the northern hemisphere which leads to a globally significant CH₄ sink (Jiang et al., 2025), and as previously mentioned (see Sect. 4.1), the sink magnitudes reported here are larger than in many previous studies. The soil CH₄ sink is still poorly constrained in global CH₄ budgets (Saunois et al., 2025) and should be considered in further studies.

In wetland environments, soil moisture has a more complex relationship to CH₄. Often the wetland types with highest soil moisture (i.e. fens and marshes) are also the highest CH₄ emitters (Kuhn et al., 2021), however within those wetland types, increased soil moisture does not necessarily increase CH₄ flux (Knox et al., 2019). Our results show that, within the different fen sites, the wettest soils were not the largest emitters (Fig. 5b). Soil moisture also did not increase CH₄ emissions in bogs or permafrost bogs. In permafrost bogs, intact permafrost bogs exhibited higher soil moisture but acted as substantially smaller CH₄ sources compared to fully degraded ones. This suggests that additional factors, such as soil pH, vegetation, or active layer thickness, may play a more important role in regulating wetland CH₄ fluxes. Similar results were reported by Turetsky et al. (2014) and Rinne et al. (2018), who found that CH₄ fluxes were smaller in fully flooded wetland soils and that higher water-table levels alone did not predict higher fluxes in boreal fens. Our findings suggest that there is an upper threshold of soil moisture for methanogenesis in wetlands, beyond which further increases in soil moisture no longer increase CH₄ fluxes. Reduced CH₄ emissions at sites with extreme soil moisture are likely the result of differences in substrate availability (i.e. less labile organic material) and temperature regime (Strack et al., 2004).

Soil temperature and land cover interaction was the strongest predictor for CH₄ fluxes (Fig. 5) and soil temperature has been identified as a key driver of CH₄ fluxes in many previous studies (Turetsky et al., 2014; Knox et al., 2019; Kuhn et al., 2021; Kübert et al., 2026). Soil temperature directly affects both methanogenesis and methanotrophy as they are microbial processes and therefore causally linked to soil temperature through the temperature-dependence of microbial activity (Moore and Dalva, 1993; Biasi et al., 2005; Deslippe et al., 2012). This also indicates that rising temperatures can intensify CH₄ dynamics in opposite directions where the net effect on emissions is dependent on the ecosystem's temperature sensitivity (Zhang et al., 2021). This effect can also be seen in our modelled CH₄ flux responses (Fig. 6a), where wetland CH₄ fluxes have a positive response to higher temperatures and dry land cover classes have a negative response, with higher temperatures leading to larger sinks (i.e. net-sink environments become stronger sinks and net-source environments emit more CH₄). Additionally, warming can affect the amount of organic material available for decomposition through increased GPP as higher mean temperatures generally increase biomass which could further enhance CH₄ fluxes (Knox et al., 2021). Thus, rising temperatures under climate change can strengthen both CH₄ sinks and CH₄ sources depending on land cover type, meaning that long-term changes in net CH₄ emissions will depend largely on shifts in land cover and ecosystem distributions and their effects on soil conditions.

Soil pH had a strong but mixed influence on CH₄ fluxes depending on land cover type (Figs. 5 and 6c). This response appeared to be related to the overall pH of these environments. Low-pH systems, such as permafrost bogs, bogs, and evergreen forests (Fig. 4d), showed a positive response to increasing pH, meaning higher CH₄ fluxes. In contrast, higher-pH environments, such as tundra wetlands and deciduous forests, showed the opposite pattern, meaning smaller fluxes with increased pH. It is worth noting the overall pH range was between 3.6–6.3 and as such, most vegetation and soil biota, including microbes, are adapted to acidic conditions (Kip et al., 2011; Yao et al., 2023). At the lower end of that pH range, many plant species and microbes can exist at the far end of their acidity tolerance (Yao et al., 2023) and even small increases in pH could increase their overall activity. This may explain the strong response seen in bogs and permafrost bogs as pH increases. At the opposite end, some of the tundra wetland and deciduous forest sites are affected by dolomite runoff from the surrounding watershed which significantly increases the overall soil pH. At these sites a higher pH decreased the CH₄ flux (Fig. 6c), which could indicate that methanogenesis is inhibited by higher pH while methanotrophs benefit from it. Conversely, it could also be a result of methanogenetic Archaea being more tolerant of acidic conditions than methanotrophs (Weslien et al., 2009). Other factors could also be of importance, specifically nutrient availability and vegetation are partly controlled by soil pH and directly affect CH₄ dynamics. Nonetheless, pH can exert a strong effect on CH₄ fluxes and should be considered in further studies.

4.3 Environmental drivers of N₂O

Similarly to CH₄ production and consumption, N₂O dynamics are also governed by microbial processes and are therefore influenced by factors such as soil temperature, moisture, pH, and vegetation. N₂O is mostly produced in the denitrification process where nitrates are reduced back into N₂ (Butterbach-Bahl et al., 2013). Other possible, but less important pathways include nitrification and nitrate/nitrite reduction to ammonium (Butterbach-Bahl et al., 2013). Soil oxygen content has a causal link to N₂O flux as nitrification generally only occurs in aerobic conditions and denitrification and reduction to ammonium mostly occur in anaerobic conditions (Shan et al., 2021). However, N₂O showed a much weaker response to the environmental variables compared to CH₄, and the explanatory power of the models was correspondingly smaller (Fig. 7). In addition, N₂O showed a significant difference between measurement years which remained unexplained by the other factors, suggesting that some key explanatory variables were missing from the model. Year 2024, which had the higher mean N₂O fluxes, was considerably warmer during the measurement period than 2025. Difference in growing season mean temperature for Pallas was 0.90 and 1.35 °C for Kilpisjärvi (Table 1). It is possible that these temperature differences contributed to the observed patterns. However, because soil temperature showed only a weak relationship with N₂O fluxes, any such effect was likely indirect, for example through the influence of air temperature on vegetation. Previous studies have identified sources of N₂O from many high-latitude environments such as unvegetated permafrost peat plateaus and palsas (Marushchak et al., 2011) and both vegetated and unvegetated permafrost tundra soils (Voigt et al., 2020) but also sinks from nutrient-poor permafrost peatlands (Triches et al., 2025) and grasslands and forests (Chapuis-Lardy et al., 2007). In our measurements, we found strong sinks from wetlands, while dry land cover classes had medians close to zero fluxes, however there was a noticeable amount of both negative and positive extreme values (Fig. 3b) and many land cover classes and single measurement plots exhibited both negative and positive fluxes. This type of flux behaviour is not uncommon (Lohila et al., 2010; Tupek et al., 2015) and highlights the spatiotemporally sporadic behaviour of N₂O, where many plots have zero or small positive N₂O emissions, yet a small number of plots can act as significant hotspots, both for sinks and sources (Marushchak et al., 2011; Triches et al., 2025).

Table 1 Long-term (1991–2020) and study-period (2024–2025) mean temperatures and precipitation during the growing season (June–September). Data were obtained from the nearest weather stations (Kilpisjärvi Kyläkeskus and Muonio Kirkonkylä).

Soil temperature and moisture generally had modest effects on N₂O fluxes and much of the variability was explained by pH and land cover (particularly wetlands). Our highest measured mean N₂O fluxes came from deciduous (mountain birch) forests, which also had comparatively higher soil pH (Fig. 4d). Cryoturbation and the absence of vegetation have been occasionally reported as having an increasing effect on N₂O fluxes, especially in permafrost soils (Repo et al., 2009; Marushchak et al., 2011; Voigt et al., 2020). However, we found no strong evidence for such an effect in our measurements from peat soils or tundra cryoturbation plots. Both bare tundra and permafrost bogs acted as small net sinks for N₂O (Fig. 3b), and vegetation cover only had a small effect in the model results.

It is possible that low soil nitrogen availability is the primary reason for the small N₂O fluxes, while also being a key missing variable in the model. Nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) are the preferred electron acceptors in denitrification, thus in the presence of excess soil nitrogen content, N₂O would not be reduced and would likely result in larger positive fluxes. Furthermore, there is evidence of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reliant pathways of methanotrophy that can produce N₂O in low oxygen conditions such as wetlands (Kits et al., 2015).

4.4 N₂O sinks

The largest N₂O sinks were found in wetlands, specifically in fen plots (Fig. 3b). This may appear counterintuitive as denitrification, which is the main pathway for N₂O production (Butterbach-Bahl et al., 2013), requires anoxic conditions which wetlands typically have. The reason is likely in the last step of the denitrification chain where N₂O is reduced. In anoxic conditions with high $\mathrm{C}/\mathrm{N}$ ratios or low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations, the denitrification process is fully completed, N₂O is produced but then reduced into N₂ gas (Kolb and Horn, 2012; Brummell et al., 2012; Pessi et al., 2022). In these low-nitrate conditions, any excess N₂O would be rapidly reduced, including atmospheric N₂O, resulting in a net sink.

Other less studied pathways for N₂O sinks and the inverse N₂O and CH₄ correlation involve anaerobic methane oxidation driven by denitrification (Valenzuela et al., 2020; Tan et al., 2025). It follows that in soils rich in CH₄, such as wetlands, substantial methanotrophic activity would occur throughout the soil profile, provided that environmental conditions permit it (Whalen et al., 1990). However, chamber measurements only capture the net CH₄ flux, so if net oxidation is smaller than net CH₄ production, the oxidation would not be detected. Anerobic methane oxidation (AOM) is a well-known process in marine environments (Hinrichs and Boetius, 2003; Reeburgh, 2007), but less so in wetlands, yet there is growing evidence for its importance (Smemo and Yavitt, 2011; Cai et al., 2021; Cheng et al., 2021). In AOM, CH₄ can be used as electron donor for reducing nitrite and N₂O, which would act as an additional sink of N₂O (Valenzuela et al., 2020; Cheng et al., 2021; Tan et al., 2025). However, the magnitude of these pathways in natural systems is unknown and somewhat hypothetical at this point and needs more research to be confirmed. Furthermore, the negative relationship between CH₄ and N₂O fluxes was restricted to wetland ecosystems (Fig. A2b). This suggests that the coupling between CH₄ and N₂O dynamics depends on environmental conditions and may be strongest under the wet and anaerobic conditions characteristic of wetlands. Regardless of the underlying mechanisms, based on our measurements, wetlands can act as notable N₂O sinks in the right conditions.

5 Uncertainties

The monthly repeated measurements and two-year design resulted in a robust overall dataset, even though some land cover classes, such as permafrost bogs, were represented by fewer plots. Generic environmental variables that can be measured from anywhere make the models easily applicable everywhere, but it also limits the number of usable explanatory variables. Modelling across diverse environments (i.e. not solely focusing on specific ecosystems such as wetlands) provides a valuable opportunity to capture a wide range of ecosystem behaviors, though it naturally results in more complex data structures. In particular, wetland GHG fluxes often show substantial variability, including occasional high fluxes far above median values (Kuhn et al., 2021; Saunois et al., 2025). These contrasts become pronounced when wetlands and non-wetland systems are analyzed together due to their inherently different flux magnitudes. Furthermore, capturing seasonal changes requires repeated measurements which limit the available statistical methods due to pseudoreplication. These complexities were effectively addressed in our statistical analyses through the use of robust and flexible approaches, including mixed multivariate models, appropriate data transformations, and non-parametric tests.

The multivariate and correlative nature of our models makes it challenging to directly disentangle causal pathways. Also, some variables, such as land cover, function as a composite descriptor of multiple co-varying environmental conditions rather than a direct mechanistic driver. Because land cover classes ultimately reflect underlying temperature, moisture, and other soil properties, these environmental factors are the more direct drivers of GHG fluxes than land cover per se.

Opaque (“dark”) chambers are used to reduce uncertainties caused by changes in chamber temperature or possible interactions between photosynthesis and GHG fluxes and are typically the default method for measuring CH₄ and N₂O fluxes (Livingston and Hutchinson, 2009; Pavelka et al., 2018; Clough et al., 2020). However, there is recent evidence of CH₄ and N₂O fluxes being affected by varying light conditions (Voigt et al., 2023; Triches et al., 2025). Thus, it is possible that GHG fluxes measured with transparent and opaque chambers are not directly comparable and the additional use of transparent chambers should be considered in future research.

6 Conclusions

Across all ecosystems, we found consistent contrasts in greenhouse gas dynamics, with tundra and forest sites acting as modest CH₄ sinks and wetlands functioning as strong CH₄ sources, whereas N₂O flux patterns were more subtle. Most notably wetlands acted as consistent N₂O sinks.

Land cover, pH, soil temperature and soil moisture were strong predictors of CH₄ fluxes and explained most of the observed variance, whereas green vegetation cover had only a marginal effect. CH₄ fluxes exhibited clear seasonal patterns, with the smallest sources and sinks occurring at the beginning and end of the growing season while peak fluxes were observed in July–August. These findings highlight the importance of longitudinal studies for accurately estimating annual CH₄ budgets rather than relying solely on peak growing-season measurements but also demonstrates that CH₄ fluxes can be modelled with relatively few key environmental variables.

In comparison to CH₄, N₂O exhibited relatively moderate responses to the measured environmental variables and considerably weaker seasonal patterns, emphasizing the inherent challenges in explaining N₂O dynamics across heterogeneous ecosystems. In a multivariate context, land cover and soil pH emerged as the most influential factors shaping N₂O flux variability, while the effects of soil temperature and moisture were negligible. The substantial proportion of unexplained variability in N₂O fluxes indicates that key controlling factors were not fully captured, pointing to the importance of additional drivers, particularly those related to soil biogeochemistry and microbial nitrogen transformation pathways. Furthermore, the high variability of N₂O emissions suggests a notable influence of spatiotemporal dynamics, underscoring the need for high resolution measurements to adequately characterize N₂O behavior. These findings highlight the complexity of N₂O dynamics and identify important avenues for future research to improve the understanding and modeling of N₂O fluxes.

All in all, our results demonstrate that high-latitude GHG fluxes arise from interactions between broad land cover patterns and fine-scale environmental variability. From an applied perspective, our findings suggest that accurate land cover mapping and in situ measurements of soil conditions are essential for robust modelling and upscaling of GHG fluxes. Because land cover classes exhibit such different GHG dynamics due to the underlying soil and vegetation characteristics, improved classification greatly enhances predictive performance across heterogenous landscapes. This is especially important in high-latitude regions, where rapid warming is expected to alter both ecosystem distributions and soil conditions with significant implications for current and future GHG budgets.

Appendix A

Table A1 CH₄ and N₂O fluxes in different land cover classes.

Table A2 Monthly CH₄ and N₂O fluxes in different land cover classes.

Table A3 Linear mixed-effects model summary for CH₄ flux. Model estimates are standardized.

Table A4 Linear mixed-effects model summary for N₂O flux. Model estimates are standardized.

Table A5 ART-ANOVA summary for monthly differences of CH₄ flux in different land cover types. Significance intervals P<0.001 ${}^{\ast \ast \ast }$, P<0.01 ${}^{\ast \ast }$, P<0.05 ^∗.

Table A6 ART-ANOVA summary table for monthly differences of N₂O flux in different land cover types. Significance intervals P<0.001 ${}^{\ast \ast \ast }$, P<0.01 ${}^{\ast \ast }$, P<0.05 ^∗.

Table A7 ART-ANOVA summary table for differences of CH₄ flux between land cover classes. Significance intervals P<0.001 ${}^{\ast \ast \ast }$, P<0.01 ${}^{\ast \ast }$, P<0.05 ^∗.

Table A8 ART-ANOVA summary table for differences of N₂O flux between land cover classes. Significance intervals P<0.001 ${}^{\ast \ast \ast }$, P<0.01 ${}^{\ast \ast }$, P<0.05 ^∗.

Figure A1 Yearly differences in N₂O fluxes between land cover classes. Colors indicate years, positive values are emissions and negative values are sinks. Note that permafrost bog category is not directly comparable between years due to additional plots being added in 2025. Boxplots have 25th and 75th quantiles, median (horizontal line), 1.5 IQR from the lower and higher quantiles (whiskers) and outliers (black dots). Y axis is log10 transformed.

Figure A2 (a) Scatterplot of mean CH₄ flux (log10) and mean N₂O flux (log10) with spearman correlation (b) Scatterplot of mean CH₄ flux and mean N₂O flux between non-wetland soils (shrub tundra, bare tundra, deciduous forest, evergreen forest) and wetland soils (bog, fen, tundra wetland, permafrost bog) with spearman correlation.