Appraisal of Tectonic‐Geomorphic Features in the Hindu Kush‐ Himalayas

Abstract

The appraisal of tectonic‐geomorphic features is the basis for the development and management of land, the selection of road routes, and the site selection and construction process of hydropower projects. However, properly evaluating tectonic‐geomorphic features and revealing the relationships among geomorphic evolution, landslides, and tectonic activity remain major challenges in geography and geomorphology. We take the core area of the planned Diamer‐Bhasha Dam as a study area. On the strength of digital elevation models (DEMs), geomorphic indices are extracted to evaluate the tectonic activity and geomorphic evolution. Landslide cataloging by field investigation is used to reveal the relationships between geomorphic evolution, landslides, and tectonic activity. We found that the tectonic activity and geomorphic evolution of the sub‐basin where the Diamer‐Bhasha Dam is planned and the sub‐basin along the Indus are in the stable and old stage, respectively. In contrast, the tectonic‐geomorphic features in the marginal sub‐basins far from the Indus are still active. The correlations among lithology, AF, Vf, HI, SL, and landslide indicate that tectonic activity can influence geomorphic evolution and induce landslides, whereas changes in lithology do not. In addition, landslides can exacerbate geomorphic evolution.

Introduction

Tectonic activity has occurred frequently in the Hindu Kush-Himalayas since the plates of Karakoram, Kohistan, and India collided with each other during the Cretaceous-Tertiary (Debon, 1995). Accordingly, landscape fragmentation (landslides, debris flows, glacier surging, etc.) caused by tectonic activity is widespread in this area (Ravikant et al., 2009). According to the needs of development, there are many engineering projects (road, hydropower, and energy engineering projects, etc.) in such mountainous areas to select the site (Sabir et al., 2017). Active tectonic activity and the continually fragmenting landscape have greatly hindered the construction of engineering projects (Hewitt, 1998; Li et al., 2005). Therefore, it is crucial to continuously assess and monitor the tectonic-geomorphic features in the region (Makkaveyev, 1972; Geach et al., 2017).

On the watershed scale, the recognition of tectonic-geomorphic features has an obvious role in predicting the risk of future earthquakes and secondary landslides (Lin et al., 2013). Some scholars believe that the relative active tectonic index (IRAT) is an effective tool for evaluating tectonic-geomorphic features and have combined it with landslide frequency to predict regional landslide risk (Chang et al., 2015; Cheng et al., 2016; Partabian et al., 2016). However, the geomorphic indices used to calculate the IRAT all contain the hypsometric integral (HI). Strahler (1952) proposed a theory in which the parameter HI can be used to represent the state of geomorphologic evolution. Geomorphic evolution is the result of tectonic activity rather than the cause of tectonic activity, as has been confirmed in many works (Goswami et al., 2009; Malik & Mohanty, 2007; Mathew et al., 2016; Micallef et al., 2014; Sarp, 2015). Therefore, the parameter HI, which represents the result of tectonic activity may not be used to calculate IRAT. As a result, when evaluating geomorphic-tectonic features, we should analyze the active tectonic features and the state of geomorphic evolution using other geomorphic indices.

Tectonic activity is the key factor that causes the uplift of rocks, accelerates river cutting, causes river diversion, and promotes geomorphic evolution (Cox, 1994; Jackson et al., 1998; Clark et al., 2004; Salvany, 2004; Schoenbohm et al., 2004). Landslides and changes in lithologic boundaries are both responses to tectonic activity and potential signals of geomorphic evolution (King, 1972; Menges, 1990; Hovius et al., 1998; Van der Beek & Braun, 1999). Recent works have demonstrated that there is a prominent relationship between tectonic uplift and geomorphic processes controlled by landslide development and sediment migration, and tectonic activity plays an important role in the development of catastrophic landslides (Strecker et al., 2003; Sajinkumar & Anbazhagan, 2015; Hu et al., 2019; Yang et al., 2021). Geomorphic evolution is the result of a combination of many factors, and lithologic boundary changes related to tectonic activity have been identified in some cases to be the dominant factor responsible for the development of the present landscape, especially in actively deforming areas (Gallousi & Koukouvelas, 2007; Carlini et al., 2018; Qiu et al., 2021). However, controversy remains regarding the relationships among tectonic activity, lithology, landslides, and geomorphic evolution, and few studies have detected the response of geomorphic evolution, lithology, and landslides to tectonic activity by exploring the relationship among them. Additionally, most studies have obtained this response relationship on the basis of subjective judgment.

