Ignoring carbon emissions from thermokarst ponds results in overestimation of tundra net carbon uptake

Arctic permafrost landscapes have functioned as a global carbon sink for millennia. These landscapes are very heterogeneous, and the omnipresent water bodies within them act as a carbon source. Yet, few studies have focused on the impact of these water bodies on the landscape carbon budget. We deepen our understanding of carbon emissions from thermokarst ponds and constrain their impact by comparing carbon dioxide and methane fluxes from these ponds to fluxes from the surrounding tundra. We use eddy covariance measurements from a tower located at the border between a large pond and semi-terrestrial tundra.

When we take the open-water areas of thermokarst ponds into account, our results show that the estimated summer carbon uptake of the polygonal tundra is 11 % lower. Further, the data show that open-water methane emissions are of a similar magnitude to polygonal tundra emissions. However, some parts of the pond's shoreline exhibit much higher emissions. This finding underlines the high spatial variability in methane emissions. We conclude that gas fluxes from thermokarst ponds can contribute significantly to the carbon budget of Arctic tundra landscapes. Consequently, changes in the water body distribution of tundra landscapes due to permafrost degradation may substantially impact the overall carbon budget of the Arctic.

1 Introduction

Water bodies make up a significant part of the Arctic lowlands with an areal coverage of about 17 % (Muster et al., 2017) and act as an important carbon source in a landscape that is an overall carbon sink (Kuhn et al., 2018). Intensified permafrost thaw in the warming Arctic will change the distribution of water bodies and thereby change their contribution (Andresen and Lougheed, 2015; Bring et al., 2016) to the landscape carbon budget (Kuhn et al., 2018) of tundra landscapes. However, data on greenhouse gas emissions from Arctic water bodies are still sparse, especially data with a high temporal resolution and from non-Yedoma regions (Vonk et al., 2015).

Our study site in the Lena River delta, Siberia, is located on an island mostly characterized by non-Yedoma polygonal tundra (Fig. 1). This landscape features many ponds; we define ponds as water bodies with an area of less than 8×10⁴ m², following Ramsar Convention Secretariat (2016) and Rehder et al. (2021). Within our area of interest, ponds cover about the same area as lakes (Abnizova et al., 2012; Muster et al., 2012). The ponds on Samoylov Island have formed almost exclusively through thermokarst processes: the soil has a high ice content, so when the ice melts, the ground subsides, and thermokarst ponds form (Ellis et al., 2008). These thermokarst ponds are often only as large as one polygon (polygonal ponds). When several polygons are inundated, this can cause larger shallow thermokarst ponds to form, which we term merged polygonal ponds (Rehder et al., 2021). Holgerson and Raymond (2016) as well as Wik et al. (2016) report that ponds emit more greenhouse gases per unit area than lakes, defined here as water bodies with an area larger than 8×10⁴ m². Thus, in our study area, they have greater potential than lakes to counterbalance the carbon uptake of the surrounding tundra (McGuire et al., 2012; Jammet et al., 2017; Kuhn et al., 2018). To better understand the impact of thermokarst ponds on the landscape carbon flux, we compare carbon dioxide (CO₂) and methane (CH₄) fluxes from thermokarst ponds to fluxes from the semi-terrestrial tundra. The semi-terrestrial tundra consists of wet and dry tundra and overgrown shallow water, which are the terrestrial land-surface types used by Muster et al. (2012) to classify Samoylov Island.

The main geophysical and biochemical processes that drive CH₄ fluxes are different to the ones that drive CO₂ fluxes. The microbial decomposition of dissolved organic carbon, which is introduced laterally into the aquatic system through rain and meltwater (Neff and Asner, 2001), dominates aquatic CO₂ production. When supersaturated with dissolved CO₂, ponds emit CO₂ into the atmosphere through diffusion. While photosynthetic CO₂ uptake has been observed in some clear Arctic water bodies (Squires and Lesack, 2003), most Arctic water bodies are net CO₂ sources (Kuhn et al., 2018). Estimates of CO₂ emissions range from close to zero (0.028 g m² d⁻¹ by Treat et al., 2018, and 0.059 g m² d⁻¹ by Jammet et al., 2017) to substantial (1.4–2.2 g m² d⁻¹ by Abnizova et al., 2012).

Within just one site, CH₄ emissions from a water body can vary by up to 5 orders of magnitude: 0.5–6432 mg m² d⁻¹ (Bouchard et al., 2015). The CH₄ that ponds emit is mostly produced in sub-aquatic soils and anoxic bottom waters (Conrad, 1999; Hedderich and Whitman, 2006; Borrel et al., 2011). Additionally, CH₄ might also be produced in the oxic water column (Bogard et al., 2014; Donis et al., 2017), though this location of methanogenesis is only significant in large water bodies (Günthel et al., 2020). Moreover, there is still ongoing debate as to whether methanogenesis occurs in oxic waters at all (Encinas Fernández et al., 2016; Peeters et al., 2019). CO₂ is also formed as a byproduct of the methanogenesis process (Hedderich and Whitman, 2006). Water bodies emit CH₄ produced in their benthic zone through diffusion, ebullition (sudden release of bubbles), or plant-mediated transport. The varying contributions of these three local methane emissions pathways lead to high spatial variability between water bodies and within a single water body (Sepulveda-Jauregui et al., 2015; Jansen et al., 2019). In particular, local seep ebullition causes high spatial variance of CH₄ emissions within one water body (Walter et al., 2006). Variability in the coverage and composition of vascular plant communities in a water body can also increase CH₄ variability because CH₄ transport efficiency can be species-specific (Knoblauch et al., 2015; Andresen et al., 2017).

To study spatial and temporal patterns of carbon emissions from thermokarst ponds, we analyzed land–atmosphere CO₂ and CH₄ flux observations from an eddy covariance (EC) tower on Samoylov Island, Lena River delta, Russia. We set up the EC tower within the polygonal tundra landscape at the border between a large merged polygonal pond and the surrounding semi-terrestrial tundra for 2 months in summer 2019. The polygonal structures were still clearly visible along the shore and underwater, and most of the pond was shallow (Rehder et al., 2021). Due to the tower's position, fluxes from the merged polygonal pond were the dominant source of the observed EC fluxes under easterly winds. From other wind directions, the observed EC fluxes were dominated by semi-terrestrial polygonal tundra with only a low influence from small polygonal ponds. This paper aims to deepen the understanding of carbon emissions from thermokarst ponds and constrain their impact on the landscape carbon balance. We (1) examine the temporal and spatial patterns of net ecosystem exchange (NEE) and the spatial pattern of CH₄ flux from semi-terrestrial tundra and thermokarst ponds and (2) investigate the influence of the thermokarst ponds on the landscape NEE of CO₂ during the months June to September 2019. To this end, we use a footprint model and model NEE of CO₂ using the footprint weights of semi-terrestrial tundra and thermokarst ponds.

