Three-Dimensional Numerical Modeling of Ground Ice Ablation in a Retrogressive Thaw Slump and Its Hydrological Ecosystem Response on the Qinghai-Tibet Plateau, China

Abstract

Retrogressive thaw slumps (RTSs), which frequently occur in permafrost regions of the Qinghai-Tibet Plateau (QTP), China, can cause significant damage to the local surface, resulting in material losses and posing a threat to infrastructure and ecosystems in the region. However, quantitative assessment of ground ice ablation and hydrological ecosystem response was limited due to a lack of understanding of the complex hydro-thermal process during RTS development. In this study, we developed a three-dimensional hydro-thermal coupled numerical model of a RTS in the permafrost terrain at the Beilu River Basin of the QTP, including ice–water phase transitions, heat exchange, mass transport, and the parameterized exchange of heat between the active layer and air. Based on the calibrated hydro-thermal model and combined with the electrical resistivity tomography survey and sample analysis results, a method for estimating the melting of ground ice was proposed. Simulation results indicate that the model effectively reflects the factual hydro-thermal regime of the RTS and can evaluate the ground ice ablation and total suspended sediment variation, represented by turbidity. Between 2011 and 2021, the maximum simulated ground ice ablation was in 2016 within the slump region, amounting to a total of 492 m³, and it induced the reciprocal evolution, especially in the headwall of the RTS. High ponding depression water turbidity values of 28 and 49 occurred in the thawing season in 2021. The simulated ground ice ablation and turbidity events were highly correlated with climatic warming and wetting. The results offer a valuable approach to assessing the effects of RTS on infrastructure and the environment, especially in the context of a changing climate.

1 Introduction

Retrogressive thaw slumps (RTSs), as a kind of typical thermokarst landslides, develop in gently sloping areas underlain by permafrost with high ice content or massive ground ice in Arctic regions and the Qinghai-Tibet Plateau (QTP) (Kanevskiy et al. 2013; French 2017; Lewkowicz and Way 2019; Mu et al. 2020; Morino et al. 2021). Due to the ongoing process of climatic change and the intensification of human activities, the phenomenon of RTSs has been observed to occur with increasing frequency and severity (Fraser et al. 2018; Luo, Niu, et al. 2019; Luo et al. 2022). Normally, the thaw-induced slope failure occurred following disturbances of the ground surface and destruction of vegetation by anthropic activities, wildfires, and anomalous climatic warming and wetting (Everett 1989). After active layer detachment, a steep and frozen back scarp or headwall develop at the end of the slope (Luo et al. 2022). This headwall comprises exposed ground ice that experiences continuous melting, thereby resulting in the downhill movement of thawed materials and mixed meltwater sediments, either through sliding or flowing (Jiang et al. 2022; Jiao et al. 2022).

The QTP is a major region for the occurrence of ice-rich permafrost in China. However, the plateau has experienced significant climatic warming and wetting in recent years (Zhao et al. 2004), which has led to the rapid development of thermokarst landslides over the past decades (Luo, Niu, et al. 2019; Luo et al. 2022). In recent years (2008/2013–2022), field surveys and remote sensing interpretations have identified 2669 active RTSs in the permafrost regions of the QTP, as depicted in Fig. 1a (Luo et al. 2022). This accelerated development has led to severe ecological degradation on the QTP, particularly in its hilly and mountainous regions, including the Hoh Xil Mountain (Fig. 1b), Hongliang River region (Fig. 1c), Fenghuo Mountain area (Fig. 1d), and Beilu River (BLR) region (Figs. 1e and f).

Fig. 1

Spatial distribution of a permafrost and active retrogressive thaw slumps (RTSs) on the QTP based on the Interpreted Active RTSs Data (Luo et al. 2022) and the QTP Science Data Center for China, National Science & Technology Infrastructure; and the typical RTSs developed in b the Hoh Xil Mountain, c the Hongliang River region, d the Fenghuo Mountain area, and e, f the Beilu River region

Therefore, the evolutionary progression of RTSs has provided substantial evidence of the significant impact of ice-rich sediments thawing on the hydrogeological ecosystems. RTSs are widely distributed and have been shown to affect various aspects of geomorphology, such as the thermo-erosion of coastlines or lakeshores (Nicu et al. 2021) and the formation of negative terrain (Liu et al. 2014). They are also responsible for the wasting of fine-grained materials, including fluvial sediment transport (Droppo et al. 2022) and surface runoff (Yang et al. 2021). Importantly, they can also impact local ecological environments, such as lake water quality (Kokelj et al. 2005), the salinity of soils (Lantz et al. 2009), and the release of organic carbon (Masyagina and Menyailo 2020; Du et al. 2023) and methane from buried frozen sediments (Droppo et al. 2022). These impacts not only contribute to climatic warming but also have significant implications for global carbon cycling.

