The Geometry of Large Tundra Lakes Observed in Historical Maps and Satellite Images

Abstract

The climate of the Arctic is warming rapidly and this is causing major changes to the cycling of carbon and the distribution of permafrost in this region. Tundra lakes are key components of the Arctic climate system because they represent a source of methane to the atmosphere. In this paper, we aim to analyze the geometry of the patterns formed by large (> 0.8 km${}^2$) tundra lakes in the Russian High Arctic. We have studied images of tundra lakes in historical maps from the State Hydrological Institute, Russia (date 1977; scale $0.21166$ km/pixel) and in Landsat satellite images derived from the Google Earth Engine (G.E.E.; date 2016; scale $0.1503$ km/pixel). The G.E.E. is a cloud-based platform for planetary-scale geospatial analysis on over four decades of Landsat data. We developed an image-processing algorithm to segment these maps and images, measure the area and perimeter of each lake, and compute the fractal dimension of the lakes in the images we have studied. Our results indicate that as lake size increases, their fractal dimension bifurcates. For lakes observed in historical maps, this bifurcation occurs among lakes larger than 100 km${}^2$ (fractal dimension $1.43$ to $1.87$). For lakes observed in satellite images this bifurcation occurs among lakes larger than ∼100 km${}^2$ (fractal dimension $1.31$ to $1.95$). Tundra lakes with a fractal dimension close to 2 have a tendency to be self-similar with respect to their area–perimeter relationships. Area–perimeter measurements indicate that lakes with a length scale greater than 70 km${}^2$ are power-law distributed. Preliminary analysis of changes in lake size over time in paired lakes (lakes that were visually matched in both the historical map and the satellite imagery) indicate that some lakes in our study region have increased in size over time, whereas others have decreased in size over time. Lake size change during this 39-year time interval can be up to half the size of the lake as recorded in the historical map.

1. Introduction

The Arctic is warming at approximately twice the rate of the rest of the globe [1,2], and it is thought that climatic change in this region is proceeding so rapidly that the Arctic climate system is approaching a tipping point [3]. This is causing profound changes to the cycling of carbon in the Arctic and altering the spatial distribution of permafrost [4]. Projected impacts of these changes range from local-scale alterations to the composition of vegetation in the Arctic, to socio-economic problems such as the reduced stability of buildings used for housing and industrial purposes [5].

Tundra lakes are of particular concern in the context of Arctic climate change. This is because thawing permafrost and the resulting decomposition of previously frozen organic carbon significantly enhances the amount of methane emitted into the atmosphere under a warming climate [6]. Tundra lakes, which change shape and increase in size as permafrost thaws, therefore represent a significant source of methane that provides a positive feedback to the atmosphere [7,8], and they are critical elements of the Arctic climate system. Due to the central importance of tundra lakes in the Arctic carbon cycle, understanding their spatio-temporal dynamics is a key on-going scientific challenge [9].

Previous work on the spatio-temporal dynamics of lake-like patterns in the Arctic has shown that changes in their shape and patterning influence the structure of the Arctic climate feedback. For example, studies conducted in Greenland and in the wider Arctic region have shown that surface melting during summer months has increased the abundance of supraglacial lakes on top of glaciers. These supraglacial lakes have lower albedo compared to ice and therefore absorb more of the sun’s energy, causing increased warming and potentially further melting on the ice sheet. Additionally, maps of the bathymetry of supraglacial lakes and streams, together with in situ measurements of their reflectance and depth, allowed Legleiter et al. [10] to to measure the transient melt water flux through streams. Arctic sea ice melt ponds represent visible pools of collected meltwater on the sea ice surface, and these are capable of lowering albedo by altering light scattering properties of the ice surface [11].

There are several potential sources of information on the spatio-temporal dynamic of tundra lakes, and each has its own advantages and disadvantages. It is possible to go into the field and measure the shapes and sizes of lakes manually. This has the potential to provide highly accurate information on lake geometry, but with limited spatial range, since it is not practical to survey hundreds of hectares in this fashion. Similarly, it is possible to analyze topographic maps that incorporate data on lake geometry. Such maps typically cover wide geographic areas, and can provide information on the character of the land surface from several decades ago, but they are not regularly updated owing to the time demands of gathering such data. Aerial photographs also provide a useful source of information on the distribution of elements on the land surface over wide geographic areas, but the cost of flying aircraft means that these data can have limited temporal resolution. Drone technology overcomes this barrier to a degree, but drones still require researchers to be present in the field in order to operate the craft. Satellite imagery provides a high volume of data on the land surface, and has provided valuable information on lake dynamics in Arctic permafrost regions [12,13,14]. A challenge of working with satellite images is that the boundaries of objects of interest can be difficult to define, and image segmentation (the process of partitioning a digital image into multiple segments) can be difficult (for example, see [15]).