2 Methods

2.1 Study site

Samoylov Island (72^∘22^′ N, 126^∘28^′ E) is located in the southern part of the Lena River delta (Fig. 1b). It is approximately 5 km² in size and consists of two geomorphologically different components. The western part of the island (∼ 2 km²) is a floodplain, which is flooded annually during the spring. The eastern part of the island (∼ 3 km²), a Late Holocene river terrace, is characterized by polygonal tundra. The partially degraded polygonal tundra at this study site features high spatial heterogeneity on a scale of a few meters in several aspects, including vegetation, water table height, and soil properties. Dry and wet vegetated parts of the semi-terrestrial tundra are interspersed with small and large thermokarst ponds (1 to 10 000 m²) and with larger lakes (up to 0.05 km²; Boike et al., 2015a; Kartoziia, 2019). The island is surrounded by the Lena River and sandy floodplains, creating additional spatial heterogeneity on a larger scale.

This study focuses on a merged polygonal pond (Figs. 1d and A1) on the eastern part of the island. This merged polygonal pond has a size of 0.024 km² with a maximum depth of 3.4 m and a mean depth of 1.2 m (Rehder et al., 2021; Boike et al., 2015a). In an aerial image of the pond, the polygonal structures are clearly visible under the water's surface (Boike et al., 2015c). The vegetated shoreline of this merged polygonal pond is dominated by Carex aquatilis, but it also features Carex chordorrhiza, Potentilla palustris, and Aulacomnium spp. These plants grow in the water near the shore, while the deeper parts of the merged polygonal pond are vegetation-free.

2.2 Instruments

We measured gas fluxes using an eddy covariance (EC) tower between 11 July and 10 September 2019. The EC tower was located on the eastern part of Samoylov Island, directly on the western shore of the merged polygonal pond (Fig. 1d). The EC instruments were mounted on a tripod at a height of 2.25 m (Fig. A1). The tower was equipped with an enclosed-path CO₂–H₂O sensor (LI-7200, LI-COR Biosciences, USA), an open-path CH₄ sensor (LI-7700, LI-COR Biosciences, USA), and a 3D ultrasonic anemometer (R3-50, Gill Instruments Limited, UK). All instruments had a sampling rate of 20 Hz. We also installed radiation-shielded temperature and humidity sensors at the EC tower (HMP155, Vaisala, Finland) and used data from a photosynthetically active radiation (PAR) sensor mounted on a tower approximately 500 m to the west of the EC tower (SKP 215, Skye Instruments, UK). Additional meteorological data for Samoylov Island were provided by Boike et al. (2019).

Figure 1 The location of the study site in Russia is shown in (a), and the location of Samoylov Island within the Lena River delta is shown in (b). Samoylov Island is shown in (c); the surrounding Lena River appears in light blue. The outline of the river terrace land-cover classification (Sect. 2.4.1) is indicated by the blue line. We focus on the polygonal tundra; however, large lakes are excluded (circled in yellow). In (d), the land-cover classification is drawn in blue (open water) and green (dark green: dry tundra; medium green: wet tundra; light green: overgrown water) shades. The merged polygonal pond studied here is outlined in red. The location of the EC tower is marked by a black cross. The cumulative footprint (see Sect. 2.4.2) is shown in gray shades. Of the flux, 30 % likely originated from within the dark gray area, 50 % from within the medium dark gray area, 70 % from within the medium light gray area, and 90 % from within the light gray area. Map data from © OpenStreetMap contributors 2020, distributed under the Open Data Commons Open Database License (ODbL) v1.0 (a, b) and modified based on Boike et al. (2012) (c, d).

2.3 Data processing

We performed the raw data processing and computation of half-hourly fluxes for open-path and enclosed-path fluxes (CO₂, CH₄, and H₂O) using EddyPro 7.0.6 (LI-COR, 2019). The convention of this software is that positive fluxes are fluxes from the surface to the atmosphere, while negative fluxes indicate a flux from the atmosphere downwards. Raw data screening included spike detection and removal according to Vickers and Mahrt (1997) (1 % maximum accepted spikes and a maximum of three consecutive outliers). Additionally, we applied statistical tests for raw data screening, including tests for amplitude resolution, skewness and kurtosis, discontinuities, angle of attack, and horizontal wind steadiness. All of these tests' parameters were set to EddyPro default values. We rotated the wind-speed axis to a zero-mean vertical wind speed using the “double rotation” method of Kaimal and Finnigan (1994). Further, we applied linear de-trending to the raw data following Gash and Culf (1996) before performing flux calculations. We compensated for time lags via automatic time lag optimization using a time lag assessment file from a previous EddyPro run. In this previous time lag assessment, the time lags for all gases were detected using covariance maximization (Fan et al., 1990), resulting in time lags between 0 and 0.4 s for CO₂ and between −0.5 and +0.5 s for CH₄. For H₂O, the time lag was humidity-dependent and was calculated for 10 humidity classes. We compensated for air-density fluctuations due to thermal expansion and contraction and varying water-vapor concentrations, following Webb et al. (1980). This correction depends on accurate measurements of the latent and sensible heat flux and was applied to the open-path data of the LI-7700. For the LI-7700 in particular, the correction term can be larger than the flux itself, but the correction was derived from the underlying physical equations. Because we used well-calibrated instruments as well as EddyPro, which uses an up-to-date implementation of the correction, we were confident that the LI-7700 would provide accurate CH₄ flux estimates. For enclosed-path data, we performed a sample-by-sample conversion into mixing ratios to account for air-density fluctuations (Ibrom et al., 2007b; Burba et al., 2012). Flux losses occurred in the low- and high-frequency spectral range due to different filtering effects. We compensated for flux losses in the low-frequency range in accordance with Moncrieff et al. (2004) and in the high-frequency range in accordance with Fratini et al. (2012). For the high-frequency range compensation method, a spectral assessment file was created using the method of Ibrom et al. (2007a). The spectral assessment resulted in cutoff frequencies of 3.05 and 1.67 Hz for CO₂ and CH₄, respectively. For H₂O, we found a humidity-dependent cutoff frequency between 1.25 Hz (relative humidity, RH, of 5 %–45 %) and 0.21 Hz (RH 75 %–95 %). We performed a quality check on each half-hourly flux following the 0–1–2 system proposed by Mauder and Foken (2004). In this quality check, flux intervals with the lowest quality received the flag “2” and were excluded from further analysis.

2.4 Data analysis

2.4.1 Land-cover classification