Despite the substantial body of research examining the phenomenon of RTSs (Mu et al. 2020; Morino et al. 2021), much of the previous literature has primarily focused on qualitative descriptions of the evolution of topography and geomorphology from the perspective of periglacial landforms. Such studies have often been restricted to describe the process of RTS development (Sun et al. 2017; Swanson 2021), thereby limiting the ability to assess ecosystem influences and hazards, and thus hampering efforts to comprehend the mechanism of RTS (Tebbens 2020). Indeed, remote sensing techniques have shown great potential in advancing our understanding of RTS and its impact on the eco-environmental system. Several recent studies have employed remote sensing data and techniques to investigate the spatial and temporal dynamics of RTS and to assess their impact on the surrounding environment. For example, Huang et al. (2021) used high-resolution remote sensing images to monitor the changes in size and shape over space and time of a RTS in the northern QTP. Their results showed that RTSs’ polygon had increased significantly over the study period, and that the affected area had expanded from 2016 to 2019. Van der Sluijs et al. (2018) employed thermal infrared unmanned aerial vehicle (UAV) data to study the changes in the temperature of the ground surface and the thickness of the active layer over time and space in northwestern Canada. Meng et al. (2022) conducted a study on an earthflow in the Zhimei pastureland of the QTP to understand how permafrost landforms are affected by climate change, leading to landslides. They used various techniques to analyze the deformation characteristics, internal structure, and driving factors behind the earthflow, and their findings can be used as a reference for understanding slope destabilization in frozen regions due to climate change. Turner et al. (2021) conducted a study using drones to understand the impact of climate change on the largest RTS along the Old Crow River in Yukon, Canada. They found that the RTS had exported over 29,000 m³ of sediment, affecting local travel and infrastructure, and estimated that it had also exported significant amounts of soil organic carbon and total nitrogen. Their results revealed that the warm environment of the active layer was significantly affected by the presence of RTS, resulting in a thickening of active layer depth or thawing front. Thus, the use of remote sensing data and techniques has great potential to contribute to the development of a more comprehensive understanding of RTS. To the best of our knowledge, the means to accurately quantify the melting of exposed massive ground ice in the development of RTSs are limited (Lewkowicz 1985; Chen et al. 2018; Farquharson et al. 2019; Ohara et al. 2022; Wang et al. 2022). Among the existing methods, a majority of them have used dowels that were embedded into holes drilled in the ground-ice slump floor (Kerfoot 1969; Burn 1983; Burn and Friele 1989). Notably, one of the methods involved the development of ablatometers designed to monitor the rate of ablation on the slump floor located in the Western Canadian Arctic (Lewkowicz 1985, 1986). A recent study by Wang et al. (2022) estimated the amount of ground ice in the permafrost layers of the source area of the Yellow River in the northeastern part of the QTP, and found that the amount of ground ice varies depending on the type of land and depth below the surface. Another study by Chen et al. (2018) used persistent scatterer interferometric synthetic aperture radar (PSI) to measure gradual surface subsidence due to permafrost thaw at the Eboling Mountain in the northeastern region of the QTP, finding that permafrost areas near thermal erosion gullies are more vulnerable to degradation and demonstrating the potential of using PSI to study permafrost thaw processes over vast areas. A two-dimensional numerical model based on an energy equation was also used to calculate the melting of ground ice at a RTS study site on the QTP (Jiao et al. 2023). However, to our knowledge, there has been few comprehensive quantitative assessment of ground ice ablation volume in RTS development area. These studies conducted thus far have been unable to provide a quantitative assessment of ground ice ablation volume on the RTS, thereby rendering it challenging to fully comprehend the associated natural hazards and broader impacts on the ecological environment.

This is particularly problematic given the difficulty in assessing the eco-environmental impacts of RTS without quantifying the ground ice ablation. The quantification also holds potential significance for the future carbon emissions and soil erosion. Therefore, accurately assessing the impacts of eco-environmental changes requires quantification of ground ice or ice-rich permafrost thawing. This approach can improve comprehension of the intricate interplay between physical processes and ecological dynamics, enabling the development of more effective mitigation strategies. Our study aimed to quantify the ground ice ablation and eco-environmental impacts within the studied RTS and provide insights into the potential quantitative assessment of ground ice ablation of the entire QTP through the methods applied in this study. To achieve this goal, we selected a representative RTS located on the northeast-facing slope of the Gerlama Mountain (34°50′45″N, 92°54′27″E), which is about 300 m west of the K1129 marker on the Qinghai-Tibet Railway, within an area of continuous permafrost. This specific RTS, illustrated in Fig. 2, is referred to as the K1129W RTS for convenience. We developed a three-dimensional (3D) full coupled thermal-hydraulic model with solute transport of the K1129W RTS considering the local topographical factors. Through this modeling, our objectives were three-fold: (1) to examine the thermal-hydrological characteristics and permafrost table fluctuations; (2) to quantitatively estimate the ground ice ablation; and (3) to analyze the hydrogeological changes resulting from the presence of the K1129W RTS in the nearby region. The findings provide insights into the eco-environmental impacts of RTS development on the QTP, and can serve as a reference for future assessments of such impacts in the region.

Fig. 2

Aerial view and details of the K1129W retrogressive thaw slump (RTS). a Aerial photograph of the K1129W RTS, b details in the headwall, c thawed materials and small puddles below the headwall, d elevation of uplift at the toe of the K1129W RTS

2 Materials and Methods

This section commences with an introduction to the study area under investigation, followed by a comprehensive exposition of the field survey and data analysis methodologies, encompassing techniques such as electrical resistivity tomography (ERT), UAV imaging, and the estimation of ground ice content.

2.1 K1129W Retrogressive Thaw Slump (RTS) Development Area

Within the continuous permafrost zone of the QTP, this research focused on a RTS situated in the BLR region (34°58′N, 92°57′E), an area characterized by a landscape shaped by river sediment deposition and flooding (Lin et al. 2020). Despite representing less than 0.01% of the total land area of the QTP, the BLR region is home to nearly 30% of all RTSs identified in this region (Luo et al. 2022). With a mean elevation of 4673 m a.s.l., this region is situated between the Fenghuo Mountain in the south and Hoh Xil Mountain in the north, and has elevations that range from 4418 to 5320 m a.s.l. (Luo, Niu, et al. 2019). The RTS features a customary curved contour and maintains an average angle of 8.6°. It spans 145–168 m from the headwall to the leading edge of the compacted region, and has a width of 98 m at its broadest point located in the center. The K1129W RTS is developed in an area with active layer depth of 1.8 to 2.5 m, overlying a 2.5 m thick ice-rich permafrost (Jiao et al. 2022).