In this study, we broadly aim to analyze the geometry of the patterns formed by tundra lakes in the Russian High Arctic. This methodological paper represents a step towards understanding the spatio-temporal dynamics of tundra lakes in the Arctic, and builds on previous work using aerial photography and high-resolution satellite data to gather information on melt pond geometry and statistics (size and shape) [11,16]. We have chosen to investigate the potential for extracting geometrical information on tundra lakes from two sources: topographical maps and satellite images. Our specific objectives are:

• To develop image-processing routines that segment historical topographic maps and recent satellite images, and allow us to visualize tundra lakes in binary (black and white) images.
• To calculate area–perimeter values for each of the individual lakes in our two sources of data.
• To measure the geometrical properties of the tundra lakes in our two sources of data by calculating their fractal dimension.

2. Materials and Methods

Our study was focussed on Western Siberia (${60}^{\circ }{00}^{\prime }{00}^{″}$ Latitude and ${75}^{\circ }{00}^{\prime }{00}^{″}$ Longitude) (Figure 1), and, in this region, a long network of tundra lakes covers the West Siberian Plain [12,17,18]. Our historical map [19] was located in the State Hydrological Institute, Russia. This map was created to display different types of wetland habitats, with the scale of 1:2,500,000, and the data underlying the map were collected from a series of field expeditions during the years preceding the map’s date of publication (1969–1973). The map was created based on aerial photographs of the West Siberian Plain and field work. Detailed information about the analyses of aerial photographs obtained during this field work can be found in [20]. The field work strategy included the setting of temporary stations in different areas of Western Siberia. Different institutions of the Academy of Sciences of the USSR collected the information from these stations during the study period and reported to the State Hydrological Institute (the field work coordinator). The summary of all field work was published in [21]. The map uses a variety of different colors and elements to display different types of wetlands, and the complexity of this color palette presents a challenge to the detection of lakes using computational vision. Our Landsat satellite images were gathered from the Google Earth Engine (G.E.E.) [22] during 2016. Each image was taken at an eye height of $151.94$ km. The Google Earth Engine is a cloud-based and high-performance computing platform for geospatial analysis contains remote sensing data from different dataset such as Landsat, MODIS, Sentinel and ASTER [23]. This is a convenient platform for scientists working across different disciplines and is considered a useful tool for mathematicians, physicists and geologists, as well as climate and environmental scientists studying methane emissions from tundra lakes. The G.E.E. has been used for mapping of population dynamics, irrigation patterns (for example, see [24]), and the spatial patterns formed by vegetation in dryland ecosystems [15].

2.1. Detection of Lakes in Historical Maps

In order to detect lakes in historical maps, the digital images of these maps were segmented based on their color intensity values. Figure 2a shows a part of that historical maps. Firstly, each image was thresholded based on dark blue intensity values (Figure 2b). Mathematically, given an input image $I_{RGB}(x,y)$, we extracted the blue channel $I_B(x,y)$ and then subtracted the gray-level $I_{gray}(x,y)$ image of the $I_B(x,y)$, which was done by

$$I_{blue}(x,y)=I_B(x,y)-I_{gray}(x,y),$$

(1)

where $I_{blue}(x,y)$ is the output image, and is the same size as the input image. This output image was converted to a binary (black and white) image in order to find the center of each lake (Figure 2c). Topographical maps are characterized by many colors and lines, which are confounding in the context of image processing. To tackle this, we created a flat linear structuring element with the pixel of interest (the pixel being processed) located at its center. This flat linear structuring element was a matrix that contains zeros and ones. The pixels with values of 1 define the neighborhood of pixels (structured as lines in the shape) that were included in the processing. The elements of the matrix with values of 0 were not included. This structuring element is symmetric with respect to the neighborhood center, and there was a distance of 10 pixels between the centers of the structuring element members at opposite ends of the line and an angle of 25 degrees from the horizontal axis, which were each chosen by trial and error. We then subtracted that linear structure from the $I_{blue}(x,y)$ image (Figure 2d). Our final image segmentation involved applying a region growing procedure to the images that were subjected to this line removal strategy. Since this region-growing step was common to both topographical maps and satellite images, we describe it separately below.

