Water chemistry and greenhouse gas concentrations in waterbodies of a thawing permafrost peatland complex in northern Norway

Thermokarst ponds in thawing permafrost landscapes play a considerable role in greenhouse gas (GHG) emissions despite their small size, yet they remain underrepresented in Earth system models. Transitions from hydrologically isolated thermokarst ponds in peat plateaus to connected wetlands can substantially alter GHG dynamics. However, the processes and GHG impacts of these shifts are not well understood, particularly in the sporadic permafrost zones of Fennoscandia. To address this, we investigated water chemistry, dissolved organic matter (DOM) processing, and GHG fluxes over two years at the Iškoras site in northern Norway, where a degrading peat plateau includes both thermokarst ponds and an adjacent wetland stream. Thermokarst ponds exhibited low pH, high organic acidity, high oversaturation of dissolved carbon dioxide (CO₂), and, especially, high concentrations of dissolved methane (CH₄). The adjacent wetland stream, however, with near-neutral pH, showed lower CH₄ concentrations and organic acidity but significantly higher CO₂ emissions, despite moderate saturations driven by turbulence and bicarbonate replenishment. By contrast, CO₂ emissions in ponds were primarily linked to DOM mineralization. Despite differences in chemistry, DOM mineralization rates were similar between ponds and the stream sites, suggesting that environmental factors like pH and microbial community differences counteract DOM lability variations. As permafrost decays and transitions from peat plateaus to wetlands, ponds as hotspots of CH₄ emissions will disappear. However, total GHG fluxes across the peatland–wetland continuum will depend on wetland emissions, where CH₄ emissions are usually considerable, and the fate of organic matter within the plateau. Lateral DOM fluxes may represent a significant loss of soil organic carbon (OC), highlighting the importance of hydrological connectivity in linking terrestrial and aquatic systems. This study emphasizes the need to account for the relationship between hydrological and chemical processes when assessing C and GHG fluxes in permafrost-impacted regions.

1 Introduction

Northern-latitude permafrost regions hold one of the largest terrestrial carbon reservoirs on the planet (Schuur et al., 2008; Schuur et al., 2015; Walter et al., 2006). Although covering only about 15 % of global soils, these regions store an estimated 1400–1600 Pg of organic carbon (OC) (Hugelius et al., 2014; Schuur et al., 2022; Strauss et al., 2025), making them a critical component of the global carbon (C) cycle. Sequestered under cold and oxygen-limited conditions, this C is increasingly vulnerable to release as permafrost thaws due to climate warming, generating significant feedbacks that complicate predictions of future climate trajectories (Hugelius et al., 2020; Schuur et al., 2008, 2015; Walter et al., 2006). As permafrost degrades, the release of greenhouse gases, particularly methane (CH₄) and carbon dioxide (CO₂), through the microbial decomposition of previously frozen organic matter (OM), can rapidly escalate the impact of this feedback (Schuur et al., 2008; Walter et al., 2008; Wik et al., 2016; Zimov et al., 2006). While the large-scale thaw of permafrost is widely recognized (Leppiniemi et al., 2023) and permafrost regions warm three to four times faster than the global average (Meredith et al., 2019), the timing, magnitude, and pathways of carbon release remain uncertain, influenced by such processes as burial, mobilization, lateral export, and mineralization (Verdonen et al., 2023; Vonk et al., 2015).

Permafrost thaw leads to irreversible landscape transformations. Peatlands in northern Norway are predominantly located in the sporadic permafrost zone, where they form distinctive landscape features, such as peat plateaus and palsas. These are peat uplands and mounds with a frozen core, elevated above the water table by the formation of segregation ice (Alewell et al., 2011; Krüger et al., 2017). As these features degrade, permafrost thaw is often abrupt and subsidence and collapse is to be expected, leading to the formation of thermokarst ponds, as excess ground ice is lost (Martin et al., 2021). More than half of the permafrost areas in the Scandinavian Peninsula are at risk of disappearing under current and projected climate conditions (Gisnås et al., 2017; Schuur et al., 2008). The areal extent of peat plateaus in this region decreased by 33 %–71 % between the 1950s and the 2010s, with rapid degradation observed during the last decade (Borge et al., 2017). This regional degradation mirrors processes observed across the northern hemisphere, including in the Canadian Arctic, European Russia, and the Kola Peninsula, highlighting the vulnerability of sporadic permafrost regions to warming climates (Krutskikh et al., 2023; Payette et al., 2004; Sannel and Kuhry, 2011). While the processes driving permafrost thaw and landscape transformations, such as thermal disturbances, vegetation shifts, and subsidence, are relatively well-studied, their consequences for greenhouse gas (GHG) fluxes and C cycling remain uncertain, limiting our ability to project future climate feedbacks (Holmes et al., 2022; Olefeldt et al., 2021; Turetsky et al., 2020).

