GLACIAL LAKE EVOLUTION BASED ON REMOTE SENSING TIME SERIES: A CASE STUDY OF TSHO ROLPA IN NEPAL

Himalayan glaciers have retreated rapidly in recent years. Resultant glacial lakes in the region pose potential catastrophic threats to downstream communities, especially under a changing climate. The potential for Glacial Lake Outburst Floods (GLOFs) has increased and studies have assessed the risks of those in Nepal and prioritised several glacial lakes for urgent and closer investigation. The risk posed by the Tsho Rolpa Glacial Lake is one of the most serious in Nepal. To investigate the feasibility of high-frequency monitoring of glacial lake evolution by remote sensing, this paper proposes a workflow for automated glacial lake boundary extraction and evolution using a time series of Sentinel optical imagery. The waterbody is segmented and vectorised using bimodal histograms from water indices. The vectorised lake boundary is validated against reference data extracted from rigorous contemporary unmanned aerial vehicle (UAV)-based photogrammetric survey. Lake boundaries were subsequently extracted at four different epochs to evaluate the evolution of the lake, especially at the glacier terminus. The final lake area was estimated at 1.61 km, significantly larger than the areal extent last formally reported. A 0.99 m/day maximum, and a 0.45 m/day average, horizontal glacier retreat rates were estimated. The reported research has demonstrated the potential of remote sensing time series to monitor glacial lake evolution, which is particularly important for lakes in remote mountain regions that are otherwise difficult to access.