2.2. Detection of Lakes in Google Earth Engine Images

In order to detect lakes in satellite images derived from the Google Earth Engine (Figure 3b), we employed a decision-making model using a support vector machine classifier (SVM) [25], which is a supervised statistical and discriminative classification technique. We trained the SVM using two training samples of input images (size $7\times 7$ pixels). The first training sample consisted of images representing lakes. The second training sample consisted of images representing the image background (not lakes). The vector information of each sample was extracted by representing each pixel p as the Euclidean distance between its corresponding coordinates on the red, green, and blue channels and the origin, which was done by

$$p(x,y)=\sqrt{{(r-0)}^2+{(g-0)}^2+{(b-0)}^2}=\sqrt{r^2+g^2+b^2},$$

(2)

where $r,g,b$ are the red, green, and blue channels of each sample ($7\times 7$ sub-image) respectively of each pixel position $p(x,y)$.

This vector information was fed to the support SVM classifier to build the SVM training model. Each input image was then divided to overlapping windows of size $7\times 7$ pixels and its vector information was extracted in order to classify each window as a lake or not a lake. Lakes were represented by white pixels and pixels that were not lakes were represented by black pixels (Figure 3d). Our final image segmentation involved applying a region growing procedure to the images that were subject to this line removal strategy.

Region growing for segmentation of historical maps and Google Earth Engine images. We applied a standard region-growing segmentation to binary images of both historical maps (Figure 3c) and satellite images (Figure 3d). The pixel at the center of each lake was thought of as the pixel of interest. The difference between the intensity value of the pixel of interest and the intensity value of the region mean was calculated. If that difference was less than or equal to a predefined threshold, a new pixel was added to the region; otherwise, the process was stopped [26,27]. We experimented with a range of threshold values and found that $0.05$ yielded the clearest segmentation of lakes. This region-growing segmentation is shown in Figure 3.

2.3. Calculating the Geometrical Properties of Lakes

We used standard connected component analysis to calculate the area and perimeter of each lake [28]. These quantities were calculated in terms of pixels, and reported in metric units by multiplying the number of pixels with their corresponding pixel size. The pixel size for the historical map is $0.21166$ km/pixel, while the pixel size for the Google Earth Engine is $0.1503$ km/pixel.

The fractal theory introduced by Mandelbrot [29] can be used as a method to study partially correlated (over many scales) spatial phenomena that are not differentiable but are continuous. This theory helps quantify complex shapes or boundaries and relate them to underlying processes that may affect pattern complexity. For simple objects like circles and polygons, the perimeter P scales as the square root of the area A. However, for complex planar regions with fractal curves as their boundaries,

$$P\sim \sqrt[D]{A},$$

(3)

where the exponent D is the fractal dimension of the boundary curve.

We employed a method for computing the fractal dimension that is an extreme value analysis based on the lower edge of the area–perimeter data points. To do this, we took the convex hull of the data points in the ($A,P$)-plane, identified the lower edge, and computed its slope. This procedure should guarantee that fractal dimension does not decrease as area increases, and explains why there are so few data points that maintain the non-decreasing property of the fractal dimension. The method may be seen as a robust alternative to other methods of computing fractal dimensions, such as the lexicographic ordering method used for calculating melt pond fractal dimension [11]. We used a scatter plot of perimeter–area obtained from the historical map and satellite images to compute fractal dimensions using MATLAB (R2015a, The MathWorks Inc., Natick, MA, USA). A summary of our method is shown in Figure 4.

3. Results

3.1. Geometrical Properties of Tundra Lakes in Historical Maps and Satellite Imagery

There were 864 tundra lakes detected in the historical topographical map. The area–perimeter values for these lakes were plotted on a log–log scale, and analysis of this plot shows that the linear trend in the data changes slope about 30 km${}^2$ in terms of area (Figure 5a). These lakes are greater than $0.0008$ km${}^2$, which are the smallest lakes that our technique can detect without any false positives. There are numerous features in the map that are smaller than $0.0008$ km${}^2$. Visual inspection of our images indicate that these small features are not lakes, but our method could falsely label them as lakes because of their small size.