Among the new landscape forms that emerge from degrading peat plateaus, thermokarst ponds and wetlands play a critical role in greenhouse gas dynamics. These small aquatic systems, formed by the thaw and collapse of permafrost, are characterized by high concentrations of dissolved organic carbon (DOC) and dissolved inorganic carbon (DIC) (Abnizova et al., 2012; Martin et al., 2021; Matveev et al., 2018). Thermokarst ponds, in particular, act as hotspots for CH₄ and CO₂ emissions due to unique biogeochemical conditions, including hydrological isolation, anoxic sediments, and high organic matter availability (in 't Zandt et al., 2020; Polishchuk et al., 2018; Vonk et al., 2015; Ward and Cory, 2015). Despite the small size of thermokarst ponds, these waterbodies can contribute significantly to regional C fluxes, with CH₄ and CO₂ supersaturation levels often surpassing those of larger lakes – whether thermokarst or not – or surrounding tundra ecosystems (Abnizova et al., 2012; Kuhn et al., 2018; Shirokova et al., 2012). However, the contributions of thermokarst ponds are often overlooked in large-scale C assessments, as they remain difficult to detect using satellite-based methods because of their small size (Holgerson and Raymond, 2016; Muster et al., 2017).

As permafrost thaw progresses, the transition of isolated thermokarst ponds to interconnected wetland systems further alters GHG dynamics. While northern permafrost wetlands currently act as a C sink, the inclusion of thaw pond emissions into broader wetland carbon budgets reveals their potential to offset the sink capacity by 39 % (Kuhn et al., 2018). Compared with thermokarst ponds, wetlands have sustained CH₄ fluxes over larger areas, due to persistent waterlogging and OM decomposition (Pirk et al., 2024; Swindles et al., 2015; Turetsky et al., 2020), thus constituting important long-term CH₄ sources (Bansal et al., 2023). The transformation from stable permafrost to thermokarst landscapes is accompanied by shifts in hydrology, OM lability, and microbial activity, which collectively shape CO₂ and CH₄ production pathways (Holmes et al., 2022; Laurion et al., 2020). Understanding the dynamics of these evolving permafrost and wetland systems is critical for assessing the broader impacts of permafrost thaw on regional C uptake and emissions, as well as global C cycles.

Northern Norway's sporadic permafrost zone, with its abundant small thermokarst ponds and emerging wetlands, provides a valuable opportunity to investigate carbon cycling in rapidly evolving subarctic landscapes. The region's degrading peat plateaus host significant C stocks, yet small aquatic systems, especially those in Fennoscandia, remain underrepresented in Earth system models (Abnizova et al., 2012; Muster et al., 2019; Muster et al., 2017). While existing studies emphasize the importance of quantifying CH₄ and CO₂ fluxes in these environments and their implications for C budgets (Abnizova et al., 2012; Matveev et al., 2018), the interactions between hydrology, vegetation, and carbon processing are not well understood. Yet such processes are central to key questions regarding how transitions between permafrost, thermokarst, and wetland systems influence C dynamics and whether these landscapes function as net C sources or sinks under changing climatic conditions (Sim et al., 2021). In particular, peatland ponds and thermokarst waterbodies exhibit unique biogeochemical cycling from lakes, driven more by internal dynamics than external watershed inputs (Arsenault et al., 2022). These differences remain poorly represented in both observational datasets and Earth system models.

This study aims to address these gaps by examining the GHG dynamics and C biogeochemistry of thermokarst ponds and wetland streams in the sporadic permafrost zone of northern Norway. Over two years, we collected a novel dataset combining biogeochemical and dissolved gas measurements with C flux data from thermokarst ponds and a wetland stream within a small permafrost peatland plateau undergoing rapid permafrost degradation. This setting captures a landscape in active transition from isolated thermokarst ponds to interconnected wetlands. We hypothesize that (1) thermokarst ponds serve as hotspots of CH₄ and CO₂ production, relative to the wetland stream; (2) the transition from isolated ponds to wetlands significantly alters GHG emission pathways, driven by shifts in hydrology and OC availability; and (3) recently mobilized OM from thawing permafrost presents a labile source of C, promoting CO₂ production in thermokarst water bodies compared with the wetland stream. By exploring these dynamics, this study provides insights into the role of small water bodies in permafrost-C feedbacks, advancing our understanding of subarctic and boreal C cycling.

2 Methods

2.1 Study area

The Iškoras field site (69.34° N, 25.29° E; 381 m a.s.l.) is a permafrost peatland plateau located in the interior of Finnmark province, northern Norway, on the Finnmarksvidda (Fig. 1). The region of the Finnmarksvidda lies between 300 and 500 m a.s.l. and is characterized by a subarctic continental climate. The topography was shaped by Pleistocene glaciations, which deposited ground moraines and glaciofluvial and glaciolacustrine sediments (Sollid et al., 1973). The depressions in the landscape are commonly filled with peatlands (Borge et al., 2017); peat plateaus underlain by permafrost are common.

Figure 1 Location map of the study area in (a) Europe and (b) the Iškoras catchment (aerial image by the Norwegian Mapping Authority, Norgeibilder, 2024), determined from the outlet station (Outlet). (c) Close-up of the wetland (Norgeibilder, 2024) with regular sampling sites and main water flow direction. (d–f) Pictures of the three main pond sites.

The Iškoras peat plateau covers an area of approximately 4 ha and is part of a 3.3 km² subarctic headwater catchment that drains into the Báhkiljohka river (91 km²). Mean annual air temperature and precipitation for the 30-year normal (1991–2020) period were −1.9 °C and 513 mm, respectively (Table 1). For our study period, 2021–2022, MAAT and MAP were −1.1 °C and 589.5 mm, respectively (SeNorge, 2023). Iškoras lies within the zone of sporadic permafrost; the peat soils extend down to about 1.5 m in the plateau areas (Kjellman et al., 2018), with active-layer depths up to 90 cm. The plateau exhibits a complex surface of intact and degrading palsas, along with thermokarst ponds, and is surrounded by wetlands and a stream to the northwest (Martin et al., 2019). Between 2019 and 2022, up to 0.8 m of subsidence of palsas was measured at localized sites (Pirk et al., 2024). The site is located about 90 km south of the nearest coastal fjord and is dominated by mountain birch (Betula pubescens) forest and tundra vegetation, including dwarf birch (B. nana). The plateau consists primarily of low heath shrubs, Ericaceae (Empetrum nigrum, Rhododendron tomentosum), lichen crusts, mosses, and cloudberry (Rubus chamaemorus) or bare ground, while the surrounding wetlands are dominated by sphagnum mosses, sedges (Carex spp.), and cotton-grass (Eriophorum spp.) (Kjellman et al., 2018; Martin et al., 2019).

Table 1 Mean and interannual standard deviation (shown as mean ± SD) of climate parameters for the Iškoras catchment for the normal period (1991–2020) and the study period (2021–2022). Summer is defined as May to September.

The study area included water bodies within a peat plateau and the adjacent wetland, selected for sampling and monitoring. Measurements and samples were taken approximately monthly in the ice-free season from May until October in 2021 and 2022 (full details in Table S1). The waterbodies consisted of three thermokarst ponds (TK-Pond 1, 2, and 3), a seasonal drainage channel (TK-Drain) connecting the peat plateau to the wetland, and the wetland's inlet and outlet streams (Inlet, Outlet), with the Outlet also marking the terminus of the Iškoras catchment. From mid-May to early November, monitoring showed that the thermokarst ponds and Inlet were ice-free for 170 ± 5 d, while the Outlet remained ice-free for 184 d (Table S2).

The thermokarst ponds varied in hydrological connectivity and permafrost influence, reflecting differences in age and physical characteristics. TK-Pond 1 (0.4 m deep) is small and located at the peat plateau–wetland transition, experiencing periodic hydrological isolation. TK-Pond 2, the largest and deepest (1.5 m deep), lies centrally on the peat plateau, surrounded by degrading palsas. TK-Pond 3 (0.6 m deep) is situated at the plateau's edge and was initially isolated by a surrounding permafrost mound.

The TK-Drain is a shallow, ephemeral drainage channel that provides the primary hydrological connection between the peat plateau and wetland. The wetland's inlet stream (0.6 m wide, 15–40 cm deep) begins approximately 200 m upstream of the plateau, flowing through birch forest and mires without permafrost before entering the wetland, where the channel becomes less defined. The outlet stream (0.8 m wide, 30–60 cm deep) re-emerges approximately 700 m downstream at the far end of the wetland, serving as the catchment outlet.

Between September 2020 and October 2022, a total of nine field campaigns were conducted for regular sampling of water chemistry, dissolved gases, CO₂ emissions, GHG production using the dark incubation method, and high-frequency monitoring of water height and temperature.

2.2 Water chemistry

Water samples for chemical analysis were collected using standardized procedures. From each site, 500 mL unfiltered water was collected in rectangular HDPE bottles (Emballator Melledrud AB, Stockholm, Sweden) after rinsing with sample waters three times, kept dark after sampling, carried out of the field, and stored within hours after sampling at 4 °C. The samples were then transported by car and plane back to the laboratory and delivered for chemical analysis; they were kept in the laboratory at 4 °C until analysis. To ensure that our sampling procedures were suitable for the determination of nutrient concentrations, we performed a comparison of different field procedures, including acidification and filtration in the field (see the Supplement). We found no significant differences between procedures (Table S4). Chemical analysis of pH, electrical conductivity (EC), and concentrations of sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), silica (SiO₂), ammonium (NH${}_{\mathrm{4}}^{+}$), nitrate (NO${}_{\mathrm{3}}^{-}$), total phosphorous (totP), total organic carbon (TOC), DOC, and particulate organic carbon (POC) in the water samples was conducted at accredited laboratories at the Norwegian Institute for Water Research (NIVA); methods for analysis and quality control are described in the ICP Waters Programme Manual (Gundersen et al., 2025). The samples were not fully digested according to standard procedures required for the determination of total nitrogen (totN), hence totN values are expected to be underestimated and are therefore not shown in the manuscript, although the values were sufficient to confirm that the dominant form of N was organic (Thrane et al., 2020)

