The limits of watershed delineation: implications of different DEMs, DEM resolutions, and area threshold values

A watershed is a natural geo-hydrological unit where water flows downhill, moves through a common outlet by a system of streams and distributes to rivers, reservoirs, or lakes (Paranjape et al. 1998; Sowmya et al. 2020). It consists of a natural water divider separating drainage basins and offering a natural boundary for planning, management, and development (LaMoreaux et al. 2008; Kumar & Dhiman 2014). Watershed delineation is a prerequisite in hydrologic and environmental evaluation of a watershed, such as runoff estimation, water quality modelling, flood assessment, and disaster risk assessment (Luo et al. 2011; Giridhar et al. 2015; Ray 2018). Before the advancement of remote sensing and Geographical Information System (GIS) technologies, watershed delineation was predominantly manual drawings of the watersheds and stream networks using visual interpretations and digitization of the topographic and contour maps (Ehsani et al. 2010; Salih & Hamid 2017). With recently advanced GIS technology, watershed delineations are done using remotely sensed terrain features processed by different digital elevation models (DEMs) at various spatial resolutions (Munoth & Goyal 2019).

Most hydrological models use DEM to account for geomorphological and hydrological variables in a watershed spatio-temporal modelling. The spatial resolution of a DEM can affect the outputs for elevation, slope, stream network, sub-watersheds, watershed area, sediment and nutrient loads (Chaubey et al. 2005; Zhang et al. 2014; Reddy & Reddy 2015; Buakhao & Kangrang 2016; Xu et al. 2016; Wu et al. 2017; Ray 2018; Munoth & Goyal 2019). A DEM with a coarser resolution can simulate larger basin areas with shorter streams, flatter slopes and insignificant differences in altitude, which reduces heterogeneity and accuracy (Chaplot 2005; Chaubey et al. 2005; Liu et al. 2010). Conversely, a finer resolution DEM can delineate a more accurate river network but less accurate watershed boundary (Liu et al. 2010; Ray 2018). Moreover, DEMs vary in horizontal and vertical outputs, primarily affecting stream network and delineation in relatively flatter topography (Li & Wong 2010; Luo et al. 2011; Al-Khafaji & Al-Sweiti 2017; Shafiq et al. 2020). The outcomes of GIS-based hydrological models are greatly influenced by the spatial resolution and accuracy of DEMs (Piwowar & LeDrew 1990; Moore et al. 1991; Wolock & Price 1994; Wolock & McCabe 1995; Cho & Lee 2001). Hence, a DEM selection is critical in the watershed delineation process and subsequent catchment solute transport study. The use of two or more DEMs was suggested for efficient hydrogeological research and watershed management of an area (Ray 2018), though the use of a higher number of DEMs has not been reported in the literature until now.

Apart from DEM resolution, Area Threshold Values (ATVs) also control the quality of hydrological model outputs (Munoth & Goyal 2019). A decrease in ATV detects minute differences within the extracted drainage parameters by different DEMs that produce a proportionately large number of drainage channels and, in turn, affects hydrological model outputs such as sub-basin count and surface runoff (Reddy et al. 2018; Gautam et al. 2019; Munoth & Goyal 2019). However, none of these researches suggested any optimum ATV for a specific range of watersheds in any hydrological modelling study. A hydrological modelling software, such as the Soil and Water Assessment Tool (SWAT), derived drainage channel delineation affected by DEM resolution and ATVs (Kalcic et al. 2015).

SWAT is a widely-used physically-based basin-scale distributed continuous hydrologic model (Arnold et al. 1998). In SWAT, the stream burning process modifies DEM elevation data to overlap the surface drainage patterns with existing stream network locations (Wang et al. 2011; Lindsay 2016). The SWAT-generated stream network depends on ATVs and allows users to divide the basin into sub-basins based on topography to incorporate spatial details (Munoth & Goyal 2019). Since a lower ATV gives a higher sub-basin count and subsequent surface runoff reduction, an ATV significantly influences basin hydrological characterization and modelling.