There were 1103 tundra lakes detected in the satellite images. These lakes were greater than $0.8362$ km${}^2$. This minimum size was set in order to make the total number of detected lakes from both sources close to each other. This allowed us to compare the fractal dimension of the detected lakes in Western Siberia using both sources of data. We are able to make the following observations concerning the detection of lakes: (1) if the minimum size of the detected lakes is set to be the same in both the satellite imagery and the historical map, the total number of detected lakes would be much higher in the dataset of satellite images, and this may compromise the comparison of fractal dimension values between the two sources of information; and (2) there are physical limitations on hand-drawn polygons. In the context of the historical map we have analyzed here, these limitations encompass issues ranging from the idea that it is difficult to hand-draw very small objects consistently to the notion that some lake regions may have been colored incorrectly by the person manually rendering the map. We suggest that these factors underlie the observation that the total number of detected lakes is lower in our historical map compared to the satellite imagery in the same region (see Figure 9 below).

The area–perimeter values for each of these lakes were also plotted on a log–log scale, and analysis of this plot shows that the linear trend in the data changes slope about 10 km${}^2$ in terms of area (Figure 5b).

To investigate this deviation in detail, we computed the fractal dimension $D\left(A\right)$ as a function of tundra lake area detected on the historical map. A plot of lake fractal dimension against lake area shows that lakes within the size range 1–70 km${}^2$ generally fall on a linear trend-line and have an average fractal dimension of $1.62$ (Figure 6a). However, the fractal dimension of lakes larger than 100 km${}^2$ displays interesting behavior: the fractal dimension of some lakes increases above the linear trend line to $1.87$, while the fractal dimension of other lakes falls below the linear trend to $1.43$ (Figure 6a). This indicates that the fractal dimension of large tundra lakes in the that historical maps we have analyzed bifurcates.

The fractal dimension of tundra lakes in our satellite images within the size range 1–50 km${}^2$ remains constant at about $1.70$ (Figure 6b). However, the fractal dimension of the lakes larger than ∼100 km${}^2$ changes dramatically. For some lakes, the fractal dimension is almost $1.95$, whereas others are just $1.31$ (Figure 6b).

As expected given our method of computation, we obtain few data points for the fractal dimension of the lakes that we have investigated. While this may not look sufficient, especially when we have generated a large data set, these plots ensure that we include only those data points that retain the non-decreasing property of the fractal dimension. To confirm the observed features of lake pattern fractal dimension that we report here, we have employed another parameter called the elasticity and defined the variance $\sigma$ of $\log \left(P\right)$. Simply stated, variance describes how far a data set is spread out, and formally this is the average of the squared differences from the mean. The elasticity covers the entire cluster of points on the $(A,P)$-plane. In contrast, the fractal dimension plots show only those data points that maintain the non-decreasing fractal dimension property. The onset of fractal dimension splitting may be identified with the beginning of elasticity cleavage as shown in Figure 7.

3.2. Statistical Properties of Tundra Lakes Distribution in Historical Maps and Satellite Imagery

These area–perimeter data allow us to define some statistical properties of lake distribution in order to use this information in future climate modeling. The distribution of small and large tundra lakes is not stable in time but varies substantially depending on the evolution of climate. What statistical distributions do tundra lakes exhibit depending on historical period? The probability distribution function (PDF) $\mathrm{prob}\left(A\right)$ in both historical map and satellite imagery can be described by a power law scaling $prob\left(A\right)$∼$A^{\zeta }$ with the scaling exponent $\zeta$.

Figure 8a shows the log–log plot of the PDF with the scaling exponent $\zeta \approx -1.45$ for tundra lakes detected from our historical map. In addition, we can see that the lakes with a length scale larger than 100 km${}^2$ are power-law distributed with a tail exponent ($\tau =2.28$). Similar exponents have also been observed for other types of lakes on Earth [30].

Figure 8b shows the log–log plot of the PDF $\mathrm{prob}\left(A\right)$ with the scaling exponent $\zeta \approx -1.80$ for tundra lakes detected in satellite imagery derived from the Google Earth Engine. In this case, lakes with a length scale larger than 70 km${}^2$ are power-law distributed with a tail exponent ($\tau =1.93$).

4. Discussion

4.1. Image Analysis Effectiveness and Changes in Lake Size over Time