Here, we take the core area of the planned Diamer-Bhasha Dam as the study area. Leveraging digital elevation models (DEMs), geomorphic indices are extracted within an ArcGIS environment to evaluate the tectonic activity (the asymmetry factor (AF) and the valley floor width to valley height ratio (Vf)) and geomorphic evolution (the HI, the stream length gradient index (SL), and the longitudinal river profile). The main objectives of this study are as follows: (1) evaluate tectonic-geomorphic features; (2) quantify the relationships among tectonic activity, landslides, and geomorphic evolution; and (3) reveal the responses of geomorphic evolution and landslides to tectonic activity. The results are expected to provide a scientific reference and theoretical basis for the construction of hydropower projects.

Geological and Geomorphic Setting

The planned Diamer-Bhasha Dam is located on the Indus in northern Pakistan, approximately 315 km upstream of the Tarbela Dam, 165 km downstream of the Northern Areas capital of Gilgit, and 40 km downstream of Chilas. The height of the planned dam is 272 m, the reservoir area is 9.25 × 10³ km², and the reservoir capacity is 7.89 × 10⁹ m³ (Amin et al., 2004). In this paper, with the location of the Diamer-Bhasha Dam as the core, the distance from the starting point of the reservoir to the dam as the radius, and the Indus as the main channel, 33 sub-basins are extracted as the research area based on a DEM (Figure 1a).

Figure 1. Overview of the study area. (a) The elevation distribution of the basin and the location of the planned Diamer-Bhasha Dam; (b) the distribution of lithologies, faults, and landslides in the basin.

The study area is located in the Himalayan seismic belt, which is the largest seismic belt in the world. Active structures that cause riverbed uplift, river incision, bedrock erosion, and geomorphic evolution are widespread and include thrust faults, normal faults, and anticlines/synclines (Figure 1b). The geological conditions are complex, and the dominant lithologies include diorite, gneiss, metamorphic rocks, sedimentary rocks, amphibolite, limestone, andesite, loess, ultramafic rock, glacier debris, anorthosite, and per-metamorphic rocks with ages ranging from Cambrian to Quaternary (Figures 1b and 2). The region has a dry and harsh climate with an annual rainfall of less than 250 mm, and the large temperature difference between day and night causes alternating freeze-thaw cycles (Salma et al., 2012). In this area, the tectonic-geomorphic features are extremely unstable, and landslides occur frequently. Through the interpretation of Google Earth images and field investigations, 203 landslides have been mapped in the study area. Eighty percent of these landslides are large-scale landslides with areas greater than 0.5 km², and only 5% are small-scale landslides (area less than 0.2 km²) (Figures 1b and 2).

Figure 2. Field investigation of landslides and lithologic markers. (a), (b), (d), and (e) Photographs taken on-site. (c) An image from Google Earth.

Materials and Methodology Data Sources

