Size Distribution, Surface Coverage, Water, Carbon, and Metal Storage of Thermokarst Lakes in the Permafrost Zone of the Western Siberia Lowland

Despite the importance of thermokarst (thaw) lakes of the subarctic zone in regulating greenhouse gas exchange with the atmosphere and the flux of metal pollutants and micro-nutrients to the ocean, the inventory of lake distribution and stock of solutes for the permafrost-affected zone are not available. We quantified the abundance of thermokarst lakes in the continuous, discontinuous, and sporadic permafrost zones of the western Siberian Lowland (WSL) using Landsat-8 scenes collected over the summers of 2013 and 2014. In a territory of 105 million ha, the total number of lakes >0.5 ha is 727,700, with a total surface area of 5.97 million ha, yielding an average lake coverage of 5.69% of the territory. Small lakes (0.5–1.0 ha) constitute about one third of the total number of lakes in the permafrost-bearing zone of WSL, yet their surface area does not exceed 2.9% of the total area of lakes in WSL. The latitudinal pattern of lake number and surface coverage follows the local topography and dominant landscape zones. The role of thermokarst lakes in dissolved organic carbon (DOC) and most trace element storage in the territory of WSL is non-negligible compared to that of rivers. The annual lake storage across the of their annual delivery by WSL rivers to the Arctic Ocean from the same territory. However, given that the concentrations of DOC and metals in the smallest lakes (<0.5 ha) are much higher than those in the medium and large lakes, the contribution of small lakes to the overall carbon and metal budget may be comparable to, or greater than, their contribution to the water storage. As such, observations at high spatial resolution (<0.5 ha) are needed to constrain the reservoirs and the mobility of carbon and metals in aquatic systems. To upscale the DOC and metal storage in lakes of the whole subarctic, the remote sensing should be coupled with hydrochemical measurements in aquatic systems of boreal plains.


