A conceptual model for evaluating the stability of high-altitude ice-rich slopes through coupled thermo-hydro-mechanical simulation 通过热-水-机械耦合模拟评估高海拔富冰斜坡稳定性的概念模型
Mingdong Wei, Limin Zhang *, Ruochen Jiang 魏明东，张利民 *，蒋若晨Department of Civil and Environmental Engineering, The Hong Kong University of Science and Technology, Hong Kong, China 中国香港，香港科技大学土木与环境工程系
A R T I C L E I N F O
Keywords: 关键词：
Ice-soil mixture 冰土混合物
Slope stability 斜坡稳定性
Thermo-hydro-mechanical model 热-水-机械模型
Climate change 气候变化
Glacial ablation 冰川消融
Abstract 摘要
A B S T R A C T The collapse of glacial, permafrost, and ice-rich moraine slopes in high-altitude mountainous areas not only threatens downstream residents and infrastructure but also displaces ice mass to lower and warmer elevations, accelerating glacial ablation to some extent. Despite being considered rare, ice-rich slopes collapse more frequently than commonly thought due to climate change. For example, two adjacent mountain glaciers in the Aru Range of the Tibetan Plateau, characterized by large volumes ( and , respectively) and low surface slope angles ( and , respectively), collapsed surprisingly in 2016 . While the mechanisms behind these collapses have garnered broad attention, the ability to quantitatively assess the instability of ice-rich slopes remains limited due to the complex interplay of multi-physical processes. Taking the Aru glacier collapses as reference cases, this paper presents a conceptual model, implemented through coupled thermohydro-mechanical simulation, to evaluate the stability of high-altitude ice-rich slopes due to climate change, rainfall and ice ablation. Results indicate that the methodology captures the effects of temperature change, rainfall and meltwater on the instability events well, demonstrating promising potential in evaluating potential collapse zones of ice-rich slopes similar to the Aru glaciers. Furthermore, the role of climate change in the wellknown Aru events is demonstrated using a state-of-the-art global climate reanalysis dataset. Findings reveal that the increase in liquid water infiltrating the Aru glaciers since 2010 was a critical factor leading to the instability events. A B S T R A C T 高海拔山区的冰川、永久冻土和富冰碛斜坡崩塌不仅威胁下游居民和基础设施，而且还会将冰体移至海拔较低、温度较高的地方，在一定程度上加速冰川消融。尽管富冰斜坡被认为是罕见的，但由于气候变化，富冰斜坡崩塌的频率比通常想象的要高。例如，青藏高原阿鲁山脉两座相邻的高山冰川分别具有体积大（ 和 ）和表面坡角小（ 和 ）的特点，却在2016年出人意料地发生了崩塌。虽然这些塌方背后的机制引起了广泛关注，但由于多种物理过程的复杂相互作用，定量评估富冰斜坡不稳定性的能力仍然有限。本文以阿鲁冰川崩塌为参考案例，提出了一个概念模型，通过热-水-机械耦合模拟来评估气候变化、降雨和冰消融导致的高海拔富冰斜坡的稳定性。结果表明，该方法能很好地捕捉温度变化、降雨和融水对不稳定事件的影响，在评估类似阿鲁冰川的富冰斜坡潜在崩塌区方面显示出巨大的潜力。此外，利用最先进的全球气候再分析数据集证明了气候变化在著名的阿鲁事件中的作用。研究结果表明，2010 年以来渗入阿鲁冰川的液态水增加是导致不稳定事件的关键因素。
1. Introduction 1.导言
The mountain cryosphere is susceptible to the impacts of climate change, where elevated temperatures expedite dissolution of glaciers and thawing of permafrost, ultimately culminating in slope failures (Rechberger and Zangerl, 2022; Rechberger et al., 2021; Guo et al., 2023; Peng et al., 2022; Sun et al., 2024). However, an increasing number of engineering projects are being undertaken in the mountain cryosphere, such as hydropower development activities in the Himalayan region (Wang et al., 2022). The abrupt instability of slopes composed of glaciers, permafrost, and moraine poses significant threats to downstream residents, infrastructure and engineering projects (Huggel et al., 2005; Jiang et al., 2021; Kääb et al., 2021; Shugar et al., 2021), as well as displaces ice masses to lower and warmer elevations, thus indirectly accelerating glacial ablation, as illustrated by an ice-rich slope instability event in Fig. 1. 山区冰冻圈容易受到气候变化的影响，温度升高会加速冰川溶解和永久冻土融化，最终导致斜坡崩塌（Rechberger 和 Zangerl，2022 年；Rechberger 等人，2021 年；Guo 等人，2023 年；Peng 等人，2022 年；Sun 等人，2024 年）。然而，越来越多的工程项目正在山区冰冻圈进行，如喜马拉雅地区的水电开发活动（Wang 等人，2022 年）。由冰川、永久冻土和冰碛组成的斜坡突然失稳，会对下游居民、基础设施和工程项目造成严重威胁（Huggel 等人，2005 年；Jiang 等人，2021 年；Kääb 等人，2021 年；Shugar 等人，2021 年），同时也会将冰块移至更低和更温暖的海拔高度，从而间接加速冰川消融，图 1 中的富冰斜坡失稳事件就说明了这一点。
A series of ice-rich landslide disasters have prompted researchers to delve into the study of ice-rich slope instability (Chen et al., 2020). In 一系列富冰滑坡灾害促使研究人员深入研究富冰斜坡的不稳定性（Chen 等人，2020 年）。在
April 2000, a massive long-distance landslide occurred in Yigong, Tibet, China. The landslide formed a dam with a volume of approximately 300 million cubic meters, blocking the Yarlung Tsangpo River's tributary, Yigong Tsangpo. Subsequently, the dam broke, and the resulting flood affected tens of thousands of people downstream (Guo et al., 2020; Zhao et al., 2023). On February 7, 2021, a high-altitude long-distance rock and ice avalanche chain disaster occurred at Chamoli, Indian Himalaya. The resulting high-altitude rock and ice avalanche, originating at an elevation of , triggered the entrainment of saturated loose mixtures of glacial debris in the valley. This led to the formation of a longdistance-moving mudflow, devastating two hydropower projects under construction downstream and killing over 200 lives (Shugar et al., 2021; Jiang et al., 2021; Zhang et al., 2023). 2000 年 4 月，中国西藏易贡发生大规模长距离山体滑坡。山体滑坡形成了一个体积约为 3 亿立方米的大坝，堵塞了雅鲁藏布江的支流易贡藏布。随后，大坝决堤，造成洪水泛滥，下游数万人受灾（Guo 等人，2020 年；Zhao 等人，2023 年）。2021 年 2 月 7 日，印度喜马拉雅山脉的 Chamoli 发生了高海拔长距离岩冰雪崩连锁灾害。由此引发的高海拔岩冰雪崩源于海拔 ，引发了山谷中饱和松散混合物冰川碎屑的夹带。这导致形成了长距离移动的泥石流，摧毁了下游 两个在建的水电项目，造成 200 多人死亡（Shugar 等人，2021 年；Jiang 等人，2021 年；Zhang 等人，2023 年）。
Generally, the instability of ice-rich slopes occurs on bedrock steeper than (Kääb et al., 2018). Somewhat surprisingly, substantial parts of some low-angle ice-rich slopes have collapsed in recent years, and this kind of ice-rich landslides appears to be more frequent than commonly thought. In 2002, the Kolka glacier collapsed (Evans et al., 2009); of ice and rock detached from the bed with an average slope of only , claiming 135 lives. This event was formerly considered an extraordinary case of low-angle glacier detachments, as the glacier was located on a dormant volcano (Haeberli et al., 2004; Huggel et al., 2005). Nevertheless, two adjacent low-angle mountain glacier detachments (Fig. 2) occurred afterward in the Aru range, Tibet, releasing approximately of mass, resulting in the deaths of nine herders and hundreds of livestock (Gilbert et al., 2018; Kääb et al., 2018; Tian et al., 2017). Subsequently, several detachments similar to the Aru glacier collapses have also been documented, albeit at smaller scales (Falaschi et al., 2019; Jacquemart et al., 2020). 一般来说，富冰斜坡的不稳定性发生在陡于 的基岩上（Kääb 等人，2018 年）。令人略感意外的是，近年来一些低角度富冰斜坡的很大一部分已经坍塌，这种富冰滑坡似乎比通常想象的更为频繁。2002 年，科尔卡冰川崩塌（Evans 等人，2009 年）； 冰和岩石从平均坡度仅为 的冰床上脱离，造成 135 人丧生。由于冰川位于休眠火山上，这一事件以前被认为是低角度冰川脱离的特殊案例（Haeberli 等人，2004 年；Huggel 等人，2005 年）。然而，西藏阿鲁山脉随后发生了两次相邻的低角度山地冰川脱离（图 2），释放了约 的质量，造成 9 名牧民和数百头牲畜死亡（Gilbert 等人，2018 年；Kääb 等人，2018 年；Tian 等人，2017 年）。随后，与阿鲁冰川塌方类似的几次脱离也被记录在案，尽管规模较小（Falaschi等人，2019年；Jacquemart等人，2020年）。
These ice-rich landslides raise the following scientific questions: (I) why did seemly stable low-angle ice-rich slopes detach catastrophically, (II) how to assess ice-rich slopes' stability quantitatively, and (III) how to identify the easiest-to-detach part of an ice-rich slope for hazard mitigation purposes? Many efforts have been made to answer Question I and have provided much instructive understanding (Faillettaz et al., 2015; Gilbert et al., 2018; Huggel, 2009; Jacquemart et al., 2020; Kääb et al., 2021). Previous studies reported that the instability of the Aru glaciers was caused by various factors, such as water infiltration due to thawing and rain, glacier surface steepening, polythermal glacier regime, and topography (Kääb et al., 2018). Gilbert et al. (2018) presented a detailed analysis based on observation and thermo-mechanical modelling to back-analyse the stress, friction and temperature changes that led to the Aru glacier collapses, but they did not attempt to show explicitly the change in the stability (such as the factor of safety) of the glaciers before their collapses. Although these investigations are constructive for understanding the details of the Aru glacier collapses, relatively easy-toimplement approaches that are practical for engineers to project icerich slope instability and potential detachment zones, are still not available. This is likely because the instability of ice-rich slopes is a multidisciplinary problem involving complex multi-physical processes, such as ice crevassing, subglacial water flows, and changes in the thermal structure of the ice-rich slopes. Therefore, further research is needed to evaluate ice-rich slope instability, particularly for low-angle slopes. 这些富冰滑坡提出了以下科学问题：(I) 为什么看似稳定的低角度富冰斜坡会发生灾难性脱离；(II) 如何定量评估富冰斜坡的稳定性；(III) 如何识别富冰斜坡最易脱离的部分以达到减灾目的？为回答问题一，人们做出了许多努力，并提供了许多具有启发性的认识（Faillettaz 等人，2015 年；Gilbert 等人，2018 年；Huggel，2009 年；Jacquemart 等人，2020 年；Kääb 等人，2021 年）。之前的研究报告称，阿鲁冰川的不稳定性是由多种因素造成的，如解冻和降雨造成的水渗透、冰川表面陡峭化、多热冰川机制和地形（Kääb 等人，2018 年）。Gilbert 等人（2018 年）基于观测和热机械模型进行了详细分析，反向分析了导致阿鲁冰川崩塌的应力、摩擦力和温度变化，但他们并未试图明确显示冰川崩塌前的稳定性变化（如安全系数）。尽管这些研究对了解阿鲁冰川崩塌的细节具有建设性意义，但对于工程师来说，仍然没有相对容易实施的方法来预测冰川斜坡的不稳定性和潜在的脱离区。这可能是因为富冰斜坡的不稳定性是一个多学科问题，涉及复杂的多物理过程，如冰裂缝、冰川下水流和富冰斜坡热结构的变化。因此，需要进一步开展研究，以评估富冰斜坡的不稳定性，尤其是低角度斜坡的不稳定性。
Taking the Aru glacier collapses as reference cases, the present study aims to develop methodologies to assess ice-rich slope stability, considering the multi-physical processes (e.g., thermo-hydromechanical) involved in an ice-rich slope system. Unlike investigations dedicated to revealing the details of the Aru glacier collapses, such as ice flow velocity, the main purpose of this article is to present a model for estimating the factor of safety (FOS) and potential collapse zone of an ice-rich slope that lacks monitoring information (e.g., flow rates) and whose stability is mainly affected by rainfall, meltwater, and temperature change instead of other factors like seismicity. Therefore, some features of the Aru glacier collapses may be overlooked in the present study, such as stress transfer from predominant basal drag toward predominant lateral shearing in the detachment areas (Gilbert et al., 2018), as they are not the main focus here. However, because meteorological information is collected to drive the coupled thermo-hydro-mechanical modelling examples (see Section 4) that are based on the Aru glaciers' geometries and the proposed model, this article also provides some insights into the Aru glacier collapses through processing and analysing a state-of-the-art global climate reanalysis dataset, i.e., the ERA5-Land hourly dataset (Muñoz-Sabater et al., 2021). The air temperature and precipitation characteristics in the Aru Range are analysed, aiming at further examining how and when climate change began to affect the stability of the glaciers significantly. Our analysis supports quantitative stability analysis and disaster risk assessments of glacial, permafrost, and ice-rich moraine slopes experiencing impacts from climate change. 本研究以阿鲁冰川崩塌为参考案例，旨在开发评估富冰斜坡稳定性的方法，同时考虑富冰斜坡系统中涉及的多种物理过程（如热-水-机械）。与专门揭示阿鲁冰川崩塌细节（如冰流速度）的研究不同，本文的主要目的是提出一个模型，用于估算富冰斜坡的安全系数（FOS）和潜在崩塌区，因为这种斜坡缺乏监测信息（如流速），其稳定性主要受降雨、融水和温度变化而非地震等其他因素的影响。因此，本研究可能会忽略阿鲁冰川崩塌的一些特征，如脱离区域的应力从主要的基底阻力向主要的横向剪切力转移（Gilbert 等人，2018 年），因为它们不是本研究的重点。然而，由于气象信息的收集是为了驱动基于阿鲁冰川几何形状和拟议模型的热-水-力学耦合建模实例（见第 4 节），本文还通过处理和分析最先进的全球气候再分析数据集，即ERA5-Land 小时数据集（Muñoz-Sabater 等人，2021 年），对阿鲁冰川崩塌提供了一些见解。我们分析了阿鲁山脉的气温和降水特征，旨在进一步研究气候变化如何以及何时开始对冰川的稳定性产生重大影响。我们的分析有助于对受到气候变化影响的冰川、永久冻土和富冰碛斜坡进行定量稳定性分析和灾害风险评估。
2. The Aru glacier collapses 2.阿鲁冰川崩塌
At 11:15 a.m. Beijing Time on 17 July 2016, the lower part (5800-5190 m a.s.l.) of a glacier (termed Aru-1) in the Aru Range of the Tibetan Plateau ( ) underwent a sudden collapse, despite a low surface slope angle of . This glacier detachment spanned a length of over , had a maximum depth of over , and encompassed a volume of about (Kääb et al., 2018). Two months later at 5:00 a.m. - 11:20 a.m. on 21 September 2016, a similar event occurred on Aru-2 Glacier (Fig. 2), with a surface slope angle of , resulting in the detachment of the glacier between 5800 and a.s.l. The released glacier ice formed two flow fronts, totaling (Kääb et al., 2018). 北京时间2016年7月17日上午11时15分，青藏高原阿鲁山脉（ ）一条冰川（称为阿鲁-1）的下部（海拔5800-5190米）在表面坡角 较低的情况下发生突然崩塌。这次冰川脱离的长度超过 ，最大深度超过 ，体积约为 （Kääb 等人，2018 年）。两个月后的2016年9月21日凌晨5:00-11:20，阿鲁-2冰川（图2）也发生了类似事件，表面坡角为 ，导致冰川在海拔5800- 之间脱离，释放的冰川冰形成了两个流动锋面，总计 （Kääb等人，2018）。
Early studies have ruled out abnormally high geothermal fluxes or earthquake activities as triggers for the collapses, inferring summer melt and rain as important factors (Gilbert et al., 2018; Kä̈b et al., 2018). Indeed, satellite images captured many transverse crevasses (at ) above the detachment zones prior to the collapses (Fig. 3). Moreover, the glacier surfaces above the crevasses remained relatively intact and were steeper than those containing crevasses (Fig. 4). These landform features enabled supraglacial meltwater and rainwater to penetrate the crevasses and reach the bases of the detachment zones, thereby reducing basal traction. Indeed, a polythermal structure, characterized by a temperate (thawed) basal part and a relatively cold overlying part, was inferred for the Aru glaciers by Gilbert et al. (2018). Additionally, a mudflow-like sediment fan was observed on the surface of Aru-1 two days before its collapse, indicating the presence of substantial subglacial liquid water at that time. Gilbert et al. (2018) also concluded that the Aru glaciers were under prevailing temperate basal conditions over the detachment areas, and changes in pore pressure led to the weakening of the underlying till and sediment in the temperate area. Ultimately, the 早期研究排除了异常高的地热通量或地震活动作为塌陷的触发因素，推断夏季融水和降雨是重要因素（Gilbert 等人，2018 年；Kä̈b 等人，2018 年）。事实上，卫星图像捕捉到了塌陷前脱离带上方的许多横向裂缝（ 处）（图 3）。此外，裂缝上方的冰川表面仍然相对完整，而且比包含裂缝的冰川表面更加陡峭（图 4）。这些地貌特征使得超冰川融水和雨水能够穿透裂缝，到达脱离带的底部，从而减少了基底牵引力。事实上，Gilbert 等人（2018 年）推断阿鲁冰川具有多热结构，其特点是基底部分为温带（解冻），上覆部分相对寒冷。此外，在阿鲁-1 崩塌前两天，在其表面观察到类似泥流的沉积扇，表明当时存在大量冰川下液态水。Gilbert 等人（2018 年）还得出结论，阿鲁冰川在脱离区域处于盛行的温带基底条件下，孔隙压力的变化导致温带地区下层畋猎和沉积物的减弱。最终
Fig. 1. The collapse of an ice-rich slope in Bomi, Tibet, which has damaged the road and resulted in a substantial amount of meltwater. 图 1.西藏波密一处富冰斜坡的坍塌，损坏了道路并导致大量融水。
Fig. 2. Locations and profiles of the Aru glaciers and detachment zones: (a) Location of the Aru Range, (b) Satellite image of the Aru glaciers captured after the collapses, and a close-up of post-collapse geomorphology of Aru-2, (c, d) Longitudinal cross sections along kinematic centerlines of Aru-1 and Aru-2, respectively. The solid red loops delineate outlines of the detachment zones. The dashed blue lines represent the kinematic centerlines of the glaciers, determined by Gilbert et al. (2018). The yellow arrows indicate the glaciers' movement directions. Alluvial fans I and II were formed by two main flows in the Aru-2 glacier collapse. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 图 2.阿鲁冰川和脱离带的位置和剖面图：（a）阿鲁山脉的位置；（b）阿鲁冰川崩塌后的卫星图像，以及阿鲁-2 崩塌后的地貌特写；（c、d）分别沿阿鲁-1 和阿鲁-2 运动学中心线的纵向剖面图。红色实心环线勾勒出脱离带的轮廓。蓝色虚线代表由 Gilbert 等人（2018 年）确定的冰川运动中心线。黄色箭头表示冰川的运动方向。冲积扇 I 和冲积扇 II 由阿鲁-2 冰川塌陷中的两股主要水流形成。(有关本图例中颜色参考文献的解释，请读者参阅本文的网络版）。
Apr 2009 2009 年 4 月
Dec 2011 2011 年 12 月
Nov 2015 2015 年 11 月
Fig. 3. The evolution of crevasses on the surface of Aru-2 Glacier (the last two images are reproduced after Kääb et al. (2018)). Few crevasses were visible on Aru-2 Glacier prior to April 2009. However, crevasses proliferated in the following years, indicating severe damage accumulation on the glacier during this period. 图 3.阿鲁-2 号冰川表面裂缝的演变（最后两张图片根据 Kääb 等人（2018 年）复制）。2009 年 4 月之前，阿鲁-2 冰川上几乎看不到裂缝。然而，在接下来的几年中，裂缝大量增加，表明在此期间冰川上积累了严重的破坏。
Fig. 4. Schematical diagram of Aru-1 Glacier and the corresponding finite element model. The glacier and the permeable layer are modelled as elastic-plastic deformable bodies with the outward displacements at their upper left ends constrained, whereas the bedrock is treated as a fixed rigid body. The triangles around the bedrock portion symbolize that the displacement of the bedrock is completely constrained. In the modelling, meltwater and rainwater are injected into the permeable soil-ice mixture layer from its left end. 图 4.阿鲁-1 冰川示意图及相应的有限元模型。冰川和渗透层被模拟为弹塑性可变形体，其左上角的向外位移受到约束，而基岩则被视为固定刚体。基岩部分周围的三角形表示基岩的位移完全受限。在模型中，融水和雨水从左端注入透水性土冰混合物层。
Fig. 5. Schematic diagram of a conceptual model for ice-rich slopes akin to the Aru glaciers. Water infiltration and thermodynamics affect the glacier's stability, with liquid water sourced from the catchment area upstream of the detachment zone. Glacier collapse is conjectured to be primarily initiated by the weakening of the basal permeable layer and the overlying glacier, with their strength properties affected by temperature and fluid pressure. 图 5.类似阿鲁冰川的富冰斜坡概念模型示意图。水的渗透和热力学影响冰川的稳定性，液态水来自脱离区上游的集水区。据推测，冰川崩塌主要是由基底渗透层和上覆冰川的减弱引起的，它们的强度特性受温度和流体压力的影响。
plastic rheology of the underlying till combined with low bedrock roughness induced the instability of the Aru glaciers. 底层冰碛的塑性流变与基岩的低粗糙度导致了阿鲁冰川的不稳定。
Referencing these characteristics, we establish a conceptual model based on coupled thermo-hydro-mechanical modelling in the next section to quantify the stability evolution of ice-rich slopes under the action of meltwater, rainwater and temperature change, so as to provide methodologies for assessing the stability of ice-rich slopes exhibiting the aforementioned features. 根据这些特点，我们在下一节建立了一个基于热-水-机械耦合模型的概念模型，以量化富冰斜坡在融水、雨水和温度变化作用下的稳定性演变，从而为评估具有上述特点的富冰斜坡的稳定性提供方法。
3. Conceptual model and formulations 3.概念模型和表述
Fig. 5 depicts the conceptual model, wherein the subglacial water flow is assumed to occur within an ice-soil mixture layer between the glacier and the bedrock, both deemed impervious for simplification (Hewitt, 2011). The permeable layer is used to approximate subglacial till, such as that underlying the Aru glaciers. In the model, the subglacial water flow is modelled as a variably saturated flow described by Richards' equation, and the hydraulic properties of the permeable layer are described using the van Genuchten equation (van Genuchten, 1980). Glacial till has been widely observed to behave with plastic rheology, with a shear strength strongly dependent on the effective normal stress (Iverson et al., 1998; Clarke, 2005; Iverson, 2010). Thus, an excessive water input in the permeable layer will increase the fluid pressure, thereby reducing the shear strength of the permeable layer and inducing potential instability of the ice-rich slope. Furthermore, the subglacial water flow can warm up the permeable layer. Indeed, temperate and polythermal glaciers are common at high latitudes and high altitudes (e. g., Rocky Mountains, Scandinavian Mountains, Alps, and Himalaya) (Aschwanden and Blatter, 2009); particularly in the context of a warming climate, temperatures of ice-rich slopes are expected to elevate further (Haeberli and Beniston, 1998; Gruber and Haeberli, 2007). Because a considerable amount of liquid water flows into the Aru glaciers through crevasses upstream of the collapsed zones, and under the influence of substantial subglacial water flow, the base of the ice is considered to be at the pressure melting point (Fig. S1 in Supplementary Material). Moreover, given the fact that these low-angle glaciers ultimately experience collapses and that these collapses are not triggered by seismic factors, the parts of the glaciers in contact with the bedrock are speculated to be similar to temperate glaciers, which are commonly assumed to be at the pressure melting point (Lliboutry, 1971). 图 5 描述了这一概念模型，其中假定冰川下水流发生在冰川与基岩之间的冰土混合物层中，为简化起见，两者都被视为不透水层（Hewitt，2011 年）。透水层被用来近似冰川下的冰碛物，例如阿鲁冰川下的冰碛物。在该模型中，冰川下水流被模拟为理查兹方程描述的可变饱和水流，而渗透层的水力特性则使用 van Genuchten 方程描述（van Genuchten，1980 年）。据广泛观察，冰川沉积物具有塑性流变特性，其剪切强度与有效法向应力密切相关（Iverson 等人，1998 年；Clarke，2005 年；Iverson，2010 年）。因此，透水层中过多的水输入会增加流体压力，从而降低透水层的剪切强度，导致富冰斜坡的潜在不稳定性。此外，冰川下水流还能使渗透层升温。事实上，温带和多热源冰川在高纬度和高海拔地区（如落基山脉、斯堪的纳维亚山脉、阿尔卑斯山和喜马拉雅山）非常常见（Aschwanden 和 Blatter，2009 年）；尤其是在气候变暖的背景下，富冰斜坡的温度预计会进一步升高（Haeberli 和 Beniston，1998 年；Gruber 和 Haeberli，2007 年）。由于大量液态水通过塌陷区上游的裂缝流入阿鲁冰川，并受到大量冰川下水流的影响，冰层底部被认为处于压力熔点（补充材料图 S1）。 此外，鉴于这些低角度冰川最终会发生崩塌，而且这些崩塌并非由地震因素引发，因此推测冰川与基岩接触的部分类似于温带冰川，而温带冰川通常被认为处于压力熔点（Lliboutry，1971 年）。
Besides subglacial water, the atmosphere also affects the thermal structure of an ice-rich slope. Heat transfer between the ice-rich slope and air occurs through convection and radiation, which are mainly influenced by wind speed and slope surface albedo, receptively, in addition to the temperature difference between the slope and air. Heat transmission within the ice-rich slope is mainly through conduction, contingent on the slope's temperature gradient and thermal conductivity. 除了冰川下水，大气也会影响富冰斜坡的热结构。富冰斜坡与空气之间的热量传递是通过对流和辐射进行的，这主要受风速和斜坡表面反照率的影响，也受斜坡与空气之间温差的影响。富冰斜坡内部的热量传输主要通过传导，取决于斜坡的温度梯度和导热性。
The hydraulic head governs the variation in pore water pressure within the permeable layer, influencing the change in effective stress and shear strength of the permeable layer. The vertical thickness of a glacier typically varies along the slope. Consequently, the location of water influx determines the extent to which shear strength may decrease due to rainfall and meltwater. However, as mentioned earlier, the topography of the Aru glaciers facilitates water infiltration from the upper left end of the collapse zone. Moreover, in subglacial water flow simulations, ice is often treated as impermeable (Hewitt, 2011). Therefore, the current conceptual model only considers infiltration from the upper left end of the collapse zone and water flow along the permeable layer. Under the combined effect of increasing temperature and subglacial water pressure, the strength of ice-soil mixtures is degraded. Thus, the slope stability can be analysed by determining the temperature, subglacial water pressure and stress state of the slope based on the link between the ice-soil mixture strength, the thermal structure and the water pressure. Herein, a strength reduction method, widely used in slope stability analysis (Hua et al., 2022; Wu et al., 2018), is utilized to assess the stability of ice-rich slopes, and this conceptual model is formulated through coupled thermo-hydro-mechanical modelling to evaluate the pre-collapse stability of the Aru glaciers, which is affected by water infiltration and temperature rise. 水头控制着渗透层内孔隙水压力的变化，影响着渗透层有效应力和剪切强度的变化。冰川的垂直厚度通常沿斜坡变化。因此，水流入的位置决定了剪切强度因降雨和融水而降低的程度。不过，如前所述，阿鲁冰川的地形有利于水从塌方区的左上方渗入。此外，在冰川下水流模拟中，冰通常被视为不可渗透的（Hewitt，2011 年）。因此，当前的概念模型只考虑了来自塌陷区左上方的渗透和沿渗透层的水流。在温度升高和冰川下水压力的共同作用下，冰土混合物的强度会降低。因此，可以根据冰土混合物强度、热结构和水压力之间的联系，通过确定温度、冰川下水压力和斜坡应力状态来分析斜坡稳定性。在此，利用广泛应用于边坡稳定性分析的强度还原法（Hua 等人，2022 年；Wu 等人，2018 年）来评估富冰边坡的稳定性，并通过热-水-机械耦合建模来制定这一概念模型，以评估阿鲁冰川受水渗透和温度上升影响的崩塌前稳定性。
3.1. Thermodynamics modelling 3.1.热力学建模
In the conceptual model, an ice-rich slope is treated as a nonNewtonian viscous fluid. Based on this treatment, the heat transfer process between the ice-rich slope and the ambient is simulated using equations established for heat transfer in fluid interfaces, which have been widely accepted for modelling thermal structures of glaciers. The energy balance equation in the thermodynamics modelling is as follows (Tannehill et al., 1997): 在概念模型中，富冰斜坡被视为非牛顿粘性流体。在此基础上，利用流体界面热传递方程模拟富冰斜坡与环境之间的热传递过程，该方程已被广泛接受用于冰川热结构建模。热力学模型中的能量平衡方程如下（Tannehill 等人，1997 年）：
where is the ice-rich slope's temperature; is the density; and are the heat fluxes by convection and radiation, respectively; and are the specific heat capacity and thermal conductivity, respectively. For the Aru glacier modelling examples, and are calculated using the following equations (Ritz, 1987): 其中， 是富冰斜坡的温度； 是密度； 和 分别是对流和辐射的热通量； 和 分别是比热容和热导率。在阿鲁冰川建模实例中， 和 是通过以下公式计算得出的（Ritz，1987 年）：
where , and are parameters used for the principle of dimensional homogeneity. 其中 和 是用于维度均匀性原则的参数。
The convective heat flux between the air and the ice-rich slope's surface is modelled using Newton's law of cooling (Eq. (4)). Regarding the Aru glacier modelling examples, Eq. (5) outlines a method for determining the heat transfer coefficient ( ) (Incropera et al., 2006). 空气与富冰斜坡表面之间的对流热通量采用牛顿冷却定律（公式 (4)）建模。关于阿鲁冰川建模实例，公式 (5) 概述了确定传热系数 ( ) 的方法（Incropera 等人，2006 年）。
where is the outward unit normal vector; denotes the air temperature; is a characteristic length; is the Reynolds number; is the Prandtl number; represents the wind speed on the the slope surface. 其中， 为向外单位法向量； 表示空气温度； 为特征长度； 为雷诺数； 为普朗特数； 表示斜面上的风速。
The radiation from the slope surface to the air is considered using the Stefan-Boltzmann law 从斜面到空气的辐射是通过斯蒂芬-玻尔兹曼定律来考虑的
where is the Stefan-Boltzmann constant, and is the surface albedo. 其中， 是斯蒂芬-玻尔兹曼常数， 是表面反照率。
To determine the pressure melting point (PMP) in each part of the Aru glaciers, the glacier is simulated as a creeping flow using the Stokes equations. The governing equations of mass and momentum conservation are as follows (Tannehill et al., 1997): 为了确定阿鲁冰川各部分的压力熔点（PMP），使用斯托克斯方程将冰川模拟为蠕动流。质量和动量守恒的支配方程如下（Tannehill 等人，1997 年）：
where is the velocity vector, is the pressure, is the density, is the gravitational acceleration, is the identity matrix, and is the dynamic viscosity, which can be described using Glen's flow law (Glen, 1955) as 其中， 为速度矢量， 为压力， 为密度， 为重力加速度， 为特征矩阵， 为动态粘度，可使用格伦流动定律（格伦，1955 年）描述为
follows: 如下所示：
where is the flow law exponent, and is the shear rate, defined as the norm of the strain rate tensor: 其中， 是流动规律指数， 是剪切速率，定义为应变速率张量的法线：
The flow rate factor is usually simulated using the Arrhenius law: 流速系数 通常使用阿伦尼乌斯定律进行模拟：
where is a flow rate constant, is the activation energy, and is the universal gas constant. 其中 是流速常数， 是活化能， 是通用气体常数。
By referring to the User's Guide of Comsol Multiphysics (2022), and are determined using the following equations in the Aru glacier modelling examples: 参考《Comsol Multiphysics 用户指南》（2022 年），在阿鲁冰川建模实例中， 和 是通过以下公式确定的：
The relationship between the PMP and pressure is given by (Wagner et al., 1994) PMP 与压力之间的关系如下（Wagner 等人，1994 年）
where and are the standard atmospheric pressure and the corresponding melting point temperature of water, respectively; is the Clausius-Clapeyron constant. 其中 和 分别是标准大气压和水的相应熔点温度； 是克劳修斯-克拉皮隆常数。
3.2. Subglacial water flow simulation 3.2.冰川下水流模拟
As mentioned above, the subglacial water flow is modelled using Richards' equation (Richards, 1931) 如上所述，冰川下水流是利用理查兹方程（理查兹，1931 年）模拟的
where is the water pressure, is the fluid velocity vector, denotes the effective saturation, represents the specific moisture capacity, is the storage coefficient, is the density of water, is the fluid source (positive) or sink (negative), is the gravitational acceleration, denotes the saturated hydraulic conductivity, and is the relative permeability. 其中 为水压， 为流体速度矢量， 表示有效饱和度， 表示比容湿度， 为储水系数、 是水的密度， 是流体源（正）或汇（负）， 是重力加速度， 表示饱和导水性， 是相对渗透性。
Hydraulic characteristics of the permeable layer are described using the van Genuchten equation (van Genuchten, 1980): 透水层的水力特征采用 van Genuchten 方程（van Genuchten，1980 年）进行描述：
where is the pressure head; represents the relative permeability; denotes the specific moisture capacity; and are the saturated and residual liquid volume fractions, respectively; , and depend on the porous media. 其中， 为压头； 表示相对渗透率； 表示比容湿度； 和 分别为饱和液体体积分数和残余液体体积分数； 和 取决于多孔介质。
In the present model, the volume of meltwater and rainwater is calculated from the hourly precipitation (rainfall plus snowfall) in the ERA5-Land reanalysis dataset (Muñoz-Sabater et al., 2021). 在本模型中，融水和雨水量是根据ERA5-Land 再分析数据集（Muñoz-Sabater 等人，2021 年）中的每小时降水量（降雨量加降雪量）计算得出的。
The meltwater volume is determined through a degree-day factor method (Kuchment and Gelfan, 1996): 融水体积是通过度日系数法确定的（Kuchment 和 Gelfan，1996 年）：
where is the daily meltwater in equivalent water depth (m); is the degree-day factor day represents the total number of grids of the target catchment; is the average air temperature above at the th grid and a function of elevation: 其中， 为等效水深的日融水（米）； 为度日因子 日 表示目标流域的网格总数； 为第 个网格上 的平均气温，是海拔高度的函数：
where is the daily average air temperature at a reference elevation is the elevation of the th grid; is the temperature lapse rate. 其中， 是参考海拔高度处的日平均气温； 是第 个网格的海拔高度； 是温度失效率。
The method for obtaining the amount of rainwater is as follows (Kääb et al., 2018): 获取雨水量的方法如下（Kääb 等人，2018 年）：
where is the daily rainfall in equivalent water depth (mm); is the daily total precipitation including rainfall and snowfall; is described as 其中， 为等效水深的日降水量（毫米）； 为包括降雨和降雪在内的日总降水量； 描述为