INTRODUCTION
Glacial lakes are constantly evolving due to glacier dynamics, which are affected by climate and environmental conditions, as well as anthropogenic activities (Richardson and Reynolds, 2000). The risk of Glacial Lake Outburst Floods (GLOFs) increase when the glacier retreat rate increases as a result of rising temperatures. A GLOF typically involves a sudden flood wave descending downstream and can cause catastrophic damages to communities, resources and infrastructure (Rana et al., 2000;Richardson and Reynolds, 2000). The glaciers in the Himalayas have retreated significantly in recent decades, indicative of the impact of global climate change. As a result, glacial lakes in the region have grown rapidly, and so has the risk of GLOFs.
Since 1960s, 23 GLOF events have been recorded for Nepal. Out of the 23 events, 10 GLOFs had a trans-boundary impact as damage occurred both inside and outside Nepal, specifically in Tibet Autonomous Region (TAR) as well as China (ICIMOD, 2011). Over the last five years, two GLOF events with catastrophic trans-boundary impacts were reported. The first event occurred on the 5 th July 2016 in Zhangzangbo valley, TAR, China, causing damage along the Bhotekoshi River in Nepal (IMHE, 2016;The Kathmandu Post, 2017). The second event occurred on the 20 th April 2017 in Barun Valley, Nepal (Byers et al., 2019). These recent events demonstrate that continuous investigations to evaluate the risks of GLOFs are a necessity.
There are more than a thousand glacial lakes in Nepal. Research has been conducted to assess the risks of glacial lakes and those with the highest risks have been prioritised for further * Corresponding author investigations. 21 lakes have been identified as critical, one of them being Tsho Rolpa (ICIMOD, 2011). Tsho Rolpa is located northeast of Kathmandu in the Rolwaling Valley, Dolakha District, Nepal, at an altitude of 4,580 m ( Figure 1). The lake is the largest moraine-dammed proglacial lake in the the Nepalese Himalayas (Rana et al., 2000). It was formed only in the last 60 years or so, and its surface area and water volume have evolved rapidly (Shrestha and Nakagawa, 2014). It is currently c. 3.45 km long and 0.50 km wide, fed by the Trakarding Glacier to the southeast (ICIMOD, 2011).
The study of Tsho Rolpa began in the early 1990s (Rana et al., 2000). The lowering of the lake water level and monitoring water fluctuation were already recommended at that time. The monitoring of the lake area, depth and the growth pattern were followed, and topographic surveys of downstream villages were carried out. Later studies suggested that buried ice was present in the end moraine and near islands located in the lake. Geotechnical analysis concluded that the moraine dam was stable if no large-scale event such as an earthquake occurred, but that it will eventually fail within a timescale of 100 years if no mitigation was taken. A mitigation programme was therefore introduced in 2000, including the construction of an open channel through the end moraine. The water level was subsequently lowered by three metres and the risk of a GLOF was reduced by an estimated 20%. Nevertheless, the Tsho Rolpa Lake is still considered to be one of the most dangerous lakes in Nepal and is continuously monitored by the Nepalese government.
Research by Kargel et al. (2016) reported apparent tension cracks on the moraine dam at Tsho Rolpa that were caused by the movement (c. 1.5 m horizontally and 0.5 m vertically) of moraine materials due to the 2015 Gorkha earthquake. Another study by Shrestha and Nakagawa (2014) analysed the outburst flood scenarios of Tsho Rolpa using numerical models to simulate the failure of moraine dam due to seepage and water overtopping. However, there is a consensus that any attempt to forecast a GLOF event is extremely difficult. There is therefore an urgent need to monitor the lake on a regular basis. As the lake is located in a remote valley and at a high altitude, regular visits can prove challenging and even dangerous. This research therefore aimed to demonstrate the capability of remote sensing to monitor the evolution of the Tsho Rolpa Glacial Lake.
Remote sensing satellite images, both radar and optical, have been widely used to measure and monitor the characteristics of glaciers, to delineate the changes of a glacier's terminus position (Racoviteanu et al., 2008), as well as to detect water content on the Earth surface (Yue and Liu, 2019). With increased spatial and temporal resolutions, it has proved feasible to measure the boundary changes of glacial lakes and the glacier terminus position through time (Racoviteanu et al., 2008;Li and Sheng, 2012;Strozzi et al., 2012;Yue and Liu, 2019). Accurate glacial lake boundary estimation by remote sensing is important for at least two reasons: 1) to monitor lake evolution on a regular basis with a high level of automation; 2) to help estimate water volume given bathymetry data, which generally have to be obtained onsite but not necessarily regularly.
Regarding Tsho Rolpa Glacial Lake, some remediation measures have been applied to reduce the risk of dam failure. Due to the outlet channel constructed at the end moraine (northwest of the lake) in 2000, the water volume was reduced and does not change significantly due to constant water discharge. It has been assumed that the lake boundary does not vary significantly, especially the boundary close to the dam. This implies that the boundaries extracted from remote sensing data in recent years should not differ much. By 2000, Landsat imagery showed that the lake area was 1.53 km 2 after the water level had been reduced (ICIMOD, 2011). Since then, the lake's growth slowed, and was mainly driven by the melting of the hanging glacier at the opposite end of the lake. Satellite images acquired in 2005 (Landsat) and 2007 (AVNIR-2) indicated the area as 1.535 km 2 and 1.538 km 2 , respectively. Bathymetric surveys have been repeated several times since 1993, with a maximum measured depth of 132 m. Through field investigation, the lake is observed to be deepening and the rate increases from northwest to southeast across the end moraine, with an average rate of 0.43 m/year. Using bathymetric data, the water volume of the lake can be estimated and the reported volume in 2011 was 85.94 million m 3 (ICIMOD, 2011).
As aforementioned, open-source low resolution remote sensing data, e.g., Landsat, has already been used to estimate lake boundaries (Li and Sheng, 2012), and has been validated by sparse ground measurements. This research utilises Copernicus Sentinel satellite remote sensing imagery to extract more accurate lake boundaries at a higher spatial resolution. The results are validated by Unmanned Aerial Vehicle (UAV)-based photogrammetry and compared with previous estimations to investigate the lake evolution. Ultimately, the technique could be expanded to other glacial lakes that are not remediated but are at equally high risk to GLOFs.