The land-cover classification covers the Late Holocene river terrace of Samoylov Island (3.0 km², area within the blue line in Fig. 1c). It is based on high-resolution near-infrared (NIR) orthomosaic aerial imagery obtained in the summer of 2008 (Boike et al., 2015b). We used a subset of the existing classification of Muster et al. (2012) as a training dataset to perform semi-supervised land-cover classification using the maximum likelihood algorithm in ArcMap Version 10.8 (Esri Inc., USA). We then applied the ArcMap majority filter tool to the new classification. The land-cover classification has a resolution of 0.17 m × 0.17 m,. It is projected onto WGS 84 UTM Zone 52N, and the land-cover classes include open water (15.7 %), overgrown water (7.0 %), dry tundra (65.1 %), and wet tundra (12.1 %), as defined by Muster et al. (2012). We summarize the classes overgrown water, dry tundra, and wet tundra in the land-cover type of semi-terrestrial tundra. The river terrace consists of this semi-terrestrial tundra, large lakes, and thermokarst ponds. Since small ponds are an integral part of the polygonal tundra, we use the term “polygonal tundra” to refer to the area of the river terrace covered by semi-terrestrial tundra and by thermokarst ponds.

2.4.2 Footprint model

In deploying an EC measurement tower, the tower's location and sensor height are crucial parameters. A lower measurement height results in a smaller footprint. The tower's footprint describes the source area of the flux within the surrounding landscape. As we installed sensors at a height of 2.25 m next to the merged polygonal pond, we expected to observe substantial flux signals from the adjacent water body as well as from the surrounding semi-terrestrial tundra. Each land-cover type's contribution to the flux signal depended on the wind direction and turbulence characteristics. We implemented the analytical footprint model proposed by Kormann and Meixner (2001) in MATLAB (2019). We combined the footprint model with land-cover classification data described in Sect. 2.4.1 to estimate the contribution of each land-cover type to each half-hourly flux (from now on referred to as the weighted footprint fraction). The model accounted for the stratification of the atmospheric boundary layer and required a height-independent crosswind distribution and horizontal homogeneity of the surface. The input data required stationarity of atmospheric conditions during the flux intervals of 30 min.

We derived the vertical power-law profiles for the eddy diffusivity and the wind speed for each 30 min flux depending on the atmospheric stratification (see Eq. 6 in Kormann and Meixner, 2001). We used an analytical approach to find the closest Monin–Obukhov (M–O) similarity profile (see Eq. 36 in Kormann and Meixner, 2001). Next, we calculated a two-dimensional probability density function of the source area for each flux (from Eqs. 9 and 21 in Kormann and Meixner, 2001). We combined each probability density function with the land-cover classification of Samoylov Island's river terrace with its four land-cover types (see Sect. 2.4.1). The resolution of the footprint model was set to the land-cover classification resolution of 0.17 m × 0.17 m. Hence, we were able to estimate how much a given grid cell contributed to each 30 min flux. We also knew each grid cell's dominant land-cover type from the land-cover classification. We combined both pieces of information for each grid cell and calculated the sum of the fraction fluxes within the source area for each of the four land-cover types (dry tundra, wet tundra, overgrown water, and open water) and determined the contribution of each land-cover type with respect to each 30 min flux (a_{dry tundra}, a_{wet tundra}, a_{overgrown water}, and a_{open water}). We refer to this contribution of each land-cover type as the weighted footprint fraction.

We also summed all 30 min two-dimensional probability density functions over the entire deployment time. This sum is referred to as the cumulative footprint (gray-shaded area in Fig. 1c–d).

2.4.3 Gap-filling the CO₂ flux

To gap fill the net ecosystem exchange (NEE) fluxes of CO₂, we used the bulk-NEE model proposed by Runkle et al. (2013). The model is specifically designed to model NEE in Arctic regions: it takes impacts of the polar day into account by allowing both respiration and photosynthesis to occur simultaneously throughout the day. The bulk-NEE model uses the sum of total ecosystem respiration (TER) and gross primary production (GPP) to describe NEE, our target variable:

$$\begin{array}{}\text{(1)}& \mathrm{NEE}=\mathrm{TER}+\mathrm{GPP},\end{array}$$

where TER and GPP are in units of µmol m⁻² s⁻¹. TER is approximated as an exponential function of air temperature T_air:

$$\begin{array}{}\text{(2)}& \mathrm{TER}=R_{\mathrm{base}}\cdot Q_{\mathrm{10}}^{\frac{T_{\mathrm{air}}-T_{\mathrm{ref}}}{\mathit{\gamma }}},\end{array}$$

where T_ref=15 ^∘C and γ=10 ^∘C are constant, independent parameters. R_base (µmol m⁻² s⁻¹) describes the basal respiration at the reference temperature T_ref, and Q₁₀ (dimensionless) describes the sensitivity of ecosystem respiration to air temperature changes.

GPP is described as a rectangular hyperbolic function of PAR (µmol m⁻² s⁻¹):

$$\begin{array}{}\text{(3)}& \mathrm{GPP}=-{\displaystyle \frac{P_{\max }\cdot \mathit{\alpha }\cdot \mathrm{PAR}}{P_{\max }+\mathit{\alpha }\cdot \mathrm{PAR}}}\end{array}$$

where α (µmol µmol⁻¹) is the initial canopy quantum use efficiency (slope of the fitted curve at PAR =0) and P_max (µmol m⁻² s⁻¹) is the maximum canopy photosynthetic potential for PAR →∞.

The parameters R_base, Q₁₀, P_max, and α were fitted simultaneously. To account for seasonal changes in plant physiology, we fitted the parameters for 5 d running windows as proposed in Holl et al. (2019).

We split the datasets into training (70 %) and validation (30 %) datasets to test model performance. We implemented the bulk-NEE model in MATLAB 2019b (MATLAB, 2019) using the fit function with the NonLinearLeastSquares fitting method. We used the coeffvalues function to estimate the four parameters (R_base, Q₁₀, P_max, and α) and the confint function to estimate their 95 % confidence bounds. All partitioned fluxes were converted into CO₂–C fluxes in units of g m⁻² d⁻¹ before data analysis.

2.4.4 Separating CO₂ fluxes from semi-terrestrial tundra and water bodies

We wanted to extract fluxes from thermokarst ponds and semi-terrestrial tundra to analyze the influence of thermokarst ponds on the carbon balance of a polygonal tundra landscape. However, due to the strong heterogeneity of the landscape and the relatively small size of the merged polygonal pond compared to the EC footprint, we measured a mixed signal from all wind directions. In other words, each flux that was measured with the EC method contained information from different land-cover types. We divided the footprint into two classes – semi-terrestrial tundra and thermokarst ponds – to assess the impact of thermokarst ponds on the carbon balance.

Similar approaches of analyzing heterogeneous eddy covariance fluxes in Arctic environments have been conducted for CO₂ and CH₄ (e.g., Rößger et al., 2019a, b; Tuovinen et al., 2019). Rößger et al. (2019a, b) extracted CO₂ and CH₄ fluxes from two different land-cover classes on a floodplain, while Tuovinen et al. (2019) separated CH₄ fluxes from nine individual land-cover classes, including water, and combined them into four source classes (with no separate class for water). All three studies differentiate between fluxes from different vegetation types. Our method is dedicated to distinguishing between fluxes from semi-terrestrial tundra and water bodies.

