1 Introduction

The Karakoram highway (KKH) traverses the confluence of the world’s three youngest mountain systems (the Himalaya, the Karakorum and the Hindu Kush). It is the highest international highway in the world, with the average altitude reaching more than 3000 m and the highest point (4733 m) at Khunjerab (Liao et al. 2012). Due to the geological structure, the tectonic movement along the KKH is highly active, with strong seismic activity, rapid crustal uplifting and significant terrain undulation (Sobel and Dumitru 1997; Li et al. 2016). In addition, there are large amounts of loose solid material, abundant ice and snow melt water runoff, and frequent rainstorms in the region (Owen and Derbyshire 1989; Shroder 1998; Ge et al. 2015). Under the joint influence of internal and external forces of the Earth’s crust, debris flows are commonly seen and active along the KKH. On 28 July 2021, the intensification of glacier/snow melt, coupled with the heavy rainfall in local mountains, caused multiple debris flows in several sections of the KKH. Traffic was disrupted in both directions, nearly 200 vehicles were blocked, and more than 400 people were stranded. Frequent debris flows are disrupting road transportation and hindering socioeconomic development in this region (Qing et al. 2020; Ali et al. 2021).

The West Kunlun Mountains, where the Kashgar-Khunjerab section of the KKH is located, is one of the regions with the most and largest glaciers in the world at low and middle latitudes (Schoenbohm et al. 2014; Zhang et al. 2016; Li et al. 2020). Due to the fragile ecosystems and sparse vegetation coverage, this area is among the most sensitive and vulnerable regions to climate change (Yao et al. 2022; Pei et al. 2023). A warming climate has the most significant impact on glaciers and ice-marginal areas and geomorphic processes such as debris flows are particularly sensitive to the response of climate change (Chiarle et al. 2007). The overall susceptibility of debris flows in the study area will increase with climate change, and the high susceptibility watersheds will reach 13.5% of the total (Zou et al. 2021). Therefore, it is crucial to assess the risk of debris flows in a warming climate and to propose strategies to deal with debris flow disasters under different warming scenarios, which can facilitate road route planning to reduce the threat of debris flows and support integrated disaster risk management. However, the risk of debris flows in the KKH in the context of climate change has not been adequately explored because of the lack of long-term historical observations and the related methods of analysis.

Previous studies of debris flow activity and disaster risk in the context of climate change mainly focused on the debris flow hazard, risk analysis of different rainfall return periods and the debris flow characteristics under future temperature and precipitation scenarios. Staffler et al. (2008) used the two-dimensional flood routing model to simulate the depth and impact range of debris flow accumulation under increased rainfall scenarios in the Tyrol region of Italy. Stoffel et al. (2014) used a long-term dataset of past debris flows and future climate scenario data to analyze the impact of climate change on the frequency and magnitude of debris flows in the Swiss Alps from 2001 to 2100. Ouyang et al. (2019) used the Massflow model to analyze the debris flow risk under different rainfall recurrence intervals in Niwan Gully, Gansu Province of China.

However, the comprehensive and systematic assessments of debris flow hazard and road vulnerability with consideration of temperature factors are rarely seen and thus it is difficult to recognize and quantify the magnitude of debris flow risk in road systems under future climate change. In this study, we collected the meteorological observation data of 1957–2017, future climate model data of 2009–2050 and debris flow activity data along the KKH to assess the debris flow risk under future climate change.

The three objectives of this study are: (1) to analyze the spatiotemporal changes of air temperature and precipitation; (2) to construct a road debris flow risk assessment model that includes an improved key index screening and weight assigning method and separately analyze rainfall debris flow and glacier-rainfall mixed debris flow; (3) to assess the debris flow risk under future climate change. We hope our findings can help improve the debris flow risk management of linear projects under future climate change and help understand the impact of climate change on global debris flow hazards. The remainder of this article is organized as follows. The overview of the study area is introduced in Sect. 2. The data and methods are presented in Sect. 3. The results are analyzed in detail in Sect. 4. Some discussions are presented in Sect. 5. Finally, Sect. 6 summarizes the whole text.