Absorption spectra of dissolved organic matter (DOM) were measured at NIVA for wavelengths between 200 and 900 nm, using 1 nm intervals, with a 5 cm cuvette length and Milli-Q water as a reference, using a Lambda 40 UV/Vis spectrophotometer (Perkin Elmer, USA) and expressed in absorbance per centimetre. In two samples, incomplete filtration caused excess scattering; these spectra were removed. The absorbance values at 254 nm (A_λ254 nm) were used to calculate specific UV absorbance, expressed as sUVa = A_λ254nm $/$ mg C L⁻¹, and the specific UV absorption ratio (SAR = $A_{\mathit{\lambda }\mathrm{254}\mathrm{nm}}/A_{\mathit{\lambda }\mathrm{400}\mathrm{nm}}$) was calculated for each sample.

2.3 Dissolved gas analysis

Dissolved gases (CO₂, CH₄) were sampled in the field using the acidified headspace technique (Åberg and Wallin, 2014). Duplicate gas samples were collected, according to Valiente et al. (2022), with 50 mL syringes. These were filled and sealed underwater without air bubbles to prevent gas loss. Excess water was expelled to retain 30 mL, and 20 mL of ambient air was drawn in to create a headspace. All samples were acidified with 0.6 mL of 3 % HCl to achieve pH<2, ensuring DIC was present as CO₂. Equilibrium was reached by shaking for 1 min, followed by a 30 s rest, thrice; then 15 mL of headspace gas was transferred to 12 mL evacuated vials. Water temperature in the syringe was measured immediately after gas transfer. Samples were stored at room temperature and flown to southern Norway for analysis. A 15 mL ambient air sample was taken daily for background correction.

Analysis was performed via automated gas chromatography (GC) at the Norwegian University of Life Sciences (NMBU), as described by Yang et al. (2015). A GC autosampler (GC-Pal, CTC, Switzerland) injected 2 mL headspace samples into an Agilent 7890A GC unit (Santa Clara, CA, USA) with a 20 m wide-bore PoraPLOT Q column at 38 °C, using He as the carrier gas to separate CH₄ and CO₂ from Ar, N₂, and O₂. For calibration, certified standards of CO₂ and CH₄ in He were used (AGA, Germany) and N₂, O₂, and Ar were calibrated using laboratory air. CH₄ was measured with a flame ionization detector (FID). A thermal conductivity detector (TCD) was used to measure all other gases.