To research the tectonic-geomorphic features and the relationships among tectonic activity, landslides, and geomorphic evolution, a geological map of northern Pakistan at a scale of 1:650,000 is used in this paper to determine the locations of lithological boundaries and faults. The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) with a resolution of 30 m is used to calculate geomorphic indices and extract water systems. Because of abnormal values in the original DEM, we need to eliminate these values before calculating the geomorphic indices and extracting the water systems. First, a TIN is generated in ArcGIS. Then, ArcScene is used to open the TIN (with contour lines and elevation points) and the original DEM. Abnormal values in the original DEM are found and eliminated. For the extracted water system, we use a high-precision (5 m) Google Earth image for spatial correction. In addition, Google Earth images with a resolution of 5 m are used for landslide interpretation. The geological map of northern Pakistan was obtained from the Center of Excellence in Geology, Peshawar University. Google Earth images with 5 m resolution were obtained from the Enterprise Edition 91 Weitu Assistant. DEM data were obtained from Geospatial Data Cloud (http://www.gscloud.cn).

Methods

Geomorphic indices are indicators capable to detect landform responses to recent deformation processes and therefore have been broadly used to characterize sectors deformed by active faults (El Hamdouni et al., 2008; Pedrera et al., 2009). Geomorphic indices can be used to quantify geomorphologic responses to tectonic activity. The objectives of this paper are to analyze the tectonic-geomorphic features of the study area and reveal the relationships among tectonic activity, landslides, and geomorphic evolution in this area. Therefore, we select geomorphic indices under the two main themes of tectonic activity and geomorphic evolution. We select two parameters, AF, representing tectonic horizontal tilt, and Vf, representing tectonic vertical uplift, to quantify the state of tectonic activity in the study area. We select three geomorphic parameters (SL, longitudinal river profiles, and HI) to evaluate the geomorphic evolution of the watershed from three dimensions: Point, Line, and Polygon (Figure 3).

Figure 3. Selection of geomorphic indices for detecting tectonic-geomorphic features. (a) The result of tectonic tilt (modified after Keller et al., 1997); (b) the result of tectonic uplift and fluvial incision. The photographs were taken during the field investigation.

Asymmetry Factor (AF)

The AF is a key index to analyze whether a basin has been tilted by active tectonics, and the value of AF can sensitively reflect tectonic activity in the horizontal direction (Hare & Gardner, 1985; Keller et al., 1997). When AF is significantly greater than or less than 50, the basin is controlled or eroded by active tectonics and is tilted. When the AF approaches 50, the basin is little affected by tectonic activity and exhibits little or no tilting. Hence, we define AF as the absolute value of an area ratio minus 50 and use an arrow to indicate the direction of the basin tilt (El Hamdouni et al., 2008). [Image Omitted. See PDF]where A_r is the area on the right side of the mainstream of the sub-basin (facing downstream) and A_t is the total area of the sub-basin. Both A_r and A_t are calculated by ArcGIS.

Table 1 AF Value and State of Each Sub-Basin

Basin no.	AF	State	Basin no.	AF	State
1	22.23	Active	18	7.82	Stable
2	9.04	Complex	19	1.22	Stable
3	14.50	Complex	20	5.03	Stable
4	1.88	Stable	21	17.41	Complex
5	14.71	Complex	22	25.14	Active
6	0.39	Stable	23	20.51	Active
7	4.40	Stable	24	7.29	Stable
8	13.49	Complex	25	4.82	Stable
9	9.07	Complex	26	12.82	Complex
10	10.29	Complex	27	27.29	Complex
11	1.17	Stable	28	14.78	Complex
12	6.76	Stable	29	24.88	Active
13	7.13	Stable	30	20.15	Active
14	0.22	Stable	31	11.14	Complex
15	1.79	Stable	32	27.39	Active
16	37.99	Active	33	14.01	Complex
17	6.09	Stable

Abbreviation: AF, Asymmetry Factor.

Valley Width to Valley Height Ratio (Vf)

Vf is a geomorphic index that reflects tectonic uplift. The value of Vf can distinguish whether a sub-basin is a V-shaped valley or a U-shaped valley (Bull & McFadden, 1977). Vf is defined as: [Image Omitted. See PDF]where Vf_w is the valley floor width, E_ld and E_rd are the elevations of the left and right sides of the valley watershed (facing downstream), and E_sc is the elevation of the valley floor (Figure 4c). A small value of Vf (Vf < 1) indicates that the basin has been strongly uplifted, resulting in a high incision rate and a V-shaped valley, whereas a larger Vf value (Vf > 1) indicates that the tectonic activity is relatively static and that the shape of the basin is a flat U-shaped valley (Keller et al., 1997).

Figure 4. The calculation principle of the valley width to valley height ratio (Vf). Figures (a) and (b) are the locations of the cross-sections. Figure (c) is a two-dimensional perspective view of the parameters of the valley floor width to valley height ratio.

Since valleys in mountain fronts are often very narrow, we need to set a definite distance from the outlet of the sub-basin to select the cross-section for calculating the Vf (Ramírez-Herrera, 1998; Silva et al., 2003). For a large-scale sub-basin, if only one cross-section is selected to calculate Vf, the result will be insufficient to represent the tectonic uplift of the basin. Thus, we need to select multiple cross-sections in different locations of these sub-basin rivers to calculate Vf and take the average value to assess tectonic uplift. Accordingly, in this paper, we set a definite distance of 2–4 km from the outlet of each river and select 64 cross-sections to calculate Vf values (Figures 4a and 4b).

Longitudinal River Profile

The longitudinal river profile is a crucial index reflecting the influence of tectonic activity on the fluvial incision (Kirby & Whipple, 2001; Whipple, 2004). Longitudinal river profiles are often used to investigate the relationship between erosion and uplift (Hack, 1973; Hovius, 1998; Keller et al., 1997). Hovius (1998) defines three types of profile morphologies for observing geological processes and tectonic activity. Concave profiles reflect a balance between uplift and erosion, S-shaped profiles show the dominance of erosion, and convex profiles reflect the dominance of uplift (active tectonics).

The longitudinal river profile of each sub-basin is extracted from the DEM by Global-Mapper software. Owing to the influence of regional geographical location, basin size, and data resolution, the extracted longitudinal river profiles are chaotic and cannot well reflect the impact of tectonic activity on the fluvial incision. However, the normalization principle of the longitudinal river profiles can resolve this issue (Demoulin, 1998). Therefore, the Y coordinate (elevation) and X coordinate (horizontal distance) used to draw each longitudinal river profile are normalized in this paper. After normalization, the abscissa and ordinate are L/Lo and H/Ho, respectively. H is the river elevation at the measuring point, Ho is the river elevation from the river mouth to the headwaters, L is the distance from the measuring point to the headwaters, and Lo is the distance from the headwaters to the river mouth (Antón et al., 2014).

Stream Length Gradient Index (SL)

AF is one of the vital indices used to measure the change in the longitudinal profile of a river (Hack, 1973). SL values reflect not only the response of channel erosion to tectonic activity but also the influence of tectonic uplift on the lithology (Font et al., 2010; Pérez-Peña et al., 2010). The formula for calculating SL is as follows: [Image Omitted. See PDF]where ΔH is the distance between contours, ΔL is the length of the segment by contour, and L_sc is the total channel length from the midpoint to the drainage watershed (Figure 5).

Figure 5. The method used to calculate the stream length gradient index (SL) (modified after Hack, 1973).

When we calculate SL, we first define the projection of WGS-1984-UTM-Zone-42N for rivers in ArcGIS. Second, the DEM is converted into contours with an interval of 30 m, the intersection function tool of ArcGIS is used to automatically segment the river, and the midpoint of each section is extracted (Figure 5). Then, the property sheet containing ΔH and ΔL is exported to calculate the SL in a Microsoft Excel worksheet (Troiani et al., 2014). Finally, the ordinary kriging interpolation of discrete values is used to generate SL maps (Cressie, 1990) (Figure 10).

Figure 6. Directions of tectonic tilting (the white arrows) and the asymmetry factor values of the sub-basins, which are divided into three classes: strongly asymmetric basins (active), moderately asymmetric basins (complex), and weakly asymmetric basins (stable).

Figure 7. The values of the valley width to valley height ratio (Vf) of the 33 sub-basins in the study area, which are divided into three classes: strong tectonic uplift (active), moderate tectonic uplift (complex), and little or no tectonic uplift (stable).

Figure 8. Longitudinal river profiles of the 33 major drainages in the research area. Note that all drainages are normalized to allow for representation at the basin scale. Boxes A, B, C, and D indicate knickpoints. T-T' is the longitudinal river profile of the Indus. Curves of the sub-basins to the north of the Indus are red, whereas curves of the sub-basins to the south of the Indus are blue.

Figure 9. A conceptual model relating the results of the longitudinal river profile and SL of the Indus. The formation of knickpoints (A, B, C, and D) is shown to be a result of tectonism (faults and earthquakes), river capture, and glacier melting (landslides). The photographs demonstrate typical landforms in these locations.

Figure 10. Stream length gradient index (SL) distribution and the result of SL interpolation.

Hypsometric Integral (HI)

The state of geomorphic evolution is usually quantified by HI values and hypsometric curves (Guarnieri & Pirrotta, 2008; Ma et al., 2021; Pérez-Peña et al., 2010; Willgoose & Hancock, 1998). A hypsometric curve can reflect the distribution of basin area and elevation, and its shape represents the phase of geomorphic evolution. Convex, concave, and S-shaped curves signify basins in the infancy stage, old stage, and mature stage of geomorphic evolution, respectively (Hack, 1973). The HI value represents the degree of geomorphic erosion (Strahler, 1952). A low HI value indicates that the basin is experiencing little erosion and is already in the old stage of geomorphic evolution, while a high HI value indicates that the basin is experiencing severe erosion and is in the infancy stage of geomorphic evolution. The drawing of hypsometric curves is based on Willgoose and Hancock (1998), and the formula for HI is defined as (Pike & Wilson, 1971): [Image Omitted. See PDF]where H_max, H_min, and H_mean are the maximum, minimum, and average elevation of a basin, respectively.

Results Evaluation of Tectonic Activity Asymmetry Factor (AF)

The AF values of the 33 sub-basins are calculated to evaluate the tectonic tilt of the basins. The arrows indicate the tilt direction of the basins. The DEM-based AF values of the 33 sub-basins in the study area range from 0.22 (sub-basin 14) to 37.99 (sub-basin 8) and can be divided into three grades by the natural breakpoint method in ArcGIS: stable (0.22–7.82), complex (7.83–17.41), and active (17.42–37.99) (Figure 6 and Table 1). The lower AF values are all located in the sub-basins 4, 6, 7, 11, 12, 13, 14, 15, 17, 19, 20, and 24 near the Indus, indicating that the tectonic activity along the Indus is largely stable and that the sub-basins exhibit little or no tilting. The AF values in sub-basins 1, 16, 22, 23, 27, 29, 30, and 32 far from Indus are higher, suggesting that these sub-basins are controlled by active tectonics and exhibit severe tilting.

Valley Width to Valley Height Ratio (Vf)

The Vf values in the study area vary from 0.20 (sub-basin 33) to 1.46 (sub-basin 4). The Vf values are divided into three classes: active (0.20–0.43), complex (0.43–1.01), and stable (1.01–1.46) by using the natural breakpoint method (Figure 7 and Table 2). Lower Vf values (0.20–0.43) are mapped in sub-basins 1, 2, 19, 13, 29, 30, 31, 32, 33, 21, 24, and 27, which are V-shaped and have high degrees of tectonic uplift and rates of fluvial incision relative to those of the other basins. Higher Vf values (1.01–1.46) occur along the main channel where the sub-basins (4, 6, 7, 11, 13, 14, 15, 17, and 18) are U-shaped and the tectonic activity is relatively stable.

The spatial distributions of the geomorphic indices (AF and Vf) used to analyze the tectonic activity in the study area show that the tectonic activities of the sub-basin where the Diamer-Bhasha Dam is planned and the sub-basin along the Indus are in a stable state. In the marginal regions far from the Indus, severe tectonic tilting and frequent tectonic uplift indicate extremely active tectonic activity in these regions.

Table 2 Vf Value and State of Each Sub-Basin

Basin no.	Vf	State	Basin no.	Vf	State
1	0.43	Active	22	0.54	Complex
2	0.25	Active	23	0.42	Active
3	1.46	Stable	24	0.31	Active
4	1.37	Stable	25	0.58	Complex
5	0.57	Complex	26	0.63	Complex
7	1.38	Stable	27	0.27	Active
8	0.98	Complex	28	1.32	Stable
9	0.57	Complex	29	0.21	Active
10	0.49	Complex	30	0.39	Active
16	0.69	Complex	31	0.21	Active
19	0.31	Active	32	0.24	Active
20	0.64	Complex	33	0.20	Active
21	0.41	Active

Note. When we calculate Vf values, the sub-basins (No. 6, 7, 11, 12, 13, 14, 15, 17, 18) in which the Indus River is located are combined into one sub-basin (No.7).

Geomorphic Evolution Analysis Longitudinal River Profiles

The longitudinal river profile analysis results are presented in Figure 8. All basin drainage channels are visually comparable to the analysis of tectonic tilt and uplift (Figures 6 and 7), with a dominance of active tectonic-geomorphic features in the marginal regions of the basin and stable conditions recorded along the Indus. This would suggest that the sub-basins are well coupled with their drainage systems and that misfit conditions (e.g., an active channel incising in a relatively stable sub-basin) do not exist. At a basin scale, with all drainage values normalized, according to the shape of the normalized longitudinal river profile, we divided the results into three types: active, stable, and complex (Figure 8). The longitudinal river profiles of sub-basins 1, 2, 22, 23, 30, and 31 are convex, which indicates that these sub-basins are subject to serious runoff erosion and that the river has a high rate of fluvial incision in the vertical and transverse directions. The longitudinal river profiles of sub-basins 6, 7, 8, 11, 12, 13, 14, 15, 17, and 18 are concave. A concave profile suggests an equilibrium between uplift and erosion, and the tectonic-geomorphic features are relatively stable.

Knickpoints, which are steep reaches in the longitudinal profile, can be caused by a resistant lithology, an increase in shear stress, or surface uplift (Bishop et al., 2005). In each sub-catchment, knickpoints mark the boundary between the steady state and adjusting landscape. These anomalies in a profile can indicate either a stream in equilibrium where the upstream retreat communicates changes in base level to an upstream valley (Bishop et al., 2005) or, in some cases, a dynamic equilibrium between fluvial processes and tectonic movements (Snow & Slingerland, 1990). It is worth noting that there are knickpoints (A, B, C, and D) at the confluence of the Indus and the rivers in sub-basins 14, 10, 29, and 5 (Figure 8). We analyze the reason for the occurrence of these knickpoints (Figure 9) based on the longitudinal river profile, SL values, and field investigations of the Indus. Two of the knickpoints (A and B) are caused by melting glaciers and landslides. Melting glaciers can erode the surface and induce landslides, which in turn can cause geomorphic changes and the emergence of knickpoints (A and B). The investigation also ascertained that knickpoint C is caused by surface water erosion, which has resulted in the removal of bedrock, an increase in the river gradient, and disruption of the original shape of the river longitudinal profile. As we can see from Figure 9, the landscape at knickpoint D changes dramatically and shifts from one landscape type to another because of the large number of landslides in this area, which suggests that landslides can accelerate geomorphic evolution.

Stream Length Gradient Index (SL)

The SL analysis results provide a more detailed appraisal of local changes in the stream length gradient than can longitudinal river profiles. At the basin-scale, it is possible to identify regions where the relative change is greatest. The calculated SL values in this paper are between 0.3 and 25,133.9. Based on the ordinary kriging and natural breakpoint methods in ArcGIS, the SL values are interpolated and divided into five classes: very low (0.3–1,048.6), low (1,048.6–2,472.5), medium (2,472.5–4,690.9), high (4,690.9–9,384.7) and very high (9,384.7–25,133.9) (Figure 10). As shown in Figure 10, the SL values in most sub-basins are generally low, and only the SL values in sub-basins 1, 2, 10, 14, 19, 26, 28, 32, and 33 are higher than 4,690.9. The low SL values (0.3–1,048.6) denote that the river valley in the sub-basin is wide and flat, whereas the high SL values (9,384.7–25,133.9) indicate that the gradient of a river valley in the sub-basin is steep and the relative changes are large. In addition, the mapping reveals that the SL interpolation results have spatial clustering, with higher SL values distributed in areas far from the Indus and along faults.

Hypsometric Integral (HI)

The HI values and hypsometric curves in the 33 sub-basins calculated based on the DEM are shown in Figure 11. A regional pattern of geomorphic evolution is notable across the study area. Old-stage sub-basins (6, 7, 11, 12, 13, 14, 15, 17, and 18) occur in the center of the basin and along the Indus, where the HI values are low (HI < 0.33) and the hypsometric curves are concave. In sub-basins 3, 4, 5, 8, 9, 19, 20, 21, 24, 25, 26, 27, 28, 32, and 33, the hypsometric curves are S-shaped, and the HI values are between 0.35 and 0.53, which indicates that these sub-basins are in a mature stage of geomorphic evolution. The geomorphic evolution in nine marginal sub-basins (1, 2, 10, 16, 22, 23, 29, 30, and 31) far from the Indus is in the young stage, the HI values are high (HI > 0.54), and the curves are convex.

Figure 11. The state of geomorphic evolution based on the hypsometric integral value (HI) and hypsometric curve. Curves of the sub-basins to the north of the Indus are red, whereas curves of the sub-basins to the south of the Indus are blue.

Considering together the results of the river longitudinal profile (Figure 8) and SL (Figure 10) in the study area, we observed that the geomorphic evolution of the sub-basin where the Diamer-Bhasha Dam is planned and the sub-basin along the Indus are in the old stage. In contrast, the geomorphic evolution in the marginal sub-basins where many landslides, faults, and glaciers are distributed (Figure 1) are still in the young stage.

Exploring the Relationships Among Geomorphic Evolution, Landslides, and Tectonic Activity

The Kruskal-Wallis test is adopted to obtain insight into the relationship between lithology, HI, and SL, and Kendall’s rank coefficient method is used to reveal the relationships among AF, Vf, and HI in R software (Figure 12). The Kruskal-Wallis test revealed no correlations between lithology and HI underlying the watersheds (Figure 12a) and weak correlations between lithology and SL (Figure 12b). Kendall’s rank coefficient analysis revealed a positive correlation between AF and HI (Figure 12c) and a negative correlation between Vf and HI (Figure 12d). Furthermore, these relationships signify that tectonic activity can influence geomorphic evolution, while changes in lithology do not. However, the HI and SL values of the basins containing glaciers are high (Figures 12a and 12b), which indicates that the melting of glaciers can cause geomorphic evolution. Recently, with global climate warming, the melting of glaciers in this region has accelerated (Prasad et al., 2009). On the one hand, the melting of glaciers can induce landslides, debris flows, and other mountain disasters, which can change the state of the landscape; on the other hand, the melting of glaciers can also form runoff that can erode the landscape.

Figure 12. Relationships among lithology (a), tectonic activity (b, c, and d), and geomorphic evolution. Relationships are computed using the R language. Lithologies: Ka: Amphibolite, Cv: Andesite, Su: Anorthosite, Kb: Diorite, Sg: Gneiss, Y: Limestone, Pz: Metamorphic rocks, CC: Ultramafic rocks, (g) Glacier. L: Loess.

Figure 13 displays the relationships among AF, Vf, and landslides and the relationships between landslides and each of HI and SL based on Kendall’s rank coefficient method. The positive correlation between AF and landslides (Figure 13a), the negative correlation between Vf and landslides (Figure 13b), and the positive correlations between landslides and both HI and SL (Figures 13c and 13d) indicate that tectonic activity can induce landslides and that landslides affect geomorphic evolution.

Figure 13. Correlations between tectonic activity and landslides (a), (b) and the relationship between landslides and geomorphic evolution (c), (d).

The lack of correlations between HI and lithology and the weak correlations between lithology and SL signify that geomorphic evolution in the region is not controlled by changes in lithology. The positive correlations of AF with HI and landslides, the negative correlations of Vf with HI and landslides, and the positive correlations of landslides with HI and SL demonstrate that tectonic activities can cause geomorphic evolution and induce landslides that affect geomorphic evolution.

Discussion A Comparison of the Methods for Calculating SL Defined by Geostatistics

In this paper, we used a DEM with a 30 m resolution to extract the water systems and corrected them with a high-resolution (5 m) Google Earth image. Then, the drainages were segmented by automatic interactive tools to calculate the SL value in the study area (Keller et al., 1997). When we take ΔH = 30 m as a setting value, the calculated SL has significant correlations with HI and landslides, and the SL interpolation results are consistent with the results of the field investigation. In fact, when calculating SL, besides treating ΔH as a setting value, we can treat ΔL as a setting value (Hack, 1973). In previous studies, most scholars used a setting ΔL to calculate SL (El Hamdouni et al., 2008; Pedrera et al., 2009; Troiani & Della Seta, 2008), and some scholars compared the results of calculating SL with a setting ΔH value versus calculating it with a setting ΔL value (Troiani et al., 2014). They generally concluded that calculating SL by adopting a setting ΔL value is more consistent with the tectonic-geomorphic features observed in field investigations. Initially, in this paper, the SL value was calculated by adopting a setting ΔL value. However, the calculated SL showed a decreasing trend from the inlet to the outlet along the flow direction of the river, whereas the field investigation found that the gradient was higher at the middle and outlet portions of the river. Consequently, in different research areas with differences in basin size, data resolution, topographic features, geological environment, and other factors, we should choose the most suitable method to calculate SL according to local conditions.

Regional Dependence of Lithology on Geomorphic Evolution

Traditional studies on tectonic-geomorphic features have interpreted the lithologic change as evidence of geomorphic evolution, which has been observed and discussed by many authors (Bowman et al., 2007; Mather, 2000; Stokes et al., 2002). Through comprehensive analysis of landscape development and geomorphic deformation in the Moravka River, Škarpich et al. (2013), from the perspective of contemporary and historical factors, considered key factors causing geomorphic evolution to be human activities (road location and dam construction), changes in land cover and land use, and lithology changes. However, different lithologies in different areas may not all play a role in geomorphic evolution (Hurst et al., 2013; Lifton & Chase, 1992). Mathew et al. (2016) quantitatively and qualitatively analyzed the tectonic-geomorphic features on the sub-catchment scale by quantizing geomorphic indices and systematically demonstrated for the first time the response of the geomorphic evolution in Negeri Sabah to tectonic activity. Similarly, Gao et al. (2013) evaluated the recent uplift of the northeastern margin of the Tibetan Plateau, China, with geomorphic indices and geodetic leveling data on different time scales obtained by their investigation and found that the geomorphologic evolution of the region was not affected by lithology. In this paper, the study area is located in the Himalayan seismic belt, which is the largest seismic belt in the world. Since the Cretaceous-Tertiary, frequent tectonic activity has caused geomorphic evolution (Debon, 1995; Ravikant et al., 2009). The positive correlations of AF with HI and landslides, the negative correlations of Vf with HI and landslides, and the positive correlations of landslides with HI and SL indicate that tectonic activities can cause geomorphic evolution and induce landslides that affect geomorphic evolution (Figures 12 and 13), and the weak correlations of lithology with HI and SL indicate that geomorphic evolution is not controlled by lithology (Figures 12a and 12b). Therefore, lithology does not ubiquitously play a role in geomorphic evolution and is regionally dependent on geomorphic evolution.

Conclusions

This study detects tectonic-geomorphic features and the correlations among geomorphic evolution, landslides, and tectonic activity for a core area of the planned Diamer-Bhasha Dam using geomorphic indices and landslides. The tectonic activity in the sub-basin along the Indus is in a stable state. In the marginal regions far from the Indus, severe tectonic tilting and frequent tectonic uplift have led to extremely active tectonic activity in these regions. In the sub-basin along the Indus and the sub-basin where the Diamer-Bhasha Dam is planned, geomorphic evolution is in the old stage. The lack of correlations between HI and lithology, the weak correlations between lithology and SL, the positive correlations of AF with HI and landslides, and the negative correlations of Vf with HI and landslides demonstrate that tectonic activities can influence geomorphic evolution and induce landslides, whereas changes in lithology do not. In addition, the positive correlations of landslides with HI and SL signify that landslides can impact geomorphic evolution. We anticipate that the conclusions of this work will be a theoretical reference for the construction of the project (Diamer-Bhasha Dam).

Acknowledgments

This work was funded by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (grant no. 2019QZKK0902), International Science & Technology Cooperation Program of China (grant no. 2018YFE0100100), Natural Science Basic Research Program of Shaanxi (grant no. 2021JC-40), National Natural Science Foundation of China (grant no. 41771539), Strategic Priority Research Program of Chinese Academy of Sciences (grant no. XDA 20030301), and International Partnership Program of Chinese Academy of Sciences (grant no. 131551KYSB20160002).

Data Availability Statement

Data used to generate plots in the paper are available in the following public domain repository: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/KAJX9B