The images of tundra lakes that we have examined in this study, especially the satellite imagery derived from the Google Earth Engine, contain complex information such as color and shape. This presents a substantial challenge in terms of image segmentation. As mentioned in Section 2.2, our image segmentation technique is based on the difference in the intensity values between the pixel of interest and the region of interest. Therefore, small intensity value differences within a lake region itself may cause a single lake to be detected as two or more lakes by our algorithm. This will provide an erroneous number of total detected lakes and will also result in area and perimeter values that are incorrect. Although we attempted to account for this by experimenting with different thresholding values during the image segmentation process, it is possible that these confounding factors mean that some of our detected lake boundaries are artificial. Some of the image segmentation challenges we have faced are shown in Figure 9.

In this paper, we do not have a time-series of satellite image or a time-series of historical maps with which to investigate how tundra lake geometry has changed through time in our study region. However, we have compared specific lakes that are recorded in both our historical map and our sample of satellite images, and this comparison is located in

Table 1. Lake area and perimeter comparison.

Historical Map				Google Earth Engine
Original Image	Binary Image	Area (km${}^{\mathbf{2}}$)	Perimeter (km)	Original Image	Binary Image	Area (km${}^{\mathbf{2}}$)	Perimeter (km)
		116.883	48.309			98.920	41.899
		95.916	35.515			76.275	32.057
		52.416	30.615			117.316	47.914
		478.553	115.059			159.262	53.789

. This preliminary analysis illustrates the nature of the variation in lake size in our study area during the time interval that separates the historical map and the satellite imagery. This time interval is around 39 years. These focused example lakes indicate that, in our study region, some lakes have increased in size over time, whereas others have decreased in size over time, and changes in lake size during this time interval can be up to half the size of the lake as recorded in the historical map.

Clearly, it is difficult to draw substantial conclusions from a comparison of this nature because the two data sources are so different, each with its own advantages and limitations. For instance, Figure 9 suggests that there are some lakes in the Google Earth Engine imagery (see Figure 9b) that do not exist in the historical map (see Figure 9a). However, visual comparison of the two datasets indicates that the lakes were in fact present in the historical map, but they were mislabeled. Specifically, they were colored green instead of dark blue, which is the color used for the vast majority of lakes in this map, and so they were missed by our algorithm during the segmentation of the historical map. We have included this comparison to highlight an emerging challenge in this research domain: the integration and comparison of different sources of information on tundra lake geometry. This challenge arises because certain types of data on historical lake geometry, such as the historical map used here, may simply not exist for the present day. Similarly, satellite data for a given region may span a decade or so, but other sources of information, such as maps or aerial photographs, could provide valuable historical context for such data.

4.2. Phase Transitions and the Complex Geometry of Tundra Lake Patterns

The fractal dimension of the tundra lake patterns in both the historical map and satellite images indicate that lake patterns with a short length scale (100 km${}^2$) have a non-changing fractal dimension of around 1.6–1.7. This is in agreement with previous research on the fractal dimension of some elements of the Siberian tundra landscape (where the fractal dimension of tundra lakes is generally 1.68) [31], and is typical of Earth’s lakes in general, which have no tendency to change fractal dimension over time and also remain close to a limit [30].

Although several natural patterns such as clouds and rivers exhibit constant fractal dimension, there is a system of melt ponds on Arctic sea ice that is characterized by dramatic geometrical variability. Melt ponds on Arctic sea ice develop due to phase transitions (sea-ice to ocean-water) as a result of increased temperature and radiative forcing. It has been shown that the fractal dimension of Arctic melt ponds shifts abruptly from 1 to 2 as ponds grow in size, with the critical transition zone around 100 m${}^2$ [11]. Fundamentally, the thawing of permafrost and development of tundra lakes is also a phase transition phenomenon, where a solid turns to a liquid, albeit on a large regional scale and over a period of time that depends on environmental forcing and other more local factors. We suggest that the "bifurcation" we observe in the fractal dimension of the tundra lakes that we have studied represents a transition phenomenon that is analogous to that seen in melt ponds [11].

