Carbon emission from thermokarst lakes in NE European tundra

Emission of greenhouse gases (GHGs) from inland waters is recognized as highly important and an understudied part of the terrestrial carbon (C) biogeochemical cycle. These emissions are still poorly quantified in subarctic regions that contain vast amounts of surface C in permafrost peatlands. This is especially true in NE European peatlands, located within sporadic to discontinuous permafrost zones which are highly vulnerable to thaw. Initial measurements of C emissions from lentic waters of the Bolshezemelskaya Tundra (BZT; 200,000 km2) demonstrated sizable CO2 and CH4 concentrations and fluxes to the atmosphere in 98 depressions, thaw ponds, and thermokarst lakes ranging from 0.5 × 106 to 5 × 106 m2 in size. CO2 fluxes decreased by an order of magnitude as waterbody size increased by > 3 orders of magnitude while CH4 fluxes showed large variability unrelated to lake size. By using a combination of Landsat‐8 and GeoEye‐1 images, we determined lakes cover 4% of BZT and thus calculated overall C emissions from lentic waters to be 3.8 ± 0.65 Tg C yr−1 (99% C‐CO2, 1% C‐CH4), which is two times higher than the lateral riverine export. Large lakes dominated GHG emissions whereas small thaw ponds had a minor contribution to overall water surface area and GHG emissions. These data suggest that, if permafrost thaw in NE Europe results in disappearance of large thermokarst lakes and formation of new small thaw ponds and depressions, GHG emissions from lentic waters in this region may decrease.