To estimate CO₂ fluxes from the merged polygonal pond (F_pond), we first fitted the bulk-NEE model to training data, excluding fluxes from the direction of the merged polygonal pond (30^∘ < WD < 150^∘, where WD denotes wind direction). We obtained a dataset consisting of information about as much semi-terrestrial tundra as possible. We performed this step since we expected little to no photosynthetic activity in the open-water part of the merged polygonal pond. This gap-filled CO₂ flux (hereinafter F_modeled,mix) represents the polygonal tundra surrounding the EC tower, meaning the flux is dominated by semi-terrestrial tundra, but also includes polygonal ponds from the wind directions of north, west, and south. In the model input, we excluded 30 min CO₂ fluxes with an absolute value of more than 4 g m⁻² d⁻¹. In 38 windows of 5 d duration, we found an R² above 0.9 between the model output and the validation set. In 18 cases, we obtained an R² between 0.8–0.9; in six instances, we obtained an R² below 0.7. The final RMSE between the model input and the gap-filled NEE had a value of 0.29 g m⁻² d⁻¹.

We assumed that the total observed flux was a linear combination of the fluxes from the land-cover types weighted by their respective contribution to the footprint. Thus, we postulated that the observed CO₂ flux (F_obs,mix, not gap-filled) was the sum of the individual land-cover type fluxes (F_modeled,mix and the merged polygonal pond F_pond) each multiplied with their weighted footprint fraction (a_mix and a_pond), with a_{open water}=a_pond, $a_{\text{mix}}=a_{\text{sum}}-a_{\text{pond}}$, and a_sum being the sum over all land-cover classes:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{F}_{\text{obs,mix}}& =a_{\text{pond}}\cdot F_{\text{pond}}+a_{\text{mix}}\cdot F_{\text{modeled,mix}}\\ \iff F_{\text{pond}}& ={\displaystyle \frac{F_{\text{obs,mix}}-a_{\text{mix}}\cdot F_{\text{modeled,mix}}}{a_{\text{pond}}}}.\end{array}\end{array}$$

To improve data quality, we excluded 30 min fluxes of F_pond when a_pond<50 %. Then, we used the median of F_pond for further calculations, and we assumed that all thermokarst ponds in the EC footprint emitted the same amount of CO₂.

As mentioned above, the observed CO₂ flux from the wind directions of north, west, and south (F_obs,mix) was influenced by polygonal ponds to a small degree. Since our aim was to assess the impact of thermokarst ponds (both polygonal ponds and merged polygonal ponds) on NEE, we needed to eliminate the influence of polygonal ponds from our NEE estimate. To extract uncontaminated CO₂ flux data from the semi-terrestrial tundra (F_{modeled,tundra}), we subtracted the previously estimated pond CO₂ flux F_pond from the observed CO₂ flux F_obs,mix:

$$\begin{array}{}\text{(5)}& F_{\text{modeled,tundra}}={\displaystyle \frac{F_{\text{obs,mix}}-a_{\text{pond}}\cdot F_{\text{pond}}}{a_{\text{mix}}}}.\end{array}$$

We then used this estimated CO₂ flux from the semi-terrestrial tundra F_{modeled,tundra} as the regressand variable for the bulk-NEE model to obtain a gap-filled dataset regarding CO₂ flux from the semi-terrestrial tundra. This gap-filling modeling of CO₂–C flux had an RMSE of 0.31 g m⁻² d⁻¹.

To evaluate the impact of thermokarst ponds on landscape CO₂ flux, we estimated a polygonal tundra landscape–CO₂ flux from the Late Holocene river terrace of Samoylov Island (F_landscape) by combining thermokarst ponds and semi-terrestrial tundra linearly:

$$F_{\text{landscape}}=A_{\text{pond}}\cdot F_{\text{pond}}+A_{\text{tundra}}\cdot F_{\text{modeled,tundra}},$$

where F_pond describes the CO₂ emissions from the open-water areas of thermokarst ponds (Eq. 4), F_{modeled,tundra} describes the modeled CO₂ flux from the semi-terrestrial tundra (Eq. 5), A_pond=0.07 is the fraction of the river terrace area of Samoylov Island that is covered by thermokarst ponds (from the land-cover classification; see Sect. 2.4.1), and $A_{\text{tundra}}=\mathrm{1}-\mathrm{0.07}$ is the fraction of the entire river terrace area that consists of other land-cover types. We did not account for larger or deeper lakes in this upscaling approach as we expected different greenhouse gas emission dynamics from these lakes and there were no lakes in our footprint and therefore within our observation range. Thus, we scaled the above numbers to $A_{\text{tundra}}+A_{\text{pond}}=\mathrm{1}$, which results in A_pond=0.076 and A_tundra=0.924.

2.4.5 CH₄ flux partitioning

The data show that the CH₄ emissions from the heterogeneous landscape around the tower were less spatially uniform than the CO₂ emissions. Therefore, we could not use a gap-filling model for CH₄ that was similar to the bulk model we used for CO₂, so we investigated CH₄ emissions in a different way. Based on preliminary results from our analysis and the aerial image of the study site, we focused on four wind sectors instead of extracting the fluxes from the land-cover types:

• Tundra. At least half of the footprint consisted of dry tundra, and the wind direction was larger than 170^∘.

• Shore 50^∘ (denoted shore${}_{\mathrm{50}{}^{\circ }}$). Less than 40 % of the footprint consisted of dry tundra, and water comprised at least 30 % of the footprint. The wind direction range was 30^∘ < WD < 65^∘.

• Pond. At least half of the footprint consisted of open water, and the wind direction range was 65^∘ < WD < 110^∘.

• Shore 120^∘ (denoted shore${}_{\mathrm{120}{}^{\circ }}$). Less than 40 % of the footprint consisted of dry tundra, and water comprised at least 30 % of the footprint. The wind direction range was 110^∘ < WD < 130^∘.

2.4.6 CH₄ permutation test

To evaluate whether the differences in flux medians between the four wind sectors were significant, we applied a permutation test (Edgington and Onghena, 2007). In this test, we randomly assigned each 30 min flux to one of two groups and calculated both groups' medians and the differences between the group's medians. We conducted six tests in total, using all possible combinations of pairs with the four wind sectors. After repeating this step 10 000 times, we plotted the resulting differences in medians in a histogram and performed a one-sample t test to evaluate whether the observed difference in medians differed significantly (p<0.01) from the randomly generated differences.

3 Results

3.1 Meteorological conditions

During the measurement period between 11 July and 10 September 2019, half-hourly air temperatures range from −0.5 to 27.6 ^∘C with a mean temperature of 8.7 ^∘C (Fig. A2a). The maximum wind speed measured at the EC tower at a height of 2.25 m is 8.9 m s⁻¹ (Fig. A2b). PAR reaches values of up to 1419 µmol m⁻² s⁻¹ with decreasing maximum values during the measurement period (Fig. A2c). Throughout the measurement period, there are 28 cloudy days, determined by identifying days with low PAR values (maximum values below ∼ 500 µmol m⁻² s⁻¹).