The upper layers predominantly consist of silty clay, reaching a maximum thickness of 8 m, while the lower layers mainly comprise of mudstone (Lin et al. 2011). The region experiences a mean annual air temperature (MAAT) of − 3.8 °C and receives an average annual precipitation of about 300 mm, mainly concentrated between May and August, based on data recorded from 1960 to 2016 (Luo, Niu, et al. 2019). The mean annual evaporation is 1316.9 mm. According to meteorological data from the past 50 years, the MAAT has increased by 0.31 °C/10a, and the precipitation increment was 20 mm/10a (Fig. 3). The K1129W RTS has a total disturbed area of 12,881 m², with a perimeter length of 604 m and a total displaced material of 13,284 m³. In the western region of the K1129W RTS, a headwall measuring between 1.4 and 2.4 m in height has developed, resulting in alpine meadow vegetation covering approximately 8–20% of the scar area (Fig. 2).

Fig. 3

Variations of mean annual air temperature and precipitation, 1950–2020

During each thawing season, tension fissures develop in the scarp of the K1129W RTS due to ice-rich permafrost thawing. These fissures grow in a circular pattern, leading to the retreat of the steep headwall and the downslope movement of the thawed soil and vegetation.

2.2 Field Survey

Fieldwork was carried out to examine the distribution of permafrost, ice content, and the depth of the permafrost table in the K1129W RTS region. The fieldwork consisted of several investigations, including geophysical survey, borehole sampling, and UAV investigation. First, to investigate the distribution of ground ice and permafrost table of the RTS development area, three ERT investigation lines (A–A′, B–B′, C–C′) were laid out in September 2021 (Fig. 4). The active layer had almost completely thawed during this period. The location of these three ERT investigation lines was shown in Fig. 4. Electrical resistivity values were measured using a set of direct-current resistivity system [Advanced Geosciences Inc. (AGI) Supersting R8] with 36 electrodes and 3 m spacing by using the Dipole–Dipole array. Second, in order to improve our understanding of the ground ice content in this region, we collected 125 samples from various sampling sites in the BLR region between 2016 and 2020 (see site locations in Fan et al. 2021, 2022). Twelve boreholes were situated at the shady slope (with an angle of 7° to 10°) of the Gerlama Mountain. Last, the UAV (DJI Matrice 300 RTK) that was equipped with a Zenmuse P1 camera captured high-resolution images of the study area. The images were used to create a digital elevation model (DEM). The UAV investigation, conducted in September 2021, provided valuable information on the topography and landscape of the study area. Overall, the data collected from these investigations were used to better understand the dynamics of permafrost in the region and to develop strategies for managing RTS.

Fig. 4

Unmanned aerial vehicle image-based DEM of the K1129W retrogressive thaw slump (RTS). Electrical resistivity tomography (ERT) survey lines A–A′, B–B′, and C–C′ are highlighted in bold for identification purposes.

2.3 Data Processing

To improve our analysis of ground ice content at the study site, we used the calorimetric method to directly measure gravimetric moisture contents (ω_M) of soil samples in the field (Lin et al. 2020). The equation for calculating the ω_M (%) is:

$$\omega_{M} = \frac{{m_{0} - m_{s} }}{{m_{s} }}$$

(1)

where m₀ is the mass of the natural core (g), and m_s is the dry weight of soil (g).

Volumetric ground ice content (v_ice) was estimated by permafrost bulk density (ρ) and ice content of the permafrost (Lin et al. 2020):

$$v_{ice} = \left( {\rho - \frac{{m_{s} }}{V} - \omega_{u} - \omega_{a} } \right)/\rho_{i} ,$$

(2)

where V is the volume of the sample (cm³), ρ_i is the density of ice (0.91 g/cm³), ω_u is the mean unfrozen water content (%), and ω_a is the air content (%). Through the computed tomography scan, ω_u and ω_a can be estimated and they are 1% and 2% in undisturbed samples (Lin et al. 2020; Fan et al. 2021).

Since the ERT survey was carried out in the K1129W RTS area, we did not have a reference value of the permafrost after the freeze-thaw erosion. To better interpret the electrical resistivity values of permafrost in the K1129W RTS area, a regressive equation between the electrical resistivity and unfrozen water content under different unfrozen water content was then defined as (Fortier et al. 2008):

$$R = e^{{ - \frac{{\omega_{u} }}{{\omega_{u0} }}}} \cdot R_{{\omega_{u0} }} ,$$

(3)

where R is the electrical resistivity value of permafrost (Ω m), $R_{{\omega_{u0} }}$ is the reference electrical resistivity value of 12,820 Ω m when the unfrozen water content ω_u0 = 5%, and ω_u is the unfrozen water content (%).

2.4 Estimating the Presence of Ground Ice in the Beilu River (BLR) Region

The moisture content of the studied soil samples was determined by analyzing 125 undisturbed samples, as reported in Fan et al. (2021, 2022). The soil moisture content was measured in the active layer of the boreholes, revealing an increasing trend with depth, as illustrated in Fig. 5. Additionally, 25 samples with gravimetric moisture content exceeding 20% were selected for further analysis. Based on the gravimetric moisture content of the upper permafrost (2 m to 5 m below the ground surface) displayed in Fig. 5, approximately 72% of the analyzed samples exhibited the ice-rich permafrost conditions (ω_M > 20%). These results provide valuable insights into the soil moisture dynamics and permafrost conditions in the studied region.

Fig. 5

Soil gravimetric water contents in the 12 boreholes. The permafrost table is marked, and the blue and green squares represent the ice-rich borehole and the ice-poor borehole, respectively.