Arctic and subarctic lakes emit sizeable amounts of greenhouse gases (GHGs) to the atmosphere (Tranvik et al. 2009;Kosten et al. 2010). High CO 2 and CH 4 fluxes are partly driven by microbiological and physicochemical processing of labile organic carbon released from thawing permafrost (van Huissteden et al. 2011;Vonk et al. 2015). Because climate warming and permafrost thaw may strongly enhance this impact Tan and Zhuang 2015;Bartosiewicz et al. 2016), numerous works have measured carbon emissions from inland waters in circumpolar territories (Rautio et al. 2011;Langer et al. 2015;Wik et al. 2016). In these regions, abundant thermokarst lakes are formed via thawing of frozen organic or mineral soils Sjöberg et al. 2013;Bouchard et al. 2017). There are also instances where certain lakes disappear entirely due to drainage (Smith et al. 2005). However, there is significant bias in the geographical coverage of available field measurements (see Fig. 1). The vast majority of seasonally resolved data on thermokarst lakes and ponds were collected in Northern Europe de Wit et al. 2018;Kuhn et al. 2018), Alaska (Sepulveda-Jauregui et al. 2015;Stackpoole et al. 2017), and Canada (Tank et al. 2009;Laurion et al. 2010;Arsenault et al. 2018;Matveev et al. 2018); only sparse measurements were conducted in Eastern Siberia (Walter et al. 2006(Walter et al. , 2007Abnizova et al. 2012;Walter Anthony et al. 2014;Tan et al. 2017). Specifically, as Matveev (2019) also note, there are a limited number of C measurements in peatland thermokarst lakes as compared to those for yedomalike permafrost.
Recently, a multiseasonal study with large geographic coverage in western Siberian peatlands demonstrated high C emissions from both rivers and thermokarst lakes (Serikova et al. 2018(Serikova et al. , 2019. While western Siberia can serve as a proxy for peatlands in northern Eurasia, the lack of data on continental planes with permafrost peatlands prevents the direct extrapolation of western Siberia lake GHG concentrations and fluxes to other regions. Among large territories of frozen peatbogs accessible all year round for sampling, the Bolshezemelskaya Tundra (BZT) of northeastern Europe (area~200,000 km 2 ) is especially interesting as it contains a vast amount of thermokarst lakes (Goldina 1972;Davydova et al. 1994), combines isolated, sporadic and discontinuous permafrost zones, and exhibits a generally milder climate and with a higher terrestrial biome productivity compared to western Siberia (Zamolodchikov and Karelin 2001). There is very limited information on GHG emissions from inland waters of the European Russian tundra. CO 2 and CH 4 emissions have been assessed for a single lake and river in the Vorkuta region (Heikkinen et al. 2004); however, seasonal variations of CO 2 emissions have been measured in three thermokarst lakes of this same region (Marushchak et al. 2013). There is recent interest in the BZT as over past decades lakes in this territory exhibited sizeable increases in summer temperature and pCO 2 , presumably due to enhanced bacterial respiration of allochthonous dissolved organic matter (DOM) from thawing permafrost (Drake et al. 2019).
In order to quantify GHG emissions from inland waters of the BZT, we measured CO 2 and CH 4 concentrations and atmospheric fluxes in thermokarst lakes, thaw ponds, and depressions of various sizes during the three main hydrological seasons (spring, summer, and autumn). We further measured the size distributions and water surface areas of lakes for the full BZT territory using remote sensing. We aimed to (1) assess the effect of waterbody size (from 1 to 10 5 m 2 ) and seasonality on C emissions and (2) quantify the overall emissions flux of CO 2 and CH 4 to the atmosphere from lakes of the BZT territory. We believe assessment of GHG pools and exchange fluxes in thermokarst waterbodies of the NE European tundra combined with evolving lake sizes, distributions, and surface areas over past decades can allow forsight into potential future C emissions to the atmosphere from lentic aquatic ecosystems in northern Eurasia.

Study site and methods
Thermokarst lakes and thaw ponds of BZT, NE Europe BZT (NE Europe) is the largest frozen peatland in Europe. It is covered by permafrost: discontinuous in the eastern part and sporadic to isolated in the western part (Brittain et al. 2009). Peat thickness in the BZT ranges from 0.5 to 5 m depending on the landscape context (Maksimova and Ospennikov 2012;Iglovsky 2016). The dominant elevation here is between 100 and 150 m; these altitudes are comprised of hills and moraine ridges which are mainly composed of sands and silts with intermixed boulders. Thermokarst lakes are located between these moraines and ridges. In the southern forest-tundra zone, the dominant soils are histosols of peat bogs and podzol-gleys. The mean annual temperature is −3.1 C with a mean annual precipitation of 503 mm. The dominant vegetation in the tundra zone is mosses, lichens, and dwarf shrubs.
Permafrost thaw in high latitude peatlands leads to formation of "thermokarst" (thaw) lakes due to frozen peat and mineral soil subsidence. The formation and development of thermokarst lakes has been investigated over the years and appears to follow a cycle (Matthews et al. 1997;Hinkel et al. 2003;Audry et al. 2011;Kirpotin et al. 2011;Pokrovsky et al. 2011). The cycle starts with frozen peat subsidence caused by warming and formation of a depression in moss/lichen cover that is filled by thaw water Pokrovsky et al. 2011). The depression grows due to ongoing subsidence of peat in the soil and instability of the edges until it reaches lake size. A mature thermokarst lake can eventually be drained into another lake or into the hydrological system (local streams and rivers). Gramineae-like vegetation colonizes along the drained lake bottom followed by mosses and lichens. Gradually, the renewed permafrost heaves the soil, creating small isolated mounds that merge into a flat uniform palsa plateau. Ultimately, the cycle can starts over again through formation of another wet depression. This sequence is an important part of the hydrological continuum in the BZT region, and it strongly influences organic carbon and GHG dynamics in surface waters (Shirokova et al. 2019).
Sampled waterbodies and basic physical and chemical characteristics and parameters of the CO 2 system are listed in Supporting Information Table S1. For convenience, all sampled waterbodies (generally referred to as "lakes") were classified according to their surface area: depressions (< 10 m 2 ), thaw ponds (10-1000 m 2 ), and thermokarst lakes (> 1000 m 2 ). Altogether, 48 thermokarst lakes, 10 thaw ponds, and 40 depressions were sampled for hydrochemical parameters, dissolved C, CO 2 , and CH 4 concentrations. GHG emissions for these bodies of water were calculated from CO 2 and CH 4 concentrations and wind speed (see "CO 2 and CH 4 flux measurements and calculations" section). One key site (~43 km SE from the town of Naryan-Mar) was selected for seasonal and annual monitoring of lakes with two additional sites (Shapkino and Khorey-Ver) to account for spatial heterogeneity over summer (July 2016); these sites are shown in Fig. 1 and Supporting Information Fig. S1a. In order to reveal seasonal and annual variability of GHG concentrations and fluxes into the atmosphere in BZT depressions, thaw ponds, and thermokarst lakes, we visited the Naryan-Mar site in spring (16-17 June 2017, sampling 7 lakes and 3 depressions), summer (17-23 July 2017, sampling 11 lakes and 13 depressions), and autumn (02-04 October 2017, sampling 11 lakes and 12 depressions). To account for interannual variability, three lakes from the Naryan-Mar site were sampled in July of 2015, 2016, 2017, and 2018. These four sampled years were rather contrasting in mean monthly and annual temperature and precipitation (see Supporting Information Table S2). The summer of 2015 was cold (−4.8 C below normal for July) while summers of 2016, 2017, and 2018 were warm (+5.0 C, +4.0 C, and +2.5 C, respectively, above normal for July). The summers of 2016 and 2018 were dry (31% and 36% of normal precipitation for July) whereas the summers of 2015 and 2017 were more humid (83% and 56% of normal precipitation for July).
The spatial variability of dissolved organic carbon (DOC) and GHG concentrations was assessed in July 2016. A total of 19 lakes, thaw ponds, and depressions were sampled in the Shapkino site. Another 33 lakes and ponds were visited in the Khorey-ver site. Both of these sites were compared with six lakes in the Naryan-Mar site sampled during the same period.
Thermokarst lakes and thaw ponds in the BZT, similar to those in western Siberia Manasypov et al. 2014;Serikova et al. 2019), are mainly shallow (large lakes are only 0.5-1.5 m deep and small thaw ponds are < 0.5 m deep) and are not stratified with respect to temperature. There are only a few exceptions where lakes are 1.5-2.0 m deep but these lakes constituted less than 1% of all studied (see Supporting Information Table S1). For simplicity, we refer to all waterbodies studied in this work as "lakes" although their size ranges from < 1 m 2 to > 7 km 2 .

Lake size and area inventory using satellite imagery
Remote sensing data utilized for this study include Landsat-8 medium resolution (MR) and GeoEye-1 ultra-high resolution (UHR) imagery; imagery were used to map distribution of lakes throughout the BZT territory (Supporting Information Fig. S1b). Landsat-8 images were taken in the late July through early August 2017 and 2018. GeoEye-1 images were collected in the summer months of 2011 through 2017. These summer periods, when the lakes are ice-free, correspond to minimum coverage of the territory by lakes and minimum seasonal change in water level. A combination of Landsat-8 images over 2 yr was necessary as observation over the course of 1 yr could not provide full coverage of the territory.
For the entire BZT territory affected by permafrost (195,000 km 2 ), we used a mosaic of 50 Landsat-8 MR images. Images were processed using standard ArcGIS 10.3 software tools using the Fmask algorithm (see Polishchuk et al. 2014Polishchuk et al. , 2017Polishchuk et al. , 2018 for method description). Given the sizable areal coverage and accessibility challenges, no ground calibration has been performed in BZT and we rely on ground truthing from former studies of western Siberia where there is a similar distribution of thermokarst lakes. Taking into account the 30 m resolution of MR imagery, the minimum size of lakes examined is limited to 5000 m 2 . This corresponds to six pixels 30 × 30 m (900 m 2 ) in size, which is sufficient for reliable separation of lakes from background noise. The uncertainty of measuring the surface area of a lake using Landsat-8 is 3.5% for lakes > 10 5 m 2 surface area and 7% for lakes > 2 × 10 4 m 2 (Bryksina and Polishchuk 2013;Kornienko 2014). UHR GeoEye-1 imagery, evenly distributed throughout the territory, were used to quantify the number and total areas of small lakes and ponds. For test sites, we chose representative permafrost peatland landscape, similar to those sampled in this work. With a spatial resolution of 0.6 m, pixel size in these images is 0.36 m 2 , providing a minimal pond size of 2 m 2 . GeoEye images were processed using ArcGIS 10.3 at 11 test sites distributed throughout the study area.
To plot the distribution of lakes in a wide range of sizes, we used 23 partial size ranges of lakes taken on arithmetic progression scale: from 2 to 5 m 2 , 5 to 10 m 2 , 10 to 20 m 2 , and so on up to 10 8 m 2 . Note that remote sensing did not allow for quantification of the proportion of smallest depressions (< 2 m 2 ). However, based on our previous work on similar permafrost peatlands in western Siberia, the overall contribution of such small waterbodies to water volume, surface area, and GHG pools is negligible (Polishchuk et al. 2018).
To combine the MR and UHR imagery, we used a technique developed for thermokarst lakes in western Siberia (Polishchuk et al. 2018). First, we constructed diagrams for the number and areal distribution of all thermokarst lakes detectable on MR images. Then, we built similar diagrams for small thermokarst lakes and very small thaw ponds using HR images of the 11 test sites across the BZT. This allowed for extrapolation of the number and area of lakes across the entire BZT territory, as determined by UHR images from the 11 test sites. Finally, we combined primary MR and UHR histograms into integral histograms for lake numbers and overall areas. These histograms included lakes and small and very small thaw ponds spread over five orders of magnitude in size. We estimated uncertainty in lake size and area distribution stemmed from the process of overlapping images up to 15%, following the methods applied for much larger territory of western Siberia (Polishchuk et al. 2018). The average depth was measured for all sampled lakes (see Supporting Information Table S1). For the regional estimation of water volume, we used an empirical approach developed for similar peatlands of western Siberia (Polishchuk et al. 2018), where we approximated the dependence of lake size-lake depth and used it for each lake size bin.

Sampling and chemical analyses
The surface waters were sampled from the shore (depression and small thaw ponds) or via small PVC boats which gently moved across the lake surface. On small ponds, we used a rope held by a person standing on shore. This was done to avoid any disturbance of lake sediments during sampling. At each location, we measured water temperature and dissolved oxygen concentration (WTW Oxi 3205 DurOx 325-3), pH and specific conductivity (WTW Multi 3320 multiparameter) as well as the water depth (Cole-Parmer), air temperature, and atmospheric pressure (Silva). Dissolved oxygen saturation was measured at 50 cm intervals until reaching the lake sediment interface. Water samples for DOC, dissolved inorganic carbon (DIC), and dissolved CO 2 and CH 4 were collected from 0.2 to 0.3 m depths and analyzed following methods of Manasypov et al. (2015) and Serikova et al. (2019). In brief, for DOC analyses, collected waters were filtered on-site in prewashed 30-mL polypropylene Nalgene ® bottles through single-use Minisart filter units (0.45 μm pore size, Sartorius, acetate cellulose filter). DOC and DIC were analyzed using a Carbon Total Analyzer (Shimadzu TOC VSCN) with an uncertainty better than 3% (see Prokushkin et al. 2011 for methodology). Note that in small acidic (pH < 5) thaw ponds with high DOC (> 30 ppm) and low DIC (< 1-2 ppm), the DIC measurements could be subjected to significant uncertainty. For GHG analyses, unfiltered water was sampled in 60-mL Serum bottles and closed without air bubbles using vinyl stoppers and aluminum caps and immediately poisoned by adding 0.2 mL of saturated HgCl 2 via a two-way needle system. The CH 4 and CO 2 concentrations were measured with 3-5% uncertainty using a Bruker GC-456 gas chromatograph equipped with flame ionization and thermal conductivity detectors. Replicate samplings from the same lake yielded an average reproducibility of 6-8%. Dissolved (< 0.45 μm) nutrients (NO − 3 , NH + 4 , PO 3− 4 ) and trace elements were measured using spectrophotometry and Inductively Coupled Plasma Mass Spectrometry (ICP_MS), respectively, as described elsewhere for humic lake waters (Manasypov et al. 2015;Chupakov et al. 2017).
CO 2 and CH 4 flux measurements and calculations CO 2 and CH 4 fluxes were calculated from wind speed and surface-water gas concentrations. This technique is based on the two-layer model of Liss and Slater (1974), and widely used for GHG flux assessment (Repo et al. 2007;Juutinen et al. 2009;Laurion et al. 2010;Elder et al. 2018). The gas transfer coefficient was taken from Cole and Caraco (1998): Zabelina et al.
C emission from thermokarst lakes where U 10 is the wind speed taken at 10 m height, without explicitly accounting for lake size and shape (i.e., Vachon and Prairie 2013). Average daily wind speed was retrieved from official data of the nearest weather station (Naryan-Mar) as published by Rosgidromet for the day of sampling. Further, 8 and 15 waterbodies were used for measuring emissions of CO 2 and CH 4 , respectively, using floating chambers in order to validate modeled wind-based fluxes.
The direct CO 2 fluxes were measured using small (30 cm diameter) floating chambers equipped with nondispersive infrared CO 2 logger (ELG, SenseAir; Alin et al. 2011;Bastviken et al. 2015). These CO 2 loggers were calibrated in the lab against pure N 2 -CO 2 mixtures before each sampling campaign. We placed 2-3 chambers per lake along the transect from the shore to the center of the lake, if the lake size allowed. For small thaw ponds (< 100 m 2 ) and permafrost depressions (< 10 m 2 ), two chambers were typically deployed. The CO 2 accumulation rate inside each chamber was recorded at 5 min interval. Chambers were ventilated every 10-24 h and CO 2 accumulation rate inside each chamber was calculated via linear regression for about 30 min accumulation time (Serikova et al. 2019). All measurements had a linear CO 2 increase over time with R 2 > 0.98. CH 4 flux measurements from the surface of several waterbodies (7 lakes and 1 permafrost depression in 2016; 3 lakes and 5 depressions in 2017) were performed using circular chambers (surface area = 0.154 m 2 ; volume = 20.0 L) following the design described in Guerin et al. (2007). In summary, within 45 min, four air samples were collected with a valve and a syringe from the chambers over 15 min intervals. Air samples for CH 4 were collected in 15 mL glass vials which contained 6 mol L −1 NaCl solution and capped with high density butyl stoppers and aluminum seals. Methane fluxes were calculated from the slope of the linear regression of gas concentration in the chamber vs. time. The fluxes were accepted when the coefficient R 2 of the linear regression was > 0.90. Emissions were considered as the mix of ebullition and diffusive when one or more abrupt increases of the floating chamber headspace methane concentration were observed and the R 2 of linear regression between methane concentration and the elapsed time was below 0.90. Only diffusive CH 4 fluxes were used to calculate the overall CH 4 flux. For example, any strong deviations from linear increase of CH 4 concentration linked to potential bubble evasion were discarded. Therefore, obtained numbers represent the lowest threshold of total methane emission. It should be noted that bubble evasion was strongly pronounced in all studied ponds and lakes; this could be the main reason for the strong difference between calculated and measured methane fluxes (see below). Each direct flux measurement was done together with a determination of the CH 4 concentration in surface water.
Although chamber measurements are the most accurate tool for obtaining direct flux in remote locations (Serikova et al. 2019), we used calculated fluxes for regional C emission upscaling given that direct chamber measurements represented only 7% and 19% of total CO 2 and CH 4 fluxes assessed. To estimate regional C emission, we used median values of CO 2 and CH 4 fluxes for each category of lake size, averaged over seasons. We exclusively considered lake size as the main controlling parameter on GHG, DOC, and DIC concentration because temporal and spatial variations of concentrations and fluxes were smaller than variations due to lake size (see "DOC, DIC, CO 2 , and CH 4 concentration and C emissions of thermokarst lakes as a function of lake area" section and "Seasonal, annual, and spatial variability of C concentrations and emission" section). For this, we used six size categories of waterbodies: < 10 m 2 , 10-10 2 m 2 , 10 2 -10 3 m 2 , 10 3 -10 4 m 2 , 10 4 -10 5 m 2 , and > 10 5 m 2 . Note that, similar to previous works in other Arctic lakes, GHG flux numbers represent a minimum estimate for the total annual CO 2 and CH 4 fluxes in Eastern European tundra thaw lakes because we did not include CH 4 ebullition (i.e., Sabrekov et al. 2017;Elder et al. 2018); however, we accounted for potential spring release of CO 2 and CH 4 that accumulated during the winter season under ice .

Statistical treatment of the data
A nonparameric Mann-Whitney U-test was taken for treating the data due to the short series of observation, nonnormal distribution of number of parameters, and some amount of outliers. Regressions were used to examine the relationships between DOC, CO 2 , CH 4 concentrations, and GHG fluxes to the atmosphere of the lake area. The normality of DOC, GHG concentration, and fluxes was assessed via the Shapiro-Wilk test. Because the data were not normally distributed, we used nonparametric statistics. Thus, the median, 1 st , and 3 rd quartiles were used to trace the dependence of DOC, CO 2 , CH 4 concentration, and GHG fluxes to the atmosphere on the area of the waterbody, season, and the year of sampling. This is consistent with other studies on GHG in lakes when dealing with non-normal data distribution (Erkkilä et al. 2018;Serikova et al. 2019). Differences in concentrations and fluxes between lake size classes, seasons, and years at the same and among different sites for the same seasons and lake size classes were tested using Mann-Whitney U-test for paired data sets with significance level at 0.05. The Spearman rank order correlation coefficient (Rs) (p < 0.05) was used to determine the relationship between GHG concentrations and fluxes and climatic variables. Further statistical treatment of a complete set of C concentration drivers in thermokarst lake waters included a stepwise regression, which allowed to test the effect of various hydrochemical and climatic parameters on behavior of DOC and GHG. All graphs and drawings were plotted using MS Excel 2016 and Statistica-8 software package (http://www.statsoft.com).

DOC, DIC, CO 2 , and CH 4 concentration and C emissions of thermokarst lakes as a function of lake area
The full set of measured DOC, CO 2 , CH 4 concentrations, and CO 2 and CH 4 fluxes into atmosphere are provided in Supporting Information Table S1. DOC concentrations typically ranged, for different lake classes, from 10 to 100 mg L −1 (median of 20-50 mg L −1 ) and decreased almost 10-fold as surface area increased from < 1 m 2 to 1 km 2 ( Fig. 2a and Supporting Information Fig. S2a). The DIC concentration was less than 10% of DOC and was invariant to lake size (Fig. 2b,  Supporting Information Fig. S2b). CO 2 concentration decreased up to 10-fold with increasing lake size (Fig. 2c). The median values of CO 2 concentration (100-500 μmol L −1 ) demonstrated a factor of 5 decrease with increasing lake size (from < 10 m 2 to > 1000 m 2 , Supporting Information Fig. S2c). The CO 2 fluxes progressively decreased from 500 to 1000 mmol m −2 d −1 in depressions (< 10 m 2 ) to about 100 mmol m −2 d −1 in lakes larger than 1000 m 2 (Fig. 2d). A threefold decrease in median CO 2 fluxes from ponds < 100 m 2 to lakes ≥ 1000 m 2 was also observed (Supporting Information Fig. S2d). A general decrease of DOC, CO 2 concentration, FCO 2 , and optical density at 254 nm (UV 254 ) with lake surface area across all sites and all seasons is confirmed by strong (p = 0.0000) power dependence (Supporting Information  Table S3). In contrast to CO 2 , CH 4 concentrations and fluxes to atmosphere exhibited strong variability, ranging from 0.1 to 100 μmol L −1 and mmol m −2 d −1 , respectively. There was a local maximum of CH 4 concentration and emission in thaw ponds of 10-1000 m 2 surface area (Fig. 2e,f, Supporting Information Fig. S2e,f, respectively). Considering all calculated fluxes to atmosphere (n = 135), the contribution of C-СН 4 to the total C emissions calculated by wind-based modeling ranged from 0.02% to 19% (with a mean of 2.6% and a median of 1.1%). Chamber-measured CH 4 fluxes in eight bodies of water constituted between 0.2% and 58% of the total CO 2 + CH 4 fluxes (Supporting Information Table S4a).
The CO 2 fluxes to atmosphere obtained by using wind-based models were on average 2.5 times (range of 0.8-18.9 times, n = 8) higher than those directly measured using chambers (Supporting Information Table S4a). It should be noted that Laurion et al. (2010) reported an average factor of 3.8 of the excess of calculated fluxes in comparison to measured CO 2 fluxes for similar peatland thaw ponds. The difference between measured and calculated CH 4 fluxes to atmosphere was significantly (p < 0.05) higher than for CO 2 (range 1.5-19.4, mean 6.8). In thermokarst lakes of western Siberian permafrost peatlands, CO 2 and CH 4 fluxes to atmosphere calculated with wind-based models (either Cole and Caraco 1998 or Vachon and Prairie 2013) were on average three and four times lower than measured fluxes (Serikova et al. 2019). These authors explained this discrepancy as due to the remoteness of the studied lakes in relation to their respective meteorological stations and the lack of local wind data. This is also true for most BZT sites, located approximately 40 km from the Naryan-Mar station. The most likely cause for CH 4 differences is bubble evasion (Casper et al. 2003) as it was not possible to sample shallow waterbodies (< 1 m depth) with zero disturbance of the organic detritus at lake bottom and in surrounding coastal waterlogged zones. Note that the agreement between the two methods was much better in lakes with mineral banks where visual bubble disturbance during sampling was minimal. Another reason for observed discrepancies may be due to the wind effect on the gas transfer coefficient being overestimated in small aquatic systems with small fetch and sheltered sampling points along the shoreline (Alin et al. 2011).
As a measure of variability associated with model parameters, we calculated the effects of variable wind speed for estimates of diffusive GHG fluxes to atmosphere using Eq. 1 for 25 lakes > 1000 m 2 sampled in July 2015-2018 at the Narian-Mar site (Supporting Information Table S4b). Typically, variation of the wind speed by a factor of 2 yielded a factor of 1.3-1.5 variation in relevant median values of GHG fluxes to atmosphere. Only drastic changes in wind speed (by a factor of 5 or more) had sizable effect (exceeding the interquartile range [IQR]) on CH 4 and CO 2 fluxes. Overall, this suggests that the offset between measured and calculated fluxes is not caused by uncertainties in wind speed alone.
Seasonal, annual, and spatial variability of C concentrations and emission GHG concentrations and fluxes to atmosphere showed larger spatial variability across lake sizes than temporal variability across seasons (June, July, and October). Differences in DOC, DIC, and GHG between seasons (2015-2018) were assessed for lakes > 1000 m 2 from the Naryan-Mar site (Supporting Information Table S5). For these lakes, the DIC and CO 2 concentrations were sizably lower in spring compared to summer (p < 0.05); however, there was no difference between summer and autumn (p > 0.05, Supporting Information Table S5 and Fig. S3b,c). The DOC showed an inverse pattern with spring < summer = autumn (p < 0.05, Supporting Information Fig. S3a). CH 4 concentrations and CO 2 and CH 4 fluxes to atmosphere followed the order "spring > summer" (p < 0.05) but there was no difference (p > 0.05) between summer and autumn (Supporting Information Fig. S3d-f). There were no seasonal differences (p > 0.05) in DIC, DOC, CO 2 , and CH 4 concentrations and fluxes for depressions (< 10 m 2 area). Interannual variations in DOC, CO 2 , CH 4 concentrations, and CO 2 and CH 4 fluxes (Supporting Information Fig. S4) did not exhibit any link to potential drivers like annual and/or seasonal temperature or precipitation for the specific year of observation (Supporting Information Table S2). In fact, despite sizeable contrast in mean June and July temperature and precipitation between 2015 and 2018 (p ≤ 0.01 Mann-Whitney U-test), no difference in DOC, DIC, CO 2 , and CH 4 was evident for these years (p > 0.05).
Using Spearman rank order correlation and multiple regression approach, we examined the impact of water temperature and cumulative precipitation and temperature during the summer, prior to the sampling campaign on the GHG, DOC and DIC concentrations, and GHG fluxes to atmosphere (Supporting Information Table S6a). Water temperature negatively impacted CO 2 and positively impacted DOC (R Spearman = −0.20 and 0.22, respectively, at p < 0.05). The Spearman rank order correlation yielded positive dependence of CH 4 and CO 2 on cumulative air temperature prior to sampling (R Spearman = 0.20-0.28) as well as on the mean air temperature during June and July (R Spearman = 0.26-0.30 in lakes; R Spearman = 0.43-0.46 in small thaw ponds) (Supporting Information Table S6a). Precipitation exerted a negative impact on GHG concentration and fluxes to atmosphere (R Spearman = −0.23 to −0.55) but a positive impact on specific UV absorbance (SUVA, R Spearman = 0.24-0.38). Overall, the impact of both temperature and precipitation was much more pronounced for small thaw ponds (< 500 m 2 ) compared to thermokarst lakes (> 500 m 2 ). The multiple regression approach demonstrated that the GHG concentration variability is partially determined by waterbody area and cumulative climatic parameters such as air temperature and precipitation (Supporting Information Table S6b). However, the multiple regression was not capable of explaining more than 20-40% of the overall variability in CO 2 , CH 4 , DOC, and DIC concentrations.
The variability of C concentration and emissions across sites was assessed for Naryan-Mar, Shapkino, and Khorey-Ver sites in July 2016. This comparison was conducted separately for thaw ponds with surface areas < 50 m 2 and 50-500 m 2 and for lakes with surface areas > 1000 m 2 (Supporting Information Fig. S5). The difference between sites was assessed via the nonparametric Mann-Whitney U-test (Supporting Information Table S7). All sites were similar in DOC, CO 2 , CH 4 concentrations, and C emissions to atmosphere except lakes in Shapkino which had three times higher DIC and CO 2 concentrations than those of Naryan-Mar and Khorey-Ver. This is consistent with the shallower position of mineral soils on clays and fluvioglacial deposits at the Shapkino site where some carbonate shells have been reported (Andreicheva 2007).
In order to assess possible control of hydrochemical parameters of waterbodies on CO 2 and CH 4 concentrations, we calculated Spearman rank correlation coefficients (Supporting Information Table S8). The concentrations of CO 2 and CH 4 negatively correlated with pH and O 2 , and positively correlated with DIC, DOC, and Total Dissolved Solids (TDS). CO 2 concentration positively correlated with PO 3 − 4 , P tot , NO − 3 , N tot , Si, DIC, DOC, Fe, Co, Ni and negatively correlated with NH + 4 . CH 4 concentration positively correlated with N tot , Si, DIC, DOC, Fe, Co, and Ni, and negatively correlated with NH + 4 and SO 2 − 4 (Supporting Information Table S8).

Regional-scale estimation of C emissions from BZT territory (based on satellite imagery)
For the full BZT territory (195,000 km 2 ), the total number of inventoried lakes ranging from 2 m 2 to 5 × 10 7 m 2 in size was 8.96 × 10 6 with a total water area of 7800 km 2 (Table 1 and Supporting Information Fig. S6a,b). Lakes < 5000 m 2 provided about 10% of overall lake water area in BZT. While small lakes (< 2000 m 2 ) provided 97% of the total lake numbers, their area contributed to only about 6% of total area coverage. Thermokarst lakes between 5000 and 10 8 m 2 in area accounted for only 1.5% of the lake count but provided > 90% of total lake area. Note that, from an optical remote sensing view point, there is no difference between thermokarst and nonthermokarst lakes because their reflections yield similar Table 1. Lake characteristics, C pools, and GHG fluxes inventory [median and (IQR)] for lakes of Eastern European (Bolshezemelskaya) Tundra (195,000 km 2 ).
The maximal contribution to both CO 2 and CH 4 pool is provided by lakes of > 10 5 m 2 surface area. Lakes larger than 10 5 m 2 provide > 43% of total CH 4 amount whereas the share from small lakes (< 1000 m 2 ) is < 20%. Note that small lakes contributed less to overall CO 2 and DOC levels than to CH 4 level.
Lake size range (m 2 ) Depth (m) Lake number  (Kravtzova and Bystrova 2009); this origin is also confirmed by their generally isometric shape as compared to elongated lakes with glacial origins. Both the lake abundance and area distribution obtained for the permafrostaffected part of BZT are in fairly good agreement (Supporting Information Fig. S7) with those reported for permafrostaffected parts of western Siberia (Polishchuk et al. 2017(Polishchuk et al. , 2018 and Swedish territory (Cael and Seekell 2016). In the BZT, the power law extrapolation (Pareto law) strongly overestimates the number of small lakes (< 1000 m 2 ) compared to empirical data (Supporting Information Fig. S7), which is in agreement with other regions of the world. The maximal contribution to overall GHG emissions from the territory of 195,000 km 2 (3.8 Tg yr −1 ) was provided by lakes > 100,000 m 2 whereas contributions from small thaw ponds (< 1000 m 2 ) were below 20% (Table 1; Fig. 3). The dispersity of the areal GHG fluxes to atmosphere is quite high achieving one order of magnitude in small thaw ponds and a factor of 2-5 for thermokarst lakes (Fig. 2). These uncertainties, visible as IQR in Supporting Information Figs. S2-S5, S8, do not originate from analytical procedures and rather reflect variability between sites and seasons. It is important to note that the contribution of small thaw ponds to overall C balance in BZT is below the actual uncertainties for total pools and emissions.

Discussion
C pools, fluxes to atmosphere, and possible mechanisms The highest DOC and CO 2 concentrations and CO 2 fluxes to atmosphere were observed in small (< 1000 m 2 ) thaw ponds. Overall, CO 2 and DOC concentrations and CO 2 emission decreased by an order of magnitude when lake size increased by > 3 orders of magnitude. Assuming that the size gradient of thermokarst lakes in permafrost peatlands reflects the age gradient Manasypov et al. 2014), high CO 2 flux to atmosphere from small waterbodies corresponding to the beginning of permafrost thaw and lake formation is likely largely sustained by relatively high input of allochthonous DOM. DOC concentrations in studied waterbodies were quite high (median 60-20 mg L −1 for different size bins, and decreasing with lake size) and exceeded values observed in soil pore waters of wetlands in the southern part of the studied territory (20-30 mg L −1 ; Avagyan et al. 2014). Also, given the likely short residence time of water in these small thaw ponds and depressions (Manasypov et al. 2015), this DOC input from the peat soil would be quickly replenished as well as rapidly bio-and photodegraded as was recently suggested for small waterbodies of eastern European tundra (Shirokova et al. 2019). The very low quality of DOC that remains after fast initial degradation of DOM in waterbodies of the BZT is confirmed by the high values of SUVA 254 (from 3 to 6 L mg −1 m −1 ) which was independent of the season and size of the waterbody (Supporting Information Fig. S9). Note that Fe concentrations in studied lakes (500--1000 μg L −1 , depending on the lake size) are too low to interfere with SUVA values (i.e., Poulin et al. 2014).
In large thermokarst lakes, both sediment source and suprapermafrost input of contemporary DOM likely support CO 2 fluxes to atmosphere. Allochthonous (terrestrial) organic carbon from peat and vegetation leachates is delivered into large lakes via shallow subsurface flow in summer and autumn and via surface overflow during spring (Manasypov et al. 2015;Ala-Aho et al. 2017). High concentrations of DOC and some macro-and micronutrients (P, N, Si, Fe, Co, Ni), positively correlating with pCO 2 in small thaw ponds, are due to intensive abrasion of peat at the shore and concomitant C and element delivery to the waterbodies via suprapermafrost flow Raudina et al. 2018). Upon maturation of lakes, and decrease in ratio of the circumference length to lake surface/volume, delivery of terrestrial DOM per surface area of lake or lake volume is decreasing. This increases the pH to above 4-4.5, thereby favoring phytoplankton and periphyton development as it is known in lakes of frozen peatlands of western Siberia (Pavlova et al. 2016). This may result in decreased CO 2 concentration and fluxes to atmosphere for these lakes. Furthermore, lake growth leads to certain steadystate conditions-when delivery of DOM matches its processing in the water column and burial in the sedimentsas was demonstrated in western Siberia (Manasypov et al. 2015). Thus, after some threshold lake area (10 3 -10 4 m 2 ), DOC and CO 2 concentrations remain constant with further increase in the lake area.
However, the lack of information on C content and the biodegradability of BZT lake sediments require in-depth studies of C processing in benthic ecosystems of the BZT. In particular, in order to distinguish the origin of CO 2 and CH 4 in thermokarst lake waters, radiocarbon dating of dissolved CO 2 and dissolved vs. ebullited CH 4 (i.e., Walter et al. 2006;Matveev et al. 2016;Elder et al. 2019) would be useful; however, such research is beyond the scope of this study. (light gray columns) and CH 4 (black columns) emissions for lakes in the BZT (200,000 km 2 ). To calculate these contributions, we used median concentrations and fluxes to atmosphere (averaged over all sites and all seasons, Supporting Information Fig. S2) and lake size distribution from GeoEye-1 and Landsat data (Supporting Information Fig. S6).

Regional yields and upscale of fluxes
The total carbon (CO 2 + CH 4 ) flux from the water surface of BZT lakes ranged from 100 to 400 mmol C m −2 d −1 , or 1.2 to 4.8 g C m −2 d −1 (with median values corresponding to lakes of different size bins). This value should be considered with caution, given a factor of 2.5 lower values obtained by chamber measurements compared to those calculated by windbased model. Given the scarcity of direct chamber measurements, the regional estimation is possible only based on calculated fluxes. Nevertheless, these numbers are comparable to C fluxes measured for western Siberian lakes (1.7 g C m −2 d −1 , range 0.6-5.3 g C m −2 d −1 ; Serikova et al. 2019) and with fluxes from small (4-150 m 2 ) thaw ponds of Northern European peatlands (3.4 g C-CO 2 m −2 d −1 , 0.1 g C-CH 4 m −2 d −1 ; Kuhn et al. 2018) as summarized in Supporting Information Table S9. At the same time, C fluxes in the BZT are higher than those of lakes in Arctic Alaska (Kling et al. 1992;Elder et al. 2018), thermokarst lakes in the Vorkuta region (Heikkinen et al. 2004;Marushchak et al. 2013), palsa peatlands of Canada (Matveev et al. 2016), and runnel ponds in continuous permafrost zone of the Canadian High Arctic (Negandhi et al. 2013). It should be noted that recent modeling of CO 2 emissions from Arctic (yedoma, thermokarst, and glacial) lakes across the globe recommended values between 0.036 and 0.08 g C-CO 2 m −2 d −1 (Tan et al. 2017) which is 10-100 times lower than CO 2 fluxes measured in BZT. Therefore, we caution against use of simulated values for all permafrost regions, especially peatlands. Western Siberia lowland (WSL), NE European tundra and, presumably, other peatbearing lowlands of Northern Eurasia might turn out to be prominent exceptions to these simulations; thus, they may require detailed spatial and seasonal analysis of C storage and emissions as their contribution to the global Arctic C cycle may be strongly underestimated. For example, even though the BZT represents 1% of overall permafrost-affected land area, its contribution alone can be comparable to the whole Arctic as estimated by Tan et al. (2017).
Upscaling the results to the full territory of BZT (195,443 km 2 ) yields present day CO 2 + CH 4 emission of 3.8 (AE 0.65) Tg C yr −1 (median AE IQR). This value is higher than the total annual CO 2 emissions from the North Slope of Alaska (0.95 Tg C yr −1 ) across a similar territory (231,157 km 2 ; Elder et al. 2018). Also, C emission flux to atmosphere from BZT lakes is 30% higher than the annual DOC + DIC export for the Pechora River (S watershed = 324,000 km 2 ) to the Arctic Ocean (1.1 × 10 6 t DIC yr −1 , 1.7 × 10 6 t DOC yr −1 ; Gordeev et al. 1996). When compared to the same land (watershed) surface area of this region, the C emission exceeds the lateral C export by a factor of 2.2. C emission from the BZT lakes represents 30% of that reported for the 1,300,000 km 2 of WSL (12 AE 2.6 Tg C yr −1 ; Serikova et al. 2019), a territory of similar wetlands and peatbogs. Normalized to total land area, the specific aquatic C yields to the atmosphere are, however, two times higher in the BZT than in the WSL (19.5 and (9.2 AE 2.0) × 10 −6 Tg km −2 yr −1 , respectively). It is important to note that the mean annual air temperature is 2-6 C higher in the BZT compared to the WSL. Precipitation for both regions is relatively similar. This suggests more active terrestrial C processing in inland waters of permafrost peatlands in European tundra compared to western Siberia. Because the permafrost in BZT is essentially discontinuous to sporadic, and the warming is more pronounced than in the WSL (Vasiliev et al. 2020), lakes of BZT could be expected to experience stronger and faster response to warming and permafrost thaw compared to those in western Siberia. Upscaling the average C emissions of European tundra and western Siberia peatlands ([14 AE 5] × 10 −6 Tg km −2 yr −1 ) to the overall area of permafrostaffected peatlands (2,864,000 km 2 ; Smith et al. 2007) yields a value of 40 AE 14 Tg yr −1 . However, this value should be considered with caution given the significant spatial heterogeneity of thermokarst lakes in permafrost landscapes (Vonk et al. 2019). Another issue is distinguishing between the GHG pattern of thermokarst and nonthermokarst ecosystems, including the winter period. Here, Synthetic Aperture Radar (i.e., Engram et al. 2013) can be extremely useful.
A critical issue is the emission from thaw ponds and thermokarst lakes compared to the overall C balance of permafrost wetlands. None of the BZT catchment has been subjected to extensive C exchange studies (i.e., by chamber and eddy covariance techniques coupled with aquatic export and emission estimation) in ways which it was performed for other permafrost-affected regions such as Northern Sweden (Lundin et al. 2016;Kuhn et al. 2018) or Alaska (Bogard et al. 2019). The C balance in Eastern European tundra (BZT) is poorly constrained due to strong seasonal and spatial variations that are not captured in existing measurements (i.e., Heikkinen et al. 2002;Marushchak et al. 2013). The upscaled value for terrestrial C uptake for this territory is −7.2 Tg yr −1 (Heikkinen et al. 2004) which strongly contrasts with the mean annual terrestrial C balance for the eastern part of the BZT (−0.8 AE 11.3 Tg yr −1 ; Zamolodchikov and Karelin 2001). Because the magnitude of these values is comparable to net C emissions from the waterbodies estimated in this study (+3.8 Tg C yr −1 ), lakes may play a crucial role in overall C balance of the European permafrost peatlands.
The ongoing climate warming and permafrost thaw in frozen Siberian wetlands leads to the drainage of large thermokarst lakes (Smith et al. 2005) and appearance of numerous new small thaw ponds (Polishchuk et al. 2014;Bryksina and Polishchuk 2015). Our results suggest that this may lead to a decrease in overall regional C emissions from inland water surfaces. At the same time, this change of regime is not universal across the Eurasian Arctic (Polishchuk et al. 2020); loss of lakes in discontinuous/sporadic permafrost zones may be compensated by increases in the number of lakes in regions of continuous permafrost (Smith et al. 2005). Moreover, interannual fluctuations may well mask the net changes in surface hydrology (Karlsson et al. 2014). Although remote sensing-based lake dynamics are not available for the BZT, by analogy with western Siberia and Alaska, which are also located in the discontinuous permafrost zone, one may expect that permafrost thaw in NE Europe will reduce the total area of lakes. Such lake drainage occurs via soil disturbance and increase of hydrological connectivity due to permafrost thaw and active layer thickness increase; phenomena observed both in Alaska (Hinkel et al. 2003;Riordan et al. 2006) and western Siberia (Kirpotin et al. 2008). Massive drainage of large lakes may have an unexpected impact on C balance, because air-exposed sediments have high CO 2 effluxes that greatly exceed those of surface waters (Martinsen et al. 2019). This highlights the need for systematic assessment of C emissions in large lakes including freshly drained lake basins.

Conclusions
This multiseasonal study of thermokarst lakes in the permafrost peatlands of NE European tundra (BZT) demonstrated the highest DOC and CO 2 concentrations and CO 2 emissions in depressions and small thaw ponds (< 1000 m 2 water surface area). However, CH 4 concentrations and fluxes to atmosphere exhibited a local maximum in ponds of 10-1000 m 2 surface area. The CO 2 and CH 4 concentrations and fluxes to atmosphere were highest in spring and the lowest in autumn. External spatial and temporal factors exhibited the following order of control on GHG parameters: lake size > season > year. Overall, identification of these factors allowed straightforward upscaling of field results obtained from several test sites in the BZT to the full territory of permafrost-affected NE European tundra zone. The estimated lake water surface emission flux to atmosphere from the BZT territory is two times higher than the lateral export of C (DOC + DIC) documented for the Pechora River to the Arctic Ocean. However, due to temporal and spatial variations of C fluxes, uncertainties on area GHG emissions are at least one order of magnitude for small thaw ponds and a factor of 3-5 for thermokarst lakes. Possible decrease in lake water area of Northern Eurasia, due to lake drainage induced by permafrost thaw, may decrease overall C emissions from the surface of lentic waters in this territory.