In Bangladesh, hydrological processes, which influence aquatic habitat management and flood control, have been modelled using SWAT-based hydrological models with geomorphologic inputs from open-source DEMs (Raihan et al. 2020). Such modelling helps adapt water resource management and monitoring systems for a river like Halda, a significant aquatic habitat in Bangladesh (Raihan et al. 2020). SRTM (Shuttle Radar Topography Mission) DEM of 30 m resolution was used in the SWAT model to study the effects of hydrological structures on the flow of the Halda and the sustainability of the aquatic habitat (Saha et al. 2019a) and susceptibility to flash floods to the Karnaphuli and the Sangu River basins (Adnan et al. 2019). Evaluation of the hydrological modelling using SRTM 90 m DEM in the flat terrain of twelve catchments in Bangladesh showed that the slope parameter affected river network delineation (Rahman et al. 2010). However, only SRTM DEM and one ATV were used for specific watershed areas in these studies, which apparently could not represent the variation in outputs using different sources and qualities of DEMs in different terrains.

Halda River catchment: Halda originated from the Badnatali Hills range in Ramgarh, a sub-district of Khagrachari District, and ended in Karnaphuli River as one of the major tributaries, about 35 km away from the Bay of Bengal (Saha et al. 2019b). The watershed of Halda is bounded by 91 °36.6′E–92 °1.2′E longitude and 22 °56.4′N–22 °20.4′N latitude (Figure 1). The river basin lies inside the political boundary of Bangladesh and is well known as the only natural carp-breeding habitat in Southeast Asia (Azadi 2005; Correspondent 2014; Saha et al. 2019a; 2019b). The Bangladesh government decided to declare the river as the Ecologically Critical Area to protect the sanctuary of the carp fishes (Hussain 2016). It is declared the ‘Bangabondu Fisheries Heritage’ site (Report 2020) that prohibits all water abstraction, fishing and sand quarrying activities on this river.

Sangu river catchment: Sangu River originates in the North Arakan Hills of Myanmar and enters Bangladesh from the east of Remarki, Thanchi, a sub-district of Bandarban district. The watershed of Sangu is bounded between 91 °48′E–92 °42′E and 22 °24.6′N–21 °14.4′N (Figure 1). Sangu supports irrigation, fisheries, and navigation in the region.

Chengi river catchment: Chengi spreads through 3 sub-districts of Khagrachari district and 2 sub-districts of Rangamati district, starting from Gandacherra, India. The watershed of Sangu is bounded between 91 °43.2′E–92 °12′E and 23 °42.4′N–22 °39′N (Figure 1). In hilly areas of the Khagrachari district, Chengi is the longest and one of the most potent water sources for ethnic local communities. People use the river water for domestic, agricultural, navigation, extraction, and transport of forest resources and recreational purposes. The water quality of this river has been deteriorating lately due to the disposal of untreated industrial effluents like rubber factories, municipal sewage, rubber dam, riverbank erosion, and overuse of fertilizer (Latifa et al. 2019).

In this study, AW3D30 (Advanced Land Observation Satellite World 3D 30 m), ASTGTM (Advanced Space-borne Thermal Emission and Reflection Radiometer Global Digital Elevation Model, Version 3), NASA (National Aeronautics and Space Administration, USA), and SRTM (Shuttle Radar Topographic Mission) DEMs of ∼30 m horizontal resolution; HydroSHEDS (Hydrological data and maps based on Shuttle Elevation Derivatives at multiple Scales), MERIT (Multi-Error-Removed Improved-Terrain), SRTM DEMs of ∼90 m resolution, and GMTED (Global Multi-resolution Terrain Elevation Data) DEM of ∼225 m resolution were used to carry out the analysis (Table 1 and SI – Fig. S1).

The observed stream network may differ from the extracted ones for a DEM due to different resolutions and data sources or modification of a river course by river training (Reddy & Reddy 2015). To define the accurate position of the main streams, a digitized and georeferenced main channel (polyline vector file) from Google Earth was used for each of the Halda, Sangu, and Chengi river basins during the stream network extraction process. We conducted a field survey and communicated with water management authorities to affirm the accuracy of the stream burn-in method and outlet points. Sub-basin creation followed the 8 DEMs for the threshold of 5, 10, 20, 40, 60, 80, 100, 120, and 160 km² areas, giving 72 experiments for the stream networks.