Figure 2 Polar plot of observed 30 min CO₂–C fluxes with respect to the wind direction. Negative values (inside of the dashed black line) represent CO₂ uptake, while positive values (outside of the dashed black line) represent CO₂ emissions. The values −4, −2, 0, and 2 indicate the magnitude of the CO₂–C flux in g m⁻² d⁻¹. The color of each point on the plot represents the percentage the point comprises of the total open-water weighted footprint fraction in each 30 min flux. The red boxes indicate the mean CO₂ flux of 5^∘ wind direction intervals during the 2-month observation period (red lines indicate the first standard deviation).

3.2 CO₂ fluxes

When inspecting the relation between observed CO₂ fluxes and wind direction (Fig. 2), we find that CO₂ fluxes exhibit high temporal variability between positive and negative CO₂ fluxes from most wind directions. In the wind sector between 60–120^∘, the flux source area is dominated by the merged polygonal pond. The CO₂–C fluxes from this pond sector show smaller absolute variability (${\mathrm{0.09}}_{-\mathrm{0.33}}^{\mathrm{0.38}}$ g m⁻² d⁻¹, median${}_{\mathrm{5}\mathrm{th}\mathrm{percentile}}^{\mathrm{95}\mathrm{th}\mathrm{percentile}}$) than the fluxes from all other wind directions ($-\mathrm{0.08}{}_{-\mathrm{1.56}}^{\mathrm{0.87}}$ g m⁻² d⁻¹, median${}_{\mathrm{5}\mathrm{th}\mathrm{percentile}}^{\mathrm{95}\mathrm{th}\mathrm{percentile}}$). Additionally, we observe a lower respiration rate from the merged polygonal pond than from the semi-terrestrial tundra. Figure 3 shows the observed nighttime CO₂ fluxes plotted against the respective weighted footprint fraction of open water. We define nighttime as having PAR < 20 µmol m⁻² s⁻¹; we expect that there would only be respiration and no photosynthesis during the nighttime. We find that the fluxes decrease as the pond area contribution increases. Thus, the strength of CO₂ respiration shows a dependence on the contribution of open water. We also find that low air temperatures are mostly associated with low respiration rates.

Another aspect of CO₂ flux variability stems from the diurnal cycle. We compare the diurnal cycle of the CO₂ fluxes from the merged polygonal pond (estimated in accordance with Eq. 4) and the semi-terrestrial tundra (Eq. 5, Fig. 4). The results show a less pronounced diurnal CO₂ cycle from the direction of the merged polygonal pond (blue) compared to the diurnal CO₂ cycle from the semi-terrestrial tundra (green). We combine all data from the merged polygonal pond (F_pond in Eq. 4), which results in a CO₂–C flux of ${\mathrm{0.13}}_{\mathrm{0.00}}^{\mathrm{0.24}}$ g m⁻² d⁻¹ (median${}_{\mathrm{25}\mathrm{th}\mathrm{percentile}}^{\mathrm{75}\mathrm{th}\mathrm{percentile}}$).

Figure 3 Scatterplot of observed CO₂ fluxes against the weighted footprint fraction of open water during each 30 min flux. The air temperature is represented through color. Only fluxes observed in the nighttime (PAR < 20 µmol m⁻² s⁻¹) are shown.

Figure 4 Diurnal cycle of modeled CO₂–C based on observations flux from the merged polygonal pond (blue, Eq. 4) and the semi-terrestrial tundra (green, Eq. 5) as violin plots for each half-hour flux. Blue and green crosses mark the mean CO₂–C flux during each half-hour flux. A violin plot shows the distribution of measurements along the y axis – the width of the curves indicates how frequently a certain y value occurred.

3.3 CH₄ fluxes

We plot the observed CH₄ fluxes against wind direction (Fig. 5). The results show that the CH₄ emissions peak at $\sim \mathrm{120}{}^{\circ }$, where fluxes from one shoreline of the merged polygonal pond contribute to the observed flux (Fig. 1d, from now on shore${}_{\mathrm{120}{}^{\circ }}$). We do not observe a similar peak of CH₄ emissions in the direction of the second shoreline towards $\sim \mathrm{50}{}^{\circ }$ (shore${}_{\mathrm{50}{}^{\circ }}$). These peaks did not correlate with a specifically large contribution of one of the land-cover classes to the footprint.

To further investigate the peak at shore${}_{\mathrm{120}{}^{\circ }}$, we compare the CH₄ emissions from the different wind sectors (shore${}_{\mathrm{120}{}^{\circ }}$, shore${}_{\mathrm{50}{}^{\circ }}$, pond, and tundra; Sect. 2.4.5). We find the following fluxes from the wind sectors: ${\mathrm{19.18}}_{\mathrm{14.26}}^{\mathrm{24.47}}$ (shore${}_{\mathrm{120}{}^{\circ }}$), ${\mathrm{12.96}}_{\mathrm{10.34}}^{\mathrm{15.11}}$ (shore${}_{\mathrm{50}{}^{\circ }}$), ${\mathrm{13.90}}_{\mathrm{11.02}}^{\mathrm{18.46}}$ (pond), and ${\mathrm{12.55}}_{\mathrm{9.65}}^{\mathrm{16.07}}$ mg m⁻² d⁻¹ (tundra, median${}_{\mathrm{25}\mathrm{th}\mathrm{percentile}}^{\mathrm{75}\mathrm{th}\mathrm{percentile}}$). Fluxes from shore${}_{\mathrm{120}{}^{\circ }}$ have a higher median than fluxes from the other three wind sectors (Fig. 6).

We investigated the impact of wind speed and air temperature on the CH₄ fluxes by excluding flux intervals with high wind speed (greater than 5 m s⁻¹) and high air temperature (warmer than 12 ^∘C). The randomization test (Sect. 2.4.6) provided evidence of a significant difference between CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$ and the other three wind sector classes at low wind speeds (top row in Fig. A4) and no significant difference between the CH₄ emissions from the classes pond–tundra and shore${}_{\mathrm{50}{}^{\circ }}$–tundra. The difference between the classes pond and shore${}_{\mathrm{50}{}^{\circ }}$ is significant; however, it is much smaller than the previously described differences (see central graph in Fig. A4). Note that the CH₄ emissions from pond and tundra have a similar magnitude under moderate wind-speed conditions. The results are very similar for moderate temperatures: we find evidence of a significant difference between the CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$ and the CH₄ emissions from the other three wind sector classes (top row in Fig. A5). The differences in medians between the pond and shore${}_{\mathrm{50}{}^{\circ }}$ and between the pond and tundra are significant. However, this difference is much smaller (second row in Fig. A5). In summary, neither high wind speed nor high temperatures act as a driver for the high CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$. In contrast, the peak at 180–190^∘ can be explained reasonably well using air temperature and friction velocity in multiple linear regression (R²=0.44). Using the same predictors results in an R² of 0.20 for the peak at shore${}_{\mathrm{120}{}^{\circ }}$.

