Small waterbodies reduce the carbon sink of a polygonal tundra landscape

Arctic permafrost landscapes have functioned as a global carbon sink for millennia. These landscapes are very heterogeneous, and the omnipresent waterbodies are a carbon source within them. Yet, only a few studies focus on the impact of these waterbodies on the landscape carbon budget. We compare carbon dioxide and methane fluxes from small waterbodies to fluxes from the surrounding tundra using eddy covariance measurements from a tower located between a large pond and semi-terrestrial vegetated tundra. 5 When taking the open-water areas of small waterbodies into account, the carbon dioxide sink strength of the landscape was reduced by 11%. While open-water methane emissions were similar to the tundra emissions, some parts of the studied pond’s shoreline exhibited much higher emissions, underlining the high spatial variability of methane emissions. We conclude that gas fluxes from small waterbodies can contribute significantly to the carbon budget of arctic tundra landscapes. Consequently, changes in arctic hydrology and the concomitant changes in the waterbody distribution may substantially impact the overall 10 carbon budget of the Arctic.


Introduction
Waterbodies make up a significant part of the arctic lowlands with an areal coverage of about 17%, (Muster et al., 2017) and considerably decrease the landscape carbon sink (Kuhn et al., 2018).The thaw of permafrost in the warming Arctic is going to change the distribution of waterbodies (Andresen and Lougheed, 2015;Bring et al., 2016) and thus also their contribution to the landscape-carbon budget (Kuhn et al., 2018).However, data on greenhouse-gas emission from arctic waterbodies are still sparse in space and time, especially with a high temporal resolution and in non-Yedoma regions (Vonk et al., 2015).
Our study site, in the Lena River Delta, Siberia, is located on an island mostly covered by non-Yedoma polygonal tundra (Fig. 1).This landscape features many ponds (defined here by an area < 8 • 10 4 m 2 , Ramsar Convention Secretariat (2016); Rehder et al. (2021)), as opposed to larger lakes, and in our area of interest, ponds cover about as much area as lakes (Abnizova et al., 2012;Muster et al., 2012).Ponds emit more greenhouse gases per area than lakes (Holgerson and Raymond, 2016 et al., 2016), thus, in our study area, they have a higher potential than lakes to reduce the carbon sink of the surrounding tundra (McGuire et al., 2012;Jammet et al., 2017;Kuhn et al., 2018).To estimate the impact of ponds, we compare landscape carbon dioxide (CO 2 ) and methane (CH 4 ) fluxes from the open-water area of ponds to the more commonly reported fluxes from the semi-terrestrial tundra (defined here as wet and dry tundra as well as overgrown water).
Different geophysical and biochemical processes drive pond emissions of CO 2 and of CH 4 .Aquatic CO 2 production is dominated by the microbial decomposition of dissolved organic carbon, which is introduced laterally into the aquatic system through rain and melt water (Neff and Asner, 2001).When supersaturated with dissolved CO 2 , ponds emit CO 2 to the atmosphere through diffusion.While photosynthetic CO 2 uptake has been observed in some clear arctic waterbodies (Squires and Lesack, 2003), most arctic waterbodies are net CO 2 sources (Kuhn et al., 2018).Estimates range from close to zero (0.028 g m 2 d −1 by Treat et al. (2018), or 0.059 g m 2 d −1 by Jammet et al. (2017)) to substantial CO 2 -C emissions (1.4 -2.2 g m 2 d −1 by Abnizova et al. ( 2012)).
CH 4 emissions have been found to vary even more (by up to five orders of magnitude within just one site: 0.5 -6432 mg m 2 d −1 , Bouchard et al. (2015)).In contrast to CO 2 , most CH 4 originates in the sub-aquatic soils.It is emitted from waterbodies not only through diffusion but also through ebullition (sudden release of bubbles) and plant-mediated transport, often leading to high spatial variability between waterbodies and within one waterbody (Sepulveda-Jauregui et al., 2015;Jansen et al., 2019).Local seep-ebullition events can cause high spatial variance of CH 4 emissions within one waterbody (Walter et al., 2006).Varying coverage of vascular plants in the shallow parts of a waterbody can also increase CH 4 variability through plant-mediated transport (Knoblauch et al., 2015;Andresen et al., 2017).
To study both spatial and temporal patterns, we analyze land-atmosphere CO 2 and CH 4 flux observations from an eddy covariance (EC) tower located on Samoylov Island, Lena River Delta, Russia.We set the EC tower up within the polygonal tundra next to a merged polygonal pond for two months in summer 2019.A merged polygonal pond is a larger pond which formed through the subsidence of several polygons.The polygonal structures are still clearly visible along the shore and under water, and these ponds tend to be shallow for their size (Rehder et al., 2021).Due to the tower's position, fluxes from the merged polygonal pond are the dominant source of the observed EC fluxes under easterly winds.The observed EC fluxes are dominated by vegetated polygonal tundra with only a low fraction from polygonal-center ponds from the other wind directions.We (1) compare the waterbody and tundra fluxes with a focus on temporal and spatial patterns, and we (2) investigate the influence of the merged polygonal pond on the landscape carbon balance.

