Temporal Behavior of Lake Size-Distribution in a Thawing Permafrost Landscape in Northwestern Siberia

Arctic warming alters regional hydrological systems, as permafrost thaw increases active layer thickness and in turn alters the pathways of water flow through the landscape. Further, permafrost thaw may change the connectivity between deeper and shallower groundwater and surface water altering the terrestrial water balance and distribution. Thermokarst lakes and wetlands in the Arctic offer a window into such changes as these landscape elements depend on permafrost and are some of the most dynamic and widespread features in Arctic lowland regions. In this study we used Landsat remotely sensed imagery to investigate potential shifts in thermokarst lake size-distributions, which may be brought about by permafrost thaw, over three distinct time periods (1973, 1987–1988, and 2007–2009) in three hydrological basins in northwestern Siberia. Results revealed fluctuations in total area and number of lakes over time, with both appearing and disappearing lakes alongside stable lakes. On the whole basin scales, there is no indication of any sustained long-term change in thermokarst lake area or lake size abundance over time. This statistical temporal consistency indicates that spatially variable change effects on local permafrost conditions have driven the individual lake changes that have indeed occurred over time. The results highlight the importance of using multi-temporal remote sensing data that can reveal complex spatiotemporal variations distinguishing fluctuations from sustained change trends, for accurate interpretation of thermokarst lake changes and their possible drivers in periods of climate and permafrost change. OPEN ACCESS Remote Sens. 2014, 6 622


Introduction
Climate change in the Arctic has caused significant changes in regional hydrological systems [1].For example, there is consensus that permafrost thaw is widespread across the Arctic and will experience further degradation in the next century due to shifts in climate [2,3].Permafrost thaw alters terrestrial hydrology as it increases active layer thickness and alters the water flow through the landscape with deeper and longer flow pathways to surface water discharge [4][5][6].Changes in streamflow dynamics have been observed, reflecting the influence of permafrost thaw and active layer thickness on hydrologic response at catchment (and landscape) scales (e.g., [7][8][9][10][11]), including a decreasing trend in intra-annual runoff variability [4,5,12,13].Permafrost thaw, with associated water flow and pathway effects, can further increase the connectivity between deeper and shallower groundwater and surface water [14][15][16][17], altering the terrestrial water balance [13] and distribution [17], and the waterborne transport of solutes [18] including that of dissolved organic and inorganic carbon [19,20] through hydrological catchments.
Thermokarst lakes and wetlands in the Arctic offer a window into potential hydrological and connectivity changes as they depend on permafrost and are some of the most dynamic and widespread hydrological features in Arctic lowland regions, covering up to 50% of the total land area of lowland landscapes in northwestern Canada, Siberia and Alaska [21,22].Thermokarst lakes form where melting of ground-ice and surface settlement initiate ponding [23] and Arctic warming may cause changes in the lake distributions.Changes in thermokarst lake size and number have been reported for ice-rich permafrost regions with net increases predominantly occurring in the continuous permafrost zone and net decreases predominantly occurring in the discontinuous permafrost zone [23][24][25][26][27].In the discontinuous permafrost zone, thermokarst lake changes are often characterized by spatial heterogeneity with shrinking lakes appearing alongside stable or growing lakes [24].Potential mechanisms underlying change in lake surface area are (a) changes in precipitation [28]; (b) changes in evapotranspiration [25]; (c) increases to lateral drainage caused by shoreline erosion due to permafrost thaw [11,23,24,27]; (d) increases to internal drainage if underlying permafrost thaws and open taliks appear [24]; or (e) terrestrialization, including encroachment of floating mat vegetation and basin infilling of organic matter and sediment [29].
Assessment of lake size-distributions may be important for understanding terrestrial hydrologic changes in connection to spatiotemporal thermokarst lake dynamics under climate change.Investigating temporal lake size-distributions may also help elucidate if climate change and associated permafrost thaw generates preferential loss of smaller lakes relative to larger lakes or vice versa.In addition, changes in size-distribution can have an impact on the regional biogeochemical cycling, as various physical, chemical and biological characteristics of lakes are related to their sizes [30].Direct observation of spatial changes in terrestrial hydrology via lake size distribution is, however, challenging at large scales in remote areas [31].This is exemplified in Arctic environments, which represent some of the most remote places on Earth and where process understanding, e.g., with regards to permafrost hydrology, may be limited.To determine abundance and surface area of lakes in Arctic regions, several recent studies have thus turned to remote sensing techniques, which have proven valuable for surface water detection [22,25,28].Remote sensing data (and techniques) are useful for monitoring open water extent due to their relatively high spatiotemporal resolution and ability to distinguish water areas from other surfaces.Use of remote sensing data is also a cost-effective alternative to collecting ground-based data.Remote sensing in combination with traditional in situ observations may further provide a more comprehensive understanding of the role of permafrost for catchment hydrology, whereas in situ observations alone are often limited to smaller study areas.As such, remote sensing data can facilitate monitoring of change in remote areas at large scales and can reveal complex spatial variations that cannot be readily elucidated by traditional in situ observations, thereby also helping to improve process understanding.
In this study, we use Landsat remotely sensed imagery to investigate temporal changes in thermokarst lake size distributions (abundance and surface area) that may be brought about by permafrost thaw in northwestern Siberia.Specifically, we determine lake size-distributions from three distinct periods (1973, 1987-1988, and 2007-2009) within the Nadym and Pur river basins and a small sub-basin of the former (denoted as the 7129 sub-basin).Nadym and Pur river basins are chosen because they have high lake density, offering good opportunity to investigate potential shifts in thermokarst lake size distributions.The smaller 7129 sub-basin is included for investigating resolution and water body-classification effects on changes to thermokarst lakes and their size-distributions.Also, for this region, Smith et al. [24] hypothesized, and Karlsson et al. [13] established through coupled lake and hydrological change investigations that sustained thermokarst lake changes may be good indicators of permafrost thawing.In addition to determining lake size-distributions and their changes, we also explore potential relationships between several landscape characteristics (e.g., elevation, latitude, closeness to river network, and presence of peatlands) and lake changes.