Variation of slope among 90 m DEMs is relatively lower for Sangu and Chengi rivers than 30 m DEMs for all the rivers (Table 2). Mean and standard deviation values are highest for AW3D30, and GMTED has the lowest values for all the river watershed parameters (Table 2). A maximum slope of 69.79 degrees is estimated for Sangu from NASA DEM (Table 2). The minimum slope derived from different DEMs is zero degrees. In contrast, slopes’ maximum, mean, and standard deviation are higher for 30 m resolution and lower for coarser resolutions (Table 2).

The hypsometric curves of the watersheds were concave-up (Figure 3). About 17% of the Halda watershed was steep sloping with more than the elevation of 83 m, and the remaining areas were relatively flat (Figure 3). Though Chengi had similar attributes, it had a steeper slope (about 32% of the area) and different curve positions for DEMs in high-altitude areas. Unlike the two watersheds, Sangu had a relatively smooth sloping pattern at the gradual increment of the area, and around 50% of the area had more than 105 m (average of all DEMs) elevation (Figure 3). Halda watershed had lower HI values (ranged between 0.07 and 0.11 with an average of 0.09) than Sangu (ranged between 0.17 and 0.20 with an average of 0.18) and Chengi (ranged between 0.13 and 0.21 with an average of 0.17) for all the DEMs.

Correlation values among the elevation data are between 0.99 and 1.00 for Halda and Sangu River and 0.98 and 1.00 for Chengi river (SI – Table S2). GMTED correspond least to other DEMs data for any river. NASA and SRTM-30; HydroSHEDS, and SRTM-90 DEMs are mostly correlated for all the rivers. All the DEMs show higher similarity to the elevation data among the studied rivers for the Sangu River. A higher correlation is observed in 90 m DEMs for the Sangu River watershed. The correlation among 30 m DEM is higher in every river watershed. MERIT DEM (90 m) correlates with all 30 m DEMs than coarser DEMs in the Chengi watershed. Additionally, all DEMs show a robust correlation (≥98%) for the elevation.

Among the different DEM slope data, correlation value varies from 0.62 to 0.95, 0.68–1.00, and 0.41 to 1.00 for Halda, Sangu, and Chengi river, respectively (SI – Table S3). HydroSHEDS and SRTM-90 are the most correlated DEMs, and GMTED has the lowest correlation value with all other DEMs for all the rivers, especially the lowest with ASTGTM in Halda and AW3D30 in Sangu and Chengi river watersheds. ASTGTM has lower correlation values with all coarser resolution (90 m and 225 m) DEMs in the Halda watershed, but AW3D30 has the same for other river watersheds. NASA and SRTM-30 have the highest correlation among all other 30 m resolution DEMs. The correlation among 90 m DEMs is higher than among finer DEMs (30 m).

The size and shape of the Halda, Sangu, and Chengi river watersheds are different as computed by the DEMs. Halda, Sangu, and Chengi rivers intersected common areas have 6.79, 12.38, and 0.03% differences with the overlapped total area, respectively, except AW3D30 DEM (53.75%) for the Chengi watershed (Figure 4). The delineated shapes of Halda, Sangu, and Chengi river watersheds for different DEMs are shown in Figure 4 (and SI – Figure S1).

Northeast and south of the Halda watershed showed minor differences at every bend of the river, which may be caused by loss of information at lower elevations by different DEMs. The variation of the boundary in the east and west parts was inferior. The shape differences for the Sangu River watershed were seen on the northwest side, mainly the coastal area with low elevation. Though a low deviation is observed in the southeast, delineated watersheds are typical for all the DEMs, except AW3D30 DEM in the Chengi river watershed. In agreement with the mountainous Gilgit watershed in Pakistan, SRTM, ASTER, and GTOPO30 DEM showed close agreement in delineating the watershed, but SRTM (30 m) was more accurate than others in ridge demarcation (Shafiq et al. 2020).