The ratio of CO₂–C to CH₄–C emissions at night (PAR < 20 µmol m⁻² s⁻¹) has a value of CH₄ $/$ CO${}_{\mathrm{2}}={\mathrm{0.060}}_{\mathrm{0.049}}^{\mathrm{0.076}}$ for fluxes with an open-water weighted footprint fraction of more than 60 %, whereas the ratio amounts to CH₄ $/$ CO${}_{\mathrm{2}}={\mathrm{0.020}}_{\mathrm{0.015}}^{\mathrm{0.024}}$ (median${}_{\mathrm{25}\mathrm{th}\mathrm{percentile}}^{\mathrm{75}\mathrm{th}\mathrm{percentile}}$) for fluxes with an open-water weighted footprint fraction of less than 20 %.

Figure 5 Polar plot of 30 min observed CH₄–C flux with respect to the wind direction at the EC tower. Positive values outside the dashed black line represent CH₄ emissions, while values inside the line represent CH₄ uptake during one half-hour period. The values 0, 20, 40, and 60 indicate the magnitude of the CH₄–C flux in mg m⁻² d⁻¹. The color of each point on the plot represents the percentage the point comprises of the total open-water weighted footprint fraction in each 30 min flux. The red boxes indicate the mean CH₄ flux of 5^∘ wind direction intervals during the 2-month observation period (red lines indicate the first standard deviation).

Figure 6 Violin plots of observed CH₄ emissions at the EC tower separated into four different wind sector classes. A violin plot shows the distribution of measurements along the y axis – the width of the shapes indicates how frequently a certain y value occurred. Medians of CH₄ emission distributions are shown as red lines, and 75th and 25th percentile are shown as black lines. On the right, the wind sectors with the eddy covariance tower in the center (black cross) are shown.

3.4 Upscaled CO₂ flux

We use the estimated open-water CO₂ flux from the merged polygonal pond and the modeled CO₂ flux from the semi-terrestrial tundra to linearly upscale the CO₂ flux for the polygonal tundra of Samoylov Island (excluding larger lakes, the method described in Sect. 2.4.4). As we have not obtained estimates for the CH₄ fluxes from tundra and pond land-cover types, we only upscale CO₂.

We estimate that when one includes the CO₂ flux from thermokarst ponds, the river terrace landscape's CO₂ uptake is ∼ 11 % lower than the uptake of semi-terrestrial tundra without ponds. The modeled CO₂–C flux from the semi-terrestrial tundra (without consideration of thermokarst pond fluxes) accumulated to −16.29 ± 0.43 g m⁻² during the observation period (60.5 d). If separated into months, the modeled CO₂–C flux from the semi-terrestrial tundra amounts to −15.01 ± 0.26, −3.56 ± 0.33, and +2.35 ± 0.11 g m⁻² in July (19.8 d), August (31 d), and September (9.7 d), respectively. When one includes the CO₂ flux from the merged polygonal pond to represent all thermokarst ponds on Samoylov Island, the resulting estimate of the landscape CO₂ flux amounts to −14.47 ± 0.40 g m⁻² (60.5 d), with monthly fluxes of −13.75 ± 0.24, −2.99 ± 0.31, and +2.27 ± 0.10 g m⁻² in July (19.8 d), August (31 d), and September (9.7 d), respectively. Thus, the results show that thermokarst ponds have the largest impact on the landscape's CO₂ flux in August. In September, accounting for thermokarst ponds leads to a 3.5 % lower estimate of landscape CO₂ emissions.

4 Discussion

4.1 CO₂ flux

Only a limited number of EC CO₂ flux studies from permafrost-affected ponds and lakes are available (studies with “EC” in Table 1). Estimates of open-water EC CO₂–C flux range from 0.059 (Jammet et al., 2017) to 0.11 (Eugster et al., 2003) to 0.22 g m⁻² d⁻¹ (Jonsson et al., 2008). Our estimate of ${\mathrm{0.12}}_{\mathrm{0.0014}}^{\mathrm{0.24}}$ g m⁻² d⁻¹ is, therefore, well within the range of open-water CO₂–C fluxes observed with the EC method. Other studies using different methods report a wider range of open-water CO₂ fluxes in Arctic regions. These fluxes range from a minor CO₂–C uptake (−0.14 g m⁻² d⁻¹; Bouchard et al., 2015) to substantial emissions of CO₂–C (up to 2.2 g m⁻² d⁻¹; Abnizova et al., 2012). A modeling study involving multiple lakes in northeast European Russia found that they produce almost zero emissions (0.028 g m⁻² d⁻¹; Treat et al., 2018).

Abnizova et al. (2012)Jammet et al. (2017)Jonsson et al. (2008)Eugster et al. (2003)Jansen et al. (2019)Bouchard et al. (2015)Sepulveda-Jauregui et al. (2015)Treat et al. (2018)Sieczko et al. (2020)Ducharme-Riel et al. (2015)Repo et al. (2007)Lundin et al. (2013)Kling et al. (1992)

Table 1 Daily mean water–atmosphere CO₂ and CH₄ fluxes from different study sites. TBL is the abbreviation for thin boundary layer model, EC for eddy covariance, CH for chamber measurement, MOD for modelled fluxes, STO for storage fluxes, and NEW for the method used in this study. All fluxes are given ± standard deviation, except fluxes from this study are given as median${}_{\mathrm{25}\mathrm{th}\mathrm{percentile}}^{\mathrm{75}\mathrm{th}\mathrm{percentile}}$.

Strikingly, our estimates of open-water CO₂ emissions are approximately 12–18 times smaller than those that have been previously reported for open-water CO₂ emissions at the same study site (Abnizova et al., 2012). One reason for the divergent results might be the different methods used. In Abnizova et al. (2012), the thin boundary layer (TBL) model, following Liss and Slater (1974), was applied to estimate CO₂ emissions from CO₂ concentrations. However, one other study found good agreement between the EC method and the TBL one (Eugster et al., 2003). In addition, in contrast to the larger merged polygonal pond we focus on, Abnizova et al. (2012) measured two polygonal ponds (they took 46 water samples in August and September 2008). These two ponds might have had exceptionally high CO₂ concentrations and might not be representative of polygonal ponds in our study area. If the polygonal ponds in the footprint of our EC measurements emitted CO₂ in the quantities suggested by Abnizova et al. (2012), we would expect to see their signal more clearly in our measurements.