Study Site
The study site Samoylov Island (72 • 22'N, 126 • 28'E) is located in the southern part of the Lena River Delta (Figure 1 Kartoziia (2019)).The island is surrounded by partially and regularly flooded branches of the Lena River and sandy floodplains, creating more spatial heterogeneity on a larger scale.We focus on a merged polygonal pond (Figure 1, d), which is located in the eastern part of the island.The merged polygonal pond in our study has a size of 0.024 km 2 with a maximum depth of 3.4 meters and a mean depth of 1.2 meters (Rehder et al., 2021;Boike et al., 2015a).On an aerial image, the polygonal structures are still clearly visible under the water surface (Boike et al., 2015c).
The vegetated shoreline of this merged polygonal pond is dominated by Carex aquatilis interspersed with Carex chordorrhiza, Potentilla palustris and Aulacomnium spp.. Some plants grow in the water of the pond close to the shore.The deeper parts of the pond are vegetation free.

Instruments
We measured gas fluxes using an eddy covariance (EC) tower between July 11 and September 10, 2019.The EC tower was located on the eastern part of Samoylov Island, directly at the western shore of the merged polygonal pond (Figure 1, d).The EC instruments were mounted on a tripod at the height of 2.25 meters.The tower was equipped with a closed-path CO 2 /H 2 O sensor (LI-7200, LI-COR Biosciences, USA), an open-path CH 4 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.
Additional meteorological data for Samoylov Island was provided by Boike et al. (2019).We also installed radiation-shielded temperature and humidity sensors at the EC tower (HMP 155,Vaisala,Finland) and used data from a photosynthetically active radiation (PAR) sensor mounted at a tower approximately 500 meters to the west (SKP 215, Skye Instruments, UK).

Data Processing
We perform the raw data processing and computation of half-hourly fluxes for open-path and closed-path fluxes using EddyPro 7.0.6 (LI-COR, 2019).Raw data screening includes spike detection and removal according to Vickers and Mahrt (1997) (1% maximum accepted spikes and a maximum of 3 consecutive outliners).Additionally, we apply statistical tests for raw data screening, including tests for amplitude resolution, skewness and kurtosis, discontinuities, angle of attack, and horizontal winds steadiness.All parameters of these tests are set to EddyPro default values.We rotate the wind-speed axis to a zero-mean vertical wind speed using the "double rotation"-method by Kaimal and Finnigan (1994).We apply linear de-trending following Gash and Culf (1996) to the raw data prior to the flux calculation.We compensate time lags by 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 are detected by covariance maximization (Fan et al., 1990) resulting in time lags of 0 -0.4 s for CO 2 and -0.5 -+0.5 s for CH 4 .For H 2 O, the time lag is humidity-dependent and is calculated for ten humidity classes.We compensate for air-density fluctuations due to thermal expansion/contraction and varying water-vapor concentrations following Webb et al. (1980).This correction is only applied to open-path data; for closed-path data, we perform a sample-by-sample conversion into mixing-ratios to account for air-density fluctuations (Ibrom et al., 2007b;Burba et al., 2012).Flux losses occur in the low-and high-frequency spectral range due to different filtering effects.In the low-frequency range, we compensate flux losses following Moncrieff et al. (2004)  and in the high-frequency range following Fratini et al. (2012).For applying the latter method, a spectral assessment file is created using the method by Ibrom et al. (2007a).The spectral assessment results in cut-off frequencies of 3.05 Hz and 1.67 Hz for CO 2 and CH 4 , respectively.For H 2 O, we find a humidity-dependent cut-off frequency between 1.25 Hz (RH 5 -45%) and 0.21 Hz (RH 75 -95%).We perform a quality check of each flux interval following the 0-1-2 system by Mauder and Foken (2004).In this quality check, flux intervals with the lowest quality receive the flag "2" and are excluded from further analysis.