The intersected areas between the defined watersheds from various DEMs are extracted. The intersected areas by different DEMs are 97.62–99.73%, 93.08–99.92%, and 67.80–99.92% for Halda, Sangu, and Chengi rivers, respectively (Figure 4 and SI – Table S4). The difference in the intersected areas is small for Halda and larger for Chengi. The ASTGTM is comparatively less intersected by other DEMs (97.62%) among the delineated watersheds for Halda. But the highest value (99.73%) for SRTM-30 when AW3D30 is intersected by it. For Sangu, the lowest intersection area (93.08%) is found between AW3D30 and NASA. It is evident from Figure 4 that for the Chengi, AW3D30 intersects the least with other DEMs and has the lowest value (67.80%) for intersecting GMTED.

SWAT follows the power function as y = ax^b, in sub-basin and drainage network extraction based on an ATV (Wu et al. 2017). The ‘a’ values range from 1,045.60 to 1,517.60, 2,252.50 to 2,699.50, and 2,344.90 to 3,731.20 for Halda, Sangu, and Chengi river watersheds, respectively, where R² is >0.9473 (Figure 5(a) and Table 3). Variations in ‘b’ value show that Halda has the lowest range (−0.97∼−1.07) values following the Sangu (−1.03∼−1.09) and Chengi has the highest values (−1.36∼−1.69). These curve fitting values are varied for different DEMs having different resolutions, but these variations are basin-specific.

Around 85.86, 77.19, and 76.83% of diminution happened in stream length derivation for all DEMs at ATV 5–160 km² for Halda, Sangu, and Chengi river watersheds, respectively (Figure 5(b)). The extracted stream had increased networking, and obscure streams were more visible when the imposed ATVs decreased, irrespective of DEM sources and resolutions (Figure 5(c)). Usually, the extracted total stream length from the finer-resolution (30 m) DEMs (except AW3D30 in Halda and Chengi) are longer than from the coarser resolutions (SI – Table S5). GMTED shows the lowest total stream length among all other DEMs for different ATVs, and AW3D30 shows the most extensive stream length (834.61 km) for Halda. The minimum and maximum sum of stream lengths for the Sangu river watershed were 291.61 km (GMTED), 1,519.85 km (AW3D30) at ATVs 160 km² and 5 km², respectively. As AW3D30 delineated the smallest watershed (1,253.98 km²) for Chengi, it had the smallest stream network and the lowest total stream length among all the DEMs (Figures 4 and 5(c) and SI – Table S5).

Considering hypsometric curves, 9, 40, and 17% of landmass elevations are above mean elevations (40.02 m, 173.76 m, and 102.34 m) for the Halda, Sangu, and Chengi rivers watersheds, respectively. These basins are at the monadnock stage, concave up with low integrals (HI < 0.35), usually formed from isolated bodies of resistant rock of major hills (monadnock) that are found above the subdued surface and old and deeply dissected landscapes (Strahler 1952; Kusre 2013). Three river basins have varied estuary connections and topography. All the lowest elevation values of the watershed are near the stream network.

ASTGTM shows an unlikely lower bound value (0 m, identical to MSL), and AW3D30 shows a meagre elevation value (−30 m or −18 m) on stream surface area for the Halda and Sangu watersheds which does not physically exist (observed during field investigation). Halda connects with Karnaphuli about >35 km before the estuary of the Bay of Bengal. Here, MERIT (90 m) shows a comparable elevation value to other DEMs elevation ranges. These elevation discrepancies are derived from the different acquisition and update time and calculation errors of DEMs (Stefanescu et al. 2012; Xu et al. 2016), resulting in high-value depressions, which do not exist. However, ≥98% correlation values among the DEMs for elevation data confirm that DEM sources and resolutions have minimal effect on elevation data for specific watershed types and topographical variation.