Our approach of combining a footprint model with land-cover classification to extract fluxes from different land-cover classes allows us to determine the thermokarst pond CO₂ flux. We report an uncertainty range with respect to the thermokarst pond CO₂ flux; however, identifying the full uncertainty in this flux is not possible using this approach due to the footprint analysis' unknown degree of uncertainty. Still, the results with respect to the thermokarst pond CO₂ flux are plausible and on the expected order of magnitude for two reasons. First, reduced diurnal variability is observed when the merged polygonal pond influences the flux signal (Fig. 4). This reduction indicates that the respiration rate from the merged polygonal pond is lower than the respiration rate from the semi-terrestrial tundra, where ample oxygen is available in the upper soil layer. Additionally, since the thermokarst ponds have a lower vegetation density than the tundra, there is less photosynthesis. Second, when focusing on nighttime fluxes, when only respiration occurs (i.e., no carbon is taken up), there is a decrease in CO₂ emissions with an increasing weighted footprint fraction of open water (Fig. 3); this also indicates that there is reduced decomposition in the merged polygonal pond. Overall, based on the data, the finding that thermokarst ponds have lower CO₂ emissions than the semi-terrestrial tundra is reasonable.

4.2 CH₄ flux

We observe large differences in CH₄ emissions from the four wind sectors. CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$ are significantly higher than from shore${}_{\mathrm{50}{}^{\circ }}$, pond, and tundra (Sect. 3.3). Notably, we tested the dependence of these higher fluxes on wind speed and air temperature. We expect high wind speeds to enhance turbulent mixing of the water column and diffusive CH₄ outgassing at the water–atmosphere interface. High wind speeds are also associated with pressure pumping, which potentially fosters the ebullition of CH₄. On the other hand, peak temperatures can lead to peak CH₄ production and emissions due to enhanced biological activity. However, the high emissions from shore${}_{\mathrm{120}{}^{\circ }}$ do not coincide with either of two key meteorological conditions, high wind speeds and high temperatures, which would especially favor high emissions. Thus, the difference in methane flux dynamics between shore${}_{\mathrm{120}{}^{\circ }}$ and shore${}_{\mathrm{50}{}^{\circ }}$ is astounding since the shorelines share many other characteristics.

Both shorelines extend radially (in a fairly straight line) from the EC tower (Fig. 1), thus contributing similarly to the EC flux. The underwater topography does not vary significantly between the two shorelines. Meters away from the shore, both shorelines have a water depth of a few centimeters and a few decimeters (see data from Boike et al., 2015a). As previously described in Sect. 2.1, both shorelines are dominated by Carex aquatilis, and from visual inspection, we could not identify differences in shoot density. We, therefore, assume that the characteristics of the emergent vegetation do not play a major role in explaining the differences between the CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$ and shore${}_{\mathrm{50}{}^{\circ }}$. We also examine the evolution of the shorelines at the merged polygonal pond to check whether erosion along the shoreline could cause the high CH₄ emissions. We compare an image from 1965 (Earth Resources Observation and Science (EROS) Center, 2018) with the current (2019) shoreline, yet we cannot identify signs of recent erosion. Furthermore, high-resolution aerial images of this pond from 2008 (Boike et al., 2015b, resolution > 0.33 m) and 2015 (Boike et al., 2015c, resolution > 0.33 m) show no signs of erosion. We therefore assume that past erosion is unlikely is unlikely to have been a factor that caused the high levels of CH₄ emissions we observed in 2019.

Local ebullition of the merged polygonal pond could lead to high CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$. We applied the method proposed by Iwata et al. (2018) to check for signs of ebullition events. This method uses the 20 Hz raw CH₄ concentration data to detect short-term peaks in CH₄ that originate from ebullition events. However, we cannot detect ebullition events in the 20 Hz raw data.

In summary, meteorological conditions (wind speed and temperature), characteristics of emergent vegetation, coastal erosion, and intense ebullition events are unlikely to be the main driving factors of the increased CH₄ emissions we observed. Another possible driver of higher CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$ is a small but steady seep ebullition hot spot close to this shoreline (such as ebullition class kotenok in Walter et al., 2006). Seep ebullition hot spots have been reported to occur heterogeneously in clusters in Alaskan lakes (Walter Anthony and Anthony, 2013). Unfortunately, seep ebullition has not previously been reported in water bodies in our study area, so we did not include measurements targeting this process in our measurement campaign. In future studies, visual inspection of trapped CH₄ bubbles in the ice column during wintertime, as proposed by Vonk et al. (2015), could reveal more information about the cause of the higher CH₄ emissions from shore${}_{\mathrm{120}{}^{\circ }}$, as could funnel or chamber measurements with high spatial coverage.

The results show that the merged polygonal pond emits a similar magnitude of CH₄ to the polygonal tundra surface under similar meteorological conditions and when excluding the high emissions from shore${}_{\mathrm{120}{}^{\circ }}$. However, substrate availability and temperature dynamics differ substantially. Additionally, in dense soils, methane diffuses slowly enough through soil layers containing oxygen for the methane to be oxidized before reaching the surface. In contrast, methane emitted in ponds can reach the surface quickly through ebullition or plant-mediated transport in addition to diffusion. Therefore, we expect to see larger differences between the CH₄ emissions from the merged polygonal pond and the polygonal tundra, more akin to the differences that have been detected in a subarctic lake and fen by Jammet et al. (2017). However, we see no significant difference in the CH₄ emissions from the open-water areas of the merged polygonal pond and the polygonal tundra surface (Figs. 6 and A4).

Since many other thermokarst ponds in our study area are smaller than the merged polygonal pond (making them unsuitable to study using the EC method) and since smaller ponds tend to be greater emitters of methane (Holgerson and Raymond, 2016; Wik et al., 2016), our measurements might provide a lower limit of overall thermokarst pond CH₄ emissions.

We estimate a CH₄–C flux of ${\mathrm{13.38}}_{\mathrm{10.55}}^{\mathrm{15.92}}$ mg m⁻² d⁻¹ (${\mathrm{median}}_{\mathrm{25}\mathrm{th}\mathrm{percentile}}^{\mathrm{75}\mathrm{th}\mathrm{percentile}}$) from the merged polygonal pond and ${\mathrm{12.96}}_{\mathrm{10.34}}^{\mathrm{15.11}}$–${\mathrm{19.18}}_{\mathrm{14.26}}^{\mathrm{24.47}}$ mg m⁻² d⁻¹ from the shores of this pond. This is higher than the fluxes measured by Jammet et al. (2017) from a subarctic lake (Table 1). The authors report a mean annual CH₄–C flux of 13.42 ± 1.64 mg m⁻² d⁻¹ and a mean ice-free-season CH₄–C flux of 7.58 ± 0.69 mg m⁻² d⁻¹. A study focusing on 32 non-Yedoma thermokarst lakes in Alaska found CH₄–C emissions similar to our results (16.80 ± 8.61 mg m⁻² d⁻¹; Sepulveda-Jauregui et al., 2015). Also, a synthesis of 149 thermokarst water bodies north of $\sim \mathrm{50}{}^{\circ }$ N reports CH₄–C emissions on the same order of magnitude (27.57 ± 14.77 mg m⁻² d⁻¹; Wik et al., 2016). However, other recent studies have reported considerably lower CH₄–C emissions of 2.95 ± 0.75 mg m⁻² d⁻¹ in northern Sweden (Sieczko et al., 2020), and, in contrast, a study found CH₄–C emissions of up to 6432 mg m⁻² d⁻¹ in northeast Canada (Bouchard et al., 2015). The wide range of water-body methane emissions militates in favor of caution when generalizing our results, even for Samoylov Island, especially since the emissions within the merged polygonal pond have been shown to be heterogeneous. Instead, after finding a hot spot in CH₄ emissions at the pond shore, we would like to highlight that the gathering of additional measurements – for example employing funnel traps or counting bubbles in ice – will help to better constrain thermokarst pond CH₄ dynamics in their full complexity. Nevertheless, our measurements provide a robust lower limit of thermokarst pond CH₄ emissions.