Land Cover Classification
The land cover classification covers the late-Holocene river terrace of Samoylov Island (Siberia, Russia).It is based on highresolution near-infrared (NIR) orthomosaic aerial imagery obtained in the summer of 2008 (Boike et al., 2015b).We use a subset of the existing classification by Muster et al. (2012) as a training dataset to perform a semi-supervised land cover classification using the maximum likelihood algorithm in ArcMap Version 10.8 (ESRI Inc, USA).We then apply the ArcMap majority filter tool to the new classification.The land cover classification has a resolution of 0.17 m x 0.17 m, it is projected onto WGS 1984 UTM Zone 52N and the classes include open water, overgrown water, dry tundra, and wet tundra, as defined by Muster et al. (2012).

Footprint Model
The tower location and sensor height are crucial parameters in the deployment of an EC measurement tower.A lower measurement height results in a smaller footprint.The footprint describes the source area of the flux from the surrounding landscape.
With our sensors installed at the height of 2.25 m next to the merged polygonal pond, we expect to observe substantial flux signals from the adjacent waterbody as well as from the surrounding polygonal tundra.Each land cover type's contribution to the flux signal depends on the wind direction and turbulence in the atmospheric boundary layer.We implement the analytical footprint model by Kormann and Meixner (2001) in Matlab 2019b(MATLAB, 2019) and combine the footprint model with land cover classification data described in section 2.4.1 to estimate the contribution of each land cover type to each flux signal (hereinafter referred to as the weighted footprint fraction).The model accounts for the stratification of the atmospheric boundary layer and requires a height-independent crosswind distribution and horizontal homogeneity of the surface.The input data require the stationarity of atmospheric conditions during the flux interval of 30 minutes.We derive the vertical power-law profiles for the eddy diffusivity and the wind speed for each 30-minute flux interval depending on the atmospheric stratification (equation 6 in Kormann and Meixner (2001)).We use an analytical approach to find the closest Monin-Obukhov (M-O) similarity profile (equation 36 in Kormann and Meixner (2001)).Next, we calculate a two-dimensional probability density function of the source area for each flux interval (from equation 9 and 21 in Kormann and Meixner (2001)) and combine each probability density function with the land cover classification of the river terrace of Samoylov Island, with its four land cover types (see section 2.4.1).The footprint model's resolution is set to the land cover classification resolution of 0.17 m x 0.17 m.Hence, for each grid cell within the source area, we can estimate the probability of the fraction of flux originating from this grid cell for each 30-min interval.We also know the dominant land cover type in each grid cell from the land cover classification.We combine both information for each grid cell and calculate 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 obtain the contribution of each land cover type to each 30-minute 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 combine the contributions of the dry tundra, wet tundra, and overgrown water to a single land cover class for the semi-terrestrial tundra a tundra = a dry tundra + a wet tundra + a overgrown water .
We also take the sum of all 30-min two-dimensional probability density functions over the whole deployment time.This sum is referred to as the cumulative footprint.The cumulative footprint is shown as a gray shaded area in Figure 1, c and d.
The light gray area's outer boundary represents the 90%, and the light gray area's inner boundary is the 70% isoline of the cumulative footprint.This means that there is a probability of 10% that fluxes observed at the EC tower originate from areas outside of the light gray area.Medium gray represents 50-70%, medium-dark gray 30-50%, and dark gray indicates that there is a probability of less than 30% that the observed flux originates from within the marked area.

Bulk Model / Gap-Filling CO 2 Flux
We use the bulk net-ecosystem exchange of CO 2 (NEE) model by Runkle et al. (2013) to gap-fill and partition our NEE flux observations.This model was specifically developed for modeling NEE in arctic regions taking into account the polar day.We estimate all model parameters for running 5-day periods to capture changing plant physiology during the measurement period.
NEE is partitioned into two components (equation 3): total ecosystem respiration TER (µmol m −2 s −1 , equation 1) and gross primary production GPP (µmol m −2 s −1 , equation 2).Parameters of both components are fit simultaneously.TER is modeled as an exponential function of air temperature T air : where T ref = 15 • C and γ = 10 • C are constant, independent parameters.R base (µmol m −2 s −1 ) describes the basal respiration at the reference temperature T ref and Q 10 (dimensionless) the sensitivity of ecosystem respiration to air temperature changes.
GPP is modeled as an rectangular hyperbolic function of PAR (µmol m −2 s −1 ): where α (µmol µmol −1 ) is the initial canopy quantum use efficiency (slope of the fitted curve at PAR= 0) and P max (µmol m −2 s −1 ) the maximum canopy photosynthetic potential for PAR → ∞.
We sum both components to estimate the modeled NEE F CO2,mod : We implement the bulk model in Matlab 2019b (MATLAB, 2019) using the fit function with the fit-method of NonLinear-LeastSquares.We use the coeffvalues-function to estimate the four parameters and the confint-function to estimate their 95% confidence bounds.All partitioned fluxes are converted into CO 2 -C fluxes in the unit g m −2 d −1 prior to the data analysis.