We observed that the slope data vary with different DEMs and resolutions which is in agreement with Buakhao & Kangrang (2016). However, coarse-resolution DEMs (90 m) show lower slope parameter variation (e.g., maximum). Hence, coarser-resolution DEMs ignore topographic information details (Liu et al. 2010; Chen 2013; Xu et al. 2016). In topographic parameter extraction, DEM resolutions have more impact than DEM sources (Xu et al. 2016), indicating that elevation and slope data are more sensitive to the resolution of DEMs (Dixon & Earls 2009; Tan et al. 2015). Though the variation in elevation is not significant, these complex DEM data are processed (filled, conditioned) to extract water flow direction and channel properties such as channel length, width, depth, slope, and watershed separation (Liu et al. 2010; Rao & Yang 2010). For these reasons, the effect of DEM resolution is minimal in streamflow simulation in SWAT (Zhang et al. 2014; Gautam et al. 2019; Munoth & Goyal 2019). On the contrary, TOPMODEL, Hydrological Simulation Program-FORTRAN (HSPF) streamflow showed sensitivity against DEM resolution compared to DEM sources and different resampling techniques (Suliman et al. 2019; Roostaee & Deng 2020) or in SWAT (Rocha et al. 2020).

In the study, Halda, Sangu, and Chengi river watersheds differed in size and shape for all the DEMs. We estimated larger areas for Halda and Chengi and lower perimeter values for all three rivers from 90 m DEMs, indicating that coarser-resolution DEMs represent lower information than finer resolutions DEMs. This is contrary to the findings of Al-Khafaji & Al-Sweiti (2017). They inferred inconspicuous trends in the watershed area with DEM resolutions. A larger delineated area (3,842.82 km²) for the Sangu river watershed is reported using SRTM-30 DEM (Adnan et al. 2019), where the outlets of the streams vary in our study. Akhter (2015) found a 1,727.06 km² area by SRTM-90 DEM using HEC-GeoHMS, and Chowdhury et al. (2020) found a 1,756.91 km² area by ASTER DEM for the Halda river. Area differences in various studies could be due to the different outlet/discharge point selections for the main channel and individual streams and the choice of DEMs (source and resolution) or software packages used in watershed delineation. Considering minor deviation, the SRTM 90 m DEM delineated a higher watershed area than other DEMs and resolutions. The boundary delineated by DEM may overlook an area delineated by another (Buakhao & Kangrang 2016). Though primarily the difference is derived from the local sink-related problem for all hydrological tools, SWAT produces a minimum error in watershed delineation than others (Ray 2018). Using high-resolution DEM in the automated watershed delineation process provides no significant benefit (Buakhao & Kangrang 2016), though it extracts more accurate streams (Li & Wong 2010). However, coarser-resolution DEMs perform more precisely than higher resolution DEMs in the delineation process (Ray 2018) due to fewer topographic variations at boundaries (Xu et al. 2016).

One of the reasons for the difference between the delineated overlapped total area and intersected common area by the DEMs could be the difference in vertical datum projection like ITRF97 - GRS80 in AW3D30 than other DEMs – WGS84 or converting them to UTM 46 N. Moreover, the MSL elevations of the reflective surfaces of any features are elevated above the bare earth by AW3D30 DSM. However, it cannot be confirmed for the Halda and Sangu River watersheds. There has been no demarcated watershed boundary for the management of these rivers. Hence, for initiating such approaches, the common intersected area, independent of the DEMs and their resolutions, can be used as a starting point to define basin boundary by keeping the total overlapped area as a buffer or boundary zones for the watersheds.

The correlation values among delineated watersheds by different DEMs and resolutions indicate a close association among DEMs. ASTGTM and NASA, AW3D30 showed lower intersected areas than all the studied DEMs (SI – Table S4), reflecting the correlation values of slope (SI – Table S3). DEMs accurately delineate watersheds surrounded by mountains (Shafiq et al. 2020). Hence, along the flood plain, coastal areas, or flat topography near the discharge points of the catchment, Figure 4 reveals a significant variation in the size and shape of the catchment. Hence DEM selection in flat areas may have more significant implications on discharge and pollution load estimation in any hydrological study.

