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 其中， 为等效水深的日降水量（毫米）； 为包括降雨和降雪在内的日总降水量； 描述为
where is a threshold temperature above which rain instead of snow occurs. 其中， 是一个临界温度，超过这个温度就会出现雨而不是雪。
3.3. Simulation of slope stability 3.3.模拟斜坡稳定性
In evaluating the stability of ice-rich slopes, which is deemed a mechanical problem involving matrix failure, ice-soil mixtures are treated as solids. The force balance is expressed as follows: 富冰斜坡的稳定性被认为是一个涉及基质破坏的力学问题，在评估富冰斜坡的稳定性时，冰土混合物被视为固体。力平衡表示如下
where denotes stress in solid substrate; is external stress, which corresponds to fluid pressure in the present study; represents gravity. In the Aru glacier modelling examples, the glacier is assumed to be linear elastic until local stress reaches the Mohr-Coulomb failure criterion. By referencing Fu et al. (2021), the cohesion and friction angle of the glacier are set to be negatively related to glacier temperature in the present study 其中， 表示固体基质中的应力； 表示外部应力，在本研究中相当于流体压力； 表示重力。在阿鲁冰川建模实例中，假定冰川为线性弹性体，直到局部应力达到莫尔-库仑破坏准则。参考 Fu 等人（2021 年）的研究，本研究将冰川的内聚力 和摩擦角 设置为与冰川温度负相关。
where , and . 其中， 和 。
Subsequently, the stability is evaluated through the strength reduction method (Terzaghi and Peck, 1967): 随后，通过强度降低法对稳定性进行评估（Terzaghi 和 Peck，1967 年）：
where is a strength reduction coefficient. 其中 是强度降低系数。
The threshold value of , above which the elastic-plastic finite element computation fails to converge, is defined as the factor of safety (FOS) of the glacier. Note that this operation is based on the following premises: (i) model settings (e.g., initial conditions, element size and shape) are reasonably selected to ensure that the numerical calculation can proceed smoothly; (ii) the loss of convergence is caused by excessive element distortion due to glacier failure rather than by artificial numerical effects. 的临界值被定义为冰川的安全系数（FOS），当超过该临界值时，弹塑性有限元计算将无法收敛。请注意，这一操作基于以下前提：(i) 合理选择模型设置（如初始条件、元素大小和形状），以确保数值计算能够顺利进行；(ii) 收敛失败是由冰川破坏导致的元素过度变形引起的，而不是人为的数值效应。
4. Modelling examples and results 4.建模实例和结果
4.1. Model settings 4.1.模型设置
As mentioned above, the geometry and climatic conditions of the Aru glaciers are referenced to demonstrate the specific implementation process of the proposed model. The glaciers are simulated in two dimensions under plane strain conditions, as depicted in Fig. 4. The model geometry is based on longitudinal cross-sections at the kinematic center lines of the glaciers, which are referenced from Gilbert et al. (2018). The glacier domain is discretized using an unstructured triangular mesh, while the permeable layer and bedrock domains are represented with structured mapped quadrilateral grids. The elements have dimensions ranging from to , except for the two layers of elements in the permeable layer, which is assumed to be thick and a mixture of ice (density: ) and sediment (density: (Singh and Ramanathan, 2018)) in a volume ratio of 1:1. In the present conceptual model, the permeable layer between glacier and bedrock is considered thin, and it is assumed that the ice and sediment each contributed half of the solid material to the permeable layer. Since the permeable layer is assumed to possess a porosity of (see Table 1), the saturated density of the permeable layer is calculated as . This value falls within the reported range of to for an ice-rich permafrost layer (Hausmann et al., 2007). It is noteworthy that the density of the permeable layer has a negligible impact on the overall glacier stability due to its limited thickness of only , which is minor compared to the thickness of the ice (with a maximum depth exceeding 如上所述，我们参考了阿鲁冰川的几何形状和气候条件，以展示拟议模型的具体实施过程。如图 4 所示，冰川是在平面应变条件下进行二维模拟的。模型几何基于冰川运动中心线的纵向截面，参考了 Gilbert 等人（2018 年）的研究。冰川域采用非结构化三角形网格离散，而渗透层和基岩域则采用结构化映射四边形网格。元素的尺寸从 到 不等，但渗透层中的两层元素除外，假设渗透层的厚度为 ，由冰（密度： ）和沉积物（密度： ）混合而成： （Singh 和 Ramanathan，2018 年））的体积比为 1:1。在本概念模型中，冰川与基岩之间的透水层被认为很薄，假定冰和沉积物各占透水层固体物质的一半。由于假设渗透层的孔隙率为 （见表 1），因此计算出渗透层的饱和密度为 。该值在报告的富冰永久冻土层 至 的范围内（Hausmann 等人，2007 年）。值得注意的是，由于渗透层的厚度有限，只有 ，与冰层厚度相比微不足道（最大深度超过 ），因此渗透层的密度对冰川整体稳定性的影响可以忽略不计。
Table 1 表 1
Parameters in the thermo-hydro-mechanical modelling. 热-水-机械模型中的参数。
Parameter 参数
Description 说明
Value 价值
Density of ice 冰的密度
Density of subglacial sediment 冰川下沉积物密度
（辛格和拉马纳坦，2018 年）
(Singh and
Ramanathan, 2018 )
Clausius-Clapeyron constant 克劳修斯-克拉皮戎常数
吉尔伯特等人，2014年
Gilbert
et al., 2014
Saturated liquid volume fraction 饱和液体体积分数
0.2（Vlahou 和 Worster，2015 年）
0.2 (Vlahou and Worster,
2015)
Residual liquid volume fraction 残留液体体积分数
0.05 (Li et al., 2022) 0.05（Li 等人，2022 年）
Effective hydraulic conductivity 有效水导率
Stress exponent 应力指数
3 (Gilbert et al., 2014) 3（吉尔伯特等人，2014 年）
van Genuchten model parameter 范-杰努赫腾模型参数
(Kool et al., 1985) （Kool 等人，1985 年）
van Genuchten model parameter 范-杰努赫腾模型参数
0.5 (Kool et al., 1985) 0.5（Kool 等人，1985 年）
van Genuchten model parameter 范-杰努赫腾模型参数
0.5（Jiménez-Martínez 等人，2009 年）
0.5 (Jiménez-Martínez et al.,
2009 )
热力学建模中的特征长度
Characteristic length in
thermodynamics modelling
Degree-day factor 度日系数
0.006 (m w.eq. day) )IPCC, 2019)
0.006 (m w.eq. day )
IPCC, 2019)
Temperature lapse rate 温度失效率
(Kääb et al., 2018) （Kääb 等人，2018 年）
冰川表面的融点温度
Melting-point temperature of the
glacier surface
273.15 (K) (Kääb et al., 2018)
Snow/rain threshold 雨雪阈值
275.15 (K) (Kääb et al., 2018)
). Sensitivity analyses, considering variations in permeable layer thickness from to (see details in Supplementary Material), indicate that the thickness of the permeable layer does not significantly affect the simulated pre-collapse FOS and collapse zones of the Aru glaciers. Prior to determining the number of elements, several sensitivity analyses are conducted (Please refer to Supplementary Material). Results demonstrate that variations in the number of elements or element size do not induce significant changes in the modelling results if the number of elements exceeds 8000 and the element sizes in the glacier domain are relatively uniform. Ultimately, a total of 26,770 and 28,840 linear Lagrange elements are employed for the Aru-1 and Aru-2 modelling examples, respectively. In the model, the glacier and the permeable layer are treated as elastic-plastic deformable bodies with constrained outward displacements at their upper left ends (Fig. 4). Simultaneously, the bedrock is modelled as a fixed rigid body, indicating that the displacement of the bedrock is entirely constrained. The corresponding mechanical parameter values and other parameter settings of the models are summarized in Table 1. ）。敏感性分析考虑了渗透层厚度从 到 的变化（详见补充材料），结果表明渗透层厚度对模拟的阿鲁冰川塌陷前FOS和塌陷区没有显著影响。在确定元素数量之前，进行了几项敏感性分析（请参阅补充材料）。结果表明，如果元素数量超过 8000 个，且冰川区域的元素大小相对一致，那么元素数量或元素大小的变化不会导致建模结果发生显著变化。最终，阿鲁-1 和阿鲁-2 模型分别采用了 26,770 和 28,840 个线性拉格朗日元素。在模型中，冰川和透水层被视为弹塑性可变形体，其左上角的向外位移受到约束（图 4）。同时，基岩被视为固定刚体，表明基岩的位移完全受限。表 1 总结了模型相应的力学参数值和其他参数设置。
As outlined above, the ERA5-Land reanalysis dataset is processed to acquire the hourly data of air temperature, wind velocity, glacier surface albedo, meltwater and rainfall for the Aru Range during 2000-2016. The processes for determining amounts of rainwater and meltwater from the catchments above detachment zones are elucidated as follows. Initially, the catchment area above the crown of each detachment zone ( a.s.1.) is identified based on a digital elevation model (resolution: 12.5 ) of the Aru Range. The determined catchment areas are delineated in Fig. 6. Subsequently, the rainwater and meltwater are quantified using the methods introduced in Section 3.2. Fig. 7 presents characteristics of the determined annual rainwater and meltwater amounts from 2000 to 2016. The volumes of meltwater and rainwater are of a similar order of magnitude. Consequently, both might be essential for the instability of the Aru glaciers. To further elucidate the characteristics of historical temperature variation in the Aru Range, Fig. 8 illustrates the alteration in the summer mean temperature at a.s.l. during the period of 2000-2016 and offers monthly mean temperature during 2010-2016. A rising trend at a rate of per year is observed, indicating a warming climate in the Aru Range. Moreover, July is the month with the highest monthly mean temperature each year, and the July mean temperature was higher in 2010, 2012 and 2013 than in 2016. Fig. 9 summarizes the monthly mean wind speed and albedo in the Aru Range in 2016; the results for the other years are analogous to those for 2016. The monthly mean wind speed ranges from to , and the wind speed, as well as its variability, is generally higher in the cold than in the warm season. The monthly mean albedo ranges from 0.5 to 0.7 , and the albedo variation throughout the year is similar to the temperature variation. 如上所述，通过处理ERA5-陆地再分析数据集，获得了2000-2016年期间阿鲁山脉的气温、风速、冰川表面反照率、融水和降雨量的小时数据。确定脱离带上方集水区雨水和融水量的过程如下。首先，根据阿鲁山脉的数字高程模型（分辨率：12.5 ）确定每个脱离带顶部（ a.s.1.）以上的集水区。确定的集水区见图 6。随后，采用第 3.2 节介绍的方法对雨水和融水进行量化。图 7 显示了 2000 年至 2016 年确定的年雨水量和融水量的特征。融水和雨水的数量级相近。因此，两者都可能是造成阿鲁冰川不稳定的关键因素。为了进一步阐明阿鲁山脉历史上的气温变化特征，图 8 展示了 2000-2016 年期间 海拔高度夏季平均气温的变化以及 2010-2016 年期间月平均气温的变化。从图中可以看出，气温以每年 的速度呈上升趋势，这表明阿鲁山脉的气候正在变暖。此外，每年 7 月的月平均气温最高，2010 年、2012 年和 2013 年的 7 月平均气温高于 2016 年。图 9 总结了 2016 年阿鲁山脉的月平均风速和反照率；其他年份的结果与 2016 年类似。月平均风速从 到 不等，寒冷季节的风速及其变化一般高于温暖季节。 月平均反照率在 0.5 至 0.7 之间，全年反照率变化与温度变化相似。
To ensure the efficiency and convergence of the finite element computation, quarterly mean data of these climate-related variables are employed in the modelling. The partial differential equations (Eqs. (1)(19) and (24)-(27)) were solved on the COMSOL Multiphysics platform. The main steps of the numerical computation are as follows: initially, geo-stress balance is realized so that the model is subject to gravity while its displacement is nullified. Then, the glacier temperature evolution from 1998 to 2016 is simulated with an initial temperature field determined through a steady-state calculation assuming a glacier surface temperature of in winter 1997. This initial temperature condition is selected through a trial-and-error process until the simulated temperature field at the glacier bottom during the period of interest (2000 to 2016) is nearly steady, and the extent of ice affected by air temperature changes during 2000-2016 was essentially constant. This result means that the initial temperature condition is compatible with the inherent temperature evolution of the glacier under the conditions assumed in the proposed model. 为确保有限元计算的效率和收敛性，在建模中采用了这些气候相关变量的季度平均数据。偏微分方程（式 (1)(19) 和 (24)-(27) ）在 COMSOL Multiphysics 平台上求解。数值计算的主要步骤如下：首先，实现地应力平衡，使模型受重力作用，而其位移为零。然后，假设 1997 年冬季冰川表面温度为 ，通过稳态计算确定初始温度场，模拟 1998 年至 2016 年的冰川温度变化。这一初始温度条件是通过试错过程选择的，直到在相关时期（2000 年至 2016 年）冰川底部的模拟温度场接近稳定，且 2000 年至 2016 年期间受气温变化影响的冰层范围基本不变。这一结果意味着初始温度条件符合拟议模型假设条件下的冰川固有温度演变。
In addition to the thermodynamic calculation, the subglacial fluid pressure is computed by injecting water at the upper left end of the permeable layer, as indicated in Fig. 4 with a red arrow. The water 除热力学计算外，还通过在透水层左上方注入水计算冰川下流体压力，如图 4 中红色箭头所示。水
Fig. 6. The catchments based on which the amounts of meltwater and rainwater are determined. The basin border is delineated based on altitude and slope. 图 6.据以确定融水和雨水量的集水区。流域边界根据海拔和坡度划定。
Fig. 7. The meltwater and rainwater determined from the catchment area above the crown of each detachment zone during 2000-2016. Rainfall and meltwater are of a similar order of magnitude, indicating that both are essential for the stability of glaciers. 图 7.2000-2016 年间从每个脱离带冰冠上方集水区测定的融水和雨水。降雨量和融水量的数量级相近，表明两者对冰川的稳定性都至关重要。
injection process is governed by the volume of liquid water (including rainwater and meltwater) rather than by pressure. In other words, the model assumes that the generated liquid water flows into the permeable layer from its upper left end. Lastly, the strength reduction method is implemented step by step on the glacier in the period of interest. Because Aru-1 Glacier collapsed in summer 2016, the FOS had to be 1.00 at that time. To this end, the saturated hydraulic conductivity of the assumed 1m-thick permeable layer is calibrated at through a trialand-error process. Iken et al. (1996) inferred a higher hydraulic conductivity of for the till beneath Gorner Glacier from boreholeresponse tests. Similarly, Stone and Clarke (1993) reported large values for Trapridge Glacier, with Stone et al. (1997) revising the estimate to 5 in the connected part of the drainage system using inversion theory. In our study, assuming permeable layers of and thickness resulted in effective hydraulic conductivities of and , respectively (please refer to Supplementary 注入过程由液态水（包括雨水和融水）的体积而不是压力决定。换句话说，模型假定生成的液态水从左上方流入渗透层。最后，强度降低法是在冰川上逐步实施的。由于阿鲁-1 冰川于 2016 年夏季坍塌，当时的 FOS 必须为 1.00。为此，通过反复试验，将假定的 1 米厚渗透层的饱和导流系数校准为 。Iken 等人（1996 年）通过钻孔响应测试，推断 Gorner 冰川下的冰碛物的导水率较高，为 。同样，Stone 和 Clarke（1993 年）报告了特拉普里奇冰川的较大数值，Stone 等人（1997 年）利用反演理论将排水系统相连部分的估计值修正为 5 。在我们的研究中，假设渗透层厚度为 和 ，则有效导流系数分别为 和 （请参阅补充资料）。
Material). While the calibrated effective hydraulic conductivities in this study surpass those of typical soil materials (e.g., to for clayey soils (Kluger et al., 2022)), they align with the range reported in the glacier studies mentioned above. Then, the calibrated hydraulic conductivity is applied to the Aru-2 glacier modelling, in which whether the FOS is still equal to 1 is a critical question to be checked. 材料）。虽然本研究中标定的有效水力传导率超过了典型土壤材料的有效水力传导率（例如，粘土的有效水力传导率为 到 （Kluger等人，2022年）），但它们与上述冰川研究中报告的范围一致。然后，将校准后的水力传导系数应用于 Aru-2 冰川建模，其中 FOS 是否仍然等于 1 是一个需要检查的关键问题。
4.2. Modelling results 4.2.建模结果
Fig. 10a illustrates the temperature evolution in the Aru-1 modelling example from 2000 to 2016 . The temperature along the profile ( ) reveals that the ice temperature near the top surface is more sensitive to air temperature. More notably, the glacier tends toward a warmer state; for instance, the temperature at Point C below the ice surface; see Fig. 10a) increased by from 2000 to 2016. Fig. 10b displays the simulated subglacial fluid pressure in summer 2016. The fluid pressure distributions of the two glaciers are somewhat similar, exhibiting lower values at the inlet and outlet and higher values in the middle. Nevertheless, unlike Aru-1, the subglacial water flow in Aru-2 experiences two transitions from a steeper bed to a flatter bed (at approximately and ), leading to two peaks in the fluid pressure distribution. 图 10a 展示了 2000 年至 2016 年 Aru-1 模拟实例的温度变化。沿 剖面（ ）的温度显示，顶面附近的冰温对气温更为敏感。更值得注意的是，冰川趋向于变暖状态；例如，从 2000 年到 2016 年，冰面以下 C 点 的温度（见图 10a）增加了 。图 10b 显示了 2016 年夏季冰川下流体压力的模拟情况。两座冰川的流体压力分布有些相似，入口和出口处的数值较低，中间的数值较高。然而，与阿鲁-1 不同的是，阿鲁-2 的冰川下水流经历了两次从较陡床面向较平坦床面的过渡（大约在 和 处），导致流体压力分布出现两个峰值。
This study focuses on the FOS values and collapse zones; the results are depicted in Fig. 10c. Interestingly, Aru-2 had a FOS of 0.988 in summer 2016, indicating a potentially unstable state of Aru-2 under the calculated fluid pressure and temperature conditions. More critically, this result suggests a similar stability status between the two glaciers at that time, given their nearly equal FOS values. It is worth mentioning that the FOS for both Aru-1 and Aru-2 would significantly increase when disregarding the influence of the basal fluid pressure, achieving 2.036 and 2.147, respectively. These results demonstrate the profound impact of water infiltration on the Aru glacier collapses. Regarding the effect of glacier temperature, the seasonal temperature rise from winter to summer only causes a slight decrease in the FOS, with the FOS changing less than in both modelling examples. 本研究的重点是 FOS 值和塌陷区；结果如图 10c 所示。有趣的是，2016 年夏季阿鲁-2 的 FOS 值为 0.988，这表明在计算的流体压力和温度条件下，阿鲁-2 可能处于不稳定状态。更关键的是，由于两座冰川的 FOS 值几乎相等，这一结果表明当时两座冰川的稳定性状态相似。值得一提的是，如果不考虑基底流体压力的影响，阿鲁-1 和阿鲁-2 冰川的 FOS 值都会显著增加，分别达到 2.036 和 2.147。这些结果表明了水的渗透对阿鲁冰川崩塌的深刻影响。至于冰川温度的影响，从冬季到夏季的季节性温度升高只导致 FOS 的轻微下降，在两个模拟实例中，FOS 的变化都小于 。
Besides the simulated plastic zones, the surface topography of the collapse zones, obtained through remote sensing (Gilbert et al., 2018), is also illustrated in Fig. 10c for comparison. The simulations match the reality well. Moreover, a certain degree of plastic strain concentration is simulated near the middle of Aru-2 Glacier, in addition to the collapse 除模拟塑性区外，图 10c 中还显示了通过遥感获得的塌陷区表面地形（Gilbert 等人，2018 年），以供比较。模拟结果与实际情况十分吻合。此外，在阿鲁-2 冰川中部附近模拟出了一定程度的塑性应变集中，除了塌陷外
years has potentially contributed to the 2016 Aru glacier collapse. 年可能是造成 2016 年阿鲁冰川崩塌的原因之一。
Fig. 9. Monthly mean (a) wind speed and (b) albedo in the Aru Range in 2016. The wind speed, as well as its variability, is generally higher in the cold season than in the warm season, while the glacier surface albedo and its albedo variation show a reverse relation with the season. 图 9.2016 年阿鲁山脉月平均（a）风速和（b）反照率。寒冷季节的风速及其变化一般高于温暖季节，而冰川表面反照率及其反照率变化则与季节呈反向关系。
position near the crown, indicating that the middle part may also be a relatively fragile "junction" in Aru-2. This weak "junction" may lead to the premature detachment of the lower glacier part, which explains the two flow fronts of Aru-2 Glacier (as shown in the two sections of the resultant alluvial fan, i.e., I and II in Fig. 2). 这表明中间部分也可能是阿鲁-2 号冰川相对脆弱的 "交界处"。这种脆弱的 "交界处 "可能会导致冰川下部过早脱离，这也是阿鲁-2 号冰川出现两个流锋的原因（如图 2 中冲积扇的两个切面所示，即图 I 和图 II）。
5. Discussions 5.讨论情况
5.1. The established methodology 5.1.既定方法
The calculated FOS values reveal that, if water infiltration through the crevasses is disregarded in the FOS computation, Aru-2 has higher stability than Aru-1 in summer 2016, suggesting that Aru-2 might be more stable than Aru-1 in the years with less water infiltration than in 计算得出的 FOS 值显示，如果在计算 FOS 值时不考虑水从裂缝中渗入的情况，Aru-2 在 2016 年夏季的稳定性要高于 Aru-1，这表明在水渗入量低于 Aru-1 的年份，Aru-2 可能比 Aru-1 更稳定。
2016. Fig. 11 also indicates that the FOS values of Aru-2 are higher than those of Aru-1 before failure. These results might be part of the reason why Aru-2 had a lower annual mean sliding velocity and fewer crevasses in the several years before the collapse (Gilbert et al., 2018), in addition to a high-friction patch under the Aru-2 glacier tongue (Gilbert et al., 2018). Moreover, the modelling provides quantitative evidence of similar stability levels between Aru-1 and Aru-2 in summer 2016, thereby explaining why Aru-1 and Aru-2 collapsed successively within only about two months. Additionally, the collapse zones of the glaciers are accurately reproduced by the strength reduction approach (SRA), which has been widely used to determine detachment zones of rock/soil slopes but, to the best of our knowledge, has never been validated for glacier detachment analysis. This study highlights the potential of using SRA to forecast detachment zones of large-scale low-angle ice-rich 2016.图 11 还表明，Aru-2 的 FOS 值高于崩塌前的 Aru-1 的 FOS 值。除了 Aru-2 冰川舌下的高摩擦斑块（Gilbert 等人，2018 年）之外，这些结果可能也是 Aru-2 在崩塌前的几年中年平均滑动速度较低、裂缝较少的部分原因（Gilbert 等人，2018 年）。此外，建模还提供了定量证据，证明 2016 年夏季阿鲁-1 和阿鲁-2 冰川的稳定性水平相似，从而解释了为什么阿鲁-1 和阿鲁-2 冰川仅在两个月左右的时间内相继崩塌。此外，强度降低方法（SRA）准确地再现了冰川的崩塌区，该方法已被广泛用于确定岩石/土壤斜坡的剥离区，但据我们所知，从未在冰川剥离分析中得到验证。本研究强调了使用强度还原法预测大尺度低角度富冰斜坡脱离带的潜力。
(b) Subglacial fluid pressure in summer 2016 (b) 2016 年夏季的冰川下流体压力
(c) Simulated FOS and collapse zones in 2016 (c) 2016 年模拟 FOS 和塌陷区
Fig. 10. (a) Temperature variations in the Aru-1 glacier. The glacier temperature shows an uptrend trend and the change in the glacier temperature lags behind that in the air temperature. (b) Basal fluid pressure distributions in Aru-1 and Aru-2 in summer 2016. (c) FOS values and collapse zones of the Aru glaciers. The collapse zones are represented by the concentration areas of equivalent plastic strain. For comparison, the post-collapse glacier surface topography is also displayed here. The winter here refers to December 2015 to February 2016. 图 10. (a) 阿鲁-1 号冰川的温度变化。冰川温度呈上升趋势，冰川温度的变化滞后于空气温度的变化。(b) 2016 年夏季阿鲁-1 和阿鲁-2 的基底流体压力分布。(c) 阿鲁冰川的 FOS 值和塌陷区。塌陷区由等效塑性应变集中区域表示。为便于比较，此处还显示了塌陷后的冰川表面地形。此处的冬季是指 2015 年 12 月至 2016 年 2 月。
slopes and to serve as a complementary tool to monitoring technologies (e.g., remote sensing) for related disaster risk management. The SRA may be particularly useful for ice-rich slopes that lack precursor information (e.g., large displacement) prior to detachment and for long-term predictions of detachment zones before noticeable crevasses emerge. The Aru modelling examples in this study provide references for stability assessments of ice-rich slopes analogous to the Aru glaciers. When applying the developed model to an ice-rich moraine slope, it is necessary to input the parameters specific to the moraine material. For instance, in thermodynamics modelling for an ice-rich slope with relatively high soil content, specific heat capacity and thermal conductivity are no longer provided by Eqs. (2) and (3) but need to be theoretically derived or obtained through experiments 在相关的灾害风险管理中作为监测技术（如遥感）的补充工具。特别是对于在崩塌前缺乏前兆信息（如大位移）的富冰斜坡，以及在出现明显裂缝前对崩塌区进行长期预测，SRA 可能特别有用。本研究中的阿鲁模型实例为类似阿鲁冰川的富冰斜坡稳定性评估提供了参考。在将开发的模型应用于富冰碛斜坡时，有必要输入冰碛材料的特定参数。例如，在对土壤含量相对较高的富冰斜坡进行热力学建模时，公式（2）和（3）不再提供比热容和导热系数，而是需要从理论上推导或通过实验获得。
Fig. 11. Comparison of FOS values of Aru-1 and Aru-2 between winter 2013 and spring 2016. The Aru-2 glacier had slightly higher stability than the Aru-1 glacier. These results may, to some extent, explain why Aru-2 had a lower annual mean sliding velocity and fewer crevasses than Aru-1 in the several years before their collapses. 图 11.2013 年冬季至 2016 年春季阿鲁-1 和阿鲁-2 的 FOS 值对比。阿鲁-2冰川的稳定性略高于阿鲁-1冰川。这些结果在一定程度上可以解释为什么阿鲁-2冰川在崩塌前几年的年平均滑动速度比阿鲁-1冰川低，裂缝也比阿鲁-1冰川少。
5.2. Role of climate change in the Aru glacier collapses 5.2.气候变化在阿鲁冰川崩塌中的作用
Efforts have been made to assess whether there is a relation between the changing climate and the Aru glacier collapses. Based on the ERAinterim reanalysis dataset (resolution: ), Kääb et al. (2018) modelled the cumulative glacier water equivalent (w.e.) mass balance of the Aru glaciers and reported an upward trend of air temperature of per year at a.s.l. over the Aru Range. To gain insights into the Aru glacier collapses, Zhang et al. (2023) analysed air temperature and precipitation data recorded from the Shiquanhe Meteorological Station ( a.s.l.) and Gerze Meteorological Station ( a.s.l.), which are, however, southwest and southeast of the Aru Co, respectively. In the present study, the meteorological information from the ERA5-Land hourly dataset (Muñoz-Sabater et al., 2021) has been processed to drive the above modelling examples; as much added value of ERA5-Land through the comparison of ERA5-Land has been demonstrated by Muñoz-Sabater et al. (2021). Compared to ERA-Interim, ERA5-Land is based on an improved data assimilation system, which has absorbed millions of extra observations (Muñoz-Sabater et al., 2021). Moreover, ERA5-Land has a higher horizontal native resolution of compared to ERA-Interim. 人们一直在努力评估气候变化与阿鲁冰川崩塌之间是否存在关系。基于ERAinterim再分析数据集（分辨率： ），Kääb等人（2018年）模拟了阿鲁冰川的累积冰川水当量（w.e.）质量平衡，并报告了阿鲁山脉 海拔高度气温每年 的上升趋势。为了深入了解阿鲁冰川崩塌的情况，Zhang 等人（2023 年）分析了石泉河气象站（ a.s.l.）和格尔木气象站（ a.s.l.）记录的气温和降水数据，但这两个气象站分别位于阿鲁科考站的 西南方和 东南方。在本研究中，ERA5-Land 小时数据集（Muñoz-Sabater 等人，2021 年）中的气象信息经过处理后用于驱动上述模拟实例；Muñoz-Sabater 等人（2021 年）通过对 ERA5-Land 的比较证明了 ERA5-Land 的附加值。与ERA-Interim相比，ERA5-Land基于改进的数据同化系统，吸收了数百万额外的观测数据（Muñoz-Sabater等人，2021年）。此外，与ERA-Interim相比，ERA5-Land的水平原生分辨率 更高。
Considering the state-of-the-art features of ERA5-Land, we analysed this dataset to conduct a detailed analysis of meteorological and hydrological data of the Aru Range. Fig. 12 depicts the quarterly mean air temperature at a.s.l. in the Aru Range and the total amount of meltwater and rainwater determined from the catchments above 考虑到 ERA5-Land 的先进功能，我们分析了该数据集，对阿鲁山脉的气象和水文数据进行了详细分析。图 12 描述了阿鲁山脉 海拔高度的季度平均气温，以及上方集水区测定的融水和雨水总量。
Fig. 12. Simulated lower bounds of FOS of Aru-1 in response to variations in meltwater, rainwater and air temperature at above sea level in the Aru Range. Due to the rising temperature and augmented meltwater and rainwater (in equivalent water depth) caused by climate change, glacier stability has likely undergone considerable alterations in the 2010s. Particularly, substantial loss of glacier stability occurred in 2010, 2012 and 2013 . The FOS values for 2000 to 2015 only represent the lower bounds of the glacier stability because they are also determined by assuming all upstream-coming liquid water reaching the bed, which is not true since the crevasses in 2000 to 2015 were not so developed and the amount of liquid water reaching the bed in each year from 2000 to 2010 was not as much as that in 2016 . 图 12.阿鲁山脉海平面以上 处，阿鲁-1的FOS随融水、雨水和气温变化的模拟下限。由于气候变化导致气温升高、融水和雨水增加（等效水深），冰川的稳定性在 2010 年代可能发生了巨大变化。特别是在 2010 年、2012 年和 2013 年，冰川稳定性大幅下降。2000 年至 2015 年的 FOS 值仅代表冰川稳定性的下限，因为它们也是在假定所有上游来的液态水都到达冰床的情况下确定的，但事实并非如此，因为 2000 年至 2015 年的裂缝并不发达，而且 2000 年至 2010 年每年到达冰床的液态水量也没有 2016 年那么多。
detachment zones. The air temperature at a.s.l. shows a rising rate of per year, surpassing that at a.s.1. per year) reported by Kääb et al. (2018). Moreover, the present study reveals that the summer mean air temperature at a.s.l. increased at a higher rate per year. Thus, the warming trend at higher elevations in the Aru Range, especially in warm season, might be more severe than commonly thought. 脱离区。 a.s.l.的气温以每年 的速度上升，超过了 a.s.1。Kääb等人（2018年）报告的 每年的上升率。）此外，本研究还发现， a.s.l.夏季平均气温的年增长率 更高。因此，阿鲁山脉海拔较高地区的变暖趋势，尤其是暖季的变暖趋势，可能比通常认为的更为严重。
Note that the FOS values for 2000-2015 in Fig. 12 only represent lower bounds of glacier stability because they are calculated based on the same assumptions and parameter values used for 2016. In other words, because crevasses in the glaciers were less developed before 2016, some liquid water from the basins above a.s.l. might flow away directly through supraglacial watercourses. Thus, the FOS values for 2000 to 2015 , obtained by assuming all upstream-coming liquid water reaching the bed, underestimate practical glacier stability and can only be deemed as lower bounds of the stability of the glaciers. In a practical application of the developed model, a more exact estimate of glacier stability can be achieved by more accurate monitoring of upstream-coming meltwater that flows into the crevasses. It is worth mentioning that, one reason why the two smaller glaciers to the east of Aru-2 in Fig. 6 did not collapse may be due to the underdeveloped crevasses on the glacier surfaces, resulting in less liquid water reaching the beds. Generally, the lower bound of FOS in Fig. 12 is inversely related to air temperature and water input. The Aru glaciers' stability probably fluctuated only slightly before 2010 . However, the present study indicates that the mean temperature in summer 2010 was relatively high and there was a tremendous amount of meltwater and rainwater supply. In the subsequent years, the summer mean air temperature and water amount frequently rose to an abnormal level, and the lower bound of FOS descended rapidly, often dropping to 1 from spring to summer. 请注意，图 12 中 2000-2015 年的 FOS 值仅代表冰川稳定性的下限，因为它们是根据 2016 年使用的相同假设和参数值计算得出的。换句话说，由于 2016 年之前冰川裂缝不太发达，来自 海拔以上盆地的一些液态水可能会通过冰川上水道直接流走。因此，假定上游流入的液态水全部到达冰床而得出的 2000 至 2015 年的 FOS 值低估了冰川的实际稳定性，只能被视为冰川稳定性的下限。在实际应用所开发的模型时，可以通过对流入裂缝的上游融水进行更精确的监测，对冰川的稳定性做出更准确的估计。值得一提的是，图 6 中阿鲁-2 以东的两个较小的冰川之所以没有坍塌，原因之一可能是冰川表面的裂缝不发达，导致到达冰床的液态水较少。一般来说，图 12 中的 FOS 下限与气温和输入水量成反比。阿鲁冰川的稳定性在 2010 年之前可能只有轻微波动。但本研究表明，2010 年夏季的平均气温相对较高，融水和雨水供应量巨大。在随后的几年中，夏季平均气温和水量经常上升到异常水平，FOS 下限迅速下降，从春季到夏季经常降至 1。
Notably, our results show that the summer mean air temperatures in 2010, 2012 and 2013 were similar to or higher than that in 2016. Moreover, the water amount was larger in 2010, 2012 and 2013 than in 2016. Consequently, the damage to the glaciers' stability in 2010,2012 and 2013 might have also significantly contributed to the 2016 Aru events. This inference is corroborated by the evolution of crevasses at the upper edges of the detachment zones. As shown in Fig. 3, few crevasses were visible on Aru-2 Glacier before April 2009. However, crevasses proliferated in the following years, indicating severe damage accumulation on the glaciers during this period. Therefore, the collapse of the Aru glaciers was strongly linked to long-term climate change. Such low-angle ice-rich landslides will likely continue under ongoing warming climate or extreme events of high air temperature or abundant precipitation. 值得注意的是，我们的研究结果表明，2010 年、2012 年和 2013 年的夏季平均气温与 2016 年相似或更高。此外，2010 年、2012 年和 2013 年的水量也大于 2016 年。因此，2010 年、2012 年和 2013 年对冰川稳定性的破坏可能也是 2016 年阿鲁事件的重要原因。脱离带上部边缘裂缝的演变也证实了这一推论。如图 3 所示，2009 年 4 月之前，阿鲁-2 冰川上几乎看不到裂缝。然而，在随后的几年里，裂缝大量增加，表明在此期间冰川上积累了严重的破坏。因此，阿鲁冰川的崩塌与长期气候变化密切相关。在气候持续变暖或气温过高或降水量过多的极端情况下，这种低角度富冰滑坡可能会继续发生。
6. Summary and conclusions 6.摘要和结论
The instability of low-angle ice-rich slopes involves complex multiphysical processes. The mechanisms of such ice-rich slope instability have yet to be fully elucidated, and there is a dearth of methodologies for quantitatively assessing the stability of such slopes. Referencing the collapse mechanisms of the Aru glaciers, we established a conceptual model and formulated to quantify the stability of ice-rich slopes through coupled thermo-hydro-mechanical simulation. The following conclusions are drawn: 低角度富冰斜坡的不稳定性涉及复杂的多重物理过程。这种富冰斜坡失稳的机理尚未完全阐明，也缺乏定量评估这种斜坡稳定性的方法。参照阿鲁冰川的崩塌机制，我们建立了一个概念模型，并通过热-水-机械耦合模拟来量化富冰斜坡的稳定性。得出以下结论：
The modelling examples show that Aru-1 and Aru-2 had similar stability (FOS ) in the face of meltwater and rainwater in summer 2016, elucidating why they collapsed successively within only about two months. Moreover, the safety factors of the Aru glaciers may have decreased by over due to the effects of meltwater and rainwater in summer 2016. This underscores the importance of water infiltration within an ice-rich slope system. 建模实例表明，2016年夏季，阿鲁-1和阿鲁-2冰川在融水和雨水面前具有相似的稳定性（FOS ），这就解释了它们为何在短短两个月左右的时间内相继崩塌。此外，由于2016年夏季融水和雨水的影响，阿鲁冰川的安全系数可能下降了 以上。这凸显了水在富冰斜坡系统中渗透的重要性。
In addition to the advancement in ice-rich slope stability assessment, the proposed approach also reproduces well the collapse zones of the two glaciers in a two-dimensional framework. The agreement between the simulation results and reality corroborates the proposed approach for evaluating the stability and potential collapse zones of ice-rich slopes. 除了在富冰斜坡稳定性评估方面取得进展外，所提出的方法还在二维框架内很好地再现了两条冰川的崩塌区。模拟结果与实际情况的一致性证实了所提出的富冰斜坡稳定性和潜在崩塌区评估方法。
The proposed approach enables the incorporation of climate-related variables into the assessment of ice-rich slope stability, encompassing elements like meltwater, precipitation, air temperature, glacier surface albedo, and wind speed. The approach holds promising potential for addressing future ice-rich landslides akin to the Aru events, which, as indicated in the present study, are closely linked to climate change. 所提出的方法可将气候相关变量纳入富冰斜坡稳定性评估，包括融水、降水、气温、冰川表面反照率和风速等要素。该方法有望解决未来类似阿鲁事件的富冰滑坡问题，正如本研究中指出的，这些问题与气候变化密切相关。
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 作者声明，他们没有任何可能会影响本文所报告工作的已知经济利益或个人关系。
Data availability 数据可用性
Data will be made available on request. 数据将应要求提供。
Acknowledgments 致谢
The authors are grateful to the financial support from the National Natural Science Foundation of China (Nos. U20A20112 and 41941017), the Project of Hetao Shenzhen-Hong Kong Science and Technology Innovation Cooperation Zone (HZQB-KCZYB-2020083), and Research Grants Council of the HKSAR Government (No. 16210122). 作者感谢国家自然科学基金（编号：U20A20112和41941017）、河套深港科技创新合作区项目（HZQB-KCZYB-2020083）和香港特别行政区政府研究资助局（编号：16210122）的资助。
Appendix A. Supplementary data 附录 A.补充数据
Supplementary data to this article can be found online at https://doi. org/10.1016/j.enggeo.2024.107514. 本文的补充数据可在线查阅：https://doi. org/10.1016/j.enggeo.2024.107514。
References 参考资料
Aschwanden, A., Blatter, H., 2009. Mathematical modeling and numerical simulation of polythermal glaciers. J. Geophys. Res. 114, F01027. https://doi.org/10.1029, 2008JF001028. Aschwanden, A., Blatter, H., 2009.Mathematical modeling and numerical simulation of polythermal glaciers.J. Geophys.114, F01027.https://doi.org/10.1029, 2008JF001028.
Evans, S.G., Tutubalina, O.V., Drobyshev, V.N., Chernomorets, S.S., Dougall, S., Petrakov, D.A., et al., 2009. Catastrophic detachment and high-velocity long-runout flow of Kolka Glacier, Caucasus Mountains, Russia in 2002. Geomorphology 105 (3-4), 314-321. https://doi.org/10.1016/j.geomorph.2008.10.008. Evans, S.G., Tutubalina, O.V., Drobyshev, V.N., Chernomorets, S.S., Dougall, S., Petrakov, D.A., et al., 2009.2002年俄罗斯高加索山脉科尔卡冰川的灾难性脱离和高速长流。Geomorphology 105 (3-4), 314-321.https://doi.org/10.1016/j.geomorph.2008.10.008.
Faillettaz, J., Funk, M., Vincent, C., 2015. Avalanching glacier instabilities: Review on processes and early warning perspectives. Rev. Geophys. 53, 203-224. https://doi. org RG000466 Faillettaz, J., Funk, M., Vincent, C., 2015.冰川崩塌不稳定性：过程回顾与预警展望。Rev. Geophys.53, 203-224.https://doi. org RG000466
Falaschi, D., Kääb, A., Paul, F., Tadono, T., Rivera, J.A., Lenzano, L.E., 2019. Brief communication: Collapse of of ice from a cirque glacier in the Central Andes of Argentina. Cryosphere 13, 997-1004. https://doi.org/10.5194/tc-13-997-2019. Falaschi, D., Kääb, A., Paul, F., Tadono, T., Rivera, J.A., Lenzano, L.E., 2019.Brief communication：阿根廷中部安第斯山脉一个峡谷冰川的 冰坍塌。冰冻圈 13, 997-1004.https://doi.org/10.5194/tc-13-997-2019.
Fu, Y., Jiang, Y., Wang, J., Liu, Z., Lu, X., 2021. Mechanical properties of frozen glacial tills due to short periods of thawing. Front. Earth Sci. 9, 799467 https://doi.org/ 10.3389 /feart.2021.799467. Fu, Y., Jiang, Y., Wang, J., Liu, Z., Lu, X., 2021.短时间解冻导致的冰川冻土力学特性。Front.9, 799467 https://doi.org/ 10.3389 /feart.2021.799467.
Gilbert, A., Gagliardini, O., Vincent, C., Wagnon, P., 2014. A 3-D thermal regime model suitable for cold accumulation zones of polythermal mountain glaciers. J. Geophys. Res. 119, 1876-1893. https://doi.org/10.1002/2014JF003199. Gilbert, A., Gagliardini, O., Vincent, C., Wagnon, P., 2014.适用于多热山地冰川冷积聚区的三维热机制模型。J. Geophys.119, 1876-1893.https://doi.org/10.1002/2014JF003199.
Gilbert, A., Leinss, S., Kargel, J., Kääb, A., Gascoin, S., Leonard, G., et al., 2018 Mechanisms leading to the 2016 giant twin glacier collapses, Aru Range, Tibet. Cryosphere 12, 2883-2900. https://doi.org/10.5194/tc-12-2883-2018. Gilbert, A., Leinss, S., Kargel, J., Kääb, A., Gascoin, S., Leonard, G., et al., 2018 Mechanisms leading to the 2016 giant twin glacier collapses, Aru Range, Tibet.Cryosphere 12, 2883-2900.https://doi.org/10.5194/tc-12-2883-2018.