Aquatic CO 2 Flux
In a heterogeneous landscape, fluxes observed using the EC method contain information from different land cover types.In this study, we extract fluxes primarily related to ponds and tundra from the mixed signals.We then combine these estimated fluxes to analyze the influence of ponds on a polygonal tundra landscape.
To estimate the CO 2 flux from the merged polygonal pond (F pond ), we first fit a bulk model to data from which we exclude fluxes from the merged polygonal pond (thus exclude fluxes >30 • & >150 • wind direction).This modeled CO 2 flux (F modeled,mix ) represents the vegetated tundra surrounding the EC tower, including small ponds to the north, west and south (with a weighted footprint fraction of open water of <30% in each flux signal) .In a second step, we make the assumption that individual contributions from different land cover types to the observed flux scale linearly with their contribution to the footprint.Thus, we can calculate the observed CO 2 flux (F obs,mix ) as the sum of the individual land cover type fluxes (F modeled,mix and F pond ) each multiplied with their weighted footprint fraction a tundra and a pond , respectively, where a open water = a pond , a tundra = a sum − a pond , and a sum being the sum over all land cover classes: To improve data quality, we exclude 30-min flux intervals of F pond when a pond < 50%.Then, we use the median of F pond for further calculations.
As mentioned above, the observed CO 2 flux from the north, west, and south (F obs,mix ) is still influenced by small ponds.To analyze in detail the CO 2 flux from the vegetated tundra (F modeled,veg ), we subtract the previously estimated pond-CO 2 flux F pond from the observed CO 2 flux F obs,mix : We then use this estimated CO 2 flux from the vegetated tundra F modeled,veg as the input variable for the bulk model to receive where F pond describes the CO 2 emission from the merged polygonal pond (equation 4), F tundra the modeled CO 2 flux from the vegetated tundra (equation 5), A pond = 0.07 the coverage of open pond water on the whole river terrace of Samoylov Island (from the land cover classification, section 2.4.1) and A tundra = 1 − 0.07 the coverage of other land cover types.We do no account for (larger, deeper) thermokarst lakes in this up-scaling approach, as we expect different greenhouse gas processes from these lakes and there are no lakes in our footprint.Thus, we scale the above numbers to A tundra + A pond = 1 which results in A pond = 0.076 and A tundra = 0.924.

Meteorological Conditions
During the measurement period between 11 July and 10 September 2019, half-hourly air temperatures ranged from -0.