DEM data with different resolutions have little influence on the extraction results if ATV lies between 5 and 100 km² (Figure 5(a)). However, at ATV ≤40 km², the DEM resolution greatly influences sub-basin counts. The area and perimeter of a catchment are independent of the ATVs. Wu et al. (2017) also reported that the values of a and b of the ATVs to subbasin count equation, y = ax^b, were varied for different watersheds but not for the resolutions of a specific or single DEM. The inverse relationship (as exponent b is negative) between ATVs and sub-basin counts (Munoth & Goyal 2020) become more insignificant for the watersheds for different DEM sources and resolutions at x < 60 km² for Chengi, x < 80 km² for Halda, and x < 100 km² for Sangu River. An efficient and rational sub-basin number is needed in hydrological modelling (Chang 2009). This is achieved here at ATV < 40 km² for the watershed area 1,200 km²–4,000 km². DEM resolution influences sub-basin count and stream length (Figure 5(b) and 5(c)). Additionally, the processing time by the software package also increases with higher resolution DEMs, which agrees with Munoth and Goyal's observation (2019). However, Reddy & Reddy (2015) reported that the accuracy of the sub-basin decreases with the lower resolution of DEM. Hence sub-basin count operation affects the processing time of DEM in delineation, and the input ATVs command it.

Figure 5(c) reveals that to detect a detail stream network, ≤40 km² ATV is needed. Though lowest ATVs cannot be suggested for an acceptable stream network extraction. The complex network of small streams, reservoirs, lakes, and depressions of the flat area would cause difficulty in the identification of braided streams and the natural routing process controlled by water structures or other anthropogenic changes in stream networks (Wilson & Gallant 2000; Turcotte et al. 2001; Luo et al. 2011). However, this threshold would be feckless for a larger watershed (>4,000 km²) because large stream density may increase the computation cost and complexity of the model. Luo et al. (2011) suggested that manual editing after automated delineation of sub-basin boundary and stream layers to adjust the location, range, and hydrologic connection could bring better results in the hydrologic model setup. The field survey in the catchment is essential while doing the manual edit for the studied river basins.

Depending on the ATV, the number of streams and stream lengths will vary in a watershed. A smaller ATV will increase stream length and a denser stream network. Hence, the drainage density, which is defined as total stream lengths (km)/basin area (km²), directly depends on total stream length and follows the same trend of sub-basin for a specific watershed. Cartosat-1 and ASTER DEM showed 0.05 km² as the optimum threshold value to obtain a refined and higher drainage density for the Neri watershed, Maharashtra, India, where the watershed area (about 31.96 km²) was small (Reddy et al. 2017). Therefore, threshold value selection is crucial and directly influences the basin characteristics and stream network (Gopinath et al. 2014) as well as the hydrologic model performance (Gautam et al. 2019).

This study offers a crucial answer to the question of available, open-ended, and space-borne DEMs, their resolutions, and user-defined ATVs’ influence on the physiographic and morphometric characteristics of different river watersheds. DEMs’ sources and resolutions influence different outputs of automatic watershed delineation because of the watershed features such as size and shape differences, stream network, watershed area, and sub-watershed size and shape. The studied watersheds also observed the influence of different ATVs in stream networks and sub-basin for a specific river watershed. This study suggests depending on the watershed size, slope, and topography, DEM-based watershed delineation process in SWAT needs optimum ATV values to capture a detailed and accurate subbasin area, shape and stream network. The detection and delineation of watershed areas in the mountainous part are relatively precise, but the flatter surface creates discrepancies, which suggest the use of the common intersected area of DEMs to confine a river's watershed management area. Results of this study indicate that the choice of input data (here DEMs or topographic features) can be a factor in defining terrain differences and optimum stream complexation covering expected watershed and sub-basin areas. We suggest the comparative study of outputs and operational differences of similar experiments using other geospatial analysis software packages. Moreover, river watersheds similar or different in size can affect the variation or preciseness of the impacts of DEMs, resolution, and ATVs. A base watershed is necessary to suggest a particular DEM source or resolution in watershed delineation. Accuracy assessment of a particular DEM source or resolution for stream network extraction by comparing the base watershed or observed ones is recommended for future study.

The first and second authors highly appreciate the financial support of the Planning and Development Division of Chittagong University for the research grant on the Halda watershed. The authors highly appreciate the reviewers’ suggestions and comments, which significantly improved the manuscript.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.