2 Study Area

The Kashgar-Khunjerab section of the KKH (the study area) is located within 74.78°E–76.10°E and 36.72°N–39.63°N. It traverses the western edge of the Tarim Basin and the Ghez River valley, and arrives at the border crossing point of Khunjerab, with a total length of 413 km (Fig. 1). The study area is of typical continental alpine climate with large diurnal temperature difference and abundant sunlight. Winters are cold and summers are cool, with a short warm season and a long cold season. In spring, there is less rainfall, but high winds and sandstorms are often seen. The peaks of rainfall and heat both appear in summer. Additionally, there is a distinct dry and wet season (winter), which is conducive to physical weathering. The precipitation process in the study area is mainly controlled by the westerly wind belt (Aizen and Aizen 1997; Yao et al. 2012), which brings abundant water vapor from the Mediterranean, the Black Sea and the Caspian Sea only to the western edge of the Pamir, causing aridification in the eastern Pamir (Fuchs et al. 2014). The observation data from the Tashkorgan national meteorological station (37.77°N, 75.23°E, altitude: 3090 m) during 1957–2017 show that the average annual precipitation (AAP) in the study area is 82.6 mm and the average annual temperature (AAT) is 4.1 °C. Precipitation is mainly concentrated in summer (June–August), accounting for 52% of the annual precipitation, and the average summer temperature is 15.8 °C. The AAP and AAT change rates in the study area are 5.70 mm/10a and 0.33 °C/10a, respectively, showing a clear warming and wetting trend. Since 1957, there has been a notable increasing trend of extreme precipitation and heat events.

Fig. 1
figure 1

Overview of the study area: Landform 1 represents alluvial plain; Landform 2 represents moderate high mountain; Landform 3 represents high mountain and very high mountain canyon; Landform 4 represents high mountain river valley; Landform 5 represents high mountain canyon

Full size image

The study area is strongly uplifted and deformed due to the continuous collision between the Indian and Eurasian plates, resulting in frequent and severe seismic activities (Peltzer et al. 1989; Wright et al. 2004; He et al. 2016). More than 1000 earthquakes have been recorded since 1889, including 29 earthquakes of magnitude 5 and above, 5 earthquakes of magnitude 6 and above, and 1 earthquake of magnitude 7 (Li 1994; Ma 2015; Wang 2020). The main water systems in the study area are the Ghez and Tashkorgan Rivers (Fig. 1). They are snowmelt-type rivers, with runoff mainly recharged by alpine ice and snow melt water. The typical flood season is from June to September, with the maximum annual runoff in July and August (Ren et al. 2018). There are a large number of active glaciers (Fig. 1). With tectonic movement, strong erosion by rivers, glacial action, weathering and denudation, the surface of the study area is highly undulated and the rock body is fractured, providing extremely favorable energy and material conditions for debris flows. Many sections of roadways are exposed to debris flows, which seriously threaten the traffic safety in the study area.

3 Data and Methods

This section introduces the meteorological, topographic, geomorphological, geological, seismological, hydrological, vegetation and soil data along the KKH and demonstrates the debris flow risk analysis method. A detailed description of meteorological data and data analysis method is given in Sect. 3.1, vulnerability and hazard assessment methods are presented in Sect. 3.2.

3.1 Data

3.1.1 Meteorological Data