Figure 6 shows the calculated volumetric ground ice content in the upper permafrost (about 2 m to 5 m in depth) of these boreholes. In the upper 0.5 m below the permafrost table, the mean volumetric ground ice content ranged from 30.8 to 114.2%. In the 2.5 m to 4 m depth, the mean volumetric ground ice content ranged from 22.3 to 67.2%. In the 4 m to 5 m depth, the mean volumetric ground ice content ranged from 21.2 to 46.8%. Based on the data above, the volumetric ground ice content at the uppermost part of the permafrost in this study area was approximately 50.4%. In summary, based on the ERT results and calculated volumetric ground ice content, the depth of ice-rich permafrost was about 2.5 m, which was used in the estimation of ground ice ablation.

Fig. 6

Box plots of the volumetric ground ice content of the 25 permafrost samples. The red dashed line is the mean value.

3 Numerical Simulation

The Code Heatflow software, developed for numerical analysis, was used in this study to address time-varying groundwater flow, heat transfer, and mass transfer. This software has been tested through benchmarks from InterFrost (Grenier et al. 2018), and has been precisely calibrated for simulation of permafrost mound dynamics (Dagenais et al. 2020) and thermokarst lake evolution (Perreault et al. 2021) in the permafrost regions. The code considers various parameters such as thermal conductance, heat retention, unfrozen water levels, relative permeability, latent and sensible heat, and modifications between ice and water phases (Molson and Frind 2020). The main modules used in this model are concisely outlined below.

3.1 Mathematical Model

The groundwater flux conservation equation (Molson and Frind 2020) is:

$$\frac{\partial }{{\partial x_{i} }}\left[ {K_{i,j} (T)\left( {\frac{\partial \psi }{{\partial x_{j} }} + \rho_{r} (T)n_{j} } \right)} \right] - \sum\limits_{k = 1}^{N} {Q_{k} (t) \cdot \delta (x_{k} ,y_{k} ,z_{k} )} = S_{s} \frac{\partial \psi }{{\partial t}},$$

(4)

where x_i represents the 3D spatial coordinates (m), K_i,j(T) denotes the hydraulic conductivity tensor that varies with temperature, T represents the temperature of the porous medium in degrees Celsius, ψ represents the freshwater head in meters, ρ_r(T) is the relative density of water that is dependent on temperature, n_j is a unit vector in the vertical Z-direction, Q_k represents the fluid volume flux in cubic meters per second for a source or sink located at (x_k, y_k, z_k), S_s represents the specific storage, and t denotes time. A summary of equivalent head concept in density-dependent flow is found in Post et al. (2007).

The equation governing heat transport with advection and dispersion is:

$$- \frac{\partial }{{\partial x_{i} }}(\theta S_{w} c_{w} \rho_{w} v_{i} T) + \frac{\partial }{{\partial x_{i} }}\left[ {\left( {\left( {\overline{\lambda } + \theta S_{w} c_{w} \rho_{w} D} \right)\frac{\partial T}{{\partial x_{i} }}} \right)} \right] + \Omega = \frac{{\partial (C_{0} T)}}{\partial t},$$

(5)

where porosity is denoted by θ, whereas unfrozen water saturation S_w is determined by multiplying the saturation of water in a fully thawed state with the temperature-dependent volume fraction of unfrozen water, W_u(T). Here, ρ_w, c_w, and v_i stand for density (kg/m³), specific heat (J/kg/°C), and average linear velocity (m/s) of ground water, respectively.

All medium temperature is represented by T (°C). The hydrodynamic dispersion coefficient is represented by D_i,j (m²/s), the term for thermal source/sink is represented by Ω (J/m³/s), the spatial coordinates are x_i,j (m), and time is represented by t (s). Lastly, C₀ is the capacity of the porous medium (J/m³/°C).

$$C_{0} = \theta S_{w} c_{w} \rho_{w} + \theta S_{i} c_{i} \rho_{i} + (1 - \theta )c_{s} \rho_{s} + \theta \rho_{i} L\left( {\frac{{\partial S_{w} }}{\partial T}} \right)$$

(6)

$$\overline{\lambda } = (\theta S_{w} \sqrt {\lambda_{w} } + \theta S_{i} \sqrt {\lambda_{i} } + (1 - \theta )\sqrt {\lambda_{s} } )^{2} ,$$

(7)

where the subscript w represents the water, and i and s represent ice and solids, correspondingly. L is the latent heat of water-ice phase change (J/kg). $\overline{\lambda }$ is the bulk thermal conductivity of the porous medium (like the permafrost).

The mass transport equation is:

$$\frac{\partial }{{\partial x_{i} }}\left[ {\left( {\frac{{D_{i,j} }}{R}} \right)\frac{\partial c}{{\partial x_{j} }}} \right] - \frac{\partial }{{\partial x_{i} }}\left( {\frac{{v_{i} }}{R}c} \right) + \sum\limits_{k = 1}^{N} {\frac{{Q_{k} (t)c_{k} (t)}}{R\theta }} \delta (x_{k}^{\prime } ,y_{k}^{\prime } ,z_{k}^{\prime } ) = \frac{\partial c}{{\partial t}},$$

(8)

where D_i,j and v_i stand for the vector of hydrodynamic dispersion tensor and average linear groundwater velocity, respectively. c stands for suspended sediment concentration (kg/m³), which is represented by turbidity, and R, as the retardation index, is calculated by Eq. 9.

$$R = 1 + \frac{{\rho_{b} k_{d} }}{\theta },$$

(9)

where ρ_b and K_d stand for bulk density of the saturated soil and distribution coefficient, respectively, and θ is the partition coefficient that controls the concentration of the dissolved particle.

Additionally, W_u and K_r stand for the unfrozen water saturation and relative permeability, respectively, which are functions of the ground temperature, and are defined by Eqs. 10 and 11, respectively (Palmer et al. 1992; Molson and Frind 2020):

$$W_{u} (T) = (1 - p)e^{{ - (T/q)^{2} }} + p,$$

(10)