METHODOLOGY
The methodological workflow consists of three main stages.
Stage 1 applies a pan-sharpening approach to the satellite image series to enhance spatial resolution individually per band and per epoch. Many remotely sensed images include a panchromatic band, which has the highest spatial resolution, and the pansharpening approach fuses this band with all other multispectral bands of significantly lower spatial resolution whilst maintaining their higher range of spectral information. There is a plethora of pan-sharpening approaches available, as described in Nikolakopoulos (2008). However, with the absence of a panchromatic band in the Sentinel-2 image series, the pansharpening approach has been substituted with other band fusion techniques (Lanaras et al., 2018). One example is the superresolution algorithm, Superres, developed by Brodu (2017). This approach extracts geometric information per pixel from all bands to refine the low bands' resolution while preserving the pixels' spectral indices. It uses the enhanced spatial resolution of visible and near infra-red bands to refine the low resolution of multispectral bands (Table 1). Superres is adopted here as a plugin of the open-source Sentinel Application Platform SNAP (2018) software package.
In Stage 2, various water indices are constructed using different combinations of the previously super resolved bands (see also Table 1 for band details). The six indices used here have been widely applied for water body extraction, as seen in Schwatke et al. (2019) AWEIsh = Blue + 2.5 * Green -1.5 * (NIR + SWIR1) -0.25 * SWIR2 (6) The NDWI index has been used for decades in remote sensing to distinguish water bodies from other categories such as land and vegetation (McFeeters, 1996). However, it has been found that NDWI provides unsuccessful classification when noisy land background pixels are mixed with water bodies (Xu, 2006). This has been overcome with the introduction of MNDWI (Xu, 2006) where the NIR band has been substituted with the SWIR1 band, as described in Equation 3. The NWI index based on Equation 2 considers more infra-red bands than the other indices. The constant multiplier equal to 100 is adopted here as in Schwatke et al. (2019) to increase the low values. Equation 4 describes a combination of the NDWI and MNDWI indices when added together. The two variants of the AWEI, AWEInsh and AWEInsh in Equations 5 and 6 respectively, have been mainly implemented on satellite images that include pixels with dark surfaces and shadows (Feyisa et al., 2014;Yang et al., 2017).
After implementing the aforementioned water indices in the SNAP software package, six raster images are produced that each of them classifies "water" and "non-water" regions. To generate a bimodal product of these two categories, the well-established Otsu's histogram thresholding method (Otsu, 1979), adopted in Matlab, is applied to each generated index. After applying the calculated threshold, six bimodal raster maps are constructed with digital values equal to 1 representing the waterbody and 0 corresponding to land (i.e. non-water region). These raster maps are then vectorised in ArcGIS to generate polyline features of the Tsho Rolpa Lake boundary (Figure 3).
The methodological pipeline is applied to a reference satellite epoch, which is temporally close to the UAV campaign conducted on 16-05-2019 at Tsho Rolpa. This was the first UAV campaign that has ever been carried out at the study site and was considered to constitute the reference dataset. An assessment process investigates an optimal index that closely fits to the lake boundary extracted from UAV-based photogrammetric process.
In order to achieve such comparison, polylines or polygons are typically divided into a set of points. Several measures have been used in previous studies to compare the similarity of a pair of polylines or polygons, such as a) the nearest Euclidean distance between a set of points from a reference polyline (Bishop-Taylor et al., 2019); b) the Hausdorff distance which corresponds to the maximum offset between two shapes (Alt et al., 2003); c) ratio of squared root of area to perimeter (Comber et al., 2003); d) sinuosity and relative sinuosity (Anderson et al., 2014).
Here, the sinuosity (S) and the relative sinuosity (RS) metrics were adopted, as expressed in Equations 7 and 8 below: where L(t) is the total length of a polyline, i.e. the cumulative length of all line segments, and L(sf) is the distance between the start and finish locations. S(ref) is the UAV-extracted polyline and S(i) the derived polylines with i from one to six, corresponding to the six indices. The two metrics here are tested for the most stable part of the lake around the outlet instead of the closed lake boundary to avoid zero sinuosity of a closed loop.
S was initially introduced as a metric in geomorphological studies for the derivation of the geometry of river channels (Chorley et al., 1984). It indicates the ratio between the meandering length of a stream and the straight line distance (McCuen, 1998). RS has been used to evaluate the degree of similarity between multiple polylines that represent line stream networks reconstructed from various datasets, such as lidar-based and airborne photogrammetry-derived point clouds (Anderson et al., 2014). Here, RS is implemented to assess the closeness of fit between reference (i.e. UAV-based) and Sentinel-derived polylines of the lake. To calculate S and RS, the Hawth's analysis open-source line metric toolbox (Hawth, 2020), embedded in ArcGIS, was adopted.
Stage 3 investigates the lake's boundary evolution and potentially predict the glacier's movement rate. An open source toolbox called Digital Shoreline Analysis System (DSAS; Himmelstoss et al. (2018)) has been initially developed to assess temporal coastal changes by the U.S. Geological Survey (Gibbs and Richmond, 2017). Recent research by Yue and Liu (2019) implemented the DSAS tool to estimate annual lake boundary variations in China using Landsat imagery. This tool, embedded in ArcGIS, can calculate temporal positional change of a boundary providing change rate statistics using a linear regression analysis. The DSAS tool can also estimate changes within a specified temporal horizon based on the historical movement rates with the aid of Kalman Filtering.