To better analyze the spatiotemporal changes of air temperature and precipitation, the data covering the entire study area for both historical and future periods must be obtained. However, the meteorological stations in the study area are scarce and they can only provide historical meteorological data, which is insufficient for future analysis. Therefore, we used the China Meteorological Forcing Dataset (1979–2018) (Yang and He 2019) and the Future Climate Projection over Northwest China Based on RegCM4.6 (2007–2099) (Pan and Zhang 2019) issued by the National Tibetan Plateau Data Center to represent the climate of the study area for historical and future periods, respectively. The China Meteorological Forcing Dataset (1979–2018) uses the fusion of surface observation, reanalysis data and satellite remote sensing data. The accuracy is between that of the surface observation and satellite remote sensing, and is better than the accuracy of internationally available reanalysis datasets (He et al. 2020). The Future Climate Projection over Northwest China Based on RegCM4.6 (2007–2099) uses the latest regional climate model RegCM4.6 for Northwest China, where the study area is located, to correct bias, and the data accuracy is higher than that of general global climate models (Pan et al. 2020). It has been noted that the representative concentration pathway (RCP) 8.5 is the most suitable scenario for studying climate change impacts until 2050 and is the closest to observation under existing global climate policies. Therefore, the RegCM4.6 model output data under the RCP8.5 scenario were selected to simulate future climate change scenarios under extreme conditions (Schwalm et al. 2020). These two datasets have the highest accuracy in both spatial and temporal resolutions (Table 1).

Table 1 Details of the two climate datasets
Full size table

To further improve the spatial resolution of climate projection data, the delta downscaling method was used to downscale the daily mean temperature and daily precipitation from the future climate projection data (with a spatial resolution of 25 km) to the weather element-driven dataset of the historical period (with a spatial resolution of 10 km). This method has been widely used in regional and global climate model downscaling studies because it retains important statistical features such as spatial correlation (Hay et al. 2000; Ramirez-Villegas and Jarvis 2010; Sarr et al. 2015). To obtain more accurate future climate projections, this method was used to analyze the overall climate change in the study area.

3.1.2 Debris Flow Data

The multi-source and multi-temporal remote sensing images (Landsat, Gaofen 1, Gaofen 2, Google Earth) were interpreted to obtain the distribution characteristics of debris flow gullies in the study area and analyze the spatial relationship between roads and debris flow gullies. In this study, debris flow gullies and non-debris flow gullies were determined by the source material and morphology of the gullies, combined with consideration of the morphology and material composition of the accumulation fans. Specifically, obvious debris flow gullies and non-debris flow gullies were distinguished through remote sensing interpretation. For gullies that could not be determined with certainty, field investigations were conducted to further explore the material and morphology of the gullies and accumulation fans. Finally, 257 debris flow gullies were identified in the study area based on remote sensing and field surveys. There are 230 debris flow basins in non-glacial areas and 57 debris flow basins in glacial areas (Fig. 2). On this basis, further field investigations were conducted to obtain more comprehensive and accurate data. During the fieldwork, the geological, geomorphological, meteorological, hydrological, vegetation, soil, glacier and other environmental background conditions along the KKH Kashgar-Khunjerab section were systematically investigated. The current status and impact of debris flow activities on the road safety were also comprehensively examined. The distribution and characteristics of different types of elements at risk were carefully recorded and surveyed as well. Therefore, complete information on the development of debris flow disasters along the entire route has been obtained.

Fig. 2
figure 2

Distribution of debris flow gullies and non-debris flow gullies in the study area

Full size image

3.2 Debris Flow Risk Analysis Method

When hazard, exposure and vulnerability are taken into account, debris flow risk can be assessed according to the following equation (Varnes 1984):

$${\text{R}} = H \times E \times V$$
(1)

where R is the debris flow risk, H is the hazard level, E is the exposure level and V is the vulnerability level. In this study, we chose the area of road as the exposure index.

3.2.1 Vulnerability Analysis Method

Vulnerability of elements at risk is affected by a combination of physical, social, economic and environmental factors. Measures of vulnerability can include the structure type or other settings that could be affected by debris flows. In this study, vulnerability indicators include exposure of affected elements (represented by distance to debris flow gully), structural characteristics of affected elements (represented by structure type and comprehensive disaster resilience) and post-disaster recovery capability (represented by construction cost and distance to city or town). Each indicator was normalized so that its value range is between 0 and 1. The analytic hierarchy process was used to assess the vulnerability because it can both structure complex problems and combine quantitative and qualitative attributes (Wedley 1990). This method is especially suitable for cases involving socioeconomic impact assessment (Ramanathan 2001). Taking 50 m of the road as the evaluation unit, the vulnerability indicators of each evaluation unit were calculated and then weighted to obtain the overall vulnerability index.

