Minor contribution of small thaw ponds to the pools of carbon and methane in the inland waters of the permafrost-affected part of the Western Siberian Lowland

Despite the potential importance of small (< 1000 m2) thaw ponds and thermokarst lakes in greenhouse gas (GHG) emissions from inland waters of high latitude and boreal regions, these features have not been fully inventoried and the volume of GHG and carbon in thermokarst lakes remains poorly constrained. This is especially true for the vast Western Siberia Lowland (WSL) which is subject to strong thermokarst activity. We assessed the number of thermokarst lakes and their size distribution for the permafrost-affected WSL territory based on a combination of medium-resolution Landsat-8 images and high-resolution Kanopus-V scenes on 78 test sites across the WSL in a wide range of lake sizes (from 20 to 2 × 108 m2). The results were in fair agreement with other published data for world lakes including those in circum-polar regions. Based on available measurements of CH4, CO2, and dissolved organic carbon (DOC) in thermokarst lakes and thaw ponds of the permafrost-affected part of the WSL, we found an inverse relationship between lake size and concentration, with concentrations of GHGs and DOC being highest in small thaw ponds. However, since these small ponds represent only a tiny fraction of the landscape (i.e. ~1.5% of the total lake area), their contribution to the total pool of GHG and DOC in inland lentic water of the permafrost-affected part of the WSL is less than 2%. As such, despite high concentrations of DOC and GHG in small ponds, their role in overall C storage can be negated. Ongoing lake drainage due to climate warming and permafrost thaw in the WSL may lead to a decrease in GHG emission potential from inland waters and DOC release from lakes to rivers.