SENTINEL AND UAV DATASETS
Sentinel-2 multispectral imagery (MSI) of Level-2A was acquired from the Copernicus Open Access Hub at four epochs 18-05-2019, 12-06-2019 17-07-2019 and 16-08-2019. The first epoch was the closest (two-day separation) to the reference UAV campaign with zero cloud coverage over Tsho Rolpa (Figure 1). To test the methodology a difference of approximately 30 days between epochs was selected. The Sentinel-2 bands, as listed in Table 1, were super-resolved yielding in a 10 m spatial resolution, after implementing the plugin Superres in SNAP, according to Stage 1 of the methodology. To reduce computational effort, a subset region of 4.2 km x 3 km out of the original 100 km x 100 km Sentinel image tile was extracted and utilised at all methodological stages. The false colour infrared composite, shown in Figure 1, was constructed based on the B8, B4 and B3 band combination.  In the UAV campaign (Figure 2a), 448 images were captured from a fixed-wing UAV (eBee; SenseFly (2019)), covering an area of approximately 7.064 km 2 . The eBee carries a Sony Cybershot DSC-WX220 lightweight digital camera with a 4.572 mm nominal focal length, and a 1/2.3" (7.76 mm) Exmor R TM CMOS sensor. The DSC-WX220 camera creates an image of 4896 x 3672 pixels corresponding to a chip size of 6.17 x 4.63 mm.
To ensure that UAV-photogrammetric output is georeferenced into a common fixed reference frame (WGS84/UTM Zone 45N), seven ground control points (GCPs) were surveyed using GNSS rapid static mode with eight-minute observations per point (Figure 2). An average 3D relative accuracy of 0.005 m was estimated after post-processing with GNSS observations from a base station located at a relatively stable terrain close to the outlet channel. Average root mean square errors (RMSEs) equal to 0.025 m, 0.082 m and 0.067 m were calculated in X, Y, Z axes respectively, using the Structure-from-Motion (SfM) photogrammetric software package Pix4D (2016). To achieve this the seven GCPs were incorporated in the Pix4D selfcalibrating bundle adjustment. The resulting digital elevation model (DEM) and orthomosaic were reconstructed with 0.120 m spatial resolution.