3.2.2 Hazard Analysis Method

To assess the hazard of debris flows, factors reflecting the energy, source material, and triggering conditions of debris flows were selected. Specifically, topographic and geomorphological (including slope and surface elevation difference), geological (including geotechnical fragility and fault distance effect), seismological (represented by earthquake intensity), hydrological (represented by river network density), vegetation (represented by surface exposure) and soil (represented by content of clay and sand) data were collected (Table 2). In addition, some meteorological data, including AAT, average annual daily maximum temperature (AADMT), AAP and average annual daily maximum precipitation (AADMP), were acquired through downscaled climate data as displayed in Sect. 3.1.

Table 2 Data source of factors along the Karakoram highway (KKH)
Full size table

Slope, surface elevation difference, and river network density were derived through the spatial analysis of ASTER PALSAR Digital Elevation Model, with a resolution of 12.5 m. Geotechnical fragility and fault distance effect were quantified by the lithological and fault data from the 1:250,000 geological map, respectively. The ultimate compressive strength of rocks proposed by Sun and Chen (2007) was used to reflect the compressive capacity of different lithologies in the study area. The greater the compressive strength of a rock is, the lower its geotechnical fragility is. In order to be consistent with other indicators, the compressive strength values of rocks were inverted for high and low values to obtain the geotechnical fragility indicator that contributes positively to debris flow hazard. Similarly, the Euclidean distance of each grid in the study area from the fault was calculated and inverted to obtain the fault distance effect indicator. Earthquake intensity was also calculated through the interpolation of historical earthquakes of magnitude 5 and above. An equation proposed by Xu et al. (2016) was utilized:

$$M=0.540{I}_{e}+0.530\mathrm{log}h+1.473$$
(2)

where \(M\) is the earthquake magnitude, \({I}_{e}\) is the earthquake intensity, and \(h\) represents the earthquake depth.

Surface exposure was calculated by “1—fractional vegetation cover (FVC)”. The FVC data were obtained from the FVC Dataset of Remote Sensing for Ecological Assets Assessment in Qinghai-Tibet Plateau (Liu 2019). Contents of clay and sand were acquired through A China Soil Characteristics Dataset (Dai 2019).

In this study, an improved random forest algorithm was adopted to analyze the hazard of debris flows. The method effectively combines indicator importance and classification error when selecting the key factors and assigning indicator weights (Tibshirani 1996; Wolpert and Macready 1999; Han et al. 2016). It optimizes the trade-off between the number of factors and classification accuracy, thereby improving the rationality of indicator selection and weight assignment. The screening of key factors is stopped when the minimum out of bag (OOB) error is achieved, and then the weight of each factor is determined by the average of Mean Decrease Accuracy and Mean Decrease Gini. The random forest algorithm is implemented in the R language platform by using the random forest package.

According to remote sensing interpretation and field investigation, there are 230 debris flow basins and 169 non-debris flow basins in non-glacial areas, and there are 57 debris flow basins and 29 non-debris flow basins in glacial areas (Fig. 2). For the non-glacial watershed units, 11 indicators were selected, except for the AAT and the AADMT. For the glacial watershed units, all 13 indicators were selected. All of these factors were processed as raster data with a resolution of 12.5 m. Each factor was normalized so that its value range is between 0 and 1 and calculated separately in each watershed. The average value of each factor in each watershed was used as the sample data of the random forest. The selection of key factors and the assignment of weights were carried out separately for rainfall debris flows (located in the non-glacial areas) and glacier-rainfall mixed debris flows (located in the glacial areas).

The hazard evaluation models of two different types of debris flows were acquired through random forest and then utilized to calculate the debris flow hazard indices in the study area for the baseline period and three future time periods. Taking 50 m of the road as the evaluation unit, we calculated the average value of debris flow hazard within 10 km of each evaluation unit to determine the debris flow hazard value of the road section—10 km was used here because it is the average length of debris flow gullies in our study area according to remote sensing interpretation and field investigation.