CH 4 Fluxes
Figure 5 shows the observed CH 4 fluxes plotted against the wind direction.The CH 4 emissions peak at ∼ 120 • , where fluxes from a shoreline of the merged polygonal pond contribute to the observed flux (see Figure 1 d, hereinafter Shore 120 • ).We do not observe a similar peak of CH 4 emissions in the direction of the second shoreline towards ∼ 50 • (Shore 50 • ).These peaks did not correlate with any of the four weighted footprint fractions.
To further investigate the peak at Shore 120 • we separate the CH 4 emissions into four sectors depending on wind direction (Shore 120 • , Shore 50 • , pond and tundra, see Figure 6).We find the following fluxes from the wind direction sectors: 19.18 ).Fluxes from Shore 120 • have a higher median than fluxes from the other three wind sectors (see Figure 6).High wind speeds enhance turbulent mixing of the water column and diffusive CH 4 outgassing at the water-atmosphere interphase.High wind speeds are also associated with stronger pressure pumping potentially fostering ebullition, and thereby CH 4 emission.Additionally, peak temperatures can also lead to peak CH 4 production and emissions.
So, we investigated these confounding factors by excluding flux intervals with high wind speed (larger than 5 m s ) and high air temperature (larger than 12 • C).To evaluate whether the differences in medians between the four wind sectors are significant, we apply a permutation test (Edgington and Onghena (2007), Figure A2 and A3).In this test, we randomly assign each 30-min flux signal to one of two groups, calculate the median of both groups and their difference.After repeating this step 10000 times, we plot the resulting differences in medians in a histogram and perform a one-sample t-test to evaluate whether the observed difference in medians differs significantly (p<0.01)from the randomly generated differences.In the randomization test we find evidence for a significant difference between the CH 4 emission from Shore 120 • and the other three classes at low wind speeds (top row in Figure A2) and no significant difference between the CH 4 emission from the other three classes (Shore 50 • , pond and tundra).
The results are very similar for moderate temperatures: We find evidence for a significant difference between the CH 4 emission from Shore 120 • and the other three classes (top row in Figure A3).The differences in medians between the pond and Shore 50 • as well as between the pond and the tundra are significant.However, this difference is much smaller (second row in Figure A3).
In summary, we find no meteorological parameter acting as a driver for the high CH 4 emission from Shore 120 • .Note that the CH 4 emissions from the pond and the tundra have a similar magnitude under moderate wind speed conditions.When comparing the ratio of CH 4 emissions to CO 2 emissions, we find that fluxes with an open-water weighted footprint fraction of more than 60% have a ratio of CH 4 /CO 2 = 0.057 0.104 −0.049 (Median 75% Percentile 25% Percentile ), while for fluxes intervals with less than 20% open-water contribution we observe a negative ratio (due to the negative CO 2 fluxes) with a larger spread of CH 4 /CO 2 = -0.0100.021 −0.028 (Median 75%Percentile 25% Percentile ).The distributions of these two ratios are significantly different (Mann-Whitney-U test, p < 0.01).When considering only night-time fluxes (PAR<20 µmol m −2 s −1 ), the ratio of fluxes with an open-water weighted footprint fraction of more than 60% is similar (CH 4 /CO 2 = 0.060 0.076 0.049 , Median 75% Percentile 25% Percentile ), whereas ratio with less than 20% open-water contribution is now positive (CH 4 /CO 2 = 0.020 0.024 0.015 , Median 75% Percentile 25% Percentile ).The distributions of these two ratios are still significantly different (Mann-Whitney-U test, p < 0.01).

Upscaled CO 2 flux
We use the estimated aquatic CO 2 flux from the merged polygonal pond and the modeled CO 2 flux from the semi-terrestrial tundra to linearly up-scale the CO 2 flux for the polygonal-tundra landscape of Samoylov Island (excluding larger thermokarst lakes, method described in section 2.4.5).
We estimate that the tundra CO 2 uptake would decrease by ∼ 11% when including the CO 2 flux from ponds compared to a completely semi-terrestrial tundra.The modeled CO 2 -C flux from the semi-terrestrial tundra (without consideration of pond fluxes) accumulated to -16.29 ± 0.43 g m −2 during the observation period (60.5 days).Separated into months, it amounts to -15.01 ± 0.26, -3.56 ± 0.33 and +2.35 ± 0.11 g m −2 in July (19.8 days), August (31 days), and September (9.7 days), respectively.When including the CO 2 flux from the merged polygonal pond as representative for all ponds on Samoylov island, the resulting estimate of the landscape CO 2 flux amounts to -14.47 ± 0.40 g m −2 (60.5 days) and to monthly fluxes of -13.75 ± 0.24, -2.99 ± 0.31, and +2.27 ± 0.10 g m −2 in July (19.8 days), August (31 days), and September (9.7 days), respectively.Thus, in August, the estimate of CO 2 uptake is reduced most.In September, the estimate of landscape emissions is decreased   et al., 2003) to 0.22 g m −2 d −1 (Jonsson et al., 2008).Our estimate of 0.12 0.24 0.0014 g m −2 d −1 is therefore well within the range of aquatic CO 2 -C fluxes observed with the EC method.Other studies using different methods report a wider range of aquatic CO 2 fluxes in arctic regions.These fluxes range from a CO 2 -C uptake (-0.14 g m −2 d −1 , Bouchard et al. (2015)) to substantial emissions of CO 2 -C (up to 2.2 g m −2 d −1 , Abnizova et al. ( 2012)).A modeling study involving multiple lakes in north-eastern European Russia found close to zero emissions (0.028 g m −2 d −1 , Treat et al. (2018)).Perhaps the most striking finding is that our estimates of aquatic CO 2 emissions are approx.12-18 times smaller than previously reported for aquatic CO 2 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 model (TBL) after Liss and Slater (1974) has been used to estimate CO 2 emissions from analyzed water samples.However, one other study found good agreement between the EC method and the TBL (Eugster et al., 2003), so we cannot conclusively explain the differences.
Our approach of combining a footprint model with a land cover classification to extract fluxes from different land cover classes allows us to determine the pond CO 2 flux.We report an uncertainty range of the pond CO 2 flux; however, we can not identify the full uncertainty of this flux in this novel approach due to the unknown uncertainty of the footprint analysis.
Still, the pond CO 2 flux results are plausible and in the correct order of magnitude for two reasons.First, a reduced diurnal variability has been observed when the pond influences the flux signal (figure 4).This indicates that the respiration rate from the 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 ponds have a lower vegetation density than the tundra, there is less photosynthesis.Second, when focusing on night-time fluxes only, when only respiration occurs and no carbon uptake, there is a decrease in CO 2 emission with an increasing weighted footprint fraction of open water (shown in figure 3), also indicating reduced decomposition in the pond.Overall, the lower emissions from the pond compared to the semi-terrestrial tundra are reasonable.

