An Effective Low-Cost Remote Sensing Approach to Reconstruct the Long-Term and Dense Time Series of Area and Storage Variations for Large Lakes

4.1.1. Generating Lookup Tables of Shoreline Sater Frequency and Corresponding Lake Area

According to JRC GSW occurrence maps (1984–2015), we selected lake shores with more significant and gradient changes in water frequency as the target area. The lake water frequency higher than or equal to 98% is classified as permanent water bodies. The rest are divided into about 20 grades by the frequency interval of 5%. In order to calculate the sensitivity of the water frequency grade determination, we also calculated the relative standard deviation under the grading interval of 2.5% and 7.5% (Figure S1 in supplementary information). Taking North Aral Sea in the can be referred in Table S1 as an example, the experimental results showed that the estimation results of the grading intervals as 2.5%, 5% or 7.5% are basically consistent. In order to facilitate the calculation, the interval of 5% was selected as the grading interval value for the following experiment analyses. We calculate the water body area of each grade and generate the area-frequency lookup tables for each lake.

4.1.2. Extracting Lakes Shorelines from Monthly GSW Dataset and OLI Imagery

Monthly GSW data (1984–2015) and Landsat 8 OLI/TIRS C1 Level-1 products (2013–2018) are combined as a 35-year (1984–2018) dataset for reconstructing long-term lake inundation area variations. The processing procedures of the lake extraction from the two datasets are introduced as follows:

Monthly surface water area data of GSW dataset from 1984 to 2015 can be downloaded from GEE. The monthly lake inundation data is further screened manually and the data with incomplete boundaries or obvious noisy shifts are removed. Although 320 months of data are available on GEE, it was found that not all months of our study area is usable in actual operations after the manual inspection.

Level-1 Landsat data images during 2013–2018, available in study areas, are all downloaded from USGS. The method of extracting water inundation areas via a two-step NDWI (Normalized Difference Water Index) threshold segmentation scheme was applied in this study [24]. In the first step, all potential water body pixels are filtered by a loose initial threshold segmentation based on calculated NDWI maps of selected multiple spectral images. The threshold is set to NDWI > –0.1, which is a conservative estimate to mask out various water bodies in different spectral properties at global scale; the second step is locally iterative segmentation for each potential lake entity, and whose goal for obtaining stable lake boundary via consecutive segmentations with the least shifted water area. The NDWI, generated by the normalized water index between the green (ρgreen) and near infrared (ρNIR) spectral bands, is proposed by McFeeters [25]. Moreover, mapping accuracy of water body boundaries is guaranteed by additional visual inspection and manual editing. Compared to the single-threshold segmentation method which uses a fixed threshold to delineate outlines of different lakes with obviously distinctive spectral features, this method can efficiently improve the accuracy of water body extent delineation [26]. Eventually about 101–127 months of lake shorelines can be used for the five lakes via two post-processing methods and the specific information can be referred in Table S2 of supplementary information.

4.1.3. Estimating Lake Areas via Shoreline Water Frequency from Lookup Tables

The completeness of the predefined lake shoreline is closely related to the accuracy of the reconstructed lake area based on lookup table. However, not all of extracted lake shorelines are complete due to imagery quality and mapping error. In order to improve the area reconstruction accuracy, we manually screened the lakeshores before retrieving their overlapped water frequency values, and further filtered the lakeshores in the areas with flat terrain and obvious changes. Then we convert the extracted lake shoreline into a continuous finite number of points and calculate the number of points that fall on the grades representing different water inundation frequencies. Finally, we obtain a mean value of all the points to represent the water frequency of lake shorelines at the specific observation date, which can be further converted to the lake area based on the statistical regression relationship between the shoreline water frequency and lake area in our established lookup table.

4.1.4. Accuracy Assessment of Reconstructed Lake Area Time Aeries

Two methods were used to evaluate the performances of the reconstructed lake areas. We consider the change in lake area is closely related to the water level variations, so we first compared the reconstructed lake areas with water levels from the altimetry satellites. The second method is to directly compare with mapped areas from Landsat imagery for the North Aral Sea with relatively small area. Monthly mean values of lake water level are used for evaluating the estimated area values. We pick the area and water levels of lakes in the same period and calculate the correlation with R² values. R² values is a coefficient of determination, which represents the ratio of the sum of squares caused by X in the sum of squares of Y. In the unary regression analysis, the closer R² values is to 1, the higher reference value of the related equation; on the contrary, the closer to 0, the lower reference value:

where $\overline{\mathrm{y}}$, $S{S}_{tot} , S{S}_{res}$ are the mean value of y, sum of Squares for total and sum of squares for error. Actual area value is from the North Aral Sea, because North Aral Sea is the smallest of the five study cases, and it is easiest to obtain its area by mapping two senses of Landsat images. Then, the relative error between the estimated area and the mapped actual area is calculated to evaluate the performance of area reconstruction. The relative error is calculated as follows:

where $\mathsf{\delta }$, $E, \mathrm{A}$ are the mean value of relative error, estimated area and actual area.

5.3.1. Lake Athabasca

Figure 4 presented the time series of Lake Athabasca showed overall decreasing trend in area during 1984–2018 with an obvious upward shift in 1996–1998 when the area was significantly higher than that in other years with the increase of 307.5 km². Figure 4 illustrates that the maximum inundation area was observed in April 1995 to 7643.58 km² and the minimum area was observed in July 1997 to 7336.08 km². The reason why there was a sudden increase in the area was explained [28,29]. The study mentioned that the rise in water level was caused by flooding and the floods occurred in the summer of 1996 were achieved by the formation of ice jams. These results indicated our reconstruction of water level and area confirm the description of the 1996 water mass surge in published literature. It must be mentioned that the change profile of Lake Athabasca does not show remarkable seasonality in area, water level and volume between 1984 and 1992 due to the lack of adequate time series of satellite observations. The seasonal variations and trends of DAHITI water level in study of Schwatke et al. [30] agreed well with our results. As shown in Figure 4, the volume variation of Lake Athabasca ranged from +21.44 km³ (increased value) to –19.35 km³ (reduced value) and the months at which the maximum and minimum occurred coincided with the water level time series in 1996–1998.

5.3.2. Caspian Sea

The time-series data from Figure 5 shows the area, water level and volume variation of Caspian Sea was generally decreasing with a nearly linear upward trend in 1987–1995 and 2001–2005, significant declines in 1995–2001 and 2005–2015 and relatively stable in 2015–2018. Worth noting that the maximum area occurred in June 1995 with a value of 381744.40 km² and the minimum area occurred in February 2016 with a value of 365,910.59 km² (see Figure 5). A possible explanation of these variations in areas could be largely attributed to hydrological condition changes of the Volga drainage basin [21]. For instance, Chen et al. [31] used climate model-predicted precipitation (P), evaporation (E), and observed river runoff (R) of Volga to reconstruct long-term Caspian Sea level (CSL) changes for 1979–2015. The conclusion rising surface temperatures had led to more evaporation in Caspian Sea than precipitation and runoff absorbed, leading to a reduction in area in recent decades was caused by the observed rapid CSL increase (about 12.74 cm/yr) and significant drop (~6.72 cm/yr) during the periods 1979–1995 and 1996–2015, respectively. At the same time, the annual cumulative evaporation rate of Caspian Sea will continue to increase in the foreseeable future as the Northern Hemisphere continues to warm up is pointed. In Ozyavas et al. [32], it is reported that the Volga River discharge in conjunction with evaporation over Caspian Sea were the two main components controlling water-level fluctuations of Caspian Sea. More specifically, Volga River discharge positively correlated with Caspian Sea water level rise and drop from 1993 to 1997 and in 1999, respectively. In 1997–2003, evaporation should dominate over the Volga River runoff and primarily controlled CSL fluctuations. As shown in Figure 5 the range of lake volume variation was shifted from +753.04 km³ (increased value) to −488.47 km³ (reduced value) which agree well with the overall downward trend in Figure 5. Another important factor attributed to the CS water-level decline is thought to be the transfer of CS water to the KBG bay through opening the dam in the KBG strait in 1995.

5.3.3. North Aral Sea

We reconstructed the time-series data of North Aral Sea for the 35 years. Clearly to see, the change of lake area was not characterized by a monotonically increasing tendency. The evolution processes can be divided into three sub-stages: the first large-scale fluctuation in 1984–2005, the last small-scale fluctuation in 2006–2018 and the mid-sudden rise of water level in 2005–2006, which brought about increase of the average water level (in Figure 6). More specifically, the average water level of last period was 1.13 m higher than the first. In 1984–2018, the maximum water level appeared in April 1999 with a high lake stage of 43.11 m and the minimum water level (39.71 m) appeared in October 2001 (seen in Figure 6). Kouraev et al. [33] pointed out construction and destruction of several dams between the two parts of the lake attributed to the level of the North Aral Sea oscillated between of 39 m and 42 m and other research data confirmed this view. For example, Létolle and Chesterikoff [34] found that dams have been established to isolate them and keep the level of North Aral Sea since the separation of North Aral Sea and the South Aral Sea in 1989. In July 1992, the water pressure washed away the first 1-m- height sand dam. Then the second, third, fourth dams were built one by one which all were collapsed by flood in 1992–1999. The peak water level of North Aral Sea reached in April 1999 and then fell by 2.50 m at September 1999 both due to these incidents. The water level in North Aral Sea in 2005–2006 had risen markedly and had remained steadily rising since then attributed to the construction of fifth dam in August 2005. Kouraev et al. [33] analyzed the evolution of North Aral Sea based on radar altimeter data from T/P, Jason-1, ENVISAT and GFO satellites for 1992–2006. The water level changes presented in this study showed a similar pattern as average water level from our study which increased from 41.04 m in 1984–2005 to 42.17 m in 2006–2018. Moreover, the range of volume variation in North Aral Sea is +2.11 km³ (increased value) and –9.40 km³ (reduced value) could be observed in Figure 6 by reconstructing time series of volume variation. These discussions above confirm our work in analysis reasons for water level change, and the novel method can be effectively applicable in other lakes whose changes are similar with Northern Aral Sea.