On the basis of vulnerability and hazard analyses, risk index was calculated for each evaluation unit. By using the natural breakpoint method, the hazard, vulnerability and risk indices of each evaluation unit were classified into four categories: very low, low, medium and high in the baseline period and three future time periods.

4 Results

In this section, we present the results obtained from the analyses outlined in the previous sections for the study area. The results mainly include the climate change characteristics and debris flow hazard and risk during the baseline period, 2025s, 2035s and 2045s. For road vulnerability to debris flows, we present the results during the baseline period.

4.1 Spatiotemporal Distributions of Temperature and Their Changes

During the baseline period, the AAT and AADMT in the study area ranged from − 18.64 to 12.91 °C and from 0.40 to 32.63 °C, respectively. Figure 3 shows that the Kashgar-Oytak section has the highest AAT. The Kongur Mountain in the Oytak-Bulunkol section has relatively low temperature due to its high altitude and extensive glacier coverage, while the rest of the Oytak-Bulunkol section has relatively high temperature. Similarly, the temperature is low in the Kongur Mountain and Muztagh Mountain in the Bulunkol-Tashkorgan section, while it is higher in the rest of the section.

Fig. 3
figure 3

Average annual temperature (AAT) in the baseline period (a), 2025s (b), 2035s (c) and 2045s (d)

Full size image

The AAT in the study area is projected to increase continuously over the next three periods, with a maximum of 1.29 °C by the 2045s (Fig. 3). The AADMT is projected to increase up to 2.01 °C by the 2045s. There is an increasing trend of the AAT and AADMT in the study area in the future periods and the spatial distributions of the increases vary. The greatest warming rate of the AAT is expected in the Kashgar-Oytak section and the southernmost region near Khunjerab, while the smallest warming rate is projected to be in the region near Tashkorgan (Fig. 3). The spatial variation characteristics of the AADMT demonstrate an opposite trend to the AAT, with the largest increase in the Bulunkol-Tashkorgan interval.

4.2 Spatiotemporal Distributions of Precipitation and Their Changes

Figure 4 shows that the AAP and AADMP in the study area during the baseline period ranged from 107.83 to 371.42 mm and from 15.54 to 57.30 mm, respectively. The AAP in the urban area of Kashgar is the smallest, followed by the open plateau surface section near Tashkorgan. The precipitation in the study area increases with altitude, and the topographic effect is significant (Sun et al. 2020). Consequently, the highest levels of AAP are found in the Kungey, Kongur and Muztagh Mountains, as well as in the southern part of the Sarykol Ridge.

Fig. 4
figure 4

Average annual precipitation (AAP) in the baseline period (a), 2025s (b), 2035s (c) and 2045s (d)

Full size image

The overall precipitation in the study area, except for the southernmost Khunjerab region, is expected to increase, with a feature of increasing first and then decreasing, reaching the maximum in the 2035s. The AAP in the study area increases by 11.16%, 14.62% and 11.38% in the 2025s, 2035s and 2045s, respectively. The AADMP in the study area shows a similar spatiotemopral changes with AAP. The greatest increases of the AAP and AADMP are found in the Oytak-Bulunkol interval and the southern Tashkorgan-Khunjerab region. In particular, the Oytak-Bulunkol section will experience the most significant increase in the region.

In summary, the temperature and precipitation in the study area show an overall increasing trend in the future periods. The temperature increases continuously with time, while the precipitation increases first and then decreases, reaching its maximum in the 2035s. For the Oytak-Bulunkol section, where the debris flow hazard is the most severe, the AAP and AADMP are all projected to increase significantly in the future.

4.3 Vulnerability to Debris Flows

Five vulnerability assessment indicators were calculated for each unit of assessment based on the data collected from the highway, debris flow gullies and major towns in the study area. Each indicator was normalized so that its value range is between 0 and 1. Analytic hierarchy process was then used to analyze the vulnerability indicators and obtain the weights of road vulnerability indicators (Table 3).

Table 3 Weights of vulnerability assessment factors
Full size table