CH 4 flux
We observe large differences in CH 4 emission from different wind sectors.CH 4 emissions from Shore 120 • are significantly higher than from Shore 50 • , pond or tundra, independently of meteorological conditions (see section 3.3).Especially the difference between Shore 120 • and Shore 50 • is astounding since the shorelines share many characteristics.Both extend radially (in a straight line) from the EC tower (see figure 1  vary much between the two shorelines.Both shorelines have a water depth between a few centimeters and a few decimeters within meters away from the shore (see data from Boike et al. (2015a)).As previously described in section 2.1, both shorelines are dominated by Carex aquatilis, and from visual inspection we can not find differences in shoot density.We, therefore, assume that the vegetation type does not play a major role in explaining the differences between the CH 4 emission from Shore 120 • and from Shore 50 • .We also examine the evolution of the shorelines at the merged polygonal pond to check whether erosion along the shoreline could drive the high CH 4 emissions.We compare a coarse image from 1965 (U.S. Geological Survey, EROS Center, 1965) with the current shoreline; yet, we can not identify signs of recent erosion.Also high resolution aerial images of the 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.Thus, we exclude erosion as a driving factor of high CH 4 emissions.
We also consider the possibility that local ebullition of the pond could lead to high CH 4 emissions from Shore 120 • .We apply the method proposed by Iwata et al. (2018) to check for signs of ebullition events.This method uses the 20 Hz raw concentration of CH 4 to detect short-term peaks in CH 4 that originate from ebullition events.However, we can not detect ebullition events in the 20 Hz raw data.
In summary, many causes, such as meteorological conditions, vegetation type, coastal erosion, and intense ebullition events, can be excluded as driving factors.Therefore, the most likely cause of the higher CH 4 emissions from Shore 120 • might be 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).So, a future visual inspection of trapped CH 4 bubbles in the ice column during wintertime, as proposed in Vonk et al. (2015), could reveal more information about the cause of the higher CH 4 emission from the Shore 120 • .
Excluding the high emissions from Shore 120 • , the CH 4 emission from the merged polygonal pond and the tundra surface have a very similar magnitude under similar meteorological conditions.In both landscape types, CH 4 is produced under anoxic conditions, but the emission pathways differ substantially.In dense soils, methane diffuses through upper soil layers and can oxidize before reaching the surface.In contrast, methane emitted in ponds can reach the surface quickly through ebullition or higher plant-mediated transport in addition to diffusion.We expected bigger differences between CH 4 emissions from the pond and the tundra due to the different emission pathways.Yet, as shown in figure 6, c) and A2, we see no significant difference in CH 4 emission from the open-water areas of the merged polygonal pond and the tundra surface.Since many other ponds are smaller than the pond (making them unsuitable for studying with the EC method), and since smaller ponds tend to be stronger emitters (Holgerson and Raymond, 2016;Wik et al., 2016), our measurements might be a lower limit of overall pond-CH 4 emissions.We estimate a CH 4 -C flux of 13.38  (Sieczko et al., 2020) and, in contrast, a study finding CH 4 -C emissions of up to 6432 mg m −2 d −1 in North-East Canada (Bouchard et al., 2015).The wide range of waterbody methane emissions cautions us to be careful when generalizing our results even for Samoylov Island, especially since the emissions within the pond are already heterogeneous.Instead, we would like to highlight the need for more measurements of CH 4 fluxes to understand the variability of pond-CH 4 emissions better.