Introduction
Aside from rivers and streams, lakes are also an important source of greenhouse gas (GHG) emission to the atmosphere (Cole et al 1994, Tranvik et al 2009, Kosten et al 2010, Raymond et al 2013. Recent field and modeling studies showed that methane emissions from Arctic thermokarst lakes are significant and could increase by two to four times due to global warming Zhuang 2015a, 2015b). This may be especially true for the vast Western Siberian Lowland (WSL) that represents a hot spot of CO 2 (Repo et al 2007) and CH 4 (Sabrekov et al 2014(Sabrekov et al , 2017 emissions from inland waters and covers an area of more than 2 million km 2 . A sizeable part of the WSL is represented by frozen peatlands, subject to active thermokarst development. The thermokarst activity on this lowland produces numerous lakes and ponds that have not been fully inventoried or measured for surface area. These small high-latitude lakes and thaw ponds may turn out to be very important sources of GHG to the atmosphere (Hamilton et al 1994, Riera et al 1999, Walter et al 2006, Laurion et al 2010, Rautio et al 2011, Marushchak et al 2013, Langer et al 2015, Wik et al 2016. Although very small lakes and ponds (< 1000 m 2 size) represent only 10% of overall lakes and reservoirs in terms of surface area, they can emit more than 15% of CO 2 and 40% of CH 4 (Holgerson and Raymond 2016). These estimations are based on the power law of lake size distribution using map-based information and satellite imagery (Lehner andDoll 2004, Downing andPrairie 2006) or geostatistical models (Messager et al 2016). McDonald et al (2012) were the first to challenge the Downing and Prairie (2006) global estimate for lakes based on the Pareto (power law) distribution. Later, it was argued that the Pareto model overestimates the number and areal coverage of lakes smaller than 10 000-100 000 m 2 (Seekell et al 2013, Verpoorter et al 2014, Paltan et al 2015, Cael and Seekell 2016, Polishchuk et al 2017. As a result, the GHG emission flux from these lakes may also be overestimated. The main cause for this is insufficient knowledge of the size-and numberdistribution of small (< 5000-10 000 m 2 ) ponds and lakes. Note that while several papers from the previous decade highlighted large CH 4 emissions from discrete sites (e.g. the work of Walter et al 2006, Walter Anthony et al 2014, newer studies have suggested that the overall impact is rather small based on the global landscape context. For example, Stackpoole et al (2017) concluded that methane emissions from lakes across Alaska represented a small fraction of the annual aquatic C cycle and that small lakes play a minor role in the GHG budget as they do not cover sufficiently large areas.
Over past decades, climate warming has brought about a sizeable acceleration of thermokarst processes in permafrost regions which has led to an increase in the number of small thaw ponds and thermokarst lakes, as well as concentrations and emissions of methane and CO 2 in these lakes (Karlsson et al 2012, Shirokova et al 2013, Walter et al 2007, Walter Anthony et al 2014, Sepulveda-Jauregui et al 2015. The two major scenarios of thermokarst lake evolution under climate warming and permafrost thaw in western Siberia are (1) drainage of large thermokarst lakes into hydrological networks, which is especially pronounced in discontinuous permafrost zones (Smith et al 2005, Riordan et al 2006, Polishchuk et al 2014 and (2) formation of new depressions, subsidences, and small thaw ponds (< 100-1000 m 2 ), which is evidenced across all permafrost zones of this region (Pokrovsky et al 2011, Shirokova et al 2013, Bryksina and Polishchuk 2015. The large lake drainage scenario was first suggested from a comparison of medium-resolution (MR) satellite scenes between 1950 and 2010 in thermokarst regions of western Siberia (Smith et al 2005, Riordan et al 2006, Polishchuk et al 2014. The evaluation of CH 4 and CO 2 reservoirs and fluxes primarily requires knowledge of lake and pond water surface area (Holgerson and Raymond 2016, Verpoorter et al 2014, Downing and Prairie 2006, Abnizova et al 2012, Polishchuk et al 2015b. Although remote sensing analysis of thermokarst lakes in Siberian and other northern regions confirms the power law of lake size distribution (Karlsson et al 2014, Grosse et al 2008, Polishchuk et al 2015a, 2016b, Viktorov et al 2017, the majority of these studies use MR (30 m) Landsat images which do not allow the identification of small lakes and ponds. On the other hand, highresolution (HR) characterization of lake fraction in the ridge-hollow-lake wetland landscapes (i.e. Terentieva et al 2016) does not investigate the detailed small ponds size distribution for permafrost-affected WSL territory. The most recent circum-Arctic Permafrost Region Pond and Lake database (Muster et al 2017) deals with a small part of the WSL territory and does not allow quantitative estimation of lake and pond areal coverage and water volume in this region.
It is known that small thaw ponds exhibit concentrations of CH 4 , CO 2 , and DOC that greatly exceed those in large lakes (Bastviken et al 2004, Kortelainen et al 2006, Juutinen et al 2009, Abnizova et al 2012, Kankaala et al 2013, Shirokova et al 2013, Holgerson 2015, Pelletier et al 2015a. In order to account for their contribution to overall carbon balance, HR (1-10 m) images have to be combined with MR Landsat data to take into account all lakes, from meter sized to tens of kilometers sized, of a large studied territory (e.g. Polishchuk et al 2016aPolishchuk et al , 2016b However, the methodology of such integration is not sufficiently developed and does not allow quantification of overall lake numbers and surface areas and GHG and C concentrations and fluxes in inland waters. To fill in this gap, the primary goal of this study was to develop an approach for integrating HR and MR image data into the assessment of lake count and surface coverage for the full extent of lake sizes encountered in western Siberia. The second goal of this study was to evaluate the total amount of CH 4 , CO 2 and carbon in thermokarst lakes of western Siberia using the lake and pond inventory and available data on C concentration in typical thermokarst lakes of the region. In continuous and discontinuous permafrost zones of western Siberia, the smallest thaw ponds (100-1000 m 2 ) and depressions (1-100 m 2 ) exhibit an order of magnitude higher concentration of CO 2 , up to two orders of magnitude higher methane concentration, and a factor of three to ten higher concentrations of DOC and metals (Pokrovsky et al 2011, 2013, Shirokova et al 2013. As such, even with their small contribution to the total lake surface area, these small bodies of water may display carbon storage and GHG flux to the atmosphere comparable to those of large and Landsat-8 images covered the full territory of permafrost-affected WSL (1.05 million km 2 ). medium lakes. We anticipate that the knowledge of GHG and DOC distribution in thermokarst water bodies of the permafrost-affected part of WSL territory combined with known trends of lake and thaw pond area evolution over past decades should allow foresight as to the potential for CH 4 and CO 2 emission into the atmosphere and for the lateral export of DOC to rivers from the lentic aquatic ecosystems of western Siberia.

Thermokarst lakes and ponds of western Siberia
The WSL is comprised of taiga, forest-tundra, and tundra biomes. Permafrost is abundant north of 60 • N, changing from sporadic to discontinuous and continuous northward (figure 1). The climate is humid semi-continental with mean annual temperatures ranging from −2.8 • C in the south of the cryolithozone (Surgut region) to −10.4 • C in the north (Gyda Peninsula). The annual precipitation ranges from 600 mm in the south to 350 mm in the north. Peat thickness ranges from 1 to 3 m and is not directly linked to the latitude (Kremenetski et al 2003). The bogs of the region were subjected to freezing during the Subboreal period (∼ 4500 y.a). After the increase in temperature and precipitation during the Subatlantic period (2500 y.a.), the thermokarst started (Ponomareva et al 2012, Pastukhov et al 2016. It produced thaw ponds and thermokarst lakes that were formed due to progressive thawing of frozen peat rich in ice (Pokrovsky et al 2011, Grosse et al 2016). On flat watershed divides, which dominate the territory (palsa peat bog and polygonal tundra), the lakes are subjected to a natural cycle of growth and draining which is well established in the WSL (Kirpotin et al 2011) but is also evidenced in other parts of the circum-Arctic zone (Hopkins 1949, French 1996. Because of the homogeneous nature of the surrounding substrate (Quaternary clays and sands of fluvio-glacial and lacustro-glacial origin underlying the surface peat deposits), the thermokarst lakes are generally round in shape. As such, all the lakes and thaw ponds inventoried in this work are considered as having a thermokarst origin. While satellite imagery cannot resolve the origin of small ponds, our terrestrial observations in addition to the very high number of small ponds suggest these ponds were formed after thermokarst subsidence rather than remaining in the basin following large lake drainage.
The range of thermokarst lake and thaw pond sizes in the permafrost-affected part of the WSL is from a few m 2 to hundreds of km 2 (Manasypov et al 2014). For convenience, we will distinguish between large lakes (>10 000 m 2 or 1 ha), small lakes (500-10 000 m 2 ), and very small thaw ponds (< 500 m 2 ). This conventional distinction is based on hydrochemical and gas regimes of lakes and ponds; there is an abrupt rise of DOC, CO 2 , and CH 4 in water bodies smaller than 500 m 2 (Pokrovsky et al , 2013, Shirokova et al 2013, Manasypov et al 2014, Polishchuk et al 2015b. This classification is also consistent with that of Holgerson and Raymond (2016) who considered lakes < 1000 m 2 as very small water bodies. Muster et al (2017) followed the distinction of Rautio et al (2011) and defined ponds as bodies of standing water with a surface area smaller than 10 000 m 2 .

Sampling and analyses
Water samples were collected from a PVC boat for large lakes or directly from the lake center for small (< 50 m diameter) water bodies (thaw ponds) from the middle of the water column. The sampling was performed during July and August from 2010(Pokrovsky et al 2013, Shirokova et al 2013, Manasypov et al 2014, Pavlova et al 2016, Loiko et al 2017 and in August 2015 and 2016 (this study). The typical depth of sampling was between 0.5 and 0.25 m depending on the lake depth (see section 3.3 below). A full list of collected lakes and relevant physicogeographical parameters together with the measured GHG and DOC is given in table S1 of the supplementary information available at stacks.iop.org/ERL/ 13/045002/mmedia. For GHG analyses, unfiltered water was sampled in 25 m L serum bottles that were closed without air bubbles using vinyl stoppers and aluminum caps and immediately poisoned by adding 0.2 ml of saturated HgCl 2 using a twoway needle system. In the laboratory, a headspace was created by displacing approximately 40% of water with N 2 (99.999%). Two 0.5 m L replicates of the equilibrated headspace were analyzed for their concentrations of CH 4 and CO 2 , using a Bruker GC-456 gas chromatograph equipped with flame ionization and thermal conductivity detectors. After every ten samples, a calibration of the detectors was performed using air liquid gas standards (i.e. CH 4 = 145 ppmv and CO 2 = 3000 ppmv). Duplicate injections of the samples showed that results were reproducible within ±5%. The specific gas solubility for CH 4 and CO 2 (Weiss 1974) were used in the calculation of the total CH 4 and CO 2 content in the vials, which was then recalculated to mol/L of the initial waters. For DOC analyses, collected waters were filtered onsite in pre-washed 30 mL polypropylene Nalgene Ⓡ bottles through single-use Minisart filter units (0.45 m pore size, Sartorius, acetate cellulose filter). The first 20 ml of filtrate were discarded. Blanks were used to control the level of contamination induced by sampling and filtration. DOC blanks of filtrate never exceeded 0.1 mg L −1 which is quite low for the organic-rich pore waters sampled in this study (i.e. 10-100 mg L −1 DOC). DOC was analyzed using a Carbon Total Analyzer (Shimadzu TOC VSCN) with an uncertainty better than 3% (see Prokushkin et al 2011 for methodology).

Source and description of satellite imagery
Satellite imagery from the Landsat-8 Operational Land Imager (30 m MR) [available at the USGS Global Visualization Viewer: http://glovis.usgs.gov] together with HR (2 m)) Kanopus-V scenes (available at the Roscosmos website: http://eng.ntsomz.ru/) were used to map the lake distribution over the full territory of permafrost-affected WSL (figure 1). We used MR Landsat-8 images collected at the end of July-beginning of August 2013-2014. HR Kanopus-V images were collected during summers 2013-2015. This free-ice period corresponds to the minimal coverage of the territory by lakes and minimal seasonal variation of lake water levels. A combination of images from the two years was necessary because a single year's observation could not provide full coverage of the territory.
For the entire permafrost-affected WSL territory of 1.05 million km 2 , we used a superposition of 134 MR images. Images were treated using standard tools of ArcGIS 10.3 software (Kennedy et al 2013). For automatic identification of lakes we used the Fmask algorithm (Zhu et al 2015). It allows resolving the lakes even with some cloud coverage and uses several spectral ranges of Landsat-8 using a set of empirical spectral indexes. As input data we used the intensity of the upper atmosphere reflectance and the brightness temperature presented as nominal units of pixel brightness. First, for the mosaic of Landsat-8 imagery, cloud masks were defined for individual images. Then the cloud masks were removed from the images and replaced by cloudfree fragments taken from the adjacent images during the same period of observation. The minimal pixel on the MR image was 30 × 30 m 2 , thus limiting the minimal lake size of 0.5 ha or 5000 m 2 . This corresponds to six pixels of 30 m × 30 m (900 m 2 ), which is sufficient for reliably distinguishing lakes from background noise.
A two-stage procedure was developed to separate thermokarst lakes and thaw ponds from other water objects. First, we masked coastal zones and rivers, based on vector layers of OpenStreetMap and State Water Register (http://textual.ru/gvr/). Second, we superposed this mask on a raster layer of recognized water objects forming a layer of lakes. The uncertainty of lake surface area measurement using Landsat-8 for lakes with sizes larger than 10 5 m 2 is 3.5% and for lakes larger than 2 × 10 4 m 2 is 7% Polishchuk 2013, Kornienko 2014).
It is important to note that thermokarst lakes and thaw ponds of western Siberia are different from glacial lakes of other boreal and subarctic regions. The latter are often developed on moraine and crystalline rocks and located within non-isometric elongated basins that were formed by glacier movement. In contrast, thermokarst lakes in palsa bogs of western Siberia occurring within frozen 1-3 m peat deposits are rarely oval and are usually round and isolated (Kirpotin et al 2009, Shirokova et al 2013, Pokrovsky et al 2014. According to our field observations and detailed lake mapping over 950 km latitudinal gradient of western Siberia, the fraction of lakes with an irregular shape is less than 10% (Pokrovsky et al , Manasypov et al 2014, Polishchuk et al 2017. For this reason, the influence of the irregularity of the lake border on the determination of the total lake area can be neglected. This is also confirmed by a study of the tortuosity of thermokarst lake border lines demonstrating that the deviation from circumference in lakes of sporadic, discontinuous and continuous permafrost zones produces less than 5% uncertainty on total lake area evaluation by Landsat-8 (Polishchuk and Polishchuk 2012).
HR Kanopus-V images were used for the quantification of small lakes and very small thaw ponds. For the reliable identification of lakes, at least five pixels are necessary. The minimal pond size inferred from the Kanopus-V images was 20 m 2 as the spatial resolution of 2 m yielded a pixel size of 4 m 2 allowing the required five pixels to be distinguished on the 20 m 2 area. The scenes were processed using ArcGIS 10.3 within 78 test sites (TSs), which were evenly distributed over sporadic, discontinuous, and continuous permafrost zones of the WSL (figure 1). The characteristics of these TSs are given in table 1. The number of lakes and ponds on each TS ranged from 97-2583. All TSs had similar sizes (30 km 2 ) and covered approximately 2000 km 2 which represents 0.002% of the overall permafrost-affected territory of the WSL.
We used a binary classification algorithm in ArcGIS 10.3 analogous to the 'density slice' method used by Muster et al (2017). This algorithm is based on the visual detection of certain threshold values of spectral brightness that allow us to distinguish two object classes: 'water' and 'not water'. The threshold value had to be selected individually for each TS due to differences in illumination, acquisition geometry, and different sensor spectroradiometry. The panchromatic imagery exhibited a strong contrast between water surface and the surrounding vegetation which allowed an adequate selection of threshold values needed for classification. The submerged moss and macrophytes were not detected and non-submerged moss and plants along the lake coasts had no influence on lake water surface area.

Building the integral diagrams of lake number and area distribution in a wide range of lake sizes
To build the diagrams of lake distribution in a wide range of lake sizes, we used 21 partial ranges of lake size taken within a logarithmic scale, i.e. 20-50 m 2 , 50-100 m 2 , 100-200 m 2 , 200-500 m 2 and so on until 200 000 000 m 2 . To combine MR and HR images we developed a three-step procedure. First, we constructed diagrams for the number and area distribution of all thermokarst lakes detectable on the MR images. Then, we built similar diagrams for small thermokarst lakes and very small thaw ponds using HR images of the 78 TSs across the permafrost-affected part of WSL. At this stage, the total number and area of lakes located within the TSs were extrapolated to the full territory of the permafrost-affected part of WSL, as described below.
Considering the values of lake numbers and areas in the lake size intervals suitable for Landsat-8 mapping as the most reliable, we required that within the lake size intervals where Landsat-8 and Kanopus-V images overlap, the lake number and area determined from HR scenes on TSs were to be equal to values extrapolated to the full studied territory. We believe that such extrapolation is valid because the TSs are evenly distributed across the permafrost-affected part of the WSL. The extrapolated (calculated) values for the total number ( Ncj ) and total lake area ( Scj ) in each jth interval of HR images-based histogram were obtained from: where N and S are the total number and area of lakes, respectively, that are experimentally determined based on HR scenes of all TSs, j is the number of intervals of the histogram (figure 2), and K is the empirical extrapolation coefficient which was averaged over all lake sizes and therefore independent in terms of size range. Finally, we combined primary MR and HR histograms into integral histograms of lake numbers and overall areas. These histograms included lakes and small and very small thaw ponds spread over seven orders of magnitude in size.   The uncertainty of extrapolating data from the 78 TSs to the whole territory of the permafrost-affected part of the WSL was evaluated via a comparison between data of HR and MR in the lake size interval where the HR and MR scenes overlap (figure 2). Using the number of lakes in the seven intervals of superposition, we found the relative difference in lake numbers using equation 3: where m and n represent the number of lakes in the ith overlapped interval, determined using HR and MR images, respectively. In this equation, n is chosen as a reference number given that HR images were treated for the full territory of the permafrost-affected part of WSL. The average value of the seven overlapped intervals was equal to 15% which is acceptable for lake size distribution estimation.

Results and discussion
3.1. Lake number and area in a wide range of lake sizes, from 20 to 2 × 10 8 m 2 . The MR images allowed resolving lakes of the minimal size of 5000 m 2 corresponding to six pixels which is necessary to distinguish the lakes from the surrounding landscape. The total number of lakes ranging from 5000 m 2 to 2 × 10 8 m 2 was 727 700 with a total water area of 60000 km 2 . The histograms of lake numbers and lake area distribution over discrete size ranges are given in the blue columns in figures 2(a) and (b), respectively. The HR (Kanopus-V) scenes were used to quantify the number and areas for small thaw ponds. The empirical extrapolation coefficient (equations 1 and 2) was calculated as an average of the ratio of lake number and area determined by MR scenes for the entire permafrost-affected part of the WSL territory to the total lake number and area from all HR images (at all TSs) in each lake size interval where the HR and MR histograms overlap. The histograms have a seven-interval wide overlap range from 5000 to 10 6 m 2 (figure 2). The extrapolation coefficient values calculated from the data on the number and area of lakes were equal to 85 and 86 respectively. Thus, K = 85.5 was defined as an arithmetic mean of these values. The calculated numbers and total areas of lakes in each interval of lake size using equations 1 and 2 yielded the diagrams of small lakes and very small thaw ponds distribution on the territory, shown as yellow columns in figures 2(a) and (b) for the lake size range of 20-10 6 m 2 . The integration (stitching) of HR and MR-based lake diagrams can be performed at the lake size of 20 000 m 2 with a typical uncertainty on lake determination at 7% (Bryksina and Polishchuk 2013). This point is shown in figure 2 as a red vertical line. This last step of integrated histograms yielded a lake size and area distribution between 20 and 2 × 10 8 m 2 so that lakes smaller than 2 × 10 4 m 2 were processed using HR scenes, and the lakes >20 000 m 2 were treated from MR images.
The total number of lakes was equal to 8.7 million and their total area was 6.4 × 10 10 m 2 . Very small thaw ponds (< 500 m 2 ) provided 44.2% of total lake numbers, but their area contributed to less than 0.1% of total area coverage. Big lakes (>10 000 m 2 ) which accounted for only 9.1% of lake numbers provided 91.1% of all lake areas.

Comparison of lake inventory in western Siberia with other studies
The number of lakes and lake area-lake size distribution diagrams obtained in this work for the permafrostaffected part of WSL are in fairly good agreement with the inventory of all lakes of the circum-Arctic permafrost region (Muster et al 2017); deviation between the two estimations does not exceed several percent (figure 3). A comparison with the world lake database of Cael and Seekell (2016) is shown in figure 4 where empirical distribution lakes in Sweden and throughout the world are shown. It can be seen that for lakes smaller than 10 4 -10 5 m 2 , at both the planetary and regional (WSL, Sweden) level, there is a sizeable difference between the results of direct measurement and calculation from the power law of lake distribution. Specifically, the power law extrapolation strongly overestimates the small lake number compared to empirical data. As a result, the use of power law can also lead to the overestimation of water and solute amounts in small lakes and very small thaw ponds. Indeed, a plot of the relative contribution of lakes of different size to the total area of lakes demonstrates sizeable differences in the fraction of small lakes calculated by Holgerson and Raymond (2016) for the whole planet and our direct measurements of very small thaw ponds in the permafrost-affected part of the WSL territory (figure 5). The two dependences have a generally similar shape and lakes having an  area of 10 5 m 2 provide maximal contribution to the total areal coverage. However, significant differences are pronounced for small water bodies (< 1000 m 2 ); the relative contribution of lakes evaluated by Holgerson and Raymond (2016) for the world is several times higher than that directly measured in western Siberia. For very small ponds of 100 m 2 , this difference is a factor of 6.5. As such, the estimation of total concentrations and emissions of CH 4 and CO 2 from small lakes on a worldwide scale, performed based on Monte Carlo numerical modeling using the Pareto power law (Holgerson and Raymond 2016), cannot be applied to very small thaw ponds of western Siberia. In the WSL direct remote sensing using HR images from 78 TSs allows for a more rigorous lake and pond inventory and thus provides a solid basis for the evaluation of GHG and carbon pools.
It is important to note that the use of Muster et al's (2017) Pond and Lake circum-Arctic database for the WSL territory is not possible. In the work by Polishchuk et al (2017), only two sites (Central Yamal and the Surgut region) were inventoried where lake coverage (15%-20%) was much higher than that of the total WSL permafrost-affected zone (5%-10%). In the southern part of the WSL cryolithozone, within the peatland 'Surgutskoe Polesye' of the northern taiga biome (61.5-63 • N), the proportion of open water ecosystems is around 19% based on HR (1-3 m) Google Earth inventory (Terentieva et al 2016). However, the average limnicity of the permafrost-affected part of the WSL assessed in this study is consistent with values reported in a similar context for two peat plateau areas of northeast European Russia (7%-13.6%, Sjöberg et al 2013) and the recent

Calculations of CH 4 , CO 2 , and DOC pools in small thermokarst lakes and very small thaw ponds of western Siberia
To estimate the amount of GHG and dissolved carbon in lakes of the permafrost-affected part of the WSL, published and new collected concentrations of DOC, CO 2 , and CH 4 in 146 thermokarst lakes north of the permafrost boundary were used (Pokrovsky et al 2013, Shirokova et al 2013, Manasypov et al 2014, Pavlova et al 2016, Loiko et al 2017 of the supplementary information). The CO 2 and CH 4 concentrations in small (< 500 m 2 ) thaw ponds range from 29-778 mol CO 2 L −1 and 0.02-47 mol CH 4 L −1 with average values equal to 172 and 9 mol L −1 , respectively. This is significantly (p < 0.05) higher than the range of CO 2 and CH 4 in thermokarst lakes >500 m 2 ; 29-332 with an average of 92 mol L −1 and 0.02-8.6 with an average of 0.98 mol L −1 , respectively, as illustrated in figures 6(a) and (b).
The DOC concentration in ponds is a factor of 1.5-2.0 higher than that in lakes (figure 6 (c)). Note that the ratio of CO 2 to CH 4 decreases with the decrease in water body area (figure 7) and that this is consistent with the global world trend (Holgerson and Raymond 2016). The concentrations of CO 2 , CH 4 , and DOC were not found to be significantly (p < 0.05) dependent on latitude across the WSL permafrost-affected zone surface waters, soil solutions, and suprapermafrost waters (Manasypov et al 2014, Raudina et al 2017). For this reason, as a first approximation, we will Figure 7. A decrease of CO 2 to CH 4 concentration with decrease of water body area. The relative role of CH 4 is one to two orders of magnitude higher in very small thaw ponds compared to thermokarst lakes. Table 2. Distribution of CH 4 , CO 2 , and DOC total amounts among lakes and ponds of different size across the WSL (1.05 million km 2 ). Uncertainty is estimated between 15% and 30%.
Lake area range, m 2 CH 4 , t CO 2 , t DOC, t Lake area, km 2 neglect any latitudinal trend in C concentrations in ponds and lakes and consider only lake area as the main controlling parameter of GHG and DOC concentration. It is important to note that the thermokarst lakes and thaw ponds of the WSL are extremely shallow (large lakes are 0.5-1.5 m deep and small thaw ponds are < 0.5 m deep) and are not stratified in temperature. We sampled the surface and bottom layer and have not found any measurable chemical stratification in DOC and CO 2 or O 2. Lakes were fully oxygenated over 0.5-1.5 m depth (Shirokova et al 2013, Manasypov et al 2015. The latter allowed us to assume the constant profile of CH 4 concentration in the water column. The depth of WSL thaw ponds and thermokarst lakes depends on the size of the water body (S) and can be approximated by the polynomial equation 4: Depth = 0.00007 × 3 − 0.0049 × 2 + 0.1141 × − 0.0124 ( 2 = 0.552) (4) which is based on field measurements of thaw ponds and lakes of 0.2 to 2 × 10 7 m 2 (Polishchuk et al 2017). This information provides a solid background for the estimation of GHG and DOC reservoirs in the lakes of western Siberia. The masses of CH 4 , CO 2 , and DOC in WSL lakes and ponds for seven lake size intervals are listed Table 3. The distribution of fraction of CH 4 , CO 2 , and DOC amounts over seven ranges of lake size in the WSL.
Lake area range, m 2 CH 4 , % CO 2 , % DOC, % Lake area, % The maximal contribution to the CH 4 pool in the permafrost-affected part of the WSL is provided by lakes between 10 5 and 10 6 m 2 surface area (figure 8). Lakes larger than 10 4 m 2 provide >90% of the total CH 4 amount whereas the share from small thaw ponds (< 1000 m 2 ) is only 2%. According to these data, small thaw ponds contribute less to overall CO 2 and DOC levels than to CH 4 level.

Global consequences of climate warming on pond and lake water GHG emission and DOC lateral export revisited
Remote sensing studies of the permafrost zone of western Siberia demonstrated that the number of newly formed small thermokarst lakes (5000-20 000 m 2 ) over the past three decades exceeds, by a factor of 20, the number of large lakes which tend to disappear during the same period (Bryksina and Polishchuk 2015). This increase in the number of small and very small thaw ponds over the past few decades was demonstrated through an examination of HR (2 m) Kanopus scenes. Assuming the conservative behavior of carbon in WSL thermokarst lakes and ponds and ignoring the short and long-term dynamics of carbon cycling, the draining of large thermokarst lakes and the appearance of small thaw ponds suggest a sizeable decrease in CO 2 , CH 4 , and DOC pools in thermokarst water bodies of the WSL. Thus, a 10%-20% decrease in the area of large (>20,000 m 2 ) lakes suggested in earlier studies (Smith et al 2005, Riordan et al 2006, Kravtzova and Bystrova 2009, Shiklomanov et al 2013 will lead to a corresponding decrease in CH 4 , CO 2 , and DOC concentrations in ponds and lakes of the permafrost-affected part of the WSL. At the same time, even a two to threefold increase in the number of small (< 10 000 m 2 ) lakes and a tenfold increase in very small (< 500 m 2 ) thaw pond number and area will not produce any sizeable (i.e. >2%) increase in GHG and DOC pools in WSL ponds and lakes. Taken together, the permafrost thaw and climate warming can only decrease the pools of GHG and organic carbon in thermokarst water bodies of western Siberia. Assuming a solely diffusive flux, which is governed by the GHG concentration gradient, this will decrease the CO 2 and CH 4 emission potential to the atmosphere. In addition, the lateral export of DOC from the lakes and ponds to the hydrological network will also decrease. In this regard, our study corroborates the previous results of van Huissteden et al (2011) andKessler et al (2012) that thermokarst lake drainage leads to lower GHG emissions.
There are other indications that climate warming and the permafrost thaw may not positively impact the GHG emissions from lentic aquatic ecosystems. For example, the landscape-scale model that included the entire life cycle of thaw lakes in northern Siberia (initiation, expansion, drainage, and eventual re-initiation) demonstrated that methane emissions from thaw lakes in eastern Siberia (Yakutia) are an order of magnitude less alarming than previously suggested (van Huissteden et al 2011). Moreover, the use of the Integrated Global System Model demonstrated that, on the whole the Arctic and boreal permafrost-scale increases in atmospheric CH 4 and its radiative forcing (which result from the fact that thawed, inundated emission sources under saturated (anaerobic) soils) are small, particularly when weighed against human emissions (Gao et al 2013). Considering the results presented in this work on aerobic (open lake and pond) areas, we believe the overall potential threat of permafrost thaw on GHG emission from Arctic and high latitude boreal permafrost-bearing regions may be overestimated. We do note, however, that detailed mapping of western Siberian wetland complexes (i.e. Peregon et al 2008, 2009, Terentieva et al 2016 is necessary to integrate available DOC and GHG concentrations and fluxes into a global ecosystem-specific model for this territory. Here, mapping of fine-scale heterogeneity of WSL's wetland landscapes, which is currently possible only for CH 4 emission (i.e. Bohn et al 2007, Glagolev et al 2011, will constitute the main challenge for both the DOC and CO 2 inventories. In this regard, quantifying the amount of water and carbon in waterlogged hollows-especially those abundant in the southern part of the permafrost-affected WSL territory It is further important to note that, despite a net lake area loss due to drainage of large lakes, thermokarst expansion can offset the decrease in GHG emissions and in some regions (i.e. Yedoma Ice Complex) small lake area increase can increase GHG emissions despite a net loss of large lake area (Walter Anthony et al 2016). Moreover, the ultimate response of GHG to the thermokarst lake cycle may be a function of the permafrost type (continuous versus discontinuous), the permafrost carbon stock in terrestrial landscape surrounding lakes (Yedoma, peatland or non-Yedoma mineral soils), the rate of thermokarst expansion, and/or the balance of gross lake area loss versus gain.
Another uncertainty involves the relative role of GHG ebullition versus diffusion which cannot be assessed solely from CO 2 and CH 4 concentrations. Although ebullition is known to dominate over emission in many boreal and Arctic ponds and lakes including those in Siberia (i.e. Walter et al 2006, Anthony et al 2010 there is no relationship between diffusive emissions and ebullition (Wik et al 2016, Sepulveda-Jauregui et al 2015. In lakes of permafrost-free regions of the WSL, ebullition strongly dominated over diffusion in the middle taiga zone but was comparable to diffusion in the southern taiga zone (Sabrekov et al 2017). As such, the response of GHG emission to frozen peat thaw and permafrost boundary change in the WSL may be more complex than predicted exclusively from GHG total concentrations in lentic waters.

Conclusions
A combination of MR (Landsat-8, 30 m) and HR (Kanopus-V, 2 m) satellite imagery of the WSL allowed the quantification of the number and areal coverage of lakes and ponds across the 1.05 million km 2 permafrost-affected territory. The total lake number was equal to 8.7 million with a total area of 6.4 × 10 10 m 2 . Small and very small thaw ponds (< 1000 m 2 ) accounted for 44% of the total lake number however their area contributed to less than 0.1% of the total areal coverage. The maximal contribution to the overall CH 4 pool is provided by lakes having 10 5 to 10 6 m 2 surface area. Thus, despite concentrations of DOC and GHG in small thaw ponds that are one to two orders of magnitude higher than those in thermokarst lakes, the share from small thaw ponds is less than 1.5% of the overall CH 4 pool and less than 1% of the CO 2 and DOC pools in lentic waters of the permafrost-affected WSL. A possible consequence of climate warming is that large thermokarst lakes in the permafrost-affected part of the WSL will drain and small thaw ponds will appear. This will clearly decrease, and not increase, the overall CH 4 and CO 2 emission potential from the WSL inland waters and will also decrease possible feeding of streams by DOC from lakes. However, the overall scenario of GHG emissions from inland waters due to permafrost thaw in the permafrost-affected part of the WSL will depend on the evolution of the landscape and changes in the ebullition contribution to total emissions.

Acknowledgments
This work was supported by the RNF (RSCF) grant No 15-17-10009 'Evolution of thermokarst lake ecosystems in the context of climate change' (analyses, interpretation, 50%) and RFBR grants No. 15-45-00075 and 17-55-16008 C NRS-a. We are grateful to three anonymous reviewers for their very constructive remarks. Chris Benker is thanked for the English editing. Organic and organo-mineral colloids of discontinuous permafrost zone Geochim. Cosmochim.Acta 188 1-20