In the fieldwork, a total of 40 bridges, 1 tunnel and 261 culverts were counted along the entire highway. In accordance with the Technical Specifications for Highway Bridge and Culvert Construction (Table 4), the bridges were divided into four categories: super major bridges, major bridges, medium bridges, and small bridges, which have been used to analyze the comprehensive disaster prevention capacity.

Table 4 Criteria and results for bridge classification
Full size table

Figure 5 shows that the very low-vulnerability section is mainly located in the Kashgar-Oytak area, which is characterized by flat terrain and a lack of debris flow gullies along the road. Furthermore, the area has strong post-disaster recovery ability for its proximity to major towns. The low-vulnerability section located near Oytak and Tashkorgan also has strong post-disaster recovery capability and is far from the debris flow gullies. The medium-vulnerability section takes up the highest percentage and is located close to debris flow gullies. This section is characterized by a variety of road structures, such as bridges and culverts, which are more densely distributed and thus more fragile. The high-vulnerability section is mainly located between the Kungey Mountain and the Kongur Mountain. This section is situated in a high mountain canyon area and is particularly dangerous. It is the closest to the debris flow gullies and it mostly runs across the debris flow accumulation fans. In order to reduce the hazard of debris flows on the road, fragile and high-cost structures such as the super major bridge, the tunnel and most densely packed culverts were constructed in this section. However, the improved road resilience cannot totally resist the catastrophes brought by extreme weather-triggered mega-hazards. Coupled with the intensive distribution of such fragile projects, the structure vulnerability in this section is at the highest level in the study area.

Fig. 5
figure 5

Vulnerability level of the Karakoram highway (KKH)

Full size image

4.4 Debris Flow Hazard

After the sample data of non-glacial watershed units and glacial watershed units were analyzed, the key factors were screened to effectively distinguish the debris flow gullies from the non-debris flow gullies (Table 5). The results show that the sand content and gully density factors are the least important for the non-glacial watershed unit and can be excluded to obtain higher classification accuracy, while the other factors can help to improve the classification accuracy and should be retained. For the glacial watershed units, the highest accuracy classification results are achieved by removing the clay content, gully density, sand content and AAT.

Table 5 Factor selection in non-glacial and glacial areas
Full size table

On this basis, the weight of each key indicator was separately calculated for two types of debris flows (rainfall-type debris flows and glacier-rainfall mixed debris flows) and the following risk evaluation models for the two types were obtained (Table 6).

Table 6 Hazard assessment model of two types of debris flows
Full size table

In the hazard evaluation models of these two types of debris flows, Hazard refers to debris flow hazard, S stands for slope, EI stands for earthquake intensity, FDE stands for fault distance effect, RED stands for relative elevation difference, GF stands for geotechnical fragility, SE stands for surface exposure and CC stands for clay content.

For the rainfall-type debris flows, the AADMP indicator has the greatest weight, implying that extreme rainfall events have the greatest influence on the risk of such debris flow gullies. For the glacier-rainfall mixed debris flows, the EI and FDE indicators have the largest weights, indicating that seismic intensity and fault distance effects have a greater influence on debris flow hazard, as they affect the recharge of solid material sources. In addition, the effect of temperature on the glacier-rainfall mixed debris flows is not as great as that of precipitation, suggesting that glacial meltwater alone in the region is not sufficient to trigger debris flows, and the coupling effect of rainfall is necessary to potentially stimulate glacier-rainfall mixed debris flows, which is consistent with the results of a previous study (Wang 1987). This is partly due to the low average temperatures at high altitudes, making glaciers in this region likely to be more sensitive to changes in precipitation than temperature changes (Derbyshire 1981; Owen and Benn 2005).

The results (Fig. 6) show that during the baseline time period, the very low-hazard areas are mainly distributed in the Kashgar-Oytak and Bulunkol-Tashkorgan sections. Since there are no debris flow gullies on either side, the probability for debris flows to occur in these sections is very low, and the hazard level is expected to remain unchanged in the future time periods. The low-hazard sections are mainly located near Oytak during the baseline time period and diminish in the future time periods as most of them will change to medium-hazard sections in the future, particularly in the 2025s.