4.3 Upscaling the CO₂ flux

We upscale the CO₂ emissions for the river terrace on Samoylov, an area for which we have access to high-resolution land-cover classification. We find that we overestimate the carbon dioxide uptake of the polygonal tundra by 11 % when we do not account for the thermokarst ponds' CO₂ emissions. A similar approach by Abnizova et al. (2012) found a potential increase of 35 %–62 % in the estimate of CO₂ emissions from the Lena River delta when including small ponds and lakes in the landscape CO₂ emission calculation. If we were to follow the upscaling approach by Abnizova et al. (2012) and consider overgrown water as part of the thermokarst ponds, the estimate of the landscape CO₂ uptake would decrease by 19 %. Kuhn et al. (2018) also found water bodies in Arctic regions to be an important source of carbon, which could outbalance the carbon dioxide uptake of the semi-terrestrial tundra in a future climate. In summary, our results demonstrate that open-water CO₂ emissions can substantially influence the summer carbon balance of the polygonal tundra. With respect to the night time emissions, we find that per gram CO₂–C thermokarst ponds emit 0.06 g CH₄–C whereas the semi-terrestrial tundra only emits 0.02 g CH₄–C. This finding underlines again that, especially when considering thermokarst ponds, CH₄ emissions are of significant interest. Even though mean CH₄ emissions from the semi-terrestrial tundra and open water are of similar magnitude, we expect that the impact of thermokarst ponds on the carbon balance would be even greater when accounting for CH₄ due to locally high emissions.

Our results suggest that future studies that aim to capture a representative landscape flux should pay extra attention to the water bodies in their footprint. The CO₂ flux from thermokarst ponds has the opposite sign (CO₂ emissions) to the semi-terrestrial tundra (CO₂ uptake) during the observation period. Consequently, thermokarst ponds should cover about as much area in the measurement as they do in the landscape area of interest. In this way, the chances of capturing CH₄ hot spots, which can be investigated more closely, are also greater.

5 Conclusions

We find that thermokarst ponds are a carbon source. At the same time, the surrounding semi-terrestrial tundra in our study area acts as a carbon sink during the summer period (July–September), which is in agreement with prior studies (Abnizova et al., 2012; Jammet et al., 2017), despite us observing much lower open-water CO₂ fluxes compared to previous work at the same study site (Abnizova et al., 2012). Using our approach to disentangle the EC fluxes from different land-cover classes, we posit that during the measurement period, we would overestimate the carbon dioxide uptake of the polygonal tundra by 11 % if thermokarst ponds were not accounted for. We expect lakes to have a similar effect on the carbon budget, though a smaller one, since lakes (a) cover a similar amount of surface area as the thermokarst ponds in our study site (Abnizova et al., 2012; Muster et al., 2012) and (b) are weaker emitters of greenhouse gases than ponds (Holgerson and Raymond, 2016; Wik et al., 2016).

In contrast to CO₂ emissions, which are spatially more homogeneous, small-scale heterogeneity in CH₄ emissions makes it difficult to find drivers of CH₄ emissions. We cannot pinpoint the drivers behind the high emissions along parts of the coastline, which we surmise were potentially caused by seep ebullition. Thus, we cannot estimate the impact of this heterogeneity on the landscape scale and, therefore, refrain from upscaling CH₄ emissions. Additionally, the open-water fluxes presented in this paper originate from a single merged polygonal pond since the other polygonal ponds surrounding the EC tower are too small to extract their fluxes using the footprint method applied here. Thus, we do not account for the spatial variability in CH₄ emissions between thermokarst ponds, which can be substantial (Rehder et al., 2021; Wik et al., 2016). However, we note that open-water fluxes were of a similar magnitude to the polygonal tundra fluxes. Consequently, the main impact that thermokarst ponds have on the landscape CH₄ budget might occur through plant-mediated transport and local ebullition.

While being ill-suited for the study of smaller ponds, we underline that the EC method is appropriate for observing greenhouse gas fluxes from thermokarst ponds as small as 0.024 km². The EC method has a higher temporal resolution than the TBL method. It does not disturb exchange processes like the chamber flux method, which eliminates the wind at the water surface. Especially when combining an EC footprint with land-cover classification, one can distinguish between the contribution of different land-cover classes effectively and also study the fluxes from thermokarst ponds.

We conclude that thermokarst ponds contribute significantly to the landscape carbon budget. Changes in Arctic hydrology and the concomitant changes in the water-body distribution in permafrost landscapes may cause these landscapes to change from being overall carbon sinks to overall carbon sources.

Appendix A

Figure A1 Picture of the eddy covariance tower with the merged polygonal pond in the background. Picture taken on 11 July 2019 by Zoé Rehder.

Figure A2 Timeline of observed meteorological conditions during the observation period (2019) with air temperature at 2 m height (a), wind speed at 3 m height (b), and photosynthetically active radiation (PAR) (c). Mean values and standard deviation of observations during the past 16 years are plotted as black lines and gray areas.

Figure A3 Time series of 30 min observed CO₂–C flux intervals (a) and CH₄–C flux (b) with a quality flag of 0 or 1. The blue color represents fluxes originating from the wind direction of the merged polygonal pond (30–150^∘ wind direction, mostly mixed signals from semi-terrestrial tundra and the surface of the merged polygonal pond), and the green color represents fluxes originating from all other wind directions.

Figure A4 Histogram of permutation tests between the medians of CH₄ emissions from different wind direction classes in Fig. 6. All medians from flux observations during moderate wind-speed conditions. The observed differences in medians between the different wind direction classes are shown in red vertical bars in each plot.

Figure A5 Histogram of permutation tests between the medians of CH₄ emissions from different wind direction classes in Fig. 6. All medians from flux observations during moderate air temperature conditions. The observed differences in medians between the different wind direction classes are shown in red vertical bars in each plot.