Upscaling the CO 2 flux
We upscale the CO 2 emissions for the island of Samoylov, the area where we have access to the high-resolution land cover classification.We find that the inclusion of pond-CO 2 emission would considerably (11%) decrease the estimate of the polygonal tundra landscape's sink function.A similar approach by Abnizova et al. ( 2012) found a potential increase of 35 -62 % in the estimate of CO 2 emission from the Lena River Delta when including small ponds and lakes into the aquatic CO 2 emission.
If we follow the upscaling approach by Abnizova et al. ( 2012) and consider overgrown water as part of the ponds, we even find a CO 2 emission reduction of 19%.Also, Kuhn et al. (2018) found waterbodies in arctic regions to be an important source of carbon, which could outbalance the tundra's sink function in a future climate.In summary, our results demonstrate that aquatic CO 2 emissions can substantially influence the carbon balance of the polygonal tundra during the growing season.We expect that the impact of ponds on the carbon balance would be even bigger when accounting for CH 4 due to the locally high emissions, and because from the pond more CH 4 gets emitted per mole of CO 2 compared to the tundra.

Conclusions
We find that waterbodies are a carbon source while the surrounding tundra is a carbon sink during the period July -September in agreement with prior studies (Abnizova et al., 2012;Jammet et al., 2017), even if we observe much lower aquatic CO 2 fluxes compared to previous work at the same study site (Abnizova et al., 2012).Using a novel approach to disentangle the EC fluxes from different land cover classes, we estimate that during our measurement period, the landscape CO 2 sink is reduced by 11% when including ponds rather than only considering semi-terrestrial vegetated tundra.We expect lakes to have a similar effect on the budget, though a smaller one, since lakes (a) tend to emit less greenhouse gases than ponds (Holgerson and Raymond, 2016;Wik et al., 2016) and (b) cover a similar area as ponds in our study site (Abnizova et al., 2012;Muster et al., 2012).
In contrast to the spatially more homogeneous CO 2 emissions, small-scale heterogeneity in CH 4 emissions make it difficult to find drivers of CH 4 emissions.We cannot pinpoint the drivers behind the high emissions along parts of the coastline, which might be caused by seep ebullition.Thus, we cannot estimate the impact of this heterogeneity on the landscape scale and, therefore, refrain from upscaling CH 4 emissions.Additionally, the open-water fluxes presented in this paper originate from a single merged polygonal pond since the other ponds surrounding the EC tower are too small to extract their fluxes with the footprint method applied here.So, we do not account for spatial variability of CH 4 emissions between 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 as the While being ill-suited for very small ponds, we want to underline that the EC method is appropriate for observing greenhousegas fluxes from ponds with an area as small as 0.024 km 2 .The EC method has a higher temporal resolution than the TBL ; Wik 1 https://doi.org/10.5194/bg-2021-212Preprint.Discussion started: 12 August 2021 c Author(s) 2021.CC BY 4.0 License.
, b).It has a size of about 5 km 2 and consists of two geomorphologically different parts.The western part (∼2 km 2 ) is a floodplain and regularly flooded during the annual spring flood.The eastern part (∼3 km 2 ), a late-Holocene river terrace, is characterized by polygonal tundra.The partially degraded polygonal tundra at this study site shows a high spatial heterogeneity within a https://doi.org/10.5194/bg-2021-212Preprint.Discussion started: 12 August 2021 c Author(s) 2021.CC BY 4.0 License.few meters.Dry and wet vegetated parts are interspersed with small and large ponds (<1 m 2 ->10000 m 2 ) and with large thermokarst lakes (up to 0.05 km 2 , Boike et al. (2015a);

Figure 1 .
Figure 1.Study site with an overview of Russia (a), the Lena River Delta (b), Samoylov Island with the surrounding Lena River in blue (c), and a close-up look at the study site (d).The EC tower is marked as a black cross with the cumulative footprint (see section 2.4.2) in gray shades surrounding the EC tower.The outline of the land cover classification from section 2.4.1 is shown in a blue line (c).In (d), the detailed land cover classification is shown in blue (open water) and green shades (dark green: dry tundra, medium green: wet tundra, and light green: overgrown water).The merged polygonal pond studied here is outlined in red.Map data from © OpenStreetMap contributors 2020, distributed under the Open Data Commons Open Database License (ODbL) v1.0 (a & b) and modified after Boike et al. (2012) (c & d).
We split the datasets into a training (70%) and a validation (30%) data set to test model performance.Additionally, we exclude CO 2 fluxes from the direction of the merged polygonal pond from the training dataset to obtain a dataset consisting of as https://doi.org/10.5194/bg-2021-212Preprint.Discussion started: 12 August 2021 c Author(s) 2021.CC BY 4.0 License.much semi-terrestrial tundra as possible since we do not expect photosynthetic activity in the non-overgrown merged polygonal pond.
a gap-filled dataset of CO 2 flux from vegetated tundra.https://doi.org/10.5194/bg-2021-212Preprint.Discussion started: 12 August 2021 c Author(s) 2021.CC BY 4.0 License.2.4.5 Up-scaled CO 2 fluxTo evaluate the impact of ponds on landscape CO 2 flux without the influence of ponds, we estimate a polygonal-tundra landscape-CO 2 flux (F Landscape ) including ponds by linearly combining the two landscape forms (ponds and vegetated tundra):