Study Area
The Nadym and Pur river basins (48,000 km 2 and 95,100 km 2 , respectively), and the 7129 sub-basin (1,200 km 2 ) (Figure 1) in the northwestern Siberian lowlands are characterized by a low relief (0 and 200 masl) with lakes and wetlands dominating the landscape.The permafrost distributions for the two basins (and Nadym's 7129 sub-basin) are diverse.In the region, permafrost generally extends to depths of 100-300 m below ground [32].With regards to permafrost distributions, the Nadym basin includes 92% discontinuous permafrost and 8% sporadic permafrost, the Pur basin includes 3% continuous permafrost, 77% discontinuous permafrost and 20% sporadic permafrost, and the 7129 sub-basin includes 100% discontinuous permafrost according to the coarse classification map available from Brown et al. [33].The whole study area comprises tundra, forest-tundra and taiga, whereof forest-tundra is the dominant vegetation type [34].The region also contains the most extensive peatlands on Earth, with peat occurring at up to 5 m depth and holding a substantial fraction of the global terrestrial carbon pool [35,36].It is also in these permafrost peatlands where the highest lake densities are found [37].Warming of the surface air temperature has caused warming and thawing of the permafrost in the region since the 1970s, with a temperature increase of up to 0.4 °C per decade at 10 m depth in the Nadym region and 0.45-0.6°C per decade at 8-10 m depth in the Pur region [38].A steady increase of active layer thickness is apparent in the Nadym region, which also corresponds to the increase in surface air temperature [39].A previous study by Karlsson et al. [13] has further shown that at least a part of the present study area has experienced thermokarst lake changes due to the permafrost thaw, which was also reflected in and detectable through associated hydrological changes (e.g., decreasing intra-annual runoff variability and decreasing water storage).Overall, the regional permafrost is considered to be unstable with temperature just below 0 °C [32] and both degradation and talik formations have been observed in the study area [38].

Imagery and Lake Classification
To estimate thermokarst lake size-distributions and determine changes in lake coverage over Nadym and Pur river basins, and the 7129 sub-basin, satellite imagery from Landsat Multispectral Scanner (MSS, with 60 m resolution) and Thematic mapper (TM) (with 30 m resolution) were selected for three different time slices (1973, 1987-1988, and 2007-2009) from USGS Global Visualization Viewer [40] Sensors operating in the visible-infrared, such as Landsat TM, have proven useful for monitoring open water extent for regional medium-resolution studies due to their relatively high spatial and temporal resolution.The selection of satellite imagery was made based on availability, season of acquisition and lack of obscuring clouds.Only imagery from the summer season was acquired to detect surface water.We further aimed to use imagery from the later part of the season (July-September) to avoid the snowmelt period and seasonally inundated areas (Table 1).
Table 1.Landsat images used for thermokarst lake change analysis.The image analysis included percent cloud cover at acquisition (CC), a threshold value used for water body classification, and minimum and maximum spectral values of that threshold obtained as the mean value ± one standard deviation (Min/Max).Landsat images emphasized with black borders were considered in the sensitivity analysis that used the latter values, as described in Section 3.3.Clouds and associated cloud shadows were a major limitation in acquiring useful satellite imagery, which resulted in multi-year time slices for the 1987-1988 and 2007-2009 time periods, and areas with clouds being removed for the 1973 time slice.All images were area-equalized to account for the cloud removal accordingly; hence same areas were extracted for the different time-slices (Figure 1).

Product
Landsat images were first calibrated using the Landsat calibration tool in ENVI 4.8 (ITTVIS) to convert digital numbers to spectral radiance.As multispectral bands are often highly correlated, Principal Component Analysis (PCA) was then used to produce uncorrelated output bands for removing noise components and enhancing the spectral contrast of water features in the imagery [41,42].The PCA was carried out using forward PC rotations, where covariance and eigenvalue statistics were used to generate an output of PC bands containing less noise.Inverse PC rotation was then used to transform PC images back to original format.Landsat TM (30 m) was resampled to the lower resolution of Landsat MSS (60 m) for multi-temporal analysis of lake distribution using nearest neighbor.The preprocessed and resampled Landsat images were classified using automatic segmentation of water bodies, where near infrared (NIR) band 7 was used to classify water from MSS scenes, and NIR band 5 was used to classify water from TM scenes.For the automatic segmentation, training areas were chosen to identify the spectral signatures characterized by water for each individual Landsat Scene.The automatic segmentation procedure allows for region-based classification instead of pixel-based classification, as individual image pixels are grouped into partitions according to two thresholds: minimum area and the identified spectral values.
Smith et al. [24] considered lakes larger than 40 ha to be stable landscape features with little influence of short-term seasonal variability, while other related studies have used a smaller lower lake area boundary for lake change detection (e.g., 0.1 ha [43], 1 ha [28], 20 ha [44]).To make our study comparable to the study by Smith et al. [24] as well as to other related studies, and to avoid seasonal variability, we used a lower lake area boundary of 10 ha (equivalent to 28 pixels) in defining lakes across the study area.However, lakes with surface areas smaller than 10 ha are often numerous in lake-rich regions of the Arctic [27,43,45].To explore the potential influence of the abundance of smaller lakes and ponds in the region on the resultant size distribution, we also considered a lower lake area boundary of 1 ha (3 pixels) for the 7129 sub-basin.This was used for direct comparison with results from the 10 ha lower lake area boundary, even though at the small 1 ha threshold lakes may be influenced by seasonal variability and be subject to misclassification due to image noise.

Post-Processing of Classification Results
Visual inspection of classified water bodies was carried out to check and where appropriate remove objects with water signatures that were not thermokarst lakes (e.g., rivers and oxbow lakes).Multiple images within a season were also compared to investigate intra-seasonal variability, which was found to be low for closed lake basins, while lakes connected to river networks exhibited larger seasonal fluctuations.Surface water bodies and lakes directly connected to river networks were therefore excluded, leaving only closed lake basins for subsequent analysis.All preprocessing of the Landsat images was carried out in ENVI 4.8 (ITTVIS).
The classified water surfaces from the preprocessed Landsat images for the Nadym and Pur river basins, and the 7129 sub-basin, were converted from raster to vector files (discrete objects) to compare lake number and area for the three different time slices.Subsequent analyses of thermokarst lake size-distribution over the three time periods were performed in ArcGIS (ESRI).We also explored if physical elements of the landscape, such as elevation, latitude, closeness to river network, and presence of peatlands influence thermokarst lake change with the use of spatial datasets representing landscape elements.Elevation data was derived from 1 km HYDRO1K data (United States Geological Survey) and elevation values were extracted using the centroid for each individual lake.The coordinates of the centroid of the lakes were also used to define each lake's latitudinal positions.River network data was derived from Digital Chart of the World [46] to analyze distance from each individual lake, using the lake polygon, to the nearest river.The West Siberian Lowland Peatland GIS Data Collection [47] was used to analyze lake number and water surface area change in relation to presence of peatland.

Sensitivity Analysis Regarding Water Body Classification
Apart from visual inspection of potential water misclassification, a sensitivity analysis was also carried out for the water body classification.Due to lack of higher resolution imagery and in situ observations that could serve as ground truthing for accuracy assessment of the classified Landsat images, we carried out the sensitivity analysis based on the classification method adopted.To enable consideration of different (or uncertain) possible lower lake area boundaries (1 or 10 ha), the sensitivity analysis was carried out for the 7129 sub-basin.For this analysis, the automatic segmentation based on the spectral signature for water was implemented using different spectral values defined as the mean spectral value ± one standard deviation across the three time slices for the classification of water bodies (bordered rows in Table 1).Results for lake area and number change based on the use of these different spectral values were compared with those based on the mean spectral value in order to quantify resulting effects of spectral threshold used in the basic lake classification.

Changes to Thermokarst Lake Abundance and Surface Area
There is a predominance of smaller lakes in the Pur and Nadym river basins.In the size ranges of most of our investigation (>10 ha), 83.8%-84.9%and 77.6%-78.3% of the lakes are smaller than 50 ha and 5.0%-5.8%and 9.2%-9.4% of the lakes are larger than 100 ha for the Pur and Nadym basins, respectively, across all three time periods (Figure 2).The 7129 sub-basin (considering a lower lake size threshold of 1 ha) showed similar patterns, even though with somewhat greater fluctuations between the time slices for the 66.7%-67.3% of the lakes that are smaller than 10 ha, and with 94.0%-95.6% of the lakes being smaller than 50 ha.Only 1.2%-1.3% of the lakes were larger than 100 ha in all three time periods in the 7129 sub-basin (Figure 2).Furthermore, lakes have both increased and decreased in area and abundance between different time slices.Between 1973Between and 1987Between -1988 the total area of thermokarst lakes (>10 ha for Nadym and Pur, and >1 ha for 7129) increased, from 1154 ha to 1253 ha (or a 8.6% increase) for Nadym, from 4317 ha to 5109 ha (or 18.3% increase) for Pur, and from 21 ha to 26 ha (or 19.2% increase) for 7129.In contrast, between 1987-1988 and 2007-2009 the total area of thermokarst lakes decreased instead, from 1253 ha to 1126 ha (or 10.1% decrease) for Nadym, from 5109 ha to 3773 ha (or 26.1% decrease) for Pur, and from 26 ha to 24 ha (or 7.7% decrease) for 7129.The total net area change between 1973 and 2007-2009 is then −2.4%, −12.6%, and +12.5% for Nadym, Pur and 7129, respectively.
Total lake number shows similar patterns as the lake area, with increasing total number of lakes over the first time period (1973 to 1987-1988) from 3312 to 3654 (or 10.3% increase) for Nadym, from 9418 to 11,079 (or 17.6% increase) for Pur, and from 150 to 168 (or 10.7% increase) for 7129.The second time period (1987-1988 to 2007-2009) exhibits instead decreasing total number of lakes from 3654 to 3342 (or 8.5% decrease) for Nadym, from 11,079 to 10,619 (or 4.2% decrease) for Pur, and from 168 to 158 (or 6% decrease) for 7129.Overall, the changes result in a net increase in total lake number from 1973 to 2007-2009 (i.e., between the first and the last time slice of remotely sensed data) with +0.9%, +12.8%, and +5.1% for Nadym, Pur and 7129, respectively, and a net decrease in total lake area for Nadym and Pur, but a net increase in total lake area for 7129 during the same period.
The increased lake number but net decreased lake area for the Nadym and Pur river basins can be explained by formation of small lakes as larger lakes become fragmented following partial drainage.This demonstrates the importance of analyzing both lake abundance and lake surface area in combination, in order to understand thermokarst lake dynamics, as lake abundance analysis alone may be misleading.These results also reveal the importance of considering more than two time slices when assessing lake change, as both lake area and number may fluctuate rather than undergo sustained unidirectional change over time, as shown here; when considering only two time slices the former (fluctuations) may be misinterpreted as the latter (sustained change).

Characterization and Rate of Thermokarst Lake Change
The majority of lakes that emerged as appearing or disappearing lakes in the Nadym and Pur river basins were actually lakes that increased above the 10 ha threshold or decreased below the 10 ha threshold, as opposed to newly formed lakes and lakes that disappeared completely (Figure 3).Further, Based on visual inspection of the imagery, the causal mechanisms behind these lake appearances may be both natural (climatic) and anthropogenic (new engineered infrastructure) in nature (Figure 4).This would be consistent with previous findings in this region, for example, by Kampula et al. [48] and Zakharova et al. [49].Complete (−100%) lake drainage was found in all three basins.There were  The observed thermokarst lake changes were spatially heterogeneous, with both completely drained lakes and new lakes appearing alongside stable neighboring lakes (Figure 4).Relating the lake change estimates from the remotely sensed data to the landscape characterizations exhibited no correlation patterns or spatial biases.Lakes (and their changes) were grouped according to higher versus lower elevations, using mean elevation for all lakes as a threshold, northern versus southern locations, using mean latitude for all lakes as a threshold, distance from each individual lake to nearest river (far versus near), and presence of peatlands.These various groupings had no influence on the lake change patterns or the resultant statistics (Appendix A1).
Analysis of the number of drainage events, using the definition from Hinkel et al. [22] with a > 25% reduction in surface area signifying a drainage event, showed that 225 and 566 lakes (2.4% and 5.4% of total lake abundance) drained between 1973 and 1987-1988, and between 1987-1988 and 2007-2009, respectively, for Pur river basin, while 242 and 240 lakes (7.3% and 7.2% of total lake abundance) drained in Nadym river basin during the same time periods, respectively.Altogether,

Size-Distribution of Thermokarst Lakes
As discussed above, the total area and number of lakes have fluctuated over time, with both appearing and disappearing lakes, and with high numbers of drainage events.However, the size-distribution of thermokarst lakes in the size range > 10 ha has not changed over the three time periods in any of the basins, implying that this distribution is relatively time invariant (Figure 2).The 7129 sub-basin reveals somewhat greater temporal differences, predominantly in the smaller lake size range of 1-10 ha, of the size-distribution (Figure 2); however, these differences represent short-term fluctuation rather than change in some direction that is sustained in the long term.Overall, the mean lake size has remained relatively stable during the investigated time periods, with values of 45-46 ha, 34-35 ha and 14-16 ha for the basins Nadym, Pur and 7129, respectively.
Abundance-size distributions for Nadym and Pur show a curved pattern on log-log scale (Figure 5).Fitted regression lines yield log-log line slopes of −1.66 and −1.60 (with R 2 = 0.95 and 0.96, respectively) for the relatively large lakes of >10 ha in the Nadym and Pur basins, and −0.94 (with R 2 = 0.92) for the larger size range of lakes > 1 ha in the 7129 sub-basin.For these lake sizes, the fitted results deviate from a power law distribution with log-log line slope of −1, which represents an effective theoretical upper limit for the exponent of power-law lake size distributions over large size ranges [50].As implied by findings of Muster et al. [43] and McDonald et al. [51], power-law fitting to lake size-distributions may be misleading without sufficient consideration of the relatively small lake sizes.Furthermore, when investigating lake size distributions within finite land areas, power-law cutoff effects will appear for larger lake sizes than some cutoff limit.At any rate, the present results show relative time-invariance in the long-term lake size distributions across all of the different investigated size ranges and resolutions (Figures 2 and 5); even in the lake size range 1-10 ha for the 7129 sub-basin, the distribution exhibits only short-term fluctuations rather than any sustained unidirectional change across the investigated time slices.

Sensitivity Analysis of Water Body Classification
The sensitivity analysis described in method Section 3.3 shows that number and area may differ (with ±8%) depending on the spectral threshold value used for water body classification.The analysis further revealed that use of a higher spectral threshold for water body classification may lead to areas with wet soil being misclassified as lakes, as was the case for some lakes classified based on the 1 ha lower lake area boundary in the 7129 sub-basin.However, the main results and the main conclusion of lake area and number fluctuating, rather than undergoing long-term change over the different time slices, remain robust to such differences and possible misleading aspects of lake number and area classification.

Conclusions
This study's multi-temporal remote-sensing analysis of thermokarst lake changes over three time periods (1973, 1987-1988, and 2007-2009) has revealed fluctuations, rather than sustained change, in total area and number of lakes within all of the investigated Nadym and Pur river basins and the smaller 7129 sub-basin in northwestern Siberia.Over the investigated time extent, lakes both appeared and disappeared and there were numerous drainage events.In spite of these changes, however, the lake size-distributions in the basins remained essentially constant over time, for thermokarst lakes > 10 ha as well as for those > 1 ha.
The found temporal robustness of lake-size distributions indicates that spatially variable changes in local subsurface permafrost conditions (rather than spatially uniform seasonal variation or long-term change over the whole basins) have driven the observed local lake changes.Even after these local lake changes, however, the spatial lake statistics over the whole basins have remained essentially the same as they were before the changes.
The present results highlight the importance of using multi-temporal remote sensing data (rather than such data for only two points in time), which can reveal complex spatiotemporal variation and particularly distinguish fluctuations from sustained change trends.Such distinction is needed for accurate interpretation of thermokarst lake changes and their possible drivers in periods of climate and permafrost change.

Figure 1 .
Figure 1.Location of Nadym and Pur river basins, and 7129 sub-basin, thermokarst lake distribution (as of 1973), permafrost distribution, peatland distribution, and areas excluded from the remote sensing study due to cloud cover.
undergone complete drainage (−100%) were much more frequent in the study area than newly formed lakes.Newly formed lakes were only found in the Pur river basin, where two lakes appeared between1973 and 1987-1988 and nine lakes appeared between 1987-1988 and 2007-2009.

Figure 3 .
Figure 3. Annual rate of lake number change for lakes with the lower lake area threshold of 10 ha, normalized by total lake area for Pur, Nadym river basins, and 7129 sub-basin.
in an increase of average drainage rate from 15 to 26 lakes per year for the Pur basin, and a decrease of average drainage rate from 16 to 10 lakes per year for the Nadym basin over the investigated time periods.

Figure 4 .
Figure 4. Examples of newly formed lakes (blue polygons representing the 2009 lake surface area) in the Pur river basin (A,B); and thermokarst lakes that have been drained (red polygons representing the 1973 lake surface area) in the Nadym river basin (C).White features in subset A and B are associated with infrastructure, e.g., road networks.

Figure 5 .
Figure 5. Log-abundance log-size plot for lakes in Nadym and Pur river basins with lower (10 ha) resolution, and 7129 sub-basin with higher (1 ha) resolution.For comparison, the figure also shows fitted power law functions (lines on log-log plot with slopes −1.66, −1.60, −0.94 and R 2 values 0.95, 0.96, 0.92 for Nadym, Pur and 7129, respectively).