Introduction
The quantification of the abundance, size distribution, and water storage of lakes and reservoirs has critical importance for the evaluation of carbon and nutrient storage, and the potential of greenhouse gas (GHG) exchange between the Earth's surface and the atmosphere.For these reasons, several detailed studies have documented the lake number and size distribution on the scale of our planet [1][2][3].Thus, the use of medium resolution Landsat-7 images has allowed the creation of a global database of lakes and water reservoirs, including all lakes larger than 0.2 ha or 45 m × 45 m [3,4].The total number of lakes was estimated as being 117 million, with an overall surface area of 5 × 10 6 km 2 (500 × 10 6 ha), which corresponds to 3.7% of the non-glaciated land area.The number of small lakes (0.2-1 ha) is around 90 million, whereas their overall area is equal to 25 × 10 6 ha, which is only 5% of the overall lake surface coverage.The small lakes exhibited a deviation from the general power dependence between the lake size and the lake number [4,5]; as a result, the extrapolation of power law to smaller lakes may overestimate the lake number [3].At the same time, the small lakes subjected to full freezing in winter or evaporation in summer, with a short residence time of water, play a crucial role in the integration of the carbon and other elements transported from the watershed [6][7][8], which is particularly important in high latitude regions, which are the most vulnerable to climate change [9,10].
The Arctic and subarctic permafrost-bearing regions exhibit the maximal changes in the terrestrial freshwater budget, although the hydrological responses to environmental changes strongly differ across the boreal and subarctic regions of the subarctic [11].In particular, in the tundra, continuous permafrost development strongly influences water fluxes and storage, whereas in boreal plains, slow surface and subsurface water movement produces extensive wetlands [11].Once the permafrost becomes discontinuous to sporadic in the south, this allows significant groundwater feeding of rivers [12] and, presumably, lakes [13].
In this regard, the boreal and tundra plains are extremely important for a lake inventory study because of the high coverage of the watershed area by these lakes (up to 70% in some western Siberian river watersheds [14,15]), and fast temporal dynamics of thermokarst lake landscapes, reflecting on-going climate change in their watersheds [16][17][18][19][20].The latter brings about a shorter residence time of lakes, whose size changes, especially at southern latitudes, due to the disappearance of sporadic and isolated permafrost.It is also worth noting the primary role of lakes in controlling greenhouse gas exchange with the atmosphere, both in permafrost-free [21,22] and permafrost-bearing regions [9,[23][24][25][26].
Several studies of thermokarst lakes in Alaska, Yukon, Scandinavia, and Siberia, were focused on monitoring the change in the lake area over the past 30-40 years, within relatively small regions [18,29,[41][42][43][44][45].Remote sensing studies of the permafrost zone of western Siberia demonstrated that the number of newly formed small thermokarst lakes (0.5-5 ha) over the past three decades exceeds, by a factor of 20, the number of large lakes which tend to disappear during the same period [46].Recently, the dynamics of the number and surface area of thermokarst lakes in the discontinuous permafrost zone of western Siberia, over the period of 1973 to 2009, has been studied within the watersheds of the Nadym and Pur rivers [47].According to these authors, the temporal evolution of large size (>10 ha) lakes, whose number constituted 78%-85% of all lakes, exhibited a variation within 10%.The size distribution of thermokarst lakes followed the power law, both in eastern [30] and western [47,48] Siberia.In particular, Polishchuk et al. [48] presented the results of the number of appearing and disappearing lakes in western Siberia between the 1970s and the present time, and the laws of statistical distribution of very small lakes (<0.5 ha) on several test sites of the western Siberian Lowland (WSL).In contrast, the present study encompasses a much larger territory of western Siberia and provides, for the first time, a full inventory of large to small lakes (>0.5 ha).
The main goal of the present study was to establish the law of the lake number and area distribution for the whole WSL territory and to bring together the hydrology and hydrochemistry using available data on lake depth and carbon and metal concentration in the lake water.Towards this goal, we classified and analyzed medium resolution Landsat-8 scenes, which provided complete coverage of the WSL (105 × 10 6 ha or 1.05 million km 2 ).The first specific objective of this study was to increase the resolution of the lake size to 0.5 ha for the whole territory of permafrost-affected WSL, in order to compare the results with the global database [3].Indeed, in view of the disproportionally high importance of small thermokarst water bodies relative to medium and large lakes in GHG emission and C storage [49][50][51], a rigorous quantification of the number and area of thermokarst lakes is very timely.The second objective of this work was to assess the water, carbon, and metal storage in thermokarst lakes.Recent progress in the quantification of depth, and area and lake size-averaged concentrations of major and trace elements in western Siberian thermokarst lakes [49,50,[52][53][54], allows a first-order evaluation of the water and element stocks in lakes and a characterization of the role of lakes in element storage, relative to rivers draining the same territory.

Study Area
The studied region is located within a tundra and forest-tundra zone of the northern part of the western Siberia lowland (1.05 × 10 6 km 2 ).In the northern part of the WSL, the sporadic, discontinuous, and continuous permafrost zones share 31.7%,29.1%, and 39.2% of the overall territory, respectively (Figure 1).The mean annual temperature (MAT) ranges from −0.5 • C in the permafrost-free region (Tomskaya region) to −9.5 • C in the north (Yamburg), and the annual precipitation ranges from 400 to 460 mm.For the period of the end of July-August in the central part of the studied zone (Novuy Urengoy), the average low daily temperature was 17.4 • C and 10.9 • C in 2013 and 2014, respectively.The average high temperature was 23.4 and 14.9, respectively.A detailed physico-geographical description, hydrology, lithology, and list of the soils can be found in earlier works [55,56] and in our recent limnological and pedological studies [50,53,54,57].The WSL has rather homogenous landscape conditions (palsa peat bogs, forest-tundra, and polygonal tundra), lithology (Pliocene sands and silts), soil cover (1-1.5 m thick peat, half of soil profile is frozen), and runoff (200-250 mm•year −1 ), across a large gradient of permafrost coverage [58][59][60].
The bioclimatic sub-zones of permafrost-affected WSL regions gradually change northward, from the northern taiga zone (38 × 10 6 ha) to forest tundra (13 × 10 6 ha), southern tundra (30 × 10 6 ha), and northern tundra (24 × 10 6 ha).A detailed GIS survey of the WSL allowed the quantification of the regional distribution of major wetland types and complexes [56].According to those authors, the sporadic permafrost zone, north of the Ob River (61 • -63 • N), is dominated by sphagnum bogs with pools and an open stand of trees with abundant forested (treed) shrubs-and moss-dominated mires.The discontinuous permafrost zone (63 • to 67 • N) of the forested tundra and southern tundra essentially comprises high palsa and flat palsa mires, mixed with palsa-hollow and pool-hollow patterned mires.Flat-palsa and hollow-pool flat-palsa bogs are also abundant in this region.Finally, further to the north, within the continuous permafrost zone of the southern and northern tundra (67 • -73 • N), the landscape is dominated by patterned (hollow and hollow-pool) flat-palsa bogs, polygonal mires combined with grass and moss-dominated mires [56].Thermokarst lakes are highly specific water bodies of the permafrost zone of the WSL: they have a shallow depth and are rarely connected to the hydrological network [61].Most lakes of the northern part of the WSL freeze solid in winter and have frozen sediments at their bottom throughout the year [54].The region is dominated by the presence of thermokarst lakes having a <100 ha surface area [42,62,63].The lakes are located within peat flat-mound bogs (ridge-hollow complex, palsas, and polygonal tundra); the bottom sediments of the lakes are dominated by peat detritus.An active thermokarst occurs due to the thawing of syngenetic and epigenetic segregation ice, ice wedges, and ice layers in the deep (>2 m) horizons, and is primarily due to the ice thawing of the active layer (<2 m).The thermokarst activity produces depressions, subsidences, and ponds, which are usually separated by flat mound peat bogs up to 2 m in height [64].The largest thermokarst lakes that are located within the peat bog are km-size, with a depth of 0.5-1.5 m [64,65].The overwhelming majority of lakes in all three permafrost-bearing zones of the WSL (sporadic, discontinuous, and continuous) have a thermokarst origin, i.e., thawing of frozen peat and clay surface horizons [42,53,61,66].As a result, most thermokarst lakes of the WSL exhibit quite a shallow depth, ≤1 m, in contrast to the deeper lakes of other Arctic regions, originating from surface disturbance, the melting of ground ice, and ice wedges (e.g., ref. [67]).

Remote Sensing Analysis
Satellite imagery from a Landsat-8 Operational Land Imager (OLI, with 30 m resolution), available at USGS Global Visualization Viewer [68], were used to map the lake distribution.We used medium-resolution Landsat-8 images collected at the end of July-August 2013-2014.This period corresponds to the minimal coverage of the territory by lakes and minimal seasonal variation of the lake water level.Besides, this is the period of total disappearance of ice coverage throughout all of the studied area.The sampling of the lake water and lake depth measurements was also performed during this period of the year.Note that, from the view point of optical remote sensing, there is no The lakes are located within peat flat-mound bogs (ridge-hollow complex, palsas, and polygonal tundra); the bottom sediments of the lakes are dominated by peat detritus.An active thermokarst occurs due to the thawing of syngenetic and epigenetic segregation ice, ice wedges, and ice layers in the deep (>2 m) horizons, and is primarily due to the ice thawing of the active layer (<2 m).The thermokarst activity produces depressions, subsidences, and ponds, which are usually separated by flat mound peat bogs up to 2 m in height [64].The largest thermokarst lakes that are located within the peat bog are km-size, with a depth of 0.5-1.5 m [64,65].The overwhelming majority of lakes in all three permafrost-bearing zones of the WSL (sporadic, discontinuous, and continuous) have a thermokarst origin, i.e., thawing of frozen peat and clay surface horizons [42,53,61,66].As a result, most thermokarst lakes of the WSL exhibit quite a shallow depth, ≤1 m, in contrast to the deeper lakes of other Arctic regions, originating from surface disturbance, the melting of ground ice, and ice wedges (e.g., ref. [67]).

Remote Sensing Analysis
Satellite imagery from a Landsat-8 Operational Land Imager (OLI, with 30 m resolution), available at USGS Global Visualization Viewer [68], were used to map the lake distribution.We used medium-resolution Landsat-8 images collected at the end of July-August 2013-2014.This period corresponds to the minimal coverage of the territory by lakes and minimal seasonal variation of the lake water level.Besides, this is the period of total disappearance of ice coverage throughout all of the studied area.The sampling of the lake water and lake depth measurements was also performed during this period of the year.Note that, from the view point of optical remote sensing, there is no difference between thermokarst and non-thermokarst lakes, because their reflection yield for Landsat-8 OLI is similar.However, according to published expert estimations of several zones of the WSL [42], the majority of lakes on the WSL are of thermokarst origin.
For the total territory of 1.05 million km 2 , we superposed 134 images (Figure 1).We used not only the images with <10% of cloud, but all of the images that were useful to fill the gaps in the coverage.Nevertheless, for the mosaic, only the cloudless parts of the images were used.In the case where several images were available for the same territory, the one which had the lowest cloud coverage in the middle of the summer was used.The mosaic consisted of images taken in 2013 and 2014, which were combined because the single year data could not provide full coverage of the territory.
The river waters were excluded from the analyses via the creation of the river mask.The data of the river location were taken from national river water cadasters and open street maps.The open ocean and marine coastal zones were also excluded from the analyses.The treatment of satellite images was performed using standard tools of ArcGIS 10.3 software [69], which included classification, vectorization, and surface area quantification.The automatic identification of lakes employed the Fmask algorithm developed for Landsat images, which allows resolving the lakes under some cloud coverage [70].First, for the mosaic of Landsat 8 imagery, the cloud masks were defined for individual images.Then, the cloud masks were removed from the images and replaced by cloud-free fragments taken from other adjacent images.The minimal lake size was chosen as 0.5 ha, based on following.The space resolution of Landsat-8 images is 30 m.Because the pixel size of the image is equal to 30 m × 30 m (900 m 2 ), in the area of 0.5 ha, one can distinguish 5.55 pixels in size.This number of pixels is sufficient for the reliable identification of lakes from the background digital noise of the image.According to the works of our group and Bryksina [18,31,41,46,48], the uncertainty of the lake area measurement using remote sensing is a few percent.Note that the thermokarst lakes and thaw ponds of western Siberia are different from the glacial lakes of other boreal and subarctic regions.The latter are often developed on the moraine till and crystalline rocks, and present highly irregular shapes (skinny or elongated along the glacier direction).In contrast, due to the homogeneity of the soil substrate in western Siberia (1-3 m thick frozen peat), the thermokarst processes in the peat bog of this region always produce the isomeric (round, and much less common, oval) isolated water bodies [14,31,50,61].According to our field and topographical map-based measurements across a sizeable latitudinal gradient of western Siberia, the share of lakes having an irregular shape is less than 10% [49,53,54].
To assess the latitudinal dependence of thermokarst lake properties, the studied territory was divided into latitudinal zones of 0.5 • wide.Such a division of the territory was consistent with the latitudinal gradient of the permafrost and landscape features of the WSL [56,71].Using ArcGIS 10.3 software [69], we first measured the area and number of lakes on each 0.5 • -zone of mosaic of the Landsat 8 imagery.First, we conducted the vectorization and then we determined the area of lakes using the standard procedure of all GIS software.This allowed the quantification of the number, surface area, and volume of lakes, the density of lakes, and the degree of land surface coverage by lakes, as described below.The total lake area (S tot ) in each 0.5-degree zone was computed from Equation (1): where S i is the surface of the i-th lake and n is the number of lakes.The lake fraction is calculated as the ratio of S tot /S o , where S o is the area of each 0.5-degree zone.The lake density was computed as the number of lakes (n) per unit of area, n/S o .In order to estimate the stock of carbon and related elements in lakes, the lake volume (V) was computed following Equation ( 2): where h i is the averaged depth of the i-th lake, which primarily depends on the lake size (see Section 2.3).We used a power dependence between the number and the surface area of lakes in the WSL (correlation coefficient r = 0.99, p < 0.001), in accordance with global distribution law [2]: where k is the relative number of lakes in the histogram intervals, s is the lake surface area, and A and B are the empirical constants that depend on permafrost and landscape context, respectively.

Lake Depth and Hydrochemistry
The depth of thermokarst lakes in the WSL depends on their size.The detailed depth mapping of ~50 lakes larger than 2 ha having a solid bottom (frozen sand or silt) was performed via a Humminbird GPS-echosounder from a PVC boat, along several transects of the lake.The minimal depth of probing was 20 cm.The depth of the PVC boat submersion was between 5 and 10 cm, and necessary corrections for the sensor position were made.In lakes shallower than 50 cm, a manual depth measurement with a calibrated stick was performed.The small lakes (0.5-2 ha), having essentially frozen peat at the bottom with a high amount of porous organic detritus, were monitored via the manual probing of the water depth, across the lake transect or in the middle of the lake.Based on available field measurements of the depth and surface areas of ~150 thermokarst lakes from the sporadic to continuous permafrost zone [14,31,49,50,[52][53][54]61,[72][73][74][75][76][77], the depth was approximated to be equal to 0.54 ± 0.25 m (2 s.d.) for lakes smaller than 2 ha, and 0.85 ± 0.25 m (2 s.d.) for lakes ≥2 ha.The average uncertainty of these values of h i is 20% for n = 150.One has to note that the two discrete numbers of lake depth used in this study for a survey of 727,700 lakes is a first-order approximation.However, all ~150 thermokarst lakes studied by our group over the last nine years in the WSL, across three permafrost zones, were extremely similar and exhibited an average depth of 1.0 ± 0.5 m.This is a particular feature of the WSL thermokarst lakes located within the polygonal tundra, the peat palsa bog, and the ridge-lake-bog complexes.
The total stock of dissolved organic carbon, and major and trace elements in the lakes of the permafrost zone of the WSL, was evaluated based on the available dependencies between the lake surface area and the dissolved components of the lake water, obtained during extensive sampling campaigns in July 2010, 2012, and July-August 2013-2014 [50,[52][53][54]75,76].For this, water samples were collected from the lake surface (0.3 to 0.5 m) in pre-cleaned polypropylene containers and filtered on-site or within 4 h after sampling through disposable acetate cellulose filter units (0.45 µm poresize, 33 mm diameter), using sterile plastic syringes and vinyl gloves.An ultraclean sampling procedure was used [78].The filtered samples were stored at 4-5 • C in the dark, before analysis.The concentrations of dissolved organic and inorganic carbon (DOC and DIC, respectively), cations, and trace elements (TEs), were measured using routine methods for analyzing boreal water samples in the GET laboratory (Toulouse) [79,80].The DOC and DIC were measured using a Shimadzu Total Organic Carbon Analyzer TOC 6000 with an uncertainty of 5%.The trace metals were measured using inductively coupled plasma-mass spectrometry (ICP-MS, Agilent 7500 CE and Element XR), with indium and rhenium as the internal standards and a precision better than ±5%.
For an estimation of the stock of DOC and metals in lakes, the area-averaged values of element concentration in the lake water of discontinuous to continuous permafrost zones [53], complemented with data from discontinuous [49,54,74] and discontinuous to sporadic permafrost regions [50,52], were used.The available databases included a sufficient number of lakes of different sizes, so that the lake size-element concentration dependencies could be obtained for each latitudinal range.Specifically, we sampled ~100 lakes in the discontinuous permafrost zone, ~50 lakes in the continuous permafrost zone, and 30 lakes in the sporadic permafrost zone.For most elements, the concentrations were weakly sensitive (p > 0.05) to the lake size for lakes >0.5 ha.The exception was dissolved organic carbon (DOC), whose concentration was approximated by the power dependence [DOC, mg/L] = 190 × S(m 2 ) −0.26 for lakes of 0.5 to 50 ha and is assumed equal to 10 mg/L for lakes >50 ha.Some elements exhibited a clear latitudinal trend in the thermokarst lakes of the WSL, from south to north (e.g., Ca, ref. [53]).This trend has been taken into account via an equation of polynomial dependence between the element concentration (C i ) and the latitude ( • N), applied to each latitudinal range in which the stock of water was evaluated: where a, b, and c are the empirical coefficients, specific for each element.

Results
The number and surface area of lakes as a function of latitude are shown in Figure 2A,B, and the lake density and relative coverage of the surface area are illustrated in Figure 3A,B respectively.Presented in these plots are the average values of the territory of each 0.5 • -wide latitudinal zones.In the region of 61 • -65 • N (sporadic to discontinuous permafrost), there are large (a factor of two to three) non-systematic variations of all physical parameters of the lakes.In the zone of 65 • to 69 • N (discontinuous to continuous permafrost), the lake number and the relative coverage decrease with the latitude, whereas north of 69 • -70 • N, the lake number and surface area decrease with the latitude increase.
Note that the irregular oscillations of lake density and area coverage, visible in the 0.5 • -wide latitudinal zones, are linked to the spatial non-homogeneity of the thermokarst lake distribution.They disappear after a smoothing procedure in wider (2 • ) latitudinal zones (red dashed line in Figure 3).The results of measured lake parameters within the full WSL permafrost-affected territory (1.05 million km 2 ) are listed in Table 1 and the map of the WSL coverage by lakes is given in Figure 4.There are 0.73 million lakes larger than 0.5 ha, with a total lake surface area of 5.97 million ha.It can be seen from Table 1 that the lake density and lake relative coverage increase by 19.3% and 13.8% from sporadic to continuous permafrost, respectively.The increase in the total lake number and their overall surface area from sporadic to continuous permafrost is equal to 42% and 48%, respectively.polynomial dependence between the element concentration (Ci) and the latitude (° N), applied to each latitudinal range in which the stock of water was evaluated: where a, b, and c are the empirical coefficients, specific for each element.

Results
The number and surface area of lakes as a function of latitude are shown in Figure 2A,B, and the lake density and relative coverage of the surface area are illustrated in Figure 3A,B respectively.Presented in these plots are the average values of the territory of each 0.5°-wide latitudinal zones.In the region of 61°-65° N (sporadic to discontinuous permafrost), there are large (a factor of two to three) non-systematic variations of all physical parameters of the lakes.In the zone of 65° to 69° N (discontinuous to continuous permafrost), the lake number and the relative coverage decrease with the latitude, whereas north of 69°-70° N, the lake number and surface area decrease with the latitude increase.
Note that the irregular oscillations of lake density and area coverage, visible in the 0.5°-wide latitudinal zones, are linked to the spatial non-homogeneity of the thermokarst lake distribution.They disappear after a smoothing procedure in wider (2°) latitudinal zones (red dashed line in Figure 3).The results of measured lake parameters within the full WSL permafrost-affected territory (1.05 million km 2 ) are listed in Table 1 and the map of the WSL coverage by lakes is given in Figure 4.There are 0.73 million lakes larger than 0.5 ha, with a total lake surface area of 5.97 million ha.It can be seen from Table 1 that the lake density and lake relative coverage increase by 19.3% and 13.8% from sporadic to continuous permafrost, respectively.The increase in the total lake number and their overall surface area from sporadic to continuous permafrost is equal to 42% and 48%, respectively.The partitioning of the lake number and surface area among different size ranges is listed in Table 2.The main contribution to the overall area and volume (about 85% and 87%, respectively) is provided by medium and large lakes (>5 ha), with the largest share of the total area and volume (15.5% and 16%, respectively) being kept by lakes whose size is between 20 and 50 ha.The number of lakes increases with a decrease in the size, but the overall surface area and water stock decrease for lakes <20 ha.Thus, the small lakes (0.5-1.0 ha) provide only 3% of the overall area, with less than 2% of the total water volume.It is therefore expected that the overall area of numerous lakes smaller than 0.5 ha will be lower than 3%, although high-resolution images are necessary to confirm this trend.Empirical dependencies of lake number as a function of lake size for the three permafrost zones of western Siberia are illustrated in Figure 5.The empirical coefficients of Equation (3) (A and B) for each permafrost zone of the WSL territory are listed in Table 3.
Taking into account the volume of the lake water in each lake size range and across the WSL territory (Table 3), the total amount of each (<0.45 µm fraction) element in all of the thermokarst lakes (>0.5 ha) of the WSL were estimated (Table 4).The typical uncertainty of these values ranges from ±20% to ±30%, with the exception of some elements (Zn, Cr, Ni, and Ba) exhibiting ±50% of the average value, due to a significant latitudinal trend and lake size dependence of element concentration in the lake water.The partitioning of the lake number and surface area among different size ranges is listed in Table 2.The main contribution to the overall area and volume (about 85% and 87%, respectively) is provided by medium and large lakes (>5 ha), with the largest share of the total area and volume (15.5% and 16%, respectively) being kept by lakes whose size is between 20 and 50 ha.The number of lakes increases with a decrease in the size, but the overall surface area and water stock decrease for lakes <20 ha.Thus, the small lakes (0.5-1.0 ha) provide only 3% of the overall area, with less than 2% of the total water volume.It is therefore expected that the overall area of numerous lakes smaller than 0.5 ha will be lower than 3%, although high-resolution images are necessary to confirm this trend.Empirical dependencies of lake number as a function of lake size for the three permafrost zones of western Siberia are illustrated in Figure 5.The empirical coefficients of Equation ( 3) (A and B) for each permafrost zone of the WSL territory are listed in Table 3. Relationship between the cumulative frequency (the number of lakes versus lake area) of lakes and the lake surface area for the whole territory of WSL (this study, red line), in comparison with lake distribution in the world (Global, dark blue line, [81]) and in Sweden (light blue line, [81]).Relationship between the cumulative frequency (the number of lakes versus lake area) of lakes and the lake surface area for the whole territory of WSL (this study, red line), in comparison with lake distribution in the world (Global, dark blue line, [81]) and in Sweden (light blue line, [81]).

Thermokarst Lake Area and Land Surface Coverage
Overall, the inventory of medium and large thermokarst lakes of the WSL demonstrates an agreement of size distribution and surface coverage of the lakes in this region, compared to the rest of the world.The latitudinal pattern of the number of lakes and their surface area is tightly linked to the topography and landscape conditions of the northern part of the WSL, located within the sporadic to continuous permafrost zone.Between 61 • and 64 • N, the northern taiga is represented by sphagnum-dominated bogs, with pools and an open stand of trees [56].The maximal latitudinal variability of lake coverage is observed within the watershed divide Sibirskie Uvaly (around 63 • N), where the number and proportion of lakes are strongly controlled by minor variations of local topography, such as the alternation of ridge-mire-lake complexes and taiga zones.Further to the north, the lake coverage remains fairly constant between 64.5 • and 71 • N, corresponding to the development of the peat palsa plateau with palsa-hollow patterned mires.The landscape here is highly homogeneous with a dominance of watershed divides of small and medium rivers, offering large flat surfaces suitable for the development of a thermokarst.Finally, a strong decrease in the lake number and area northward of 71 • N may be linked to the dominance of polygonal-roller and polygonal-fissure mires, combined with grass and moss-dominated mires [56].Presumably, the thermokarst processes in the polygonal mires of continuous permafrost zone are less developed than those in the peat palsas plateau, dominating the discontinuous permafrost zone.
Unlike the database that comprises all lakes of the Earth's surface [3], the present study addresses the distribution of thermokarst lakes (>0.5 ha) of the full WSL permafrost-affected territory.A consideration of very small thaw ponds (0.005 to 0.02 ha) in thermokarst-affected regions of the WSL increases the relative surface coverage by lakes to 10%-40%, with an average value of 20%, as shown using Canopus-V data on 18 test sites from 400 to 4000 ha each [82].However, the decreasing of the minimum lake size to less than 0.1 ha over the whole area of the WSL goes beyond the goals of the present study.It is important to note that the distribution of these very small thaw ponds may deviate from the power dependence (Equation ( 3)), as reported in global databases [3,5].The similarity of the B value (Equation ( 3)) among all three permafrost zones suggests a relatively weak variation of thermokarst lake size distribution patterns across the permafrost gradient in the WSL.
Noteworthy is the dramatic difference between the lake coverage of the WSL permafrost-affected territory estimated in this study (5% to 6% of the area) and the proportion of wet zones in the WSL river watersheds, assessed by ENVISAT radar altimetry (40% to 60% of the watershed area during open water period of the year [15]).These authors defined wet zones as various objects that are either constant in time (rivers, lakes, wetlands) or have seasonal variability (floodplains).It follows that the actual coverage of the WSL river watersheds by shallow (<0.1-0.5 m depth) surface water may be significantly higher than the "net" lake area.However, the estimations of the effect of flooding on land coverage by water and the lake abundance (i.e., see ref. [83] for review), or the water level fluctuation in lakes induced by evapotranspiration variation [84], were beyond the scope of this study.

Stock of DOC and Metals in Thermokarst Lakes of the WSL
The specificity of thermokarst lakes of the WSL is their low depth (≤1 m), which allowed, for the first time, a reasonable inventory of the water volume and thus an evaluation of the stocks of dissolved components (Table 4).The typical range of water residence time in the thermokarst lakes of western Siberia is between 0.5 and 1.5 years [54].The overwhelming majority of these lakes are not connected to the rivers, being isolated water bodies, protected by an impermeable permafrost layer both from the bottom (frozen sand and silt), and from the border (frozen peat).Probably for these reasons, the on-going dynamics of thermokarst lake abundances and surface areas are not yet reflected in the hydrological balance of large rivers in Western Siberia [40,47].The stock of dissolved components in lakes on the permafrost-affected WSL territory can be compared to that delivered by all rivers of the WSL from the same territory to the Arctic Ocean.For this, watershed size-averaged, year-round fluxes of carbon, and major and trace elements assessed in previous works [79,80], can be used.A diagram of element stock in thermokarst lakes, relative to that in rivers of the WSL, is presented in Figure 6.Three groups of elements can be distinguished: (i) major and trace elements, whose storage in lakes is less than 10% of that in rivers (DIC, Mg, Ca, Sr, Ba, K, Si, B, Fe, Co, Ni, Mn, and Ce); (ii) elements presenting non-negligible storage in lakes (20% ± 10% of that in rivers): DOC, Rb, Zn, Cu, V, Mo, Zr, As, Nd, and Th; and finally (iii) elements having significant, 30 to 70% storage in lakes, relative to rivers: Cd, Pb, La, Cr, and Al.It can be seen from this classification that major cations, DIC, B, Si, and metals subjected to significant redox transformations (Mn, Fe, and Ce), are essentially present in the rivers because they are delivered by groundwater feeding or shallow subsurface flux [79,80].The groundwater and subsoil feeding are very low in lakes which have frozen peat on the border and frozen sediments at the bottom, throughout the year [54].The elements exhibiting strong affinity to organic matter (Al, Cr, and rare earth elements (REEs)) and metal toxicants (Cd and Pb) enriched in moss, exhibit sizeable storage in lakes of the WSL territory.This is consistent with the surface and suprapermafrost flow that deliver the solutes to the lakes.The thermokarst lakes are typically 1-2 pH units more acidic than the surrounding rivers [50,80], and this may enhance the solubility and mobility of many low-soluble trivalent hydrolysates (Al, Fe 3+ , rare earth elements), from the lake sediment to the water column.Another important source of solutes to the lakes is surface flow from surrounding peat bogs covered by mosses and lichens, consistent with the essentially allochthonous source of DOM in thermokarst lakes [10,50].A high concentration of DOC, combined with an enrichment in Pb and Cd of the surrounding moss cover [57], may be responsible for the sizeable proportion of metal toxicants in lakes compared to rivers.Given that the DOC, Fe, and Al concentrations in the smallest (<100-1000 m 2 ) thermokarst depressions and permafrost subsidences are 3-10 times higher than that in lakes >0.5 ha inventoried in this work [50], and that most trace elements including metal micronutrients are present in the form of organic and organo-mineral colloids [61], the role of small thermokarst water bodies in element stock in surface waters and potential delivery to the hydrological network, may be particularly important and are currently strongly underestimated.For this, coupled land/water observations at a very high spatial resolution [85] may help to constrain the reservoirs and the mobility of carbon, metals, and greenhouse gases in adjacent aquatic and terrestrial biomes.
The role of small (<1000 m 2 ) thermokarst lakes is especially important for the regulation of DOC and greenhouse gas exchanges with surrounding reservoirs (hydrological network and atmosphere).According to available observations of the discontinuous to sporadic permafrost zone of western Siberia, the smallest thaw ponds (10 m × 10 m to 33 m × 33 m) and depressions (1-100 m 2 ) exhibit an order of magnitude higher concentration of CO2, up to two orders of magnitude higher methane concentration, and a factor of three to ten higher concentrations of DOC and related metals [49,50,52].As such, even with their contribution to the total lake surface area of 1%, these small bodies of water may display carbon storage and GHG flux to the atmosphere, which will be comparable to those of large and medium lakes.This hypothesis is verified in the non-permafrost European wetlands: the peatland open water pools are known to act as a net source of CO2 to the atmosphere [86,87].The importance of small (100-200 m 2 ) water bodies for CO2 emission has recently been reported in the polygonal tundra of the Lena Delta observatory [51].
The upscaling of our estimation of the DOC and metal storage in lakes, relative to the river input, requires detailed knowledge of other lakes of the subarctic, since the riverine fluxes DOC, DIC, and most major elements of the circumpolar region, are fairly well defined [88][89][90].The extrapolation of results from well-studied lakes of the Mackenzie Delta region of Canada [9,[91][92][93][94], the Yakutia alasses and yedoma lakes [20,95,96], to much larger territories of boreal plains such as the WSL peatlands, and North-Siberia and Yana-Indigirka lowlands, remains unwarranted.The lakes of these lowlands may stand apart from other studied lakes of the subarctic, in view of their high peat context, low pH, shallow depth, and very low salt content.At present, a large-scale The thermokarst lakes are typically 1-2 pH units more acidic than the surrounding rivers [50,80], and this may enhance the solubility and mobility of many low-soluble trivalent hydrolysates (Al, Fe 3+ , rare earth elements), from the lake sediment to the water column.Another important source of solutes to the lakes is surface flow from surrounding peat bogs covered by mosses and lichens, consistent with the essentially allochthonous source of DOM in thermokarst lakes [10,50].A high concentration of DOC, combined with an enrichment in Pb and Cd of the surrounding moss cover [57], may be responsible for the sizeable proportion of metal toxicants in lakes compared to rivers.Given that the DOC, Fe, and Al concentrations in the smallest (<100-1000 m 2 ) thermokarst depressions and permafrost subsidences are 3-10 times higher than that in lakes >0.5 ha inventoried in this work [50], and that most trace elements including metal micronutrients are present in the form of organic and organo-mineral colloids [61], the role of small thermokarst water bodies in element stock in surface waters and potential delivery to the hydrological network, may be particularly important and are currently strongly underestimated.For this, coupled land/water observations at a very high spatial resolution [85] may help to constrain the reservoirs and the mobility of carbon, metals, and greenhouse gases in adjacent aquatic and terrestrial biomes.
The role of small (<1000 m 2 ) thermokarst lakes is especially important for the regulation of DOC and greenhouse gas exchanges with surrounding reservoirs (hydrological network and atmosphere).According to available observations of the discontinuous to sporadic permafrost zone of western Siberia, the smallest thaw ponds (10 m × 10 m to 33 m × 33 m) 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 related metals [49,50,52].As such, even with their contribution to the total lake surface area of 1%, these small bodies of water may display carbon storage and GHG flux to the atmosphere, which will be comparable to those of large and medium lakes.This hypothesis is verified in the non-permafrost European wetlands: the peatland open water pools are known to act as a net source of CO 2 to the atmosphere [86,87].The importance of small (100-200 m 2 ) water bodies for CO 2 emission has recently been reported in the polygonal tundra of the Lena Delta observatory [51].
The upscaling of our estimation of the DOC and metal storage in lakes, relative to the river input, requires detailed knowledge of other lakes of the subarctic, since the riverine fluxes DOC, DIC, and most major elements of the circumpolar region, are fairly well defined [88][89][90].The extrapolation of results from well-studied lakes of the Mackenzie Delta region of Canada [9,[91][92][93][94], the Yakutia alasses and yedoma lakes [20,95,96], to much larger territories of boreal plains such as the WSL peatlands, and North-Siberia and Yana-Indigirka lowlands, remains unwarranted.The lakes of these lowlands may stand apart from other studied lakes of the subarctic, in view of their high peat context, low pH, shallow depth, and very low salt content.At present, a large-scale comparison of carbon and metal storage in thermokarst lakes and riverine fluxes of elements in the subarctic can only be provided for the western Siberia lowland.

Conclusions
A remote sensing analysis of thermokarst lakes (>0.5 ha) in the sporadic, discontinuous, and continuous permafrost zone of the western Siberia lowland demonstrated that the number of lakes smaller than 1 ha exceeds 33% of the total lake number, whereas their total surface area is only 2.9% of the total surface of WSL lakes.Within the full range of studied lake sizes and areas, a power dependence between the number of lakes and their surface area, consistent with the world-wide trend, is observed.The dependence of the lake number and surface coverage on the latitude exhibits: (i) a highly variable pattern (strong oscillations) between 61 • and 63 • N, within the watershed divide Sibirskie Uvaly, due to the variable topography of ridge-lake-bog complexes within the sporadic permafrost zone; (ii) stable values of lake fraction between 64 • and 71 • N of the peat palsa plateau and the discontinuous permafrost context; and finally (iii) abruptly decreasing the lake fraction northward, north of 71 • N, within the continuous permafrost zone of the polygonal tundra.The obtained laws of lake number and surface area distribution allow the calculation of the total surface area and volume of water.This yielded the dissolved metal and carbon stocks in surface aquatic systems of the permafrost-affected zone of the WSL.
The stock of C and most metals in thermokarst lakes of the WSL does not exceed 10%-20% of the riverine flux of the territory.However, the role of lakes in the storage of Al, Cr, Cd, and Pb is comparable to, or even higher than, the transport of these elements by rivers.A low pH and high DOC in WSL thermokarst lakes compared to other regions of the subarctic may be responsible for such an important role of the WSL lakes in toxic metal storage.A high-resolution (0.01-0.1 ha) inventory of small thermokarst lakes, most susceptible to permafrost thaw in key representative zones of the WSL, will aid in accounting for short-term changes in water, carbon, and metal stocks, under climate warming scenarios.The extrapolation of obtained results to the whole circumpolar region is hampered by the lack of information on other thermokarst lakes from large (million km 2 -scale) territories, notably the boreal plains.

Figure 1 .
Figure 1.Location of three permafrost zones in the study territory of western Siberia.The position of the permafrost provenances in the western Siberian Lowland (WSL) is based on extensive geocryological work in this region [59].

Figure 1 .
Figure 1.Location of three permafrost zones in the study territory of western Siberia.The position of the permafrost provenances in the western Siberian Lowland (WSL) is based on extensive geocryological work in this region [59].

Figure 2 .
Figure 2. The total number of lakes (A) and their surface area (B) along the latitude.The dashed lines mark sporadic (S), discontinuous (D), and continuous (C) permafrost zones.

Figure 2 .
Figure 2. The total number of lakes (A) and their surface area (B) along the latitude.The dashed lines mark sporadic (S), discontinuous (D), and continuous (C) permafrost zones.

Figure 3 .
Figure 3. Dependence of the lake density (A) and relative coverage of the surface (B) on the latitude.The dashed black lines mark sporadic (S), discontinuous (D), and continuous (C) permafrost zones.The dashed red lines represent the smoothing of lake number and coverage using a 2-degree grid.

Figure 4 .
Figure 4. Synthetic map of the WSL coverage by lakes based on Landsat-8 scenes.The cell size is 0.25 degree in latitude and longitude.

Figure 3 .
Figure 3. Dependence of the lake density (A) and relative coverage of the surface (B) on the latitude.The dashed black lines mark sporadic (S), discontinuous (D), and continuous (C) permafrost zones.The dashed red lines represent the smoothing of lake number and coverage using a 2-degree grid.

Figure 3 .
Figure 3. Dependence of the lake density (A) and relative coverage of the surface (B) on the latitude.The dashed black lines mark sporadic (S), discontinuous (D), and continuous (C) permafrost zones.The dashed red lines represent the smoothing of lake number and coverage using a 2-degree grid.

Figure 4 .
Figure 4. Synthetic map of the WSL coverage by lakes based on Landsat-8 scenes.The cell size is 0.25 degree in latitude and longitude.

Figure 4 .
Figure 4. Synthetic map of the WSL coverage by lakes based on Landsat-8 scenes.The cell size is 0.25 degree in latitude and longitude.

Figure 5 .
Figure 5. Relationship between the cumulative frequency (the number of lakes versus lake area) of lakes and the lake surface area for the whole territory of WSL (this study, red line), in comparison with lake distribution in the world (Global, dark blue line,[81]) and in Sweden (light blue line,[81]).

Figure 5 .
Figure 5. Relationship between the cumulative frequency (the number of lakes versus lake area) of lakes and the lake surface area for the whole territory of WSL (this study, red line), in comparison with lake distribution in the world (Global, dark blue line,[81]) and in Sweden (light blue line,[81]).

Water 2017, 9 , 228 12 of 18 Figure 6 .
Figure 6.The mass ratio of element stocks in all thermokarst lakes covering the territory of 105 million ha to annual element delivery to the Arctic ocean by WSL rivers from the same territory.

Figure 6 .
Figure 6.The mass ratio of element stocks in all thermokarst lakes covering the territory of 105 million ha to annual element delivery to the Arctic ocean by WSL rivers from the same territory.

Table 1 .
Thermokarst field parameters in different permafrost zones.

Table 1 .
Thermokarst field parameters in different permafrost zones.

Table 2 .
Lake number, lake area, and volume for different size ranges.

Table 3 .
Parameters of Equation (3) for three permafrost zones of the WSL.

Table 4 .
Dissolved organic and inorganic carbon (DOC and DIC, respectively), and major and trace element stocks in the thermokarst lakes of the Western Siberia Lowland (105 million ha).The major and trace elements are listed in the order of increasing atomic number (periodic table).

Table 2 .
Lake number, lake area, and volume for different size ranges.