Dissolved gas concentrations were calculated from headspace concentrations corrected for background air, applying temperature-adjusted Henry's law constants (Wilhelm et al., 1977) based on the recorded water temperature. At pH>4, a non-negligible amount of DIC is in the form of (bi)carbonates (HCO${}_{\mathrm{3}}^{-}$, CO${}_{\mathrm{3}}^{\mathrm{2}-}$). The bicarbonate concentrations were calculated based on pH, total dissolved CO₂ (after acidification), and the temperature-adjusted first dissociation constant (pK₁=6.41 at 25 °C; Stumm and Morgan, 1995) of the carbonic acid equilibrium. Dissolved CO₂ was calculated as DIC minus bicarbonate. To facilitate comparisons with existing studies reporting dissolved gases in microatmospheres, we converted dissolved gas concentrations to CO₂ or CH₄ saturation indices (GHG_SI), assuming atmospheric partial pressures of CO₂ and CH₄ of 400 and 1.9 µatm, respectively:

$$\begin{array}{}\text{(1)}& {\mathrm{GHG}}_{\mathrm{SI}}={\displaystyle \frac{\left[\mathrm{GHG}\right]}{{\left[\mathrm{GHG}\right]}_{\mathrm{saturation}}}},\end{array}$$

where [GHG] is the measured dissolved CO₂ or CH₄ concentration and [GHG]_saturation is the concentration of dissolved CO₂ or CH₄ at equilibrium with the corresponding atmospheric partial pressure.

2.4 Diffusive CO₂ fluxes from water to atmosphere

Measurements of CO₂ fluxes from water to atmosphere (diffusive CO₂ fluxes) were measured at each site for 30–60 min using custom-made opaque flux chambers, as described by Bastviken et al. (2015), at the water–air interface. The chamber consists of a Senseair K30 sensor (Senseair AB, Delsbo, Sweden) housed within a plastic bucket that records pCO₂, temperature, and relative humidity every 30 s. Fluxes are calculated from the linear increase in pCO₂, corrected for ambient temperature and humidity, in the chamber (Bastviken et al., 2015), considering the internal air volume and the water surface area covered by the chamber. Single measurements with a linear increase in pCO₂ with time associated with a coefficient of determination (R²) lower than 0.9 were discarded.

2.5 Dark incubation

Water samples were collected for short-term dark incubations started directly in the field and lasting between 18 and 30 h to estimate DOM mineralization and GHG production rates. Serum flasks (120 mL) were filled with 80 mL of water with a 50 mL syringe equipped with a long tube. The syringe was filled and closed under water and the water was gently pushed to the bottom of the serum flask to prevent gas loss. The remaining 40 mL were left with ambient air as headspace. The flasks were crimp-sealed with gas-tight butyl-rubber septa, sealed, covered with aluminium foil, kept at field temperature (for a maximum of 6 h), and transported back from the field to be stored at room temperature (18–20 °C). The day following the sampling (18 to 30 h after sampling), the incubations were stopped by adding 1.6 mL 3 % HCl to reach a final pH below 2, after which gas samples were taken, following the protocols described previously. Results from the dark incubation were expressed as rates of DIC production over the course of the incubation period by comparison with initial DIC concentrations and reported in units of µM h⁻¹:

$$\begin{array}{}\text{(2)}& {\mathrm{DIC}}_{\mathrm{rate}}={\displaystyle \frac{{\left[\mathrm{DIC}\right]}_{\mathrm{f}}-{\left[\mathrm{DIC}\right]}_{\mathrm{0}}}{h}},\end{array}$$

where [DIC]_f is the final solute concentration in the dark incubation, [DIC]₀ is the initial solute concentration taken in the field (µM, see Sect. 2.3), and h is the incubation duration in hours. In addition, we normalized the DIC production rate by DOC concentration to estimate DOM mineralization rates (per time unit). Also, we calculated the first-order DOM decay rate (yr⁻¹) using the exponential decay rate model (Mostovaya et al., 2017). The exponential decay model, based on early studies on sediment diagenesis (Boudreau and Ruddick, 1991; Westrich and Berner, 1984), is often the best model to describe decay rates from bioassays in closed systems (Vähätalo et al., 2010) and has been widely used to describe DOM degradation reactions. Under the exponential decay model, the decay constant (k_DOM; yr⁻¹) can be expressed as

$$\begin{array}{}\text{(3)}& k_{\mathrm{DOM}}=\ln \left({\displaystyle \frac{\mathrm{DOC}}{\mathrm{DOC}-\left({\left[\mathrm{DIC}\right]}_{\mathrm{f}}-{\left[\mathrm{DIC}\right]}_{\mathrm{0}}\right)M_{\mathrm{C}}}}\right)\times {\displaystyle \frac{\mathrm{8766}}{h}},\end{array}$$

where DOC is the DOC concentration (µg L⁻¹), M_C is the molecular mass of C (g mol⁻¹), and 8766 is the number of hours in a year. Where [DIC]_f was equal to or below [DIC]₀, we removed the values from the dataset, assuming that the temperature correction of [DIC]₀ was not precise enough (3 of 39 samples) to allow quantification of CO₂ processing rates. These occurred in September 2020 and October 2021, under cold field conditions, when [DIC]₀ was overestimated because the sample temperature in the field was unknown.

2.6 Statistical methods

Statistical analyses were conducted to evaluate differences between sites for various measured parameters. One-way analysis of variance (ANOVA) was employed to test for differences among groups. Pairwise comparisons of group means were performed using Tukey's t test using JMP 18.0.11 (2024 JMP Statistical Discovery LLC). For data that did not conform to normal distribution assumptions, non-parametric methods were applied, specifically the Wilcoxon rank-sum test, to ensure robust comparisons across sites. Results are displayed in the form of connecting letters reported within the tables. Sites with the same letter (e.g. “A” or “B”) indicate no statistically significant difference in the measured parameter between those groups at the p<0.05 significance level. Groups with different letters (e.g. “A” vs. “B”) are significantly different. When overlapping letters (e.g. “AB”) are reported, those groups are statistically similar to others with at least one shared letter but may differ from groups with entirely distinct letters. Figures were created using the ggplot2 package (Wickham, 2016) using R software (R Core Team, 2021).