RESULTS AND ANALYSIS
After conducting Stage 1 and 2 of the methodological workflow, six polylines of the lake boundary were constructed at the first Sentinel epoch, as shown in Figure 3. To evaluate the six water indices, only the area depicted in Figure 3 was used for this test, instead of the full lake boundary. This area surrounds the outlet channel (end moraine; northwest) with a distinctive separation between lake's water and land where no mixed noise is formed by snow or icy regions. Icy areas were observed at the glacier terminus in the UAV orthomosaic, which cannot clearly be identified on Sentinel images. Moreover, the boundaries around the island shown in Figure 3 were not taken into consideration here, in order to simplify the polyline assessment task. It should be noted that the freely-available Sentinel-2 Level-2A products are georeferenced in WGS84/UTM Zone 45N and are both atmospherically and radiometrically corrected. This allowed direct application of the water indices without a further reflectance correction procedure.
By visual inspection, it is evident that some indices are relatively similar to each other, such as the AWEInsh and AWEIsh, as well as the NDWIplusMNDWI with the NDWI. The step-wise pattern in all polylines is caused by the vectorisation procedure. The size of the step is equal to 10 m spatial resolution of the superresolved Sentinel bands after applying the Superres algorithm. The S and RS indices were calculated without simplifying the polylines further, and reported in Table 2.
When RS is closer to unity the polyline fits better to the UAV reference polyline. The smallest RS value of 0.81 indicates that NWI provides the worst fit among all, also shown in Figure 3. This is because NWI raster map did not provide a sharp bimodal histogram and Otsu's method did not provide successful results as it did for the other indices. According to Nakmuenwai et al. (2017), Otsu's method is limited when histograms are unimodal rather than bimodal. The NDWIplusMNDWI provided one of the clearest bimodal histograms with sharp peaks and a deep valley between them. For that reason, among those indices with RS values greater than 0.90, the NDWIplusMNDWI was chosen for use in lake boundary extraction from the Sentinel-2 image series.  An example of a sharp bimodal histogram of the NDWIplusMNDWI index is shown in Figure 4 for the final Sentinel epoch on 16-08-2019. Threshold from Otsu's method was calculated equal to 0.4922. This indicates that pixels with values greater than the threshold are part of the lake while any other pixel belongs to land. Reducing the original size of Sentinel image tiles ensured that only Tsho Rolpa Lake and its surroundings were included in the workflow with pixels indicating water, rocky surface and sediments.
The aforementioned workflow was repeated for the four Sentinel-2 image-series, as seen in Figure 5a. The vectorised polylines, representing the lake's boundaries, were generated with the NDWIplusMNDWI index and are depicted by different colours in Figure 5a. The lake's boundary on 18-05-2019 (shown in yellow in Figure 5a and b) coincides with the boundary seen on the UAV-constructed orthomosaic. It should be noted that UAV imagery was captured early in the morning to prevent sun reflections in the imagery. Because of this, some parts of the lake, mainly those close to the glacier terminus, were glacier ice covered by fresh snow, as seen in Figure 5a. However, the snow did not affect the boundary extraction from the Sentinel-2 image. The lake's boundaries in June, July and August (shown in Figure  5a in red, black and magenta, respectively), follow the cracks observed on the glacier terminus from the UAV orthomosaic. To evaluate the variations of the lake's boundaries using the DSAS toolbox, a baseline was manually generated (see Figure 5b in green). DSAS calculates the distances between the baseline and each polyline intersection point along transects (dashed line in Figure 5b). These transects were constructed every 20 m apart. DSAS applies a linear regression analysis to estimate a change rate per transect. For instance, Figure 6 demonstrates the analysis for the transect AB shown in Figure 5b. A linear regression line is determined by fitting a least squares regression to the distances shown in y-axis in Figure 6. As a result of the linear regression, the slope of the estimated line equals -0.99 and corresponds to the change rate at this particular location. The negative sign implies erosion, i.e. glacier retreat. It was estimated that the maximum horizontal glacier retreat of 0.99 m/day occurred in that location (transect AB). DSAS estimated an average horizontal glacier retreat of 0.45 m/day across all transects, as shown in Figure 5b.
The surface area and perimeter of the lake were calculated for three of the epochs and reported in Table 3. Due to cloud coverage it was not possible to extract the entire lake boundary at the final epoch in August 2019. The lake had a surface area of 1.538 km 2 in 2007 as reported in ICIMOD (2011). Since then the lake has grown with a difference of 0.074 km 2 . The monthly variations of the calculated area are possibly attributed to the glacier retreat during the monsoon period. However, a decrease in area with an increase in perimeter between May and June 2019 was observed (see Table 3). This is caused by the fact that the 17-07-2019 lake boundary was underestimated as few parts in NE of the lake were categorised as water instead of land. A comprehensive investigation to further minimise errors due to the sub-pixel segmentation and vectorisation procedure (e.g. Bishop-Taylor et al. (2019)) is essential. The aforementioned regression analysis estimated the monthly variations of the lake's boundaries, as seen in Figure 5b. As expected, the lake boundary changes mainly over the glacier terminus (SE) and not over the end moraine (NW) (Figure 1). This is probably due to the constructed dam in 2000 around the natural water outlet at the end moraine. However, to support this statement a further investigation of the lake's boundary from historical satellite images would be necessary.   Table 3. Lake surface area and perimeter per epoch.

CONCLUSIONS
The presented research proposes a workflow for waterbody boundary extraction and evolution in an automated fashion using Sentinel optical imagery at Tsho Rolpa Glacial Lake, Nepal. The research demonstrates the potential of remote sensing time series to monitor glacial lake evolution, which is particularly important for those remote lakes which are difficult to access. The research can be expanded to analyse monthly and annually boundary change rates from Sentinel image series beyond the year 2019 starting from the beginning of investigations into the Tsho Rolpa Glacial Lake in the early 1990s. In this case data sources other than Sentinel satellites should be considered. Moreover, the research has the potential to inform forecasting of glacier retreat rate into the near future.
Even though Sentinel Level-2A products are already georeferenced, further analysis is necessary to ensure that no systematic error will propagate through the workflow in the satellite time series. The water/land segmentation and vectorisation procedure will be further improved, with the aid of edge detectors, convolution kernels and smoothing algorithms (e.g. Liu and Jezek (2004)), to enhance the bimodal histograms and automatically generate a smooth lake boundary. Further considerations to identify mixed pixels of snow, ice and sediments around the lake's boundary would also enhance the segmentation process.
In parallel with the sinuosity metric, alternative approaches to compare the reference lake boundary extracted from UAV-based photogrammetry against the boundaries extracted from various Sentinel-2 indices are also being investigated (e.g. nearest Euclidean distance, Hausdorff distance, etc.). Whilst the workflow currently estimates the temporal planimetric variations of the lake's boundary, ultimately the work could be extended to monitor water volume variations by incorporating lake level observations from gauges and bathymetric measurements.