Fig. 6
figure 6

Debris flow hazard level of the whole Karakoram Highway in the baseline period (a), 2025s (b), 2035s (c) and 2045s (d)

Full size image

Medium-hazard sections are distributed throughout the highway, with the highest concentration in the Bulunkol-Tashkorgan interval. This region shows a slight increase of hazard in the 2025s, followed by a continuous decrease in the 2035s and 2045s, with little overall change. The high-hazard sections are mainly located in the upper reaches of the Ghez River and the Tashkorgan-Khunjerab interval and the proportion keeps increasing over the three future time periods (Table 7). Therefore, highly hazardous debris flows will have the largest impact on the highway in the 2045s.

Table 7 Different hazard level proportion in the baseline and future periods
Full size table

4.5 Debris Flow Risk

The results (Fig. 7) demonstrate that the debris flow risk varies among different sections under future climate change. The risk level of debris flows in the Kashgar-Oytak section is among the lowest in the entire line, owing to the lowest level of debris flow vulnerability and hazard in both the baseline and future time periods. In the Oytak-Bulunkol section, the overall risk is low in the mid-alpine area downstream of the Ghez River, due to the medium vulnerability and low hazard. However, some sections may become high-risk sections in the future. The high mountain valley area upstream of the Ghez River has an overall high risk level due to its high hazard and vulnerability levels, with a significant increase in high-risk sections in the future. The overall risk is low in the Bulunkol-Tashkorgan section as debris flow gullies are few there. However, most of the low-risk sections are expected to become medium-risk sections in the near future due to an increase in debris flow gully hazards. Most parts of the Tashkorgan-Khunjerab section are predicted to become medium-risk sections, and the percentage of high-risk sections is expected to increase significantly in the future.

Fig. 7
figure 7

Debris flow risk level of the whole Karakoram highway in the baseline period (a), 2025s (b), 2035s (c), 2045s (d)

Full size image

Table 8 shows that the highest overall risk of the highway is in the 2025s, with the smallest proportion of low-risk sections, the largest proportion of medium-risk sections, and a relatively high proportion of high-risk sections. It is because in the 2025s, the AADMP and the AAP increase in the downstream of the Ghez River and the Tashkorgan-Khunjerab interval, where rainfall debris flows are densely distributed, are at the highest level for the whole route and the whole time period, resulting in a higher overall debris flow risk level for the whole route of the highway in this time period.

Table 8 Different risk level proportion in the baseline and future periods
Full size table

5 Discussion

The results indicate that both debris flow hazard and risk are likely to increase under future climate change, with considerable spatiotemporal variability. This variability is primarily determined by different mechanisms through which the climate change impacts debris flows in non-glacial and glacial areas. In the non-glacial areas, climate change will affect debris flow activity mainly through an increase in the frequency and magnitude of extreme precipitation events (IPCC 2013, 2021; Westra et al. 2014; Donat et al. 2016). The increasing trend of AAP and extreme precipitation in the study area in the future will bring more water source for debris flows and stimulate larger and more frequent debris flows. In addition, precipitation has a significant impact on the stability of slopes. Rainwater may infiltrate into the sliding surface or shear zone, and alter the soil properties of the shear zone, thus destabilizing the slope soil (Rahardjo et al. 2007). The soil in the study area has a low clay content and a high sand content, which is highly conducive to rainfall infiltration. Consequently, an increase in future rainfall intensity in the area will deepen the infiltration and increase soil erosion intensity, and thus resulting in an increase in the initial sources of material and kinetic energy (Cui et al. 2005; Iverson 2014).