Figure 2 .
Figure2shows the observed CO 2 fluxes plotted against the wind direction.The CO 2 flux exhibits a high temporal variability between positive and negative CO 2 fluxes from most wind directions.In the wind direction sector between 60 • -120 • , the flux is dominated by the merged polygonal pond.The flux signal from this sector has a smaller variability (standard deviation of 0.073 g m −2 d −1 ) than the fluxes from the other wind direction sectors (standard deviation of 0.56 g m −2 d −1 ).Additionally, we observe a lower respiration rate from the pond than from the semi-terrestrial tundra.Figure3shows the observed CO 2 fluxes plotted against the weighted footprint fraction of open-water in each flux when only considering night-time fluxes (PAR<20 µmol m −2 s −1 , thus only respiration).The fluxes decrease with an increasing open-water contribution.Another part of the CO 2 variability stems from the diurnal cycle.We compare the diurnal cycle of the CO 2 fluxes from the merged polygonal pond (estimated following equation 4) and the semi-terrestrial tundra (equation 5) in Figure 4. We see a less pronounced diurnal CO 2 cycle from the direction of the merged polygonal pond (blue) compared to the diurnal CO 2 cycle from the tundra (green).All data from the merged polygonal pond combined result in a CO 2 -C flux of 0.13 0.24 0.00 g m −2 d −1 (Median 75% Percentile 25% Percentile ).

Figure 3 .
Figure 3. Scatter plot of observed CO2 fluxes against the weighted footprint fraction of open water in each flux interval with temperature as color.Only flux intervals at night time (PAR<20 µmol m −2 s −1 ) are shown.

Figure 4 .
Figure 4. Diurnal cycle of modeled CO2-C flux from the merged polygonal pond (blue, eq.4) and the tundra (green, eq. 5) as violin plots for each half-hour flux interval.Blue and green crosses mark the mean CO2-C flux during each half-hour flux interval.

Figure 5 .
Figure 5. Polar plot of CH4-C flux with respect to the wind direction at the EC tower.Positive values outside of the dotted black line represent CH4 emission, and inside of the line, CH4 uptake during one half-hour period.The values 0, 20, 40, and 60 indicate the magnitude of the CH4-C flux in mg m −2 d −1 .The color represents the percentage of open water weighted footprint fraction in each flux interval.The red boxes indicate the mean CH4 flux of 5°wind direction intervals (red lines indicate the first standard deviation).

Figure 6 .
Figure 6.Violin plots of CH4 emissions at the EC tower separated into four different wind direction classes.Medians of CH4 emission distributions are shown as red lines, and 75 th & 25 th percentile are shown as black lines.
), thus contribute similarly to the EC flux.The underwater topography does not https://doi.org/10.5194/bg-2021-212Preprint.Discussion started: 12 August 2021 c Author(s) 2021.CC BY 4.0 License.Table 1.Daily mean water-atmosphere CO2 & CH4 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 proposed in this study.All fluxes are given ± standard deviation, except of fluxes from this study are given as Median 75% https://doi.org/10.5194/bg-2021-212Preprint.Discussion started: 12 August 2021 c Author(s) 2021.CC BY 4.0 License.tundra fluxes.Consequently, the main impact of ponds on the landscape CH 4 budget might be through plant-mediated transport and local ebullition.

Figure A2 .
Figure A2.Histogram of permutation tests between the medians of CH4 emissions from different wind direction classes in figure 6, b).All medians from flux observations during moderate wind speed conditions.

Figure A3 .
Figure A3.Histogram of permutation tests between the medians of CH4 emissions from different wind direction classes in figure 6, c).All medians from flux observations during moderate air temperature conditions.