Automatic watershed delineation in the Tibetan endorheic basin: A lake-oriented approach based on digital elevation models

Geomorphology 358 (2020) 107127

Contents lists available at ScienceDirect

Geomorphology journal homepage: www.elsevier.com/locate/geomorph

Automatic watershed delineation in the Tibetan endorheic basin: A lakeoriented approach based on digital elevation models Kai Liu a, Chunqiao Song a,⁎, Linghong Ke b, Ling Jiang c,d, Ronghua Ma a a

Key Laboratory of Watershed Geographic Sciences, Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, Nanjing 210008, China School of Earth Sciences and Engineering, Hohai University, Nanjing 211100, China Anhui Center for Collaborative Innovation in Geographical Information Integration and Application, Chuzhou University, Chuzhou 239000, China d School of Architecture, Building and Civil Engineering, Loughborough University, Loughborough LE11 3TU, UK b c

a r t i c l e

i n f o

Article history: Received 2 March 2019 Received in revised form 22 September 2019 Accepted 26 February 2020 Available online 27 February 2020 Keywords: Watershed delineation Endorheic basin Lake classification Digital elevation models

a b s t r a c t Watershed delineation from digital elevation models (DEMs) is an essential prerequisite for land-surface hydrology studies. Although existing automatic methods perform well in exorheic basins, finding a suitable watershed delineation approach in endorheic basins is still challenging. This study develops a lake-oriented approach given that lakes are important hydrological components in the Tibetan endorheic basin. The designed workflow of the proposed approach contains four steps: flow direction calculation, hydrological structure features extractions, classification of lake types, and watershed labeling. The experiments are conducted on a sub-region of the inner Tibetan Plateau, with an area of 119,055 km2. Extraction results prove that the proposed strategies, i.e., elevation-profile judgment and overflow-point-control correction can support the accurate delineation of watershed boundaries. In addition, the proposed approach can achieve two-level watershed boundaries in which the independences of all the endorheic basins can be maintained, meanwhile, further delineation within extracted exorheic basins satisfies various application requirements at different spatial scales. Unlike manual processing and existing datasets (e.g., HydroSHEDS), the proposed approach can fully explore the advantages of high-quality DEMs and can be expanded to much larger regions to handle scientific and engineering issues. © 2020 Elsevier B.V. All rights reserved.