In addition to the two precipitation factors mentioned above, the glacier-rainfall mixed debris flow is also influenced by the AADMT. The overall debris flow hazard level reaches its maximum in the 2045s when the AADMT is the highest of all time periods. Global warming has accelerated the retreat of glaciers, releasing a large amount of liquid water and increasing glacier runoff (Huss and Hock 2018). It has been demonstrated that climate warming has resulted in a sustained rise in glacial runoff in the Tarim Basin, where the study area is situated, over the last few decades (Yao et al. 2004; Hao et al. 2008; Xu et al. 2011). The analyses indicate that the study area is projected to experience an increase in both average temperature and extreme temperature in the future. Therefore, the glacial runoff in the area is expected to increase, providing more water source for debris flows. In addition, higher temperature will cause a rise of the 0 °C isotherm, leading to more precipitation in the form of liquid water, which can result in debris flows (Beniston 2005).

Global warming has also caused significant retreat of glaciers and degradation of permafrost, exposing a large amount of moraines. With drastic changes in the hydrogeological and soil mechanical properties of the moraines, the soil cohesion and strength have been reduced. Under well-watered conditions, these moraines can quickly reach saturation and initiate destabilization, becoming the main recharge source of debris flows (Larsson 1982). In addition, the degradation of permafrost has exacerbated the retrogressive erosion of debris flows and infiltration of water sources into loose moraines, providing more materials for debris flows and further amplifying their scale (Zimmermann and Haeberli 1992). Moreover, glacial movement has caused erosion and deep cuts in the slope body under the joint influence of glacier ablation and thinning, resulting in slope destabilization due to the loss of support. It has led to ice avalanches, rock avalanches, and landslides, which not only provide water source for debris flows but also bring abundant material sources (Evans and Clague 1994). Most of these material sources are located at higher elevations. Once they are activated, more potential energy will be converted into greater kinetic energy, thereby forming larger-scale debris flows.

In this study, a comprehensive consideration of climate change impacts on debris flow risk along the Karakoram Highway has been made by selecting different factors for different types of debris flows. A systematic debris flow risk assessment method that includes an improved key indicator screening and weight assigning method has been proposed and can be applied in other climate-sensitive alpine regions.

This work presents a regional assessment of debris flow risk to an alpine highway where selected factors were used to analyze the hazard of debris flows. This kind of large-scale assessment does not reflect the physical process of debris flow activities, so there is no refined information such as inundation extent and depth. Further studies are needed to numerically simulate details of debris flow activities under climate change based on new techniques such as big data and super computing. In this way, debris flow risk can be assessed more accurately and provide more effective information for policy makers and related government departments.

6 Conclusion

In this study, the debris flow gullies along the KKH were analyzed through field survey and remote sensing interpretation to assess the debris flow risk in a warming climate and to propose strategies to deal with debris flow disasters under different warming scenarios. Historical meteorological data and future climate model data were used to analyze the future climate change of the study area and its impact on debris flow activities for the baseline time period (2009–2018) and three future time periods (2025s, 2035s and 2045s) under the RCP8.5 scenario. The results show that:

  1. (1)

    The overall AAT and AADMT in the study area show continued growth in the next three periods with adverse spatial variation trends; the AAP and AADMP are all found to be generally increasing with a characteristic of first increasing then decreasing, reaching the maximum in the 2035s. With similar spatial variations, these two precipitation factors show greatest increases in the Oytak-Bulunkol and the Tashkorgan-Khunjerab intervals. The increasing temperature and precipitation in the study area will bring more water and material sources for debris flows, and increase potential energy of material sources, consequently having a direct impact on debris flow activity in the future time periods.

  2. (2)

    The overall debris flow hazard and risk along the highway will increase in future periods, the debris flow hazard is expected to be the highest in the 2045s while the debris flow risk reaches the peak in the 2025s. Particularly, the debris flow risk and road vulnerability of the high mountain valley section in the upper reaches of the Ghez River are the highest on the entire line, warranting priority attention from the relevant government departments.

To ensure the safe and effective operation of the KKH, prevention and adaptation measures to cope with climate change must be formulated. The debris flow mitigation along the KKH under future climate change conditions should be tailored to various risk levels. For low-risk sections, the overwater pavement should be prioritized and the adequate dredging and restoration following debris flows should be ensured. For high-risk sections, the focus should be on the construction of a monitoring and early warning system, which includes rainfall, temperature, soil moisture content, mud level and other monitoring equipment.