$$k_{r} (T) = \max [10^{{ - \Psi \cdot \theta [1 - W_{u} (T)]}} ,10^{ - 6} ],$$

(11)

where p represents the water content in permafrost at a temperature range of − 1 to 0 °C, while q represents the shape factor of water content, and ψ represents the impedance factor (Jame and Norum 1980).

Table 1 lists the parameters used in Eqs. 10 and 11, along with the thermal characteristics of the solid, water, and ground ice components, and the convergence standards for the nonlinear coupling. For the benchmark InterFrost cases, we suggest using the specific temperature-dependent functions W_u and K_r.

Table 1 Physical parameters of the K1129W retrogressive thaw slump (RTS)

3.2 Quantitative Analysis of Ground Ice Melting Volume from the Retrogressive Thaw Slump (RTS)

To quantitatively analyze the melting volume of ground ice during the development of RTSs, we developed a systematic means based on the Stefan equation (Farquharson et al. 2019; Ohara et al. 2022). We first applied this equation to a single node of ice-rich permafrost to calculate its melting rate, which we then extended to the entire RTS development area. To account for the complex interactions between the ice-rich permafrost, water, and heat, we introduced the coupled equations of water and heat into our method. Through this approach and the estimated ground ice content in this region, we established a novel method for quantifying the volume of ground ice melting during the development of RTSs, which has significant implications for the sustainable development of ice-rich permafrost regions.

The schematic diagram of the method of estimating ground ice thawing was shown in Fig. 7. It is a numerical estimation method highly correlated with laboratory sample analysis and field survey. We defined the melting volume of permafrost through the variation of the calculated volume in different thawing seasons. Note that the volume of permafrost (V_p) was calculated by volume integral between the 0 °C isosurfaces and the depth of permafrost. The melting of ground ice, as simulated by the volume enclosed by the 0 °C isosurfaces of consecutive thawing seasons, is likely to have been overestimated as a result of the heterogeneity in unfrozen water moisture. So the melting amount was estimated using Eq. 12.

$$V_{a} = (V_{p,i + 1} - V_{p,i} ) \cdot v_{ice} \cdot (1 - W_{u} ),$$

(12)

where V_a is the melting of the ground ice (m³), V_p,i is the volume of the ice-rich permafrost in the i thawing season, and W_u is the volumetric unfrozen water content.

Fig. 7

Numerical estimation of the melting of ground ice of RTS

3.3 Numerical Modeling and Grid Generation

The geospatial information for the study area, which was acquired using an UAV, is presented in Fig. 4. The rectangular area delineates the plane domain of the 3D K1129W RTS numerical model. A full 3D physical model of the RTS was created using the modeling area, which has an overall slope angle of 8.6°. The DEM indicates that the elevation in the study area varies between 4590 and 4645 m a.s.l. However, for simulation purposes and to ensure compatibility with the running code, the elevation of the numerical model was scaled down to a range from 0 to 54 m. The model has dimensions of 218 m and 160 m on the X and Y axis, respectively, as shown in Fig. 8a.

Fig. 8

Numerical modeling of the K1129W retrogressive thaw slump (RTS). a Numerical meshing of the model, b 3D numerical model. The terrain of the domain is categorized as undisturbed area, headwall, scar area, and sliding mass. Four observation locations, namely points No. 1 through No. 4 in the RTS development zone, as well as longitudinal sections L1 to L3 (C–C′, A–A′, and B–B′) in the XZ direction, have been designated for analysis.

The model is composed of five layers, including mudstone, ice-poor soil, ice-rich soil, silt clay, and fine sands, arranged from bottom to top. The mesh is composed of 91 × 51 deformed hexahedral elements along the X-axis and Y-axis, respectively. Near the development site of RTS, finely tuned grids were devised with a 2 m × 2 m precision in the central region, spanning between 30 m and 200 m of the X-axis, and from 30 m to 130 m in Y-axis. And the Fig. 8b shows the meshes of the 3D model. The 3D model contains a total of 278,460 brick elements and 291,824 nodes. Laboratory experiments were performed to investigate the physical characteristics of each layer, such as its porosity, density, and thermal and hydraulic conductivity in the foundational layers. The outcomes of these tests are provided in Table 2.

Table 2 Parameters in the model

3.4 Initial and Boundary Conditions

Initial and boundary conditions for the hydro-thermal model are depicted in Fig. 9. The flow boundary and initial condition are shown in Fig. 9a. To simulate the groundwater flow-through in the active layer and recharge from ground ice meltwater, a specified hydraulic head following the topographic elevation was assigned on the top surface of the model (Kurylyk et al. 1969). Meanwhile, there was no run-off along the slope. Precipitation infiltration and replenishment were assigned to fluctuate based on the four seasons and vary spatially based on the vegetation type and surface coverage throughout the model area. The northern and southern boundaries of the model, located at a distance from the RTS area, were therefore deemed impermeable. The underside of the model was designated as impermeable too, due to its frozen state. The initial hydraulic heads were specified by the local terrain and assumed to be consistent with depth.

Fig. 9

Heat transfer boundary conditions imposed in the domain. a Flow boundary and initial condition, b boundary of heat exchange.

In regards to the boundary of heat exchange (Fig. 9b, where model depth is included for clarity), a heat flow of 0.032 W/m² was enforced at the lowermost of the model at the depth of 0 m (true elevation 4573 m a.s.l.). The left and right boundary conditions were specified as the adiabatic boundary, which is no heat exchange along the corresponding impermeable boundaries (Perreault et al. 2021). The model surface, representing ground surface, was categorized as natural ground, headwall, collapse scar, and sliding mass.

Through the linear regression equations shown in Fig. 10, the four kinds of surface temperatures were defined as the given temperature boundary conditions (Yin et al. 2021). The correlation between the temperature of the air and ground surface is:

$$T_{{\text{s}}} = {\text{F(}}T_{{{\text{air}}}} {)},$$

(13)

where T_s is ground surface temperature (°C), T_air is air temperature (°C), and F is the function of the different development terrains. The functions were obtained using a regression equation between the ground surface and air temperatures obtained from the field monitoring. It was postulated that the initial heat transfer conditions approximated the expected ground temperature and permafrost distribution, which were inferred from the geological background of the BLR region (Luo, Lin, et al. 2019).

Fig. 10

Regression functions of daily mean ground surface temperature and daily mean air temperature for different areas: a the natural area, b the recently disturbed area, c the scar area, and d the sliding mass. The linear regression coefficients and R² are shown in the figure. The solid line with numbers from 1 to 12 shows a function between ground and air temperatures.

The hydro-thermal results obtained at the conclusion of the equilibration period in 2011 served as the initial conditions for estimating ground ice ablation and assessing the hydrological influences of permafrost degradation in the RTS area from 2011 to 2021. The maximum simulation time steps for each scenario were set to 1 day. In the numerical simulations, to account for temperature-dependent physical properties, the hydraulic head and flow velocity solutions were updated at each time step. To achieve convergence, the convergence criteria for temperature and hydraulic head was set as 0.01 °C and 0.001 m, respectively.

4 Results and Discussion

This section initially presents the simulation outcomes of hydro-thermal state in the RTS development area, and comparison with the results of ERT survey. Subsequently, a quantitative evaluation is undertaken concerning the ground ice ablation and the consequential impact on the hydrological ecosystem during the process of RTS development.

4.1 Thermal State Analysis

The simulated temporal fluctuations of the thermal regime in the natural undisturbed area and the RTS area, as illustrated in Fig. 10, offer insights into the dynamic changes occurring in the active layer and distribution of warm permafrost. The maximum ground temperature differences in the thawing season were about 1.5 °C between the collapsed scar area and undisturbed terrain, and 0.8 °C between headwall and undisturbed terrain. From the 0 °C isosurface, the calculated permafrost table depth at the headwall varies from 1.7 m to 2.5 m between 2012 and 2018 (see Fig. 11). The cause of this variation may be the formation and expansion of tensional cracks behind the headwall. In 2018, the thaw penetration depth attained 202 cm in the headwall and 185 cm in the collapse scar area and it was even deeper than the adjacent natural area, while the simulated spatial distribution characteristic also clearly shows the lateral thermal erosion process of the RTS.

Fig. 11

Simulated two-dimensional ground temperature profiles of the K1129W retrogressive thaw slump (RTS) area, including a the natural area and affected areas of b the lower part of the headwall, c collapse scars, and d sliding masses. The indicated dates are the time when the maximum thawing depth occurred in the year

4.2 Moisture State Analysis

The formation and ablation of the ground ice is closely related to the soil moisture migration in the active layer, and the migration of unfrozen water is closely related to the ground temperature gradient (Cheng 1983). The variation of moisture content is a significant factor to influence the development and evolution of RTS. It changes the heat exchange condition between the atmosphere and surface and destabilizes the thermal state of the underlying permafrost (Niu et al. 2012).

The simulated unfrozen (volumetric) water contents from July to September and January to March are shown in Fig. 12, both in the natural undisturbed area (Fig. 12a) and the RTS area (Fig. 12b and c). For the undisturbed terrain, an obvious general increase of the moisture content of the active layer is displayed in Fig. 12a. In summer, the moisture content was about 15% at the top 10 cm of the active layer and gradually increased to 32% at the permafrost table. In winter, the moisture content was about 7% at the top 11 cm of the active layer and gradually increased to 13.5% at the permafrost table. The soil water migrated to the permafrost table in the natural area. The rate of downward freezing from the ground surface is faster than the upward freezing rate from the permafrost table, causing this phenomenon (Cheng 1983; French 2017).

Fig. 12

Simulated volumetric moisture content in the active and permafrost layer in the natural area and the retrogressive thaw slump (RTS) area. a The natural area, b the lower part of the headwall, and c the collapse scar area.

In the RTS area, the soil water saturation in the top layer of the collapse scar region exceeded 30%, which was notably higher than that in the natural area (see Fig. 12b and c). Analysis of the unfrozen water content at different depths and during various seasons revealed that the moisture content was lower at a depth of 0.7 m and higher at depths of 0.35 m and 0.95 m in both the headwall and the RTS area. Although water content that remained unfrozen in the upper 3 m permafrost exhibited a mostly stable pattern, with a minor increase from 2 to 3 m and a decrease from 3 to 5 m, the unfrozen water content near the permafrost table in the RTS area was significantly lower than that in the natural area due to the thawing of ground ice in the collapse scar area. Overall, there was no discernible pattern in the unfrozen water content in the RTS area. The results presented in Fig. 12b indicate that the scar area exhibited a superior capacity for downward migration of unfrozen water during the warm season. This phenomenon may be due to the different hydraulic conductivity. The crushed soil structure of the collapse scar area was easier for water migration. Therefore, since the surface conditions changed and tensional and shear crack expanded due to the RTS development, the rate of downward thawing was far faster than under the natural condition. Furthermore, depth of permafrost table in the RTS area was deeper than before and ground ice continued to melt.

4.3 Comparison Between the Simulated Ground Temperatures and the Electrical Resistivity Tomography (ERT) Survey Results

To further evaluate the actual active layer depth and the spatial distribution of the massive ground ice of the K1129W RTS area, three inversion results of the ERT survey are shown in Figs. 13a, 14a, and 15a. The geophysical investigation lines were from the bottom of the headwall to the sliding mass. The ERT inversion results show a heterogeneous distribution of electrical resistivity values in the RTS area. Due to the freeze-thaw erosion, the superficial layer of the RTS area has a large amount of mudflow materials with the ruined alpine meadow and several water puddles. Most of these electrical resistivities in the shallow strata were below 100 Ω-m.