3 Results

3.1 Water chemistry

The thermokarst water bodies were more acidic, richer in DOC and total P, and lower in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and SiO₂, compared with the wetland stream sites (all differences statistically significant; Table 2 and Fig. 2). The low pH of the ponds is consistent with their high DOC, and thus high organic acidity. The water bodies aligned along the inverse DOC–pH relationship, with TK-Pond 3 exhibiting the highest DOC and lowest pH, followed by TK-Pond 2 and TK-Pond 1. TK-Drain usually held an intermediate position between the thermokarst ponds and the Inlet and Outlet, which were found at the high-pH low-DOC end of the DOC–pH relationship. Similar patterns were found for the DOC–SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and DOC–SiO₂ relationships (Fig. 2). Particulate OC concentrations were significantly higher and more variable in thermokarst ponds (1.2–3.4 mg L⁻¹), compared with wetland stream sites (0.4–0.6 mg L⁻¹), with greater variability observed at the Outlet than the Inlet (Table 2).

Figure 2 Relationships between dissolved organic carbon (DOC) and various water quality parameters across different sites: scatter plots show correlations with (a) pH, (b) sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), (c) silica (SiO₂), and (d) total organic phosphorus (totP).

Table 2 Water chemistry parameters for thermokarst ponds and wetland sites during nine sampling campaigns. Median values with standard deviations are shown for all water chemistry variables, except for pH, which is shown as the median with minimum and maximum values. EC: electrical conductivity; SO${}_{\mathrm{4}}^{\mathrm{2}-}$: sulfate; SiO₂: silica; DOC: dissolved organic carbon; sUVa: specific UV absorbance; SAR: specific UV absorption ratio; TOC: total organic carbon; NH${}_{\mathrm{4}}^{+}$: ammonium; NO${}_{\mathrm{3}}^{-}$: nitrate; totP: total organic phosphorous; POC: particulate organic carbon (% of TOC). Letters indicate significant differences between sites for each variable (Tukey's t test; pairwise comparisons, p<0.05; see Sect. 2.6).

All water bodies had NO${}_{\mathrm{3}}^{-}$ concentrations at, or close to, the detection limit, while the thermokarst water bodies had considerable levels of NH${}_{\mathrm{4}}^{+}$, contrary to the wetland stream (Table 2). Total P was highest, and most variable, in the ponds, to some extent mirroring the pattern in DOC, understandably, given that in these nutrient-poor sites most P would be in an organic form, just like N.

The DOM quality indicator SAR was highest in the thermokarst ponds (p<0.03). SAR was positively strongly correlated with DOC concentration (positive, R²=0.57, p<0.0001), implying that lowest SAR was found in the wetland stream. The DOM quality indicator sUVa, a proxy for aromaticity, was slightly higher in the wetland stream than in the thermokarst ponds, although the difference was not significant.

The inverse relationship between DOC and pH points towards organic acidity as a strong driver of pH. Additionally, the near-to-neutral pH in the wetland stream is consistent with groundwater influences from the catchment, as well as the elevated SiO₂ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations. A limited set of water samples was analysed for base cations (Table S4), confirming that these were highest in the wetland stream.

The water chemical composition of the ponds mirrored the impact of thawing permafrost: the TK-Pond 3 is hydrologically most isolated, with the lowest pH, highest conductivity, and highest DOC. TK-Pond 1, located at the transition from peat plateau to wetland, had a higher pH and lower EC, DOC, and NH${}_{\mathrm{4}}^{+}$ than the other ponds; this is consistent with some hydrological influences from the wetland and hence less permafrost impact. TK-Pond 2 is located in the middle of the peat plateau and is by far the largest pond and, under wet conditions, hydrologically connected to neighbouring ponds. The water chemistry of TK-Drain was usually most similar to that of TK-Pond 1. An example of pH and EC gradients from the peat plateau into the wetland is consistent with the influence of thermokarst waterbodies gradually becoming less dominant in the transition from the peat plateau complex to the wetland (Fig. S1).

3.2 Dissolved gases and gas evasion

All water bodies were oxygenated and dissolved O₂ concentrations were, on average, 61 % to 81 % of water O₂ saturation (Table 3). The ponds are shallow, allowing for wind mixing, and they host sphagnum, suggesting active O₂ production through photosynthesis. All water bodies were oversaturated with CH₄ and CO₂. Dissolved CH₄ concentrations were 2000–5000 and ≈30 times higher than atmospheric equilibrium, in the ponds and in the wetland stream, respectively, indicating that all water bodies – thermokarst ponds in particular – are net sources of CH₄ to the atmosphere. The lower CH₄ oversaturation in the stream sites, compared with ponds, is likely to be related to higher CH₄ losses caused by stream turbulence and/or higher production rates of CH₄ in the thermokarst ponds (Fig. 3).

Figure 3 Variations in (a)–(c) dissolved carbon dioxide (CO₂) and (d)–(f) methane (CH₄) concentrations in relation to pH, DOC, and total organic phosphorus (totP) across sampling sites.

CO₂ saturation indices in the thermokarst waterbodies reached 5 to 20, while in the stream sites they ranged between 3 and 4 (Table 3). By contrast, the CO₂ evasion from the stream sites (e.g. 1–2 g C m⁻² d⁻¹) was higher than that from the ponds (0.3–0.5 g C m⁻² d⁻¹), consistent with the higher turbulence in the stream and the replenishment of CO₂ from bicarbonates from groundwater in the stream, which is a geological rather than a recent source of CO₂. Bicarbonates contributed about 60 %–70 % to DIC in the wetland stream (with pH between 6.0 and 7.4), while bicarbonates in thermokarst water bodies were almost negligible (with pH below 4.5); this is consistent with equilibrium between bicarbonates and CO₂ over these pH ranges. TK-Pond 3 had the lowest O₂ concentrations and the highest CH₄ and CO₂ concentrations of all thermokarst water bodies. Concentrations of CH₄ and CO₂ were positively related in the thermokarst water bodies (R²=0.48, p<0.0001, F test) but not in the wetland stream (Fig. 4). Note that DIC in the stream sites is also not correlated with CH₄ (Fig. S2).

The CO₂ emissions from the stream sites (Inlet and Outlet) were substantially larger than those from the thermokarst waterbodies (Table 3). The mean CO₂ emission at the Inlet site was 1.12 ± 0.46 g C m⁻² d⁻¹, and that at the Outlet site was 2.20 ± 1.15 g C m⁻² d⁻¹. These values are 3 to 7 times higher than the fluxes observed from the thermokarst waterbodies, which ranged from 0.30 ± 0.22 g C m⁻² d⁻¹ (TK-Pond 2) to 0.51 ± 0.28 g C m⁻² d⁻¹ (TK-Pond 3).

Annual CO₂ fluxes for the ice-free period, assuming negligible flux during the ice-covered months, ranged between 51 and 87 g C m⁻² yr⁻¹ for the thermokarst ponds, and those for the stream ranged between 190 g C m⁻² yr⁻¹ (Inlet) and 405 g C m⁻² yr⁻¹ (Outlet).

3.3 DOM processing rates

Average DIC production rates in the different water bodies were highly variable (7.8–62.5 µM d⁻¹, Table 4), but tended to be highest in the thermokarst ponds, compared with the wetland stream, while TK-Drain had the lowest rates (Tukey's t test, p<0.05). The non-parametric Wilcoxon tests supported these trends, confirming minimal site-specific effects overall, with TK-Drain showing lower activity. These results reflect the in situ processing of DOM in both thermokarst ponds and stream sites. The DOM mineralization rate did not vary significantly between sites and neither did the exponential decay rate k_DOM (Table 4). The value of k_DOM in the thermokarst ponds ranged from 4.4 to 9.0 yr⁻¹ (Table 4), while the Outlet showed higher values of k_DOM than the Inlet (8.8 and 6.3 yr⁻¹, respectively). The TK-Drain site was substantially lower (2.5 yr⁻¹).

4 Discussion

Understanding the GHG source–sink function of degrading permafrost landscapes benefits from an integrated study of water chemistry and GHG fluxes, as hydrological and biogeochemical processes are closely linked (Frey and McClelland, 2009; Vonk et al., 2015). In particular, shifts in OM mobilization, acidity, and nutrient dynamics across different thermokarst pond stages influence C cycling and GHG production. In our study area, isolated thermokarst ponds and the more hydrologically connected wetland stream represent contrasting hydrochemical environments, with ponds reflecting strong permafrost thaw inputs and the stream influenced by the catchment.

4.1 Water chemical contrasts between thermokarst ponds and water chemistry

Thermokarst ponds and the wetland stream exhibit strong contrasts in DOC concentration and acidity, with ponds showing high DOC and low pH. The inverse relationship between DOC and pH ($r=-\mathrm{0.82}$, p<0.001) suggests that organic acidity is a dominant driver of pH in thermokarst ponds. The acidity is probably driven by the leaching of DOM from recently destabilized permafrost, since the ponds are hydrologically isolated from the surrounding wetland. Elevated DOM and leaching from surrounding permafrost are observed in other thawing permafrost landscapes (Holmes et al., 2022; Ward and Cory, 2015). This is consistent with the finding that TK-Pond 3 – hydrologically the most isolated – has, on average, the highest DOC (35 mg L⁻¹) and the lowest pH (4.1), indicating strong organic acidity effects from destabilized permafrost and minimal exchange with the surrounding wetland.

In contrast, the wetland stream exhibits near-neutral pH values, which we attribute to groundwater influence. We found systematically higher concentrations of HCO${}_{\mathrm{3}}^{-}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, and SiO₂ in the wetland stream than in the ponds, and also higher concentrations of base cations (Table S3). This geochemical signature is characteristic of carbonate and silicate mineral weathering occurring along subsurface flow paths, which was found in several catchments at Iškorasfjellet, including the Báhkiljohka catchment where our study site is located (Lehmann et al., 2023). Lehmann et al. (2023) documented groundwater-driven alkalinity generation linked to both carbonate vein dissolution and silicate weathering and suggested that carbonate weathering should be considered as a potential CO₂ source in the catchment. Groundwater effects on stream waters have also been found elsewhere in permafrost landscapes (Turetsky et al., 2020; Vonk et al., 2015), suggesting that CO₂ emissions from high-pH streams can be replenished by geogenic rather than biogenic sources; this is important to account for in GHG budgets from aquatic ecosystems.

Figure 4 Relationships between CH₄ and CO₂ (a) between thermokarst waterbodies and (b) for Inlet and Outlet sites, including linear regression lines and corresponding R² and p statistics. Note the scale differences for CH₄ between thermokarst waterbodies and the wetland stream.

Table 3 Mean concentrations and standard deviations of CO₂ (µM) and CH₄ (µM), with their respective saturation ratios, along with CO₂ emission flux (g C m⁻² d⁻¹), DIC (µM), and oxygen concentrations (µM) with percentage saturation. The saturation ratio is defined as the concentration divided by the equilibrium concentration between the atmosphere and water at the given temperature. For this study, DIC is considered to be the sum of dissolved CO₂ and bicarbonate. Letters indicate significant differences between sites for each variable.

Levels of totP mirrored the pattern of DOC enrichment in thermokarst ponds (30–70 µg L⁻¹), suggesting that P in these nutrient-poor ponds is primarily organically bound (Frey et al., 2007). In addition to DOC, permafrost thaw releases organic forms of P that can affect downstream nutrient dynamics and carbon cycling (in 't Zandt et al., 2020). In contrast, the lower DOC and totP found in the wetland stream reflected the influence from the catchment upstream, including groundwater inputs. The higher POC concentrations in the thermokarst ponds (1.2–3.4 mg L⁻¹), compared with wetland stream sites (0.4–0.6 mg L⁻¹), further supports that thermokarst ponds are hotspots for OM destabilization, whereas wetland streams are influenced more by lateral transport and groundwater (Olefeldt and Roulet, 2014).

The hydrochemical contrasts between thermokarst ponds and the wetland stream sites at Iškoras, shaped by differences in DOC concentration, acidity, and groundwater influence, are key drivers of spatial variation in GHG production and emissions across the landscape. Given the strong controls of DOC and pH on CH₄ dynamics (Shirokova et al., 2012; Segers, 1998), these patterns provide important context for understanding permafrost-C feedbacks.

4.2 Thermokarst ponds as hotspots of methane emissions

Thermokarst ponds in Iškoras display CH₄ saturation indices of 2300 to 5000 (Table 3), which are among the highest values reported in the literature for natural waterbodies, particularly in northern permafrost regions. These findings align with Shirokova et al. (2012) and Matveev et al. (2018), who documented saturation indices of 50 to 5000 in Siberian thermokarst depressions and of 5 to 50 in subarctic lithalsa lakes (thermokarst lakes in mineral soils). Such high CH₄ concentrations in these poorly connected, small, and relatively protected water bodies are consistent with the established inverse relationship between CH₄ concentration and water body size, hydrological connectivity, and turbulence exposure (Abnizova et al., 2012; Kankaala et al., 2013; Polishchuk et al., 2018).

At Iškoras, the smallest pond, TK-Pond 3, exhibited the highest CH₄ oversaturation, in combination with the highest DOC and totP values and lowest pH, as mentioned earlier. This could be explained by the creation of anaerobic sediments and C-rich conditions, and low pH, by permafrost thaw and limited hydrological connectivity, that, in concert, enhance CH₄ production. The high DOM, originating from destabilized permafrost (Turetsky et al., 2020), is usually associated with low pH and elevated totP (Holmes et al., 2022; Ward and Cory, 2015). Thus, CH₄ oversaturation in thermokarst ponds could be related to particularly high DOC availability and to acidic conditions that limit CH₄ oxidation (Wik et al., 2016). Despite the presence of dissolved O₂ in thermokarst ponds, they remained highly oversaturated in CH₄. CH₄ production is known to occur mainly in anoxic sediments (Bastviken et al., 2004; Clayer et al., 2016; Wik et al., 2016), from where CH₄ is subsequently transported to overlying water. The microbial activity responsible for CH₄ production may be enhanced by the release of previously frozen OM from ongoing thermokarst development (Crevecoeur et al., 2017), as recently observed in laboratory incubations with inundated peat from the Iškoras site (Kjær et al., 2024).

As thermokarst ponds evolve into wetlands, CH₄ emission patterns may shift, due to changes in hydrological connectivity and biogeochemical cycling. Pirk et al. (2024) used a space-for-time substitution to highlight that the transition from thermokarst ponds to wetlands at Iškoras involves significant changes in GHG fluxes, finding that the degradation of palsas to thermokarst ponds led to a 17-fold increase in local GHG forcing, primarily driven by CH₄ and CO₂ emissions. This is partly because thermokarst ponds, being spatially isolated, create localized CH₄ emission hotspots (Elder et al., 2021). Wetlands, with their larger spatial extent and greater hydrological connectivity, promote slower organic matter mineralization, which can reduce CO₂ emissions and increase carbon uptake (Pirk, et al., 2024). Although wetlands continue to emit CH₄, these emissions become more diffuse, rather than concentrated at hotspots. At the same time, wetlands can act as net carbon sinks, as CO₂ uptake through plant productivity and organic matter accumulation may offset greenhouse gas emissions (Heiskanen et al., 2023; de Wit et al., 2015). This transition represents a fundamental shift from localized CH₄-dominated GHG emissions in thermokarst ponds to more spatially diffuse CH₄ emissions in wetlands, where increased CO₂ uptake and OM accumulation can contribute to a net C sink (Turetsky et al., 2020).

4.3 Carbon dioxide dynamics

Both thermokarst ponds and the stream in Iškoras are oversaturated with CO₂ (Table 3), a common feature of Arctic and subarctic aquatic systems (Allesson et al., 2022; Bastviken et al., 2004). However, the mechanisms driving CO₂ fluxes differ between ponds and streams. In ponds, despite high CO₂ concentrations, CO₂ release is lower than in streams and is likely to be limited by the lack of turbulence. In contrast, streams exhibit enhanced CO₂ fluxes, probably due to high turbulence and carbonate inputs from groundwater, as described in Sect. 4.1. Others have also found that groundwater influences can contribute both bicarbonate and dissolved CO₂ to sustain stream CO₂ fluxes (Lehmann et al., 2023; Duvert et al., 2018; Winterdahl et al., 2016). This highlights that groundwater contributions of bicarbonate and dissolved CO₂ must be accounted for when interpreting stream CO₂ fluxes, as ignoring these inputs could lead to overestimating the importance of recently mineralized DOM or newly produced CO₂ from thermokarst ponds.

Quantitatively, CO₂ efflux from the stream at Iškoras averages 0.4 g C m⁻² d⁻¹ (Table 3), aligning closely with values observed in Siberian permafrost streams (0.3–0.5 g C m⁻² d⁻¹; Shirokova et al., 2012). This flux is consistent with the combined influence of turbulent flow and groundwater inputs enriched in DIC, consistent with earlier findings of mineral weathering contributions to DIC at Iškorasfjellet (Lehmann et al., 2023). These findings complement observations from other regions, where turbulence and bicarbonate supply are key drivers of CO₂ release in streams (Lundin et al., 2013; Raymond et al., 2013). Streams probably benefit from continuous replenishment of CO₂ from bicarbonates, which account for 60 %–70 % of the DIC pool in boreal to Arctic streams and rivers (Wallin et al., 2018; Zolkos and Tank, 2020).

In contrast, CO₂ emissions from ponds at Iškoras are notably low (Table 3), similar to findings of Campeau and Del Giorgio (2014), attributed to the ponds' high DOC-to-bicarbonate ratios, which restrict bicarbonate formation and subsequent CO₂ production (Abnizova et al., 2012; Bastviken et al., 2004). The limited water mixing in ponds further diminishes CO₂ flux, due to low gas exchange rates compared with streams. However, DIC production rates in ponds (26.7–35.7 µM d⁻¹) remain sufficient to sustain CO₂ effluxes of 5–10 mol C m⁻² yr⁻¹, consistent with findings by Shirokova et al. (2012), who linked DIC production directly to DOC concentrations in Arctic aquatic systems. The relatively high values of k_DOM that we observed (4.4–9.0 yr⁻¹; Table 4) indicate faster DOM degradation than the global lake average (2.5±4.0 yr⁻¹), based on Catalán et al. (2016), and are consistent with elevated DIC production rates.

Table 4 Rates of production and decay. DIC rates reflect DIC production. The ratio DIC rate $/$ DOC indicates the relative efficiency of converting DOC to DIC, while k_DOM indicates the exponential decay rate of DOM, showing how quickly DOM is decomposed over time. Letters indicate significant differences between sites for each variable (p<0.05).

High CO₂ evasion from lower-order streams is caused by high flow velocities and associated turbulence causing high gas exchanges (Schelker et al., 2016). These drivers probably explain the relatively higher CO₂ efflux from the Iškoras stream, compared with ponds. We suggest that mineralization of DOM, influenced by DOM lability, pH, and nutrients, plays a critical role for CO₂ emissions of the thermokarst ponds. High-molecular-weight DOM, as indicated by low sUVa and SAR values, can be more resistant to microbial processing than low-molecular-weight DOM, thereby slowing DOM decay rates (Shirokova et al., 2019). Ward and Cory (2015) noted that DOM from thawing permafrost, while less aromatic and more labile, compared with active-layer DOM, may become limited by environmental factors, such as pH and nutrient availability, resulting in lower mineralization rates. The acidic conditions in ponds could shift microbial communities and affect activity (Vigneron et al., 2019), a factor that can further limit CO₂ fluxes from ponds, in addition to the low gas exchange rate.

In summary, the higher CO₂ fluxes observed in the stream sites at Iškoras are probably driven by the combined effects of turbulent flow, groundwater-derived DIC, and mineral weathering inputs, whereas lower emissions from ponds are shaped by organic acidity, limited hydrological connectivity, and limited surface exchange. These dynamics emphasize the need to account for the interaction of hydrological and chemical factors when assessing the fate of destabilized OM in water bodies in permafrost-impacted regions.

4.4 Climate feedback implications

The transition from thermokarst ponds to wetlands modifies the overall GHG footprint of the peatland–wetland continuum, balancing the loss of localized CH₄ emission hotspots with the emergence of sustained long-term CH₄ emissions from wetlands, while the fate of organic matter currently stored in permafrost remains uncertain. At the same time, CO₂ fluxes from the streams and rivers may increase, due to enhanced hydrological connectivity and increased organic matter input (Zolkos et al., 2019), in agreement with the results of our study. At Iškoras, the small spatial extent of the permafrost area limits its overall impact on the catchment-scale GHG source–sink function. However, in landscapes where peat plateaus occupy a larger area, such transitions may have more significant consequences at the regional scale. These findings reflect the complex interplay of ecological and hydrological factors shaping GHG emissions in permafrost landscapes. Turetsky et al. (2020) and Pirk et al. (2024) both emphasized the need for further research on the spatiotemporal variability of these factors, particularly during thaw cycles, as shifts in hydrological connectivity, OM transport, and microbial activity can significantly influence the GHG emissions and permafrost-C feedbacks. Improving our understanding of these dynamics is essential for refining predictions of permafrost-C feedbacks in a changing climate.

5 Conclusions

This study highlights the distinct biogeochemical roles of thermokarst ponds and wetland streams in a landscape of sporadic permafrost in subarctic Norway. Thermokarst ponds at the Iškoras site, characterized by low pH, high organic acidity, and elevated DOC concentrations, are currently hotspots for CH₄ emissions, with stable DOM lability driving sustained carbon processing. In contrast, the wetland stream exhibits higher CO₂ fluxes, largely driven by turbulence and bicarbonate replenishment from groundwater. Despite similarities in DOM mineralization rates between ponds and the stream sites, environmental constraints, such as pH, microbial community composition, and hydrodynamic mixing, are likely controls of the observed differences in GHG fluxes. As thermokarst ponds transition into wetlands, they will no longer function as hotspots for CH₄ emissions. Instead, CH₄ emissions are likely to increase across the entire landscape, as sustained waterlogging promotes elevated CH₄ production. These ecological shifts, coupled with lateral DOC losses from peat plateaus, highlight the importance of hydrological connectivity in linking terrestrial and aquatic C dynamics. Such transitions emphasize the need for integrated C budget models that account for the evolving contributions of small aquatic systems to regional and global C cycles.

Future research should prioritize direct measurements of CH₄ fluxes, microbial community contributions to DOM decomposition under varying environmental constraints, and the temporal variability of gas production and emissions. Additionally, exploring seasonal dynamics, lateral carbon transport, and hydrological processes will provide critical insights into C cycling. Investigating catchment-scale signals, such as DOC concentrations across entire river systems and their links to permafrost contributions, can further advance our understanding of landscape-level processes. By addressing these questions, we can better predict the trajectory of permafrost-impacted landscapes and their feedbacks to the global C cycle in a warming climate.