5.3.4. Lake Balkhash

As shown in Figure 7, Lake Balkhash presented a unique change pattern of area, water level and volume variation during 1989–2018. During the three decades, the tendency can be divided into four sub-periods (1994–1995, 1999–2006, 2010–2013 and 2016–2018) when lake levels showed upward shifts and four sub-periods (1989–1993, 1995–1999, 2006–2010 and 2013–2016) when lake levels showed downward shifts. Figure 7 showed the maximum area of Lake Balkhash was 17,062.79 km² in March 2006 and the minimum area was 16,522.57 km² in September 1998. Moreover, the volume variation of Lake Balkhash ranged from +19.21 km³ (increased value) to −25.60 km³ (reduced value) in Figure 7. For this, Bai et al. [35] proposed reservoirs and dams were two human factors that affect lake total runoff and lake area changes based on analyzing the surface area changes of nine major lakes from 1975 to 2007. As we all know, Lake Balkhash is mainly supplied by the Ili River, which covers 73% of the total inflow of water. Since 1974, artificial runoff had been controlled through the filling of the Kappa Gai Dam, and the annual flow of Ili River into Barhash has fallen to 77% of the mean before 1969 [36] which led water level and area to drop dramatically. Chao et al. [37] mentioned that there was a Nino Phenomenon from 1997 to 1998, which may cause the steep increase of water level in Lake Balkhash from 1998 to 2007. Propastin [38], who agreed with views from Chao et al. [37], made a point that the overall increase in precipitation and temperature in the catchment basin resulted in the water level of Lake Balkhash increased during 1998–2007. The temperature changes above Lake Balkhash could also support this view. The annual mean temperature over Lake Balkhash varied from 6.0 °C in 1976–1992 to 6.7 °C in 1998–2007. The temperature of this change directly led to a rise in water levels after 1998. These prior research findings can be well revealed in our reconstructed lake change series that area, water level and volume variations all have sharp rises after 1998.

5.3.5. Lake Issyk-Kul

Interestingly, the water level curves for both Lake Issyk-Kul and Lake Balkhash indicated that not only the general trends derived from two time series were similar but that the abrupt changes and intra-annual variations of area, water level and volume variation were basically consistent in Figure 8. For instance, the time series of Lake Issyk-Kul showed overall upward shifts in area, water level and volume variation in 1989–2018. Moreover, water level of Lake Issyk-Kul declining in 1994–1998 and rising in 1998–2006 also agreed well with it in Lake Balkhash. We suspect this is cause they both in Central Asia with Lake Issyk-Kul is locate in south of Lake Balkhash with the nearest distance less than 300 km. Moreover, we speculate the Nino phenomenon may affect a wider area of this region. These geographical characteristics led them to share the same climate change which dominates the lake variations. Specific changes of Lake Issyk-Kul showed that the maximum area was 6234.32 km² in August 2006 and the minimum area was 6192.24 km² in March 1998 (see Figure 8). Besides, the volume variation ranged from +3.64 km³ (increased value) to −5.48 km³ (reduced value) could be seen from Figure 8. Song et al. [22] pointed out that the annual mean water level change has a strong correlation with precipitation and evaporation, but there may be no obvious correlation with the melting of glaciers caused by global warming. They also pointed out that water level and inundation area of Lake Issyk-Kul had piece-wise fluctuations, which largely agree with the periodical variations of precipitation [39]. In Klerx and Imanackunov [40], the water level of Lake Issyk-Kul had been dropping at a rate of about 4 cm/yr before 1998 and the water level has increased significantly at a rate of 10 cm/yr in the past eight years. Thus, they concluded this increase perhaps may be also associated with the recent acceleration of melting of glaciers associated with regional warming.