Fig. 13

Spatial distribution of ice-rich permafrost in the K1129W retrogressive thaw slump (RTS) area revealed by the ERT investigation line A–A′ and corresponding simulation of longitudinal-section L2 (see Fig. 4 for detailed location). a The inversion results of the electrical resistivity tomography (ERT) survey line A–A′ and b the ground temperature field in section L2 on 20 September 2021 when the ERT survey was carried out. The ice-rich permafrost table is depicted by the 0 °C isotherm at this time.

Fig. 14

Spatial distribution of ice-rich permafrost in the K1129W retrogressive thaw slump (RTS) area revealed by the ERT investigation line B–B′ and corresponding simulation of longitudinal-section L3 (see Fig. 4 for detailed location). a The inversion results of the electrical resistivity tomography (ERT) survey line B–B′, and b the ground temperature field in section L3 on 20 September 2021 when the ERT survey was carried out. The ice-rich permafrost table is depicted by the 0 °C isotherm at this time.

Fig. 15

Spatial distribution of ice-rich permafrost in the K1129W retrogressive thaw slump (RTS) area revealed by the ERT investigation line C–C′ and corresponding simulation of longitudinal-section L1 (see Fig. 4 for detailed location). a The inversion results of the electrical resistivity tomography (ERT) survey line C–C′, and b the ground temperature field in section L1 on 20 September 2021 when the ERT survey was carried out. The ice-rich permafrost table is depicted by the 0 °C isotherm at this time.

As illustrated along the ERT investigation lines A–A′, B–B′, and C–C′, the thickness of the active layer was about 1.8 m to 2.3 m. The ice-rich permafrost beneath the active layer at a depth of 1.5 to 3 m within the RTS area was revealed through the use of Eq. 4, which identified its high electrical resistivity values of over 600 Ω-m even when the unfrozen moisture content was as high as 18%. The ERT model base was interpreted as mudstone with 200 to 220 Ω-m. The low electrical resistivity zone was found at the depth of 5 m to 15 m, interpreted as saline ice-poor permafrost.

To further determine the spatial distribution of ice-rich permafrost, the simulated geotemperature distribution corresponding to the ERT survey of the three XZ longitudinal-sections L1 to L3 (C–C′, A–A′, and B–B′) was extracted from the 3D model (Figs. 13b, 14b, and 15b). The location of L1 to L3 is indicated in Fig. 4. The simulated permafrost table is delineated by the isosurface of 0 °C in late September, which shows a reasonable match with the ERT results.

Based on the comparison, the massive ground ice was mainly buried under the headwall. The ice-rich permafrost within the displaced mass area was noticeably thicker than that within the collapse scar area, potentially as a result of a thickened layer in the former region. The 3D numerical result can reflect the permafrost table spatial distribution and the influence of topography on ice-rich permafrost. In total, the spatial distributions of the ice-rich permafrost provided the basis of the geometric information for the estimation of the ground ice ablation.

4.4 Estimation of Ground Ice Ablation

To assess the amount of ground ice that may have melted during the development of RTS, the corresponding parameter values, including thawing rate and depth of permafrost table, unfrozen water content, and ice content, need to be extracted from the simulation results of the annual thawing season. The simulated permafrost table was most readily discernible during the latter part of September and the beginning of October, coinciding with the complete thawing of the active layer and the subsequent emergence of a 0 °C isosurface around the distribution of ground ice. Compared to the simulated permafrost layer, the downward progression of thawing in the K1129W RTS region exhibited a temporal trend. Specifically, the RTS development of the region resulted in a progressive escalation of the average rate of permafrost thawing during the thawing season, and the 0 °C isosurface in the development area and the maximum depth of thawing have been observed to vary across different years. In the K1129W RTS area, ice-rich permafrost or massive ground ice thaws at an average rate of 45 cm (lower part of the headwall), 22 cm (the collapse scar area), and 5.5 cm (the displaced mass area) during one thawing season. The varying thawing rates across the RTS area cause irregularities in the permafrost table. The observed variations in the depth of the 0 °C isosurface and maximum depth of thawing provide insights into the changes in ice-rich permafrost or massive ground ice in the RTS area. The faster thawing rate in the lower part of the headwall compared to the collapse scar area and the displaced mass area suggests the nonuniformity of permafrost thawing. Overall, there is a strong agreement between the average thawing trend observed in the headwall and the trend predicted by the Stefan equation at a singular node, which can be characterized by an approximate negative power-law function (Fig. 16a). Furthermore, the distribution of ice-rich permafrost in the slope direction might be linked to topography-related hydrology and RTS dynamics. Additionally, the dynamics of RTS might be related to the distribution of ice-rich permafrost.

Fig. 16

Simulated depth of permafrost table in the headwall, estimated melting of ground ice, thawing index, and precipitation from 2011 to 2021. a Depth of permafrost table, b melting of ground ice, c thawing index, and d cumulative precipitation

The detailed ground ice ablation amount at the K1129W RTS area from 2016 to 2021 is shown in Fig. 16b. In the period from 2011 to 2021, the year with the highest melting amount was 2016, during which thawing index peaked at 807 °C day (as shown in Fig. 16c), and precipitation amounted to 421 mm (as indicated in Fig. 16d). According to Luo, Niu, et al. (2019) the significant triggering factor of the RTS, as illustrated in Fig. 16, was the elevated air temperature along with precipitation experienced during the thawing season of 2016. Meanwhile, calculation on ground ice ablation also indicated that, in each thawing season, the meltwater of ice-rich permafrost recharged on the permafrost table, which exacerbated the ground ice ablation. This reciprocal feedback process triggered means that the RTS is not a one-time landslide.