However, it should be emphasized that, in the case of the tundra lakes we have studied, the fractal dimension “bifurcates”, whereas, in the case of Arctic melt ponds, the fractal dimension only increases in a nonlinear fashion [11]. This may imply that the geometrical evolution of tundra lakes is more complex than Arctic melt ponds. We suggest that the bifurcation of fractal dimension in the tundra lakes we have investigated may be explained by the drainage of large lakes to the wider regional hydrological network. In particular, larger lakes become fragmented following partial drainage [9] and this leads to a decline in lake area, and thus to a lake that itself has a lower fractal dimension. This process is characteristic of a lake transforming to peatland [32].

Tundra lakes with a fractal dimension higher than average and close to 2 have a tendency to be self-similar with respect to their area–perimeter relationships. We are able to recognize a portion of a lake that is connected to a complex lake system, such that the perimeter to area ratio of the entire lake is approximately the same as that of the portion. For example, by computing the perimeter to area ratios for a whole lake and a portion of that lake (Figure 9b), we find $P/A\approx 1.2$ km${}^{-1}$ for both the lake and its portion. The measurement of complexity properties such as this may allow tundra lakes to be classified on the basis of their structure and fractal dimension.

Our statistical analyses of the size-frequency distribution of tundra lakes on the Arctic land surface also indicate that these lakes follow a power law. This highlights that a first-order property of the tundra landscape is the high abundance of very small lakes together with considerably smaller numbers of larger lakes. Such knowledge may be used in stochastic models of methane emissions from boreal permafrost [33,34]. For example, in [34], the authors assume that the lake size distribution follows the Pareto distribution (a power law with specific exponents). In such cases, the scaling exponents of the power law could be found definitely from observations using the image analyses of tundra lakes that we have developed here.

Finally, our results indicate that images derived from the free access platform Google Earth Engine can be used to detail and diagnose the fractal dimension of tundra lakes in different regions of the Arctic. By virtue of its reliance on an open-access source of imagery, our approach may represent a convenient blueprint of large-scale surveys of tundra lake geometry in other regions of the world. While multi-spectral Landsat data are widely used for analysis of high resolution tundra lake dynamics [14], we suggest that our approach of using open access imagery is also appropriate to investigate spatiotemporal dynamics of tundra lakes. At the very least, it provides a means to track lake size change through time that will be useful for numerical climate model parametrization [35,36].

5. Conclusions

Our study of tundra lake geometry is motivated by an overarching desire to understand the spatiotemporal dynamics of tundra lake patterns under climatic change. This is important because tundra lakes represent a significant source of methane to the atmosphere, and they change shape and increase in size as permafrost thaws. In this methodological paper, we have taken a step towards this goal by analyzing the geometrical properties of tundra lakes in the Russian High Arctic. These lakes were recorded by historical maps and by satellite imagery derived from the Google Earth Engine. Our specific conclusions are as follows:

(1)

An image processing strategy that combines color thresholding and region growing allows us to segment historical maps and satellite images. This isolates lakes from other information in these images. In the case of historical maps, this other information includes the distribution of wetlands, and in the case of satellite images includes vegetation, soil and other components of the land surface.

(2)

The fractal dimension of lakes within the size range 1–70 km${}^2$ observed in historical maps is on average $1.62$ (Figure 6a). However, the fractal dimension of lakes larger than 100 km${}^2$ displays bifurcating behavior: the fractal dimension of some of these large lakes is $1.87$ and the fractal dimension of some of these large lakes is $1.43$ (Figure 6a).

(3)

Similar behavior in the fractal dimension of lakes was observed in satellite imagery. The fractal dimension of lakes sized 1–50 km${}^2$ is around $1.70$ (Figure 6b), but the fractal dimension of the lakes larger than ∼100 km${}^2$ ranges from $1.31$ to $1.95$ (Figure 6b).

(4)

Area–perimeter measurements for each of the individual lakes in our analyses show that, for lakes in our historical maps, those with a length scale larger than 100 km${}^2$ are power-law distributed with a tail exponent ($\tau =2.28$). For lakes in the satellite images we have analyzed, those with a length scale larger than 70 km${}^2$ are also power-law distributed with a tail exponent ($\tau =1.93$). Similar exponents have been observed for other lakes on the Earth surface [30].

(5)

We have undertaken a preliminary analysis of the changes in lake size that have taken place in the 39 year interval between the production of the historical map we have used and the satellite images we have examined. We find that some lakes in our study region have increased in size over time, whereas others have decreased in size over time. Changes in lake size during this time interval can be up to half the size of the lake as recorded in the historical map.