1. Introduction Watershed delineation is a fundamental task serving for multidomain studies, such as hydrological modeling (Liu et al., 2014), geomorphic mapping (Evans, 2012; Xiong et al., 2018; Liu et al., 2018), soil erosion assessment (de Vente et al., 2013), and water resource management (Zhang et al., 2016). With the growing ability to access digital elevation models (DEMs) from various sources, numerous methods have been proposed to extract the drainage networks and watershed boundaries automatically (O'Callaghan and Mark, 1984; Tarboton et al., 1991; Freitas et al., 2016). Accurate watershed delineation largely depends on the quality of DEM data and the effectiveness of the methodology being used (Wilson, 2012). Recently, remote sensing techniques, such as, satellite and airborne photogrammetry, airborne/terrestrial Light Detection and Ranging, and airborne/satellite radar systems that allow the collection of abundant terrain information at different scales, have developed rapidly, thereby enabling in-depth investigations of the Earth's surface processes (Rabus et al., 2003; Liu, 2008; Tarolli, 2014; Liu et al., 2016). A considerable number of airborne or spaceborne DEMs, including the ⁎ Corresponding author at: 73 Beijing East Street, Nanjing, China. E-mail address: [email protected] (C. Song).

https://doi.org/10.1016/j.geomorph.2020.107127 0169-555X/© 2020 Elsevier B.V. All rights reserved.

Shuttle Radar Topography Mission (SRTM) DEM, Advanced Spaceborne Thermal Emission and Reflectance Radiometer (ASTER) Global DEM, and Advanced Land Observing Satellite (ALOS) World 3D 30 m resolution (AW3D30) DEM, are freely available at near-global scale. These DEMs are highly capable and have been used widely by the geoscience community (Yue et al., 2017; Yamazaki et al., 2017; Brun et al., 2017; Yao et al., 2018). In addition to the acquisition of DEM data in high quality, various improved methods of extracting stream networks and watershed boundaries have been intensively studied. At present, the most popular watershed delineation method is based on the classical algorithm of D8 flow routing, which has been used in commercial software, e.g., ArcGIS. Its crucial step is the removal of local pits or sinks caused by data errors or DEM resolution. Although there is a controversy about the pit-filling procedure, which removes true landscape depressions (Grimaldi et al., 2007), this step is still adopted widely to ensure the proper delineation of drainage networks. After obtaining the nonpits DEMs, the flow direction can be calculated. In the D8 algorithm, each cell can deliver eight possible flow directions to adjacent cells based on of the maximum gradient. Then, the accumulated flow is stored in an output raster that represents the accumulated number of cells flowing into each down-slope cell (Jiang et al., 2013). After the output raster of flow accumulation is generated, the watershed and

2

K. Liu et al. / Geomorphology 358 (2020) 107127

sub-watershed pour points can be defined, and the watershed boundary can be finally determined. Although existing automatic algorithms, including D8, have proven effective in the watershed delineation of many regions worldwide, these algorithms cannot be directly applied to endorheic basins where filling depressions can produce large uncertainties in extracted networks and watershed boundaries; thus, they are not recommended (Khan et al., 2014). The inner Tibetan Plateau (TP) is a typical endorheic basin with extensive distribution of lakes, and the internal drainage systems have no hydrological connections with marine environment (Dorsaz et al., 2013). The existing workflow fails to generate reasonable watershed boundaries for this region because of difficulties in identifying natural and invalid sinks. Thus, many misinterpretations have been made with regard to river networks and drainage relations in endorheic basins for previous global-scale hydrological drainage datasets, e.g., the HydroSHEDS (https://hydrosheds.cr.usgs.gov/index.php) and the Global Drainage Basin Database (GDBD) (http://www.cger.nies.go.jp/ db/gdbd/gdbd_index_e.html) (Phan et al., 2013; Ai et al., 2017; Gao et al., 2018). Several recent efforts, such as Phan et al. (2013) and Gao et al. (2018), have improved drainage networks in the inner TP on the basis of a HydroSHEDS dataset by utilizing historical surveying topographic maps or referring to the linkage relationship between lake and supplying glaciers to correct the misclassified drainage basins. To broaden the applicability for endorheic basins at a low time cost, an automatic approach that is independent of manual experience and local auxiliary materials needs to be developed. This study aims to develop an automatic watershed delineation approach for endorheic basins, which is still lacking and challenging. Given that lakes are the key hydrological components of endorheic drainage systems, the hydrological connections between lakes were analyzed to provide important references for watershed delineation. The proposed lake-oriented approach can overcome the limitations of the existing watershed extraction method, which is unsuitable for endorheic basins. The remainder of this study is organized into the following sections. Section 3 describes the study area and data. Section 4 introduces the general framework of the proposed method. The results and discussions are presented in Sections 5 and 6, respectively. Finally, several remarks are summarized in Section 6.

2. Study area and dataset 2.1. Study area The inner TP is located in the hinterland of TP, covering an area of about 708,000 km2 (Fig. 1). The climate types in this area range from temperate arid to semi-humid, and are dominated by dry westerlies and are limitedly influenced by the Indian monsoon (Song et al., 2013). Although the TP is known as the “the roof of the world”, the inner TP is broad and relatively less rugged compared with exorheic basins. The exorheic drainage basins can be divided into eight zones, including the Amu Darya, Yangtze, Mekong, Salween, Brahmaputra, Ganges, and Indus Rivers. Flat-topped areas and topographic barriers form a cluster of endorheic basins. As shown in Fig. 1a, the inner TP is the endorheic drainage area that has independent hydrological networks. The lakes in the inner TP are important hydrological components of alpine environments, which are sensitive to climate change (Yang et al., 2011). The TP contains approximately 1200 lakes larger than 1.0 km2, and more than 60% of these lakes are located in the inner TP. These lakes are a unique indicator of climate change due to few disturbances from human activities (Zhang et al., 2017a; Song et al., 2014). To develop an automatic watershed delineation approach in endorheic basins, a sub-region in the southeastern part of inner TP with varied terrains (elevations range from 4492 m to 7040 m) and lake sizes were selected as the test area, covering an area of 119,055 km2. The two largest lakes of Tibet, namely, Lake Siling Co and Lake Nam Co, are included in the study area (Fig. 1c). 2.2. Study data 2.2.1. DEMs DEMs are the basic data for watershed delineation and provide the topographic information with varying resolutions. Two freely available DEMs, SRTM-1 DEM, and AW3D30 DEM were selected. SRTM mission collected the land elevations by using airborne radar measurements on February 2000 and provided near global-scale DEM products at 1″ and 3″ resolution. After void-filling processing, the global SRTM-3 V4.1 DEM with 3″ resolution (approximately 90 m) and the 1″ SRTM-1

Fig. 1. Maps of study area. (a) The hydrological zone division of the Tibetan Plateau; (b) the location of the study area within the Inner TP and (c) topographic map of the study area.

K. Liu et al. / Geomorphology 358 (2020) 107127

DEM (approximately 30 m) have already been open. AW3D30 was generated using optical stereo matching technologies. Images were acquired by the Panchromatic Remote-sensing Instrument for Stereo Mapping onboard the ALOS and publicly released by the Japan Aerospace Exploration Agency in 2015 (Tadono et al., 2015). 2.2.2. Lake extent mapping results Lakes serve as the important hydrological components within their drainage basins. Lake mapping is relatively simple on an individual basin but is difficult when performed for numerous lakes across large areas and temporal scales (Pekel et al., 2016). In this study, lake water extent mapping was conducted using multi-temporal Landsat-7/8 imagery (Sheng et al., 2016). We generated the composite circa-2000 and circa-2015 lake maps covering the study area. The circa-2000 lake map and the SRTM-1 DEM were used to produce the drainage basins around 2000, and the AW3D30 DEM combined with the circa-2015 lake map were used to produce the circa-2015 drainage basins. 2.2.3. HydroSHEDS data HydroSHEDS provides global-scale hydrographical features, and is widely used in previous studies (Phan et al., 2013; Gao et al., 2018). Existing data improvement methods and newly developed algorithms have been applied in its workflow, including void-filling, filtering, stream burning, and upscaling techniques. Manual corrections were also conducted where necessary (Lehner et al., 2013). In this study, the HydroSHEDS river network (as_riv_15s) and the drainage basin (as_bas_15s) with a spatial resolution of approximately 410 m were retrieved from the USGS website (http://hydrosheds.cr.usgs.gov) and used as a reference data for comparison with the extracted results. 3. Method 3.1. Basic idea Watersheds, which are also known as drainage basins or catchment, are defined as the unit area where all the surface water from rain runoff, glacier, or snow meltwater converges to a single point or outlet to join another water body, such as a lake, river, or ocean (DeBarry, 2004). Each watershed is separated topographically from adjacent basins by a drainage division. Generally, watersheds in exorheic drainage systems drain into other watersheds through low-elevation outlets in a hierarchical pattern. However, the hydrological connection in an endorheic basin is different, with surface flows being trapped in a lake or depression without outlets. Fig. 2 demonstrates a typical endorheic basin of the inner TP which has no outlet in the basin boundary. Thus all the surface water flows within the basin cannot drain into an adjacent basin. Instead, the surface water

3

flows converge into lakes, which are classified into two types: terminal lake which has inlets but no outlets, and upstream lake which has outlets linking with terminal lake or another upstream lake. A terminal lake is usually located at the lowest elevation of the endorheic basin and is the terminus of the surface flows (Dorsaz et al., 2013). Hence, an endorheic basin may have numerous upstream lakes but only one terminal lake. Geographical interpretation of the differences between exorheic and endorheic drainage systems is the basis for developing delineation algorithms. The algorithm for extracting hydrologic networks in an exorheic drainage area should satisfy two requirements: (1) every grid cell within a DEM must have a definite flow direction and (2) the edge of the DEM can be reached by following the flow directions (Barnes et al., 2014a). As for the endorheic drainage area, e.g., the inner TP, not all cells can drain into the edge of the DEM due to terrain barriers. The traditional methods of flow direction and watershed delineation are not suitable because they cannot ensure the independence of the endorheic basins. This study proposes a lake-oriented approach for watershed delineation in the inner TP. In this approach, watershed delineation is divided into a twolevel process. The first-level delineation aims to identify the catchment areas for each terminal lake. On this basis, second-level delineation is conducted for endorheic basins with areas larger than the research needs that warrant further delineation. The workflow of the proposed algorithm includes four steps, as illustrated in Fig. 3: flow direction calculation, extraction of hydrological structure features, classification of lake types, and final watershed labeling. 3.2. The procedures of algorithm workflow 3.2.1. Flow direction calculation Flow direction is a basic hydrologic parameter and a prerequisite for watershed delineation. The widely used method is to generate the modified DEM by filling all the depressions to the elevation of their outlets (Jenson and Domingue, 1988). However, the depression-filling algorithm increases the sizes and numbers of flats in a DEM which is not suitable for processing endorheic basins. As for the commonly used D8, the large distribution of flat areas caused by the depression filling process and truly flat landscapes may result in the parallel flow problem. To ensure that the flow direction assignment is realistic and topographically consistent, some improved algorithms have been proposed (Garbrecht and Martz, 1997; Martz and Garbrecht, 1998; Zhang and Huang, 2009; Wang and Liu, 2006; Zhang et al., 2017b). In this study, we propose a strategy that integrates several flow direction algorithms (including Priority-Flood, G&M, and D8) to process different terrain units at a large scale. As for the depression area, the Priority-Flood algorithm was adopted which treated the surface depressions as individual ponds or reservoirs and to simulate how they are flooded (Wang and

Fig. 2. A typical endorheic basin in the inner TP. (a) The imagery acquired from Google Earth and (b) topographic maps representing the drainage modes.

4

K. Liu et al. / Geomorphology 358 (2020) 107127

Fig. 3. Workflow of the proposed algorithm.

Liu, 2006; Metz et al., 2010; Barnes et al., 2014a; Xiong et al., 2019). As for the flat areas, the directions were determined based on the gradient away from higher terrain and the gradient towards lower terrain, as noted in the G&M algorithm (Garbrecht and Martz, 1997; Barnes et al., 2014b). The remaining cells that were not located in depression, and flat areas can be determined by the D8 algorithm. 3.2.2. Extraction of hydrological structure features Several hydrological structure features, including source and connection points, drainage networks, and initial watershed delineation results, are required to determine the accurate watershed boundaries of each terminal lake. Flow accumulation and flow direction calculations are the prerequisites for the extraction of hydrologic structure features. After the catchment area for each cell is obtained, the drainage networks can be easily generated on the basis of a given drainage threshold value which is defined as the minimum upstream area that constitutes a stream. In this study, a low drainage threshold is recommended to ensure that all lakes can be connected with drainage networks. If a large threshold value is used, some lakes may not be connected by drainage networks. In this case, whether the unlinked lake is an upstream lake or the terminal lake cannot be determined. The subsequent step is the creation of an initial watershed delineation that corresponds with the extracted drainage. The areas of the initial watersheds are very small, which are suitable for further processing to generate accurate watershed boundaries. To determine the lake types and investigate the spatial relationships between lakes and watersheds, stream source points and connection points are needed. Stream source points are the starting cells of drainage networks that have a sufficiently large flow contribution area. On the basis of the flow directions and drainage networks, the cell is labeled as the source point if it is a part of streams while all its upstream cells are not. The connection point is associated with the

lake body extent. If two lakes are connected by the drainage networks, then the points located at both ends of the connection path are defined as the connection points. 3.2.3. Classification of lake types In accordance with the definitions of terminal and upstream lakes, the drainage networks are the important reference data for judging the lake types. As mentioned in Section 3.2.1, the drainage networks are based on the flow direction matrix. Compared with D8, the integrated algorithm can overcome the parallel flow problem in flat areas. However, this algorithm is still not specially designed for endorheic basins, in which each cell should be drained into the edge of a DEM. In this case, the generated flow direction may be inconsistent with terrain relief in some cells. The connected lakes through extracted drainage networks should be further investigated to avoid spurious connections. In this study, we developed an elevation-profile strategy to judge whether the hydrological connection between two lakes is true or not. The strategy is based on the assumption that the elevation profile along the connected flow paths would decrease monotonically if two lakes are naturally connected. If the elevation profile does not approach to the monotonic pattern, the connection relationships between the two adjacent drainage basins could be spurious. More specifically, the elevation-profile judgment was conducted as follows. First, if a cell is the connection point and outlet of the lake, then tracking analysis is conducted to record the elevation values along the connection path until the paired connection point is reached. However, our testing analyses showed that the original DEM value for each cell may be not suitable for constructing the elevation profile because DEM elevations may be influenced by local noisy errors. In this study, we used the median elevation value

K. Liu et al. / Geomorphology 358 (2020) 107127

within a 5 × 5 window size to filter out DEM outliers. Second, the Gaussian filter method was applied to smoothen processing and obtain a fitting curve of elevation profile. The derivative value of the fitting curve was calculated to identify the rise segments. The trend is increasing when the derivative value is negative and decreasing otherwise. Hence, the turning point of derivative value from positive to negative is considered a local maximum. By contrast, the turning point from negative to positive is considered a local minimum. The points between local minimum and local maximum all belong to the rising segments. After all the rising segments for one elevation profile are determined, the final step is to analyze whether the inter-lake profile is monotonic or not. In consideration of the influence of the DEM error, not all the rise segments were caused by the existing drainage divides. In this proposed approach, the standard deviation of the DEM vertical error was recommended as the threshold value for monotonic judgment. The rise segments is not considered valid unless the magnitudes of rise value is greater than the height threshold value. According to existing vertical accuracy assessment of open-access DEMs in the inner TP (Liu et al., 2019), the standard deviations of SRTM-1 DEM and AW3D3D DEM were 8.67 and 8.53 m, respectively. Hence, the two values were used as the threshold values to screen out invalid connections. Three typical examples are shown in Fig. 4 to illustrate the elevationprofile judgment. In the first example, the segments can be easily excluded, given its low rising amplitude despite the initially observed rising segments. Hence, these two lakes are connected (Fig. 4a). A disconnected example is shown in Fig. 4b. The elevation profile exhibits a rising trend with a large value. In this case, the interconnection between two lakes can be considered invalid because of the special

5

process of the flow direction algorithm in the flat area. A complex example is shown in Fig. 4c, where three rising segments were detected from the elevation profile. In this case, the maximal rise values were recorded for comparison with the threshold value. Because the maximum rise value is smaller than the height threshold value, two lakes are considered connected. As for the invalid connection, an overflow-point-control correction was designed to correct the flow direction values, which are necessary for the following watershed labeling. As shown in Fig. 5a and b, the connection path, in which the pixel with local maxima cloud be detected, is invalid. Such a point can be regarded as the overflow point of two neighbored watersheds. It is acknowledged that surface water should flow away from the local maxima instead of draining into it. Hence, the pixel draining into the overflow point is labeled as a starting point from which reverse processing of flow direction is conducted until the lake boundary is reached (Fig. 5c and d). After obtaining the corrected flow directions, the corresponding drainage networks and stream source points can be generated (Fig. 5e and f), which are applied in further processing. The elevation-profile judgment can effectively identify the “lake-tolake” relationship. However, if two lakes are too close to generate a meaningful elevation profile, then a direct strategy that compares lake surface elevations is applied. As shown in Fig. 6, the connection paths are all excessively short. Given this characteristic, these two lakes may be connected. However, the flow relationship should be determined. In this study, the medium value of the lake surface elevation is calculated. The lake with a high medium value is considered upstream one. Hence, the left lake in Fig. 6a and the right lake inFig. 6b are determined as upstream lakes.

Fig. 4. The schematic diagram for determining the drainage modes based on tracking elevation profiles. (a) An example for the connected lake pairs; (b) An example for the unconnected lake pairs and (c) an example for the uncertain situation.

6

K. Liu et al. / Geomorphology 358 (2020) 107127

Fig. 5. The schematic diagram of flow direction corrections for invalid connections. (a) The connection path is generated based on the original flow direction value. (b) The elevation profile indicates that it is an invalid connection and the overflow point located in drainage divide can be detected. (c) The original flow direction value along the connection path. (d) The corrected flow direction values. (e) The drainage networks generated from the corrected flow direction values. (f) Enlarged map of the overflow point and the drainage divide.

3.2.4. Watershed labeling Watershed delineation with the lake-oriented algorithm was conducted by a two-level process. In the first level, the catchment area for each terminal lake was identified. The independence of each endorheic basin should be maintained at this level. This study develops a flowpath tracking method for watershed labeling on the basis of lake types. The tracking analysis was conducted from the selected source point along with the drainage network until the lake was reached (Fig. 7a). If the lake is a terminal lake, then all the passing watersheds are labeled as the lake ID. If the lake is an upstream lake, then the setting value is replaced by the ID of its corresponding terminal lake. After all the source points were processed, the catchment area for each terminal lake was obtained (Fig. 7b) and the final watershed boundaries can be determined (Fig. 7c). The first-level watershed delineation process may yield several large basins, especially for some large lakes. Hence, further watershed delineation should be conducted independently within each endorheic basin to provide watershed boundaries at different spatial scales. The second-level watershed delineation process is similar to the traditional watershed delineation method, which mainly depends on given

drainage threshold values. Generally, the drainage threshold value should be less than the area of the extracted endorheic basin to ensure a successful delineation. Otherwise, the basin will not be merged with any other adjacent basin. 4. Results 4.1. Hydrological structure features Flow direction calculation is the basis for generating hydrologic structure features. To compare flow directions generated from the integrated algorithm and the traditional D8 algorithm, the extracted drainage networks, and the corresponding watersheds in the selected region are shown in Fig. 8. Obvious differences can be observed mainly in the flat region. The results generated by the D8 algorithm show severe parallel flow problems, which can be overcome by using the integrated algorithm. The reasonable hydrological networks in the flat region provide the basis for the subsequent steps that aims to utilize hydrological structure features. Similarly, the initial watersheds extracted based on the D8 algorithm are evidently unreasonable in terms of shape and

Fig. 6. The schematic diagram for determining the hydrological connection based on lake levels. (a) According to the medium value of lake surface elevations, the left lake with a higher surface elevation should drain into the right lake. (b) The right lake is determined as the upstream lake due to its higher lake level.

K. Liu et al. / Geomorphology 358 (2020) 107127

7

Fig. 7. Steps of watershed assignments in the proposed approach. (a) Tracking analysis is conducted from the source points and all passing watershed is set to 26. (b) The catchment area can be determined after processing all the source points. (c) The final extracted results.

size (Fig. 8c). The improved watersheds are suitable for further processing on the basis of the hydrological connection (Fig. 8d). The hydrological networks covering the study area are shown in Fig. 9a, and an enlarged area is selected to represent more detailed information, including source points and connection points (Fig. 9b, c). The overall features of the extracted objects based on two datasets (SRTM1, and AW3D30 DEMs) are similar at a large scale. As such, the two datasets are compared only in terms of the enlarged area. Although the two DEMs are derived from different remote sensing sources and processing methods, the same resolution level maintains great similarity in the extracted features. The distribution differences of drainage

networks and stream source points are caused by the different local topographic relief, which has little influence on the overall pattern. The lake extent change also leads to the position differences of the connection points. 4.2. Lake classification results Distinguishing terminal lakes from upstream lakes is the key step in the proposed approach. The lake type classification results are significant for understanding the hydrological characteristic of lakes in the inner TP. In general, the total lake area in the study area rapidly increased from

Fig. 8. Comparison between the traditional D8 algorithm and the integrated flow direction algorithm. (a) The calculated hydrological networks based on D8 algorithm; (b) The calculated hydrological networks based on the integrated flow direction algorithm; (c) The initial watershed delineation results based on D8 algorithm; (d) The initial watershed delineation results based on the integrated flow direction algorithm.

8

K. Liu et al. / Geomorphology 358 (2020) 107127

Fig. 9. (a) The hydrological networks extracted from SRTM-1 DEM and the Comparison in an enlarge area between (b) SRTM-1 DRM and (c) AW3D30 DEM.

8407 km2 to 9524 km2 during the period from 2000 to 2015. The accelerated lake growth in the 2000s was driven primarily by changes in precipitation and evapotranspiration (Song et al., 2014). Through an analysis of the elevation profile of connected flow paths, 47 terminal and 94 upstream lakes were identified in 2000, and 46 terminal and 85 upstream lakes were detected in 2015 (Fig. 10). Lake number decreased from 141 to 131 because some lakes were merged into a large one with the lake expansion. Although the general tendency of lake extent is expansion, several lakes shrank during this period due to different supply sources. The two largest lakes, i.e., Siling Co and Nam Co, are labeled as terminal lakes. Notably, no relationship between area and the type of lake is observed. Some small lakes are identified as terminal lakes because of local depressions within the large endorheic basin. To evaluate the accuracy of the lake classification results, validation assessment is conducted in reference to the dataset produced by Gao et al. (2018) and the dataset obtained via manual identification. To evaluate the reliability of our proposed algorithm, the lake types were also checked manually from high-resolution Google Earth images and DEM data; thus, the manual dataset can be considered as the correct results. It should be noted that the number of extracted lakes and that of the reference dataset is not the same due to different data sources and data acquisition time. Hence, the comparative analysis is conducted on the lakes that appear

in both datasets. In total, 124 lakes are selected for validation assessment. The proposed approach has better classification accuracy than the reference approach (Table 1). Five upstream lakes and four terminal lakes were misidentified in the reference dataset. Two kinds of improvements are shown in Fig. 11, the connection path between Lakes 42 and 44 has a ridge in accordance with the elevation profile. On this basis, Lake 42 cannot drain into Lake 44. Hence, Lakes 42 and 44 are two terminal, unconnected lakes. Another corrected example is shown in Fig. 11b. The declining trend of the elevation profile indicates that Lake 96 is an upstream lake that drains into Lake 92. 4.3. Watershed delineation results Following the principle of the proposed method, the delineation of watersheds in the inner TP is a two-level delineation process. The first-level delineation determines the catchment area for each terminal lake. Two different DEMs were used to generate the circa-2000 and circa-2015 watershed boundaries. Consistent with the number of terminal lakes, 47 watersheds were identified in 2000, and 46 watersheds were detected in 2015. Siling Co is the largest watershed, with an area of more than 45,000 km2, whereas the smallest watershed occupies an area of only 25 km2 (Fig. 12a, b). It is clear that the two sets of extracted

Fig. 10. The lake classification results based on circa-2000 lake product (over 1 km2).

K. Liu et al. / Geomorphology 358 (2020) 107127 Table 1 The accuracy assessments of the extracted dataset compared with the reference dataset provided by Gao et al., (2018) and the manually identified dataset. Type

Reference dataset Extracted dataset Manual dataset

Upstream lakes

Terminal lakes

Total

Correct

Wrong

Total

Correct

Wrong

80 79 79

75 79

5 0

44 45 45

40 45

4 0

watersheds represent great similarity due to the same spatial resolution and high quality. However, the uncertainty in watershed delineation caused by DEM quality should be considered before its application, especially for small-scale studies. Lake dynamic change is another factor that has a potential impact on the watershed delineation. Generally, the watershed boundary is controlled by the landforms, which are difficult to change. However, as the water level of most Tibetan lakes were rising evidently during the past two decades, the lake surface would be likely to exceed the elevation of the drainage divide; in this case, the adjacent two watersheds are merged into a large one. Although the influence of lake change is not evident in this study, this factor should be emphasized when the analysis is conducted in a large region or over a long term. To reflect the improvements of our proposed method, the HydroSHEDS data are shown in Fig. 12c. Compared with HydroSHEDS data, the performance of our method is improved in three ways. First, our method correctly judges the hydrological connection between lakes, such as the Siling Co basin; the catchment area from HydroSHEDS

9

is 28,978.40 km2, which is much smaller than the results of the proposed approach (SRTM-1: 45,516.24 km2, AW3D30: 45,454.52 km2). The inaccurate delineation divides the Siling Co basin into two parts. However, Qiagui Co, which is in west of Siling Co, is connected with Siling Co according to the elevation profile of connected flow paths. The high-resolution satellite images and the field trip to this region during the summer of 2019 also confirmed the interconnection. In this case, the outlined boundary in Fig. 12c should be a false edge. The second improvement is that our method correctly allocates some sub-watersheds. One terminal lake basin in north of Siling Co is selected as an example. As shown in Fig. 13a, the watershed boundary extracted from AW3D30 is matched well with that extracted from SRTM-1. By contrast, the watershed boundary from HydroSHEDS has some minor differences. According to the HydroSHEDS data, the subwatershed located in the east of the enlarged area was judged as a part of its left watershed. However, according to the extracted drainage networks, the surface water of this sub-region should drain into the right basin. Hence, the allocation of this sub-watershed is wrong in the HydroSHEDS data. The third improvement is likely due to the higher data resolution for delineating watersheds based on 30-m SRTM-1 and AW3D30 DEM relative to the approximately 410 m of HydroSHEDS data. The influence of resolution is much more evident in the high relief area than in other areas. As shown in Fig. 13c, the extracted boundaries fit well with the ridge line. However, the boundary from HydroSHEDS is coarser, thereby causing an obvious offset from the drainage divide. After extracting the watershed boundaries for each terminal lakes, the extracted Nam Co basin, which has a relatively large catchment area, is selected to represent the second-level delineation by using

Fig. 11. Comparisons between the lakes type classification results and the reference dataset. (a) Lake No. 42 is judged as the upstream lake in reference dataset, however the elevation profile proves that Lake No. 42 and 44 are unconnected. (b) Lake No. 96 is regarded as terminal lake in reference dataset while it is believed to be upstream lake by the proposed method.

10

K. Liu et al. / Geomorphology 358 (2020) 107127

Fig. 12. Watershed delineation results considering terminal lakes using (a) SRTM-1 DEM and (b) AW3D30 DEM. (c) The watershed boundaries from HydroSHEDS data.

different threshold values. As shown in Fig. 14, this step is the same as the traditional approach of watershed delineation in the exterior drainage. The extracted sub-basin represents the hierarchical structure that can be applied at different spatial scales. 5. Discussions 5.1. Influences of drainage and height threshold values Drainage and height threshold values involved in the proposed method, which may influence its applicability, should be further discussed. It has been discussed that the adapted drainage threshold value should ensure that all lakes can be connected by the extracted streams. Otherwise, whether the unconnected lake is the terminal lake cannot be determined. In this study, a threshold value of 3 km2 is adopted because the minimum area of the selected lakes is 1 km2,

which is relatively small. In this case, all the lakes are connected by the drainage networks. The drainage threshold value of the catchment area also influences the final watershed boundaries directly. The initial watershed boundaries may fail to delineate the drainage divides if the drainage threshold value is relatively high, especially in low relief regions. This process is similar to segmentation and classification in object-based image analysis. With regard to parameter selection, oversegmentation on a small scale is better than under-segmentation because later merging is possible (Martha et al., 2010). For the same reason, a large threshold value in initial watershed delineation may cause the under-segmentation, making the delineated watershed boundaries in the low relief regions have obvious offsets with the real drainage divides. The height threshold value is another key parameter in the proposed method, which decides whether the intra-lake connection is true or not. The adopted standard deviated value can be regarded as a good choice

Fig. 13. Comparisons between the extracted watershed boundaries and the HydroSHEDS data. (a) Overview of the differences in the selected region. (b) The sub-watershed was misclassified as a part of the selected basin. (b) The watershed boundaries from HydroSHEDS data own obvious offsets due to the low spatial resolution.

K. Liu et al. / Geomorphology 358 (2020) 107127

11

Fig. 14. The second-level watershed delineation in Nam Co using different threshold values: (a) Catchment area is 2000 km2; (b) Catchment area is 500 km2; and (c) Catchment area is 100 km2.

which depends on the accuracy of the adopted DEM instead of an empirical value. However, the height threshold value based on SRTM-1 or AW3D30 DEM is still a litter higher in processing the upcoming overflow lake, when the elevation difference between the lake surface and the overflow point is less than 8 m. In this case, the DEM with higher accuracy is needed, which can provide more reliable topographic information at a local scale. Currently, the TanDEM-X DEM with 12 m resolution can be a choice regardless of the budget.

overflow will likely happen when the lake level exceeds the drainage divide. Although watershed boundary in the study area did not change from 2000 to 2015 in this study area, the lake expansion trend has been reflected from the comparison results of closure degree in two different years. As the possibility of the watershed evolution with lake changes change in the inner TP is great, its potential impacts on lake hydrological process, environmental change, and watershed geomorphology deserve special attention.

5.2. Possibility of the watershed evolution with lake changes

6. Conclusions

Watershed boundaries are generally assumed to be constant because geological process and geomorphologic evolution occur over long time scales. However, this inertial thought should be adjusted when a study focuses on an endorheic region that is distributed with lakes with significant dynamics. As the lake level rises, lake surface water may approach the drainage divide. Once the water crosses the watershed boundary, the closure of the watershed will be broken. Recently, accelerated lake expansion in the inner TP has been confirmed due to the climate change (Song et al., 2014; Lutz et al., 2014). To assess the possibility of overflow for each watershed, the degree of closure is measured by calculating the elevation difference between lake level and its lowest overflow point along the drainage divide. As for the watershed with more than one lake, the minimum value is adopted. As shown in Fig. 15, extracted watershed are classified into five types on the basis of the elevation differences. Four watersheds as labeled as extremely low closure with values less than 10 m in 2015. In this case,

In this study, we propose a novel method of automatic watershed extraction for endorheic basins, which is different from the traditional watershed delineation workflow. In accordance with the geomorphic features of the endorheic basin, watershed delineation is conducted through a two-level process, thereby maintaining the independence of all the endorheic basins. Elevation-profile judgment is proposed to determine the drain modes of all the connected lakes and can be effectively used to differentiate terminal and upstream lakes. The tracking analysis is conducted from the source points for watershed labeling. The proposed method achieved evident improvements in lake classification and watershed delineation in the inner TP compared with HydroSHEDS and other existing datasets. Potential users of the proposed method include geomorphologists, hydrologists, and researchers involved cryosphere studies. Experiments in the study area prove that the automatic approach can fully explore the advantages of high-quality DEMs and can be expanded to a large

Fig. 15. Maps of the degree of closure for endorheic basins in (a) 2000 and (b) 2015.

12

K. Liu et al. / Geomorphology 358 (2020) 107127

region. Many high-resolution DEMs will become available in the future, and large-scale studies in High Mountain Asia are attracting considerable attention. Thus, the proposed approach can be an effective tool for handling the scientific and engineering issues.

Acknowledgments This work was partly funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA23100102), the Second Tibetan Plateau Scientific Expedition and Research (STEP) (Grant No. 2019QZKK0202), National Natural Science Foundation of China (Grant No. 41801321, 41971403, 41930102), the National Key Research and Development Program of China (Grant No. 2018YFD0900804, 2018YFD1100101), the Thousand Young Talents Program in China (Grant No. Y7QR011001), and Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences (NIGLAS2016TD01, NIGLAS2017QD15). References Ai, B.B., Qin, C.Z., Ye, Q., Zhu, A.X., Cogley, G., 2017. An approach to extracting surface supply relationships between glaciers and lakes on the Tibetan Plateau. Int. J. Digit. EARTH. 1–15. Barnes, R., Lehman, C., Mulla, D., 2014a. Priority-flood: an optimal depression-filling and watershed-labeling algorithm for digital elevation models. Comput. Geosci. 62, 117–127. Barnes, R., Lehman, C., Mulla, D., 2014b. An efficient assignment of drainage direction over flat surfaces in raster digital elevation models. Comput. Geosci. 62, 128–135. Brun, F., Berthier, E., Wagnon, P., Kääb, A., Treichler, D., 2017. A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000–2016. Nat. Geosci. 10, 668. DeBarry, P.A., 2004. Watersheds: Processes. John Wiley & Sons, Inc., Assessment and Management. Dorsaz, J.M., Gironás, J., Escauriaza, C., Rinaldo, A., 2013. The geomorphometry of endorheic drainage basins: implications for interpreting and modelling their evolution. Earth Surf. Process. Landf. 38, 1881–1896. Evans, I.S., 2012. Geomorphometry and landform mapping: what is a landform? Geomorphology 137, 0–106. Freitas, H.R.A., Freitas, C.D.C., Rosim, S., Oliveira, J.R.F., 2016. Drainage networks and watersheds delineation derived from TIN-based digital elevation models. Comput. Geosci. 92, 21–37. Gao, Y., Wang, W., Yao, T., Lu, N., Lu, A., 2018. Hydrological network and classification of lakes on the Third Pole. J. Hydrol. 560, 582–594. Garbrecht, J., Martz, L.W., 1997. The assignment of drainage direction over flat surfaces in raster digital elevation models. J. Hydrol. 193, 204–213. Grimaldi, S., Nardi, F., Benedetto, F.D., Istanbulluoglu, E., Bras, R.L., 2007. A physicallybased method for removing pits in digital elevation models. Adv. Water Resour. 30, 2151–2158. Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogramm. Eng. Rem. S. 54, 1593–1600. Jiang, L., Tang, G.A., Liu, X.J., Song, X.D., Yang, J.Y., Liu, K., 2013. Parallel contributing area calculation with granularity control on massive grid terrain datasets. Comput. Geosci. 60, 70–80. Khan, A., Richards, K.S., Parker, G.T., McRobie, A., Mukhopadhyay, B., 2014. How large is the Upper Indus Basin? The pitfalls of auto-delineation using DEMs. J. Hydrol. 509, 442–453. Lehner, B., Verdin, K., Jarvis, A., 2013. New global hydrography derived from spaceborne elevation data. EOS Trans. Am. Geophys. Union 89, 93–94. Liu, X., 2008. Airbqorne LiDAR for DEM generation: some critical issues. Prog. Phys. Geogr. 32, 31–49. Liu, J., Zhu, A.X., Liu, Y., Zhu, T.X., Qin, C.Z., 2014. A layered approach to parallel computing for spatially distributed hydrological modeling. Environ. Modell. Softw. 51, 221–227. Liu, K., Ding, H., Tang, G.A., Na, J.M., Huang, X.L., Xue, Z.G., Yang, X., Li, F.Y., 2016. Detection of catchment-scale gully-affected areas using unmanned aerial vehicle (UAV) on the Chinese Loess Plateau. ISPRS Int. J. Geo-Inf. 5, 238. Liu, K., Ding, H., Tang, G.A., Song, C.Q., Jiang, L., Liu, Y.W., Zhao, B.Y., Gao, Y.F., Ma, R.H., 2018. Large-scale mapping of gully-affected areas: an approach integrating Google Earth images and terrain skeleton information. Geomorphology 314, 13–26.

Liu, K., Song, C., Ke, L., Jiang, L., Pan, Y.Y., Ma, R.H., 2019. Global open-access DEM performances in Earth’s most rugged region High Mountain Asia: a multi-level assessment. Geomorphology 338, 16–26. Lutz, A.F., Immerzeel, W.W., Shrestha, A.B., Bierkens, M.F.P., 2014. Consistent increase in High Asia’s runoff due to increasing glacier melt and precipitation. Nat. Clim. Chang. 4, 587–592. Martha, T.R., Kerle, N., Jetten, V., 2010. Characterizing spectral, spatial and morphometric properties of landslides for semi-automatic detection using object-oriented methods. Geomorphology 116, 24–36. Martz, L.W., Garbrecht, J., 1998. The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models. Hydrol. Process. 12, 843–855. Metz, M., Mitasova, H., Harmon, R.S., 2010. Accurate stream extraction from large, radarbased elevation models. Hydrol. Earth Syst. Sci. Discuss. 7, 3213–3235. O’Callaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks from digital elevation data. Comput. Vision, Graph., Image Process. 28, 323–344. Pekel, J.F., Cottam, A., Gorelick, N., Belward, A.S., 2016. High-resolution mapping of global surface water and its long-term changes. Nature 540, 418–422. Phan, V.H., Lindenbergh, R., Menenti, M., 2013. Geometric dependency of Tibetan lakes on glacial runoff. Hydrol. Earth Syst. Sci. 17, 4061–4077. Rabus, B., Eineder, M., Roth, A., Bamler, R., 2003. The shuttle radar topography mission-a new class of digital elevation models acquired by spaceborne radar. ISPRS Int. J. GeoInf. 57, 241–262. Sheng, Y., Song, C., Wang, J., Lyons, E.A., Knox, B.R., Cox, J.S., Gao, F., 2016. Representative lake water extent mapping at continental scales using multi-temporal Landsat-8 imagery. Remote Sens. Environ. 185, 129–141. Song, C., Huang, B., Ke, L., 2013. Modeling and analysis of lake water storage changes on the Tibetan Plateau using multi-mission satellite data. Remote Sens. Environ. 135, 25–35. Song, C., Huang, B., Richards, K., Ke, L.H., Phan, V., 2014. Accelerated lake expansion on the Tibetan Plateau in the 2000s: induced by glacial melting or other processes? Water Resour. Res. 50, 3170–3186. Tadono, T., Takaku, J., Tsutsui, K., Oda, F., Nagai, H., 2015. Status of ALOS World 3D (AW3D) global DSM generation. 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), pp. 3822–3825. Tarboton, D.G., Bras, R.L., Rodriguez-Iturbe, I., 1991. On the extraction of channel networks from digital elevation data. Hydrol. Process. 5, 81–100. Tarolli, P., 2014. High-resolution topography for understanding earth surface processes: opportunities and challenges. Geomorphology 216, 295–312. de Vente, J., Poesen, J., Verstraeten, G., Govers, G., Vanmaercke, M., Van Rompaey, A., Arabkhedri, M., Boix-Fayos, C., 2013. Predicting soil erosion and sediment yield at regional scales: where do we stand? Earth Sci. Rev. 127, 16–29. Wang, L., Liu, H., 2006. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. Int. J. Geogr. Inf. Sci. 20, 193–213. Wilson, J.P., 2012. Digital terrain modeling. Geomorphology 137, 107–121. Xiong, L.Y., Zhu, A., Zhang, L., Tang, G.A., 2018. Drainage basin object-based method for regional-scale landform classification: a case study of loess area in China. Phys. Geogr. 1–19. Xiong, L.Y., Jiang, R.Q., Lu, Q.H., Yang, B.S., Li, F.Y., Tang, G.A., 2019. Improved Priority-Flood method for depression filling by redundant calculation optimization in local microrelief areas. Trans. GIS 23, 259–274. Yamazaki, D., Ikeshima, D., Tawatari, R., Yamaguchi, T., O’Loughlin, F., Neal, J.C., Sampson, C.C., Kanae, S., Bates, P.D., 2017. A high accuracy map of global terrain elevations. Geophys. Res. Lett. 322, 109–127. Yang, K., Ye, B., Zhou, D., Wu, B.Y., Foken, T., Qin, J., Zhou, Z.Y., 2011. Response of hydrological cycle to recent climate changes in the Tibetan Plateau. Clim. Chang. 109, 517–534. Yao, F., Wang, J., Yang, K., Wang, C., Walter, B.A., Crétaux, J.F., 2018. Lake storage variation on the endorheic Tibetan Plateau and its attribution to climate change since the new millennium. Environ. Res. Lett. 13, 064011. Yue, L., Shen, H., Zhang, L., Zhang, X.W., Zhang, F., Yuan, Q.Q., 2017. High-quality seamless DEM generation blending SRTM-1, ASTER GDEM v2 and ICESat/GLAS observations. ISPRS. J. Photogramm. 123, 20–34. Zhang, H., Huang, G., 2009. Building channel networks for flat regions in digital elevation models. Hydrol. Process. 23, 2879–2887. Zhang, S., Foerster, S., Medeiros, P., Araújo, J.C., Motagh, M., Waske, B., 2016. Bathymetric survey of water reservoirs in north-eastern Brazil based on TanDEM-X satellite data. Sci. Total Environ. 571, 575–593. Zhang, G., Yao, T., Shum, C.K., Yi, S., Yang, K., Xie, H.J., Feng, W., Bolch, T., Wang, L., Behrangi, A., Zhang, H.B., Wang, W.C., Xiang, Y., Yu, J.Y., 2017a. Lake volume and groundwater storage variations in Tibetan Plateau's endorheic basin. Geophys. Res. Lett. 5550–5560. Zhang, H., Yao, Z., Yang, Q., Li, S.Q., Baartman, J.E.M., Gai, L.T., Yao, M.T., Yang, X.M., Ritsema, C. J., Geissen, V., 2017. An integrated algorithm to evaluate flow direction and flow accumulation in flat regions of hydrologically corrected DEMs. Catena. 151, 174–181.