4.5 Hydrological Environment Impact

The development of RTS has been observed to intensify soil freeze-thaw erosion, which is further exacerbated by the collapse and retrogression on gentle slopes. The variation in groundwater turbidity or the concentration of total suspended materials is an indicator of the impact of the hydrogeological ecosystem surrounding the RTS area. The research conducted in this study aimed to examine the effects of RTS on soil stability and surface ponding depression water turbidity by analyzing the fine-grained migration linked to the extensive ground ice ablation on the K1129W RTS area.

Field investigations revealed that water color changes occurred in the runoff and stream or lake around the BLR region, which may be due to altered turbidity or suspended sediment concentration in the groundwater flow affected by RTS sediments accumulated in the surrounding floor plain (Lewkowicz and Way 2019). Several turbid small puddles and ponding depression were observed in the K1129W RTS area (Fig. 2). The hydro-thermal condition change caused by RTS evolution produces more meltwater formed by ground ice ablation. The thawed sediment and meltwater slide or flow down a slope, which is the source area of the sediments in the surrounding stream or lake. The melting of ground ice caused the transport of suspended sediments. To quantify the impact on a hydrosystem, a variable called relative concentration of total suspended sediments (RTSS, c/c₀) was introduced in our model. c₀ is the initial sediment concentration in the ground ice, and c is the simulated sediment concentration at different times. Thaw-induced mass wasting was transported in the active layer with meltwater, and the process was driven by the ground ice ablation under the headwall. The RCSS in the RTS area was characterized by different yield values with different seepage flow velocities (Fig. 17). The two patterns of the breakthrough curves reflected the heterogeneity of layers and the dynamic of sediment migration. These suspended particles of the icy permafrost also carried sorbed phase materials such as viruses, dissolved organic carbon (DOC), and salt ions trapped in ground ice.

Fig. 17

Breakthrough curves of suspended sediments for fine sand and clay medium at the typical flow velocities. The number marked on the curve is the corresponding flow velocity.

Based on the relationship between Nephelometric Turbidity Units (NTU) and RCSS (Felix et al. 2018; Liu et al. 2021), the temporal fluctuation of turbidity in ponding depression in one thawing season period presented that the value remained elevated in the dry period (Fig. 18). It is mainly due to the migration of suspended sediments triggered by the ground ice melt. The permafrost thawing rate was at a relatively high level in early June, as displayed in Fig. 18, and it was the main driving force of the sediment transport. Moreover, the turbidity fluctuation was correlated with the meteorological data in August. The jumping turbidity value events appeared on 29 June and 10 August, following a high air temperature of 10.1 °C and daily precipitation of 19.1 mm. Recharge from the meltwater below the headwall increased the hydraulic connectivity that facilitated the soil erosion in the RTS area. Furthermore, degraded vegetation cover provided conditions for mass wasting including suspended sediment transport, which was a nonnegligible impact on the surrounding hydrological ecosystem.

Fig. 18

Simulated a thawing depth in the retrogressive thaw slump (RTS) area and b turbidity in the ponding depression of the sliding mass area, c daily air temperature, and d daily precipitation in the Beilu River (BLR) region in 2021

4.6 Deficiencies in the Research

The current method for quantitatively assessing ground ice ablation volume in a RTS development area has several deficiencies. First, the method relies heavily on sample analysis results and ERT surveys, particularly for the complex 3D numerical model. This means that the accuracy of the estimates is limited by the quality and quantity of the data collected through these methods. Second, the method is time-consuming and complex. The use of a 3D numerical model requires significant computational power and resources, making it challenging to apply to large-scale studies. As a result, it is not efficient to estimate the amount of ground ice melting caused by the numerous active RTSs over the Tibetan Plateau, which numbers at 2669 (Luo et al. 2022).

Despite these limitations, this method presents a novel idea for estimating the melting volume of ground ice in a RTS development region. It provides a useful framework for studying the estimation of ground ice ablation from the developmental mechanism of RTSs. However, the current approach is not practical for large-scale studies, which require more efficient and accurate methods. To address these deficiencies, future research should focus on developing a more practical theoretical model of ground ice ablation estimation based on the developmental mechanism of RTSs and remote sensing data. By leveraging the vast amounts of available remote sensing data, researchers can develop more efficient and accurate methods for estimating ground ice ablation volumes in RTS development areas. This would significantly improve our understanding of the impact of climate change on the Tibetan Plateau and other similar regions, where RTS development is a significant environmental concern.

5 Conclusion

Retrogressive thaw slumps are widely developed in the permafrost regions and their development and evolution markedly influence the cryo-hydrogeological ecosystem. To quantitatively estimate ground ice ablation volume and evaluate the impact on the hydrological ecosystem of RTS, we developed a 3D hydro-thermal coupled numerical model of a RTS on the permafrost terrain at the BLR region of the QTP. Using a hydro-thermal model that was calibrated and supplemented by data from an ERT survey and sample analysis, we developed a new approach for quantitatively assessing the ground ice ablation. Our findings can be summarized as follows:

• (1) The variation of permafrost table was affected by the uneven moisture migration in the active layer. The simulated rate of thawing was found to be significantly higher in the RTS area compared to the undisturbed area. This disparity led to consequential changes in the hydro-thermal regime in close proximity to the permafrost table.

• (2) The maximum melting of ground ice occurred in 2016, with a volume of 492 m³. This confirms that climatic warming and wetting are closely linked to the ground ice ablation, specifically in the regions located beneath the headwall, which can trigger RTS development.

• (3) The freeze-thaw erosion triggered suspended sediment transport and the hydrological system impact of the surrounding area. High ponding depression water turbidity values of 28 and 49 NTU occurred on the front edge of the RTS in the thawing season of 2021 when the increment of air temperature and precipitation was up to the threshold amount.