A method for predicting thaw subsidence in permafrost

Abstract

This invention discloses a method for predicting the melting settlement of ice-rich permafrost, relating to the fields of geotechnical engineering and disaster prevention and mitigation in cold regions. The method includes: constructing a layered parameter model of the permafrost; vectorizing the temperature field based on Fourier's law and the sensible heat capacity method; and automatically identifying the soil type based on an initial void ratio threshold. For ordinary permafrost, a large-strain consolidation model is used to calculate settlement; for ice-rich layers, a volumetric collapse model is used to directly convert the melted ice volume into settlement. The method combines Lagrange dynamic mesh technology to update node coordinates and physical parameters in real time, achieving bidirectional coupling of thermodynamic parameters, and introduces adaptive time step control to improve computational efficiency. This invention, through the coupling of consolidation and collapse mechanisms, solves the problems of inaccurate settlement prediction and low computational efficiency in traditional methods for high-ice-content permafrost. It can accurately simulate the settlement deformation of complex strata such as ice wedges and ice lenses under long-term thermal action, and is suitable for stability assessment and disaster prevention in cold region engineering.

Description

Technical Field

[0001] This application relates to the fields of geotechnical engineering and disaster prevention and mitigation in cold regions, and in particular to a method for predicting the melting and settlement of ice-rich permafrost. Background Technology

[0002] In geotechnical engineering construction and disaster prevention and mitigation technologies in cold regions, permafrost, as a geological medium extremely sensitive to temperature, directly affects the safety of infrastructure due to its thermal stability. With global warming and intensified human engineering activities, permafrost has undergone significant degradation, primarily manifested in rising ground temperatures, increased active layer thickness, and the melting of underground ice. Particularly for ice-rich permafrost widely distributed in regions such as the Qinghai-Tibet Plateau, its melting is often accompanied by severe surface subsidence and deformation, leading to serious engineering disasters such as roadbed cracking, suspended pipelines, and tilted buildings. Therefore, scientifically and accurately predicting the subsidence and deformation of permafrost during the thawing process is of paramount engineering significance for the design, maintenance, and disaster prevention and mitigation of projects in cold regions.

[0003] Currently, the latest existing technology for predicting frozen soil settlement typically employs a thermo-hydraulic-mechanical (THM) coupled numerical simulation method based on thaw consolidation theory. This technique first calculates the temperature field distribution of the soil using the Fourier heat conduction equation to determine the locations of the freezing and thawing fronts. Subsequently, for the thawed soil region, it is treated as a saturated porous medium, assuming that the soil skeleton always exists and is continuous. Then, using consolidation theories (such as Terzaghi or Biot consolidation theories), the dissipation process of excess pore water pressure in the thawed soil and the effective stress increase of the soil skeleton under its own weight or external loads are calculated. Finally, the consolidation deformation of the soil layer is calculated using the soil's volumetric compressibility coefficient, and the surface settlement is obtained by summing these values.

[0004] However, before melting, the volume of pure ice layers or soils with extremely high ice content is mainly supported by ice particles. Once a phase change occurs and the ice melts, it transforms into water and is expelled. The original solid support structure experiences instantaneous volume loss and structural collapse, rather than simply compression of the soil skeleton. This leads to the current technology still forcibly applying consolidation theory to calculate compression deformation when facing high-ice-content frozen soil (ice-rich frozen soil) or pure ice layers (such as underground ice wedges or buried ice), ignoring the direct volume loss caused by phase change. As a result, the predicted settlement is far less than the actual collapse (i.e., severely distorted), failing to truly reflect the catastrophic characteristics of ice-rich strata, thus posing significant safety hazards to engineering designs. Summary of the Invention

[0005] Therefore, it is necessary to provide a method for predicting the melting and settlement of ice-rich permafrost to address the aforementioned technical problems.

[0006] The present invention adopts the following technical solution:

[0007] This invention provides a method for predicting the thawing settlement of ice-rich permafrost, comprising:

[0008] A one-dimensional stratified model of the ground in the cold region was constructed, and the one-dimensional stratified model was discretized into several grid cells along the depth direction. Each grid cell was then divided into ordinary permafrost region and ice-rich region.

[0009] Based on a one-dimensional layered model, the temperature value of each grid cell is determined; based on the temperature value of each grid cell, the grid cells that have transitioned from a frozen state to a thawing state are identified; based on the region type of the grid cells that have transitioned from a frozen state to a thawing state, the settlement calculation mode corresponding to the region type of the grid cell is determined; the settlement calculation mode includes: if the grid cell belongs to a common frozen soil region, a thawing consolidation mode is adopted, the overlying self-weight stress of the grid cell is determined based on the soil gravity of all grid cells above the grid cell, and converted into the excess pore water pressure of the grid cell and substituted into the consolidation equation for solution to obtain the consolidation settlement of the grid cell; if the grid cell belongs to an ice-rich region, a volume collapse mode is adopted, and the thickness of the melted ice within the grid cell is determined as the collapse settlement of the grid cell.

[0010] The total ground settlement in the cold region is obtained by summing the consolidation settlement and collapse settlement of all grid cells.

[0011] Preferably, the one-dimensional layered model is constructed based on geological survey data of the cold region surface; the geological survey data includes: dry density, weight water content, thermal conductivity, volume compressibility coefficient and permeability coefficient of the cold region surface soil layer; the dry density and weight water content are used to calculate the initial void ratio and the overlying self-weight stress.

[0012] Preferably, each grid cell is divided into a common permafrost region and an ice-rich region, specifically including:

[0013] Obtain the initial porosity of each grid cell;

[0014] When the initial porosity of the grid cells is greater than the preset pure ice threshold, the grid node region is determined as an ice-rich region.

[0015] When the initial void ratio of the grid cells is less than or equal to the preset pure ice threshold, the grid node region is defined as a normal permafrost region.

[0016] Preferably, the temperature value of each grid cell is determined based on a one-dimensional hierarchical model, specifically including:

[0017] The porosity, liquid water content, and specific heat capacity and thermal conductivity of soil particles, ice, liquid water and air of each grid cell are obtained, and the effective thermal conductivity and sensible heat capacity of each grid cell are calculated; wherein, the sensible heat capacity includes: the sensible heat capacity of the soil skeleton and the latent heat contribution within the phase change temperature range.

[0018] A vectorized explicit finite difference equation is constructed. Using the temperature vector of the full profile at the current moment, combined with the effective thermal conductivity and sensible heat capacity, the temperature value of each grid cell of the full profile at the next moment is calculated in one go.

[0019] Preferably, based on the temperature value of each grid cell, the grid cells that transition from a frozen state to a thawed state are determined, specifically including:

[0020] A phase transition temperature range is set; the phase transition temperature range is between the solidus temperature and the liquidus temperature of the target site;

[0021] Based on the temperature value of each grid cell, determine the liquid water content fraction of each grid cell, including:

[0022] If the temperature value is greater than the liquidus temperature, then the liquid water content fraction of the grid cell is determined to be 1;

[0023] If the temperature value is less than the solidus temperature, the liquid water content fraction of the grid cell is determined to be 0.

[0024] If the temperature value is within the phase transition temperature range, the liquid water content fraction of the grid cell changes linearly with the temperature value between 0 and 1.

[0025] Compare the liquid water content fraction of the current time step with that of the previous time step, and determine the grid cells whose liquid water content fraction changes from 0 to greater than 0 as the grid cells that have changed from frozen to thawed in the current time step.

[0026] Preferably, if the grid cell belongs to a common frozen soil region, a thawing consolidation mode is adopted. The overlying self-weight stress of the grid cell is determined based on the soil gravity of all grid cells above it, and then converted into the excess pore water pressure of the grid cell. This excess pore water pressure is substituted into the consolidation equation for solution, yielding the consolidation settlement of the grid cell. Specifically, this includes:

[0027] Calculate the total weight of the soil layer in all grid cells above the grid cell that has melted, and use it as the overlying self-weight stress of the grid cell;

[0028] The superposition of the overlying self-weight stress onto the pore water pressure of the current grid cell yields the excess pore water pressure.

[0029] The consolidation coefficient is calculated based on the permeability and volume compressibility coefficient of the grid cells. Based on the consolidation coefficient and excess pore water pressure, the one-dimensional consolidation equation is solved using the finite difference method to obtain the amount of excess pore water pressure dissipation in the grid cells within the current time step.

[0030] Multiplying the dissipation of excess pore water pressure in a grid cell by its volumetric compressibility coefficient yields the volumetric strain of the grid cell, which is the consolidation settlement of the grid cell within the current time step.

[0031] Preferably, if the grid cell belongs to an ice-rich region, a volumetric collapse mode is adopted, and the thickness of the melted ice within the grid cell is determined as the collapse settlement of the grid cell, specifically including:

[0032] The monitoring area type is the liquid water content fraction of grid cells in ice-rich areas;

[0033] If the liquid water content is greater than 0, it is determined that the ice in the grid cell has undergone phase change and melting.

[0034] Calculate the ice melt thickness of the grid cell within the current time step, and use the ice melt thickness as the collapse settlement of the grid cell;

[0035] Preferably, after the collapse settlement of the grid cell is determined, the method further includes: resetting the pore water pressure of the grid cell to zero to prevent the pore water pressure of the grid cell in the current time step from entering the calculation of the melting and consolidation mode in the next time step.

[0036] Preferably, the total ground subsidence in cold regions is calculated using an iterative method; the process of determining the time step for each iteration in the iterative calculation specifically includes:

[0037] Traverse all grid cells of the one-dimensional layered model to obtain the thermal diffusivity and consolidation coefficient of each grid cell;

[0038] By selecting the mesh element with the maximum thermal diffusivity and combining it with the vertical dimension of the mesh element, the thermal critical time step that satisfies the stability of the explicit differential solution of the heat conduction equation is calculated.

[0039] By selecting the mesh element of the maximum consolidation equation and combining the vertical dimension of the mesh element, the mechanical critical time step that satisfies the stability of the explicit differential solution of the consolidation equation is calculated.

[0040] The thermal critical time step is compared with the mechanical critical time step, and the minimum value is selected as the simulation time step for the next iteration.

[0041] Preferably, the method further includes: updating the vertical coordinates and physical property parameters of each grid cell in real time using the Lagrange method based on the consolidation settlement and collapse settlement corresponding to each grid cell, specifically including:

[0042] Using the Lagrange method, based on the consolidation settlement and collapse settlement of the corresponding grid cells, the vertical thickness of the corresponding grid cells is reduced, and the depth coordinates of each grid cell are redefined.

[0043] The porosity of each grid cell is updated based on the volumetric strain generated by the settlement of each grid cell; wherein the porosity decreases as the settlement increases.

[0044] Based on the updated porosity, the permeability and thermal conductivity of each grid cell are recalculated using the empirical formula Kozeny-Carman equation or the geometric mean method, and the updated physical property parameters are substituted into the thermo-mechanical coupling calculation for the next time step.

[0045] The above-mentioned at least one technical solution adopted in this invention can achieve the following beneficial effects:

[0046] This invention employs differentiated settlement calculation models for different regional types. For ordinary permafrost regions, the thaw consolidation model is used, converting the overlying self-weight stress into excess pore water pressure and substituting it into the consolidation equation to solve for consolidation settlement, thus describing the progressive deformation dominated by soil skeleton compression. For ice-rich regions, a volume collapse model is adopted, directly using the thickness of the melted ice as the collapse settlement to characterize the instantaneous volume loss caused by the disappearance of the ice-supported structure. This overcomes the limitations of existing thermo-hydraulic-mechanical coupling methods that forcibly apply consolidation theory and ignore the phase change collapse of high-ice-content strata. Through regional differentiation calculation driven by physical mechanisms, it retains the theoretical rationality of ordinary permafrost consolidation settlement while accurately capturing the catastrophic characteristics of structural collapse during the melting of ice-rich or pure ice layers. This achieves high-precision prediction of the total surface settlement during the entire permafrost thawing process, effectively avoiding engineering design safety hazards caused by inaccurate settlement predictions, and significantly improving the reliability and scientific nature of risk assessment for cold-region engineering sites.

[0047] In addition, this judgment mechanism, combined with adaptive time step control and Lagrange grid update, significantly improves the efficiency of long-duration calculations while ensuring simulation accuracy. It enables high-precision and high-efficiency prediction of melting and subsidence of complex permafrost strata in cold-region engineering, and is especially suitable for long-term stability assessment and disaster prevention under geological conditions containing ice wedges, ice lenses and other geological conditions. Attached Figure Description

[0048] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings:

[0049] Figure 1 A flowchart illustrating a method for predicting the melting settlement of ice-rich permafrost provided by this invention;

[0050] Figure 2 The execution flowchart of a method for predicting the melting settlement of ice-rich permafrost provided by the present invention;

[0051] Figure 3 A schematic diagram of the physical model of pure ice layer melting collapse and ordinary frozen soil consolidation in a method for predicting melting settlement of ice-rich permafrost provided by the present invention.

[0052] Figure 4 A comparison and verification diagram of the prediction results and measured data of the method for predicting the melting settlement of ice-rich permafrost provided by the present invention;

[0053] Figure 5 A comparison of prediction results from different stratigraphic models for a method to predict the melting and settlement of ice-rich permafrost provided by this invention. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments in the specification without creative effort are within the scope of protection of this application.

[0055] The technical solutions provided by the various embodiments of this application are described in detail below with reference to the accompanying drawings.

[0056] Figure 1 This is a schematic diagram of a method for predicting the melting and settlement of ice-rich permafrost in this invention, which specifically includes the following steps:

[0057] S101: Construct a one-dimensional stratified model of the ground in the cold region, and discretize the one-dimensional stratified model into several grid units along the depth direction, and divide each grid unit into ordinary permafrost region and ice-rich region.

[0058] Optionally, the one-dimensional layered model is constructed based on geological survey data of cold-region ground; the geological survey data includes: dry density, weight water content, thermal conductivity, volume compressibility and permeability of the soil layer; the dry density and weight water content are used to calculate the initial void ratio and the overlying self-weight stress.

[0059] Specifically, geological survey data refers to the set of basic parameters characterizing the vertical soil layer distribution and physical state of the target site, obtained through on-site drilling, sampling, and laboratory tests; the "one-dimensional layered model" simplifies the target site into a columnar structure model with attribute changes only along the depth direction, and its layering is based on the soil property abrupt interface (such as the boundary between clay and gravel layers); "discretization into several grid units along the depth direction" means dividing the one-dimensional model into multiple computational units of equal or varying thickness according to a preset vertical dimension, with each grid unit having independent and homogeneous initial physical parameters; this discretization process provides a spatial carrier for subsequent thermo-mechanical coupling calculations, and its grid scale must meet the basic requirements of the Courant-Friedrichs-Lewy (CFL) stability condition for spatial resolution, but no specific numerical limitations are introduced.

[0060] "Dry density" refers to the mass of solid particles per unit volume of soil. It is a fundamental physical parameter reflecting the compactness of soil and is used in this embodiment for the joint calculation of initial void ratio and overlying self-weight stress. "Gravity water content" refers to the ratio of the mass of water to the mass of solid particles in the soil. It is a key indicator characterizing the initial moisture state of the soil and, in this embodiment, together with dry density, forms the basis for calculating the initial void ratio and further supports the quantitative expression of overlying self-weight stress. "Thermal conductivity" refers to the ability of a material to transfer heat under a unit temperature gradient. It is a core parameter in Fourier's law of heat conduction that determines the rate of heat diffusion along the depth direction and is used in this embodiment to construct a temperature field evolution model. "Volume compressibility coefficient" refers to the rate of change of volumetric strain in soil under a unit effective stress increment. Compressibility is a key conversion parameter in large-strain consolidation theory that transforms excess pore water pressure into consolidation settlement. The permeability coefficient refers to the ability of soil to allow water to pass through under a unit hydraulic gradient. It is a core mechanical parameter that controls the dissipation rate of excess pore water pressure and affects the time scale of the consolidation process. This application simultaneously collects and inputs the above five types of parameters, so that the parameter assignments of the one-dimensional layered model in the two physical processes of heat conduction and consolidation deformation have clear sources and physical traceability, avoiding the uncertainty brought about by empirical values.

[0061] The initial void ratio (PVR) is the ratio of pore volume to solid particle volume in the initial state of soil. It is a core criterion for determining whether a mesh element belongs to an ice-rich region. In this embodiment, the PVR is derived from dry density and gravimetric water content using the following relationship: First, the solid volume ratio is calculated based on dry density and soil particle density. Then, the gravimetric water content is converted to volumetric water content, ultimately yielding the PVR ratio. This derivation process does not rely on the assumption of a saturated state and is suitable for modeling the initial configuration of unsaturated frozen soil. "Overburden effective stress" refers to the effective vertical stress generated vertically downwards by the weight of all soil layers above a certain depth at that depth surface. Stress is the direct load source for the excess pore water pressure formed in the melting consolidation mode. In this embodiment, the overlying self-weight stress is determined by the unit volume density of each layer based on its thickness, dry density, and weight water content, and then accumulated layer by layer. The stress calculation result is directly used as the input basis for the initial boundary conditions in the subsequent consolidation equation. This application realizes the dual function reuse of geological input parameters in thermo-mechanical coupling modeling by simultaneously using dry density and weight water content for initial porosity determination and overlying self-weight stress calculation, thereby enhancing the internal consistency and physical rationality of the model parameter system.

[0062] Specifically, the initial void ratio can be calculated as follows: based on the dry density ρ_d and the soil particle density ρ_s, first calculate the solid volume fraction V_s = ρ_d / ρ_s, then use the weight water content w to convert the water mass into water volume, combine the total mass conservation relationship of the soil to deduce the pore volume V_v, and finally obtain the initial void ratio according to e_0 = V_v / V_s;

[0063] The calculation method for overlying self-weight stress can be as follows: For all soil layers above the depth of each grid cell, calculate the unit weight of each layer in sequence as γ_i = ρ_di × g × (1 + w_i), where g is the gravitational acceleration, ρ_di is the dry density of the i-th grid cell, and w_i is the solid mass of the i-th grid cell. Then multiply the unit weight of each layer by its thickness h_i and sum them up to obtain the overlying self-weight stress σ_z = Σ(γ_i × h_i) at that depth.

[0064] This application uses dry density and weight water content as common input sources for the initial void ratio and overlying self-weight stress, so that the criteria for classifying ordinary permafrost areas and the mechanical driving conditions for thawing and consolidation settlement are consistent at the parameter source. On this basis, the differentiated calculation paths for the initial void ratio and overlying self-weight stress of different implementation methods not only meet the needs of multi-source data adaptation in engineering sites, but also ensure the robustness of the model in scenarios with missing or uncertain parameters.

[0065] Optionally, based on a preset initial void ratio threshold, the grid node regions in the one-dimensional layered model of the target site are divided into ordinary permafrost regions and ice-rich regions. Specifically, this includes: obtaining the initial void ratio of each grid unit; when the initial void ratio of the grid unit is greater than the preset pure ice threshold, the grid node region is determined as an ice-rich region; when the initial void ratio of the grid unit is less than or equal to the preset pure ice threshold, the grid node region is determined as an ordinary permafrost region.

[0066] Specifically, the initial void ratio is a dimensionless parameter characterizing the initial looseness and density of soil, calculated from dry density and weight water content. Its value directly reflects the spatial potential of the soil to accommodate ice. The preset initial void ratio threshold is an empirical critical value used to determine whether the soil has ice-rich characteristics. When the initial void ratio of the grid cell is higher than this threshold, it indicates that it has the typical physical properties of ice-rich layers with high ice storage and low skeleton support strength. This classification result forms the premise for subsequent differentiated settlement modeling—that is, under the same temperature change, the two types of regions will trigger completely different mechanical response mechanisms, rather than being uniformly described by consolidation theory.

[0067] In this context, ordinary permafrost regions refer to permafrost areas with low initial porosity, good continuity of soil particle skeletons, and ice mainly existing in cemented or infill forms. These regions typically possess considerable compressive strength and drainage consolidation capacity within their respective technical fields. In this embodiment, the determination result is used to activate the thawing-consolidation mode, ensuring reasonable modeling of the evolution of excess pore water pressure and consolidation deformation processes induced by overlying self-weight stress. The initial porosity threshold is used to define whether the soil possesses significant "ice-like" structural characteristics, namely high porosity, low skeleton support capacity, and strong phase change volumetric response. In this embodiment, the pure ice threshold serves as a critical point for classification. When the initial porosity of a grid cell exceeds this value, it indicates that its ice content is close to the level of pure ice, and the soil particle skeleton is largely encased or interrupted by ice. After melting, structural collapse is likely to occur, thus classifying it as an ice-rich region to trigger the subsequent settlement calculation path for the volumetric collapse mode.

[0068] S102: Based on a one-dimensional layered model, determine the temperature value of each grid cell; based on the temperature value of each grid cell, determine the grid cells that have changed from a frozen state to a thawing state; based on the region type of the grid cells that have changed from a frozen state to a thawing state, determine the settlement calculation mode corresponding to the region type of the grid cell; the settlement calculation mode includes: if the grid cell belongs to a normal frozen soil region, the thawing consolidation mode is adopted, the overlying self-weight stress of the grid cell is determined based on the soil gravity of all grid cells above the grid cell, and converted into the excess pore water pressure of the grid cell and substituted into the consolidation equation for solution to obtain the consolidation settlement of the grid cell; if the grid cell belongs to an ice-rich region, the volume collapse mode is adopted, and the thickness value of the melted ice in the grid cell is determined as the collapse settlement of the grid cell.

[0069] Optionally, based on a one-dimensional layered model, the temperature value of each grid cell is determined, specifically including: obtaining the porosity, liquid water content, and thermal parameters of each component of each grid cell; calculating the effective thermal conductivity and sensible heat capacity of each grid cell; wherein, the sensible heat capacity includes: the sensible heat capacity of the soil skeleton and the latent heat contribution within the phase change temperature range; constructing a vectorized explicit finite difference equation; using the temperature vector of the full profile at the current moment, combined with the effective thermal conductivity and sensible heat capacity, calculating the temperature value of each grid cell of the full profile at the next moment in one go.

[0070] Specifically, Fourier's law of heat conduction is a fundamental physical law describing the migration of heat in a medium by conduction, and its differential form is: ,in Where is the thermal diffusivity, Mesh temperature function T Depth coordinates of mesh cells z temperature gradient, t Let be the heat conduction time in the medium, where T Depth coordinates z and time t bivariate function abbreviated as T The sensible heat capacity method is a numerical processing strategy that incorporates the latent heat contribution of phase change into the heat capacity term, enabling the temperature field solution to reasonably reflect the energy absorption / release characteristics during the ice-water phase change process without introducing additional phase change interface tracking. The combination of these two methods forms the core of the thermal sub-model, whose output is the temperature value of each grid cell in the full profile at the next time step. The grid cell that changes from the frozen state to the molten state refers to the cell that completes the phase change initiation process within the current time step. The criterion for this is the state change caused by the liquid water content fraction jumping from 0 to greater than 0 when the temperature value crosses the phase change interval. This state change is the driving signal that triggers the subsequent sedimentation calculation mode switching.

[0071] Porosity refers to the proportion of pores in a unit volume of soil, used to characterize the spatial capacity of soil to hold water and air; liquid water content refers to the mass or volume of liquid water in a unit volume of soil, used to reflect the state of unfrozen water in frozen soil during heating; thermal parameters of each component include the specific heat capacity and thermal conductivity of soil particles, ice, liquid water, and air, used to support the modeling of the thermal properties of multiphase media; effective thermal conductivity is a physical quantity that comprehensively considers the coexistence state and spatial distribution characteristics of the solid phase (soil particles, ice), liquid phase (water), and gas phase (air) in the soil, and is equivalent to characterize the overall thermal conductivity of the grid cell, its value dynamically changes with the phase composition ratio; sensible heat capacity is the temperature increase of a unit mass of soil by 1 The heat required for K to absorb consists of two parts: first, the sensible heat capacity of the soil skeleton (containing mineral particles and untransformed ice crystals), which represents the thermal energy storage capacity under linear temperature changes; second, the latent heat contribution equivalent to the ice-water phase transition within the phase transition temperature range. This part transforms the latent heat of phase transition, which originally needed to be handled at the moving boundary, into a continuously distributed increase in specific heat capacity within the temperature domain, thereby avoiding explicit tracking of the phase transition front. In this embodiment, the effective thermal conductivity is used to construct the discrete expression of the heat conduction term, directly affecting the heat flux density driven by the temperature gradient. Calculation accuracy; sensible heat capacity is used as a thermal inertia parameter in the construction of the time term of the energy balance equation. Its equivalent embedding of the latent heat of phase change makes the same grid cell exhibit a significantly increased equivalent specific heat capacity in the temperature range from the solidus to the liquidus, thus truly reflecting the inhibitory effect of the phase change process on the heating rate; based on the Voigt-Reuss mixing model, the effective thermal conductivity is calculated according to the weighted average of the volume components of each phase; and the equivalent specific heat capacity method is used to divide the latent heat of phase change by the width of the phase change temperature range and then superimpose it on the sensible heat capacity of the soil skeleton to obtain the sensible heat capacity.

[0072] "Vectorization" refers to organizing the temperature variables of all grid cells across the entire profile into column vectors, constructing sparse coefficient matrices for the heat conduction and heat capacity terms, and expressing the entire heat conduction control equation in matrix-vector form. "Explicit finite difference equations" refers to a discrete format constructed using forward differencing for the time derivative and central differencing for the spatial derivative. Its solution does not require iteratively solving a system of linear equations; it can directly deduce the temperature field at the next moment based solely on the current temperature field. In this embodiment, the explicit vectorized format relies on a defined thermal critical time step constraint to ensure numerical stability. Its inputs are the current full-profile temperature vector, the effective thermal conductivity vector of each grid cell, and the sensible heat capacity vector; the output is the next full-profile temperature vector. This processing method skips the traditional point-by-point loop update logic, supporting a single call to the underlying BLAS / LAPACK library functions to complete the entire column temperature update, significantly improving computational throughput efficiency. The construction and solution method involves converting the one-dimensional heat conduction partial differential equation... Discretize into Form, and thus solve ,in, For grid cell density, For effective specific heat capacity, Let C be the effective thermal conductivity, C be the diagonal sensible heat capacity mass matrix, and K be the tridiagonal effective thermal stiffness matrix. For time derivative, It is the spatial derivative.

[0073] The discretization process of the one-dimensional heat conduction partial differential equation is as follows:

[0074] The target region is spatially discretized into N grid cells, with a discrete grid step size of... Δz spatial derivative Using central difference, for the i-th grid node, the spatial derivative term is discretized as:

[0075] ;

[0076] Effective thermal conductivity based on discrete grid and grid step size Δz Obtain the effective thermal conductivity matrix of the tridiagonal K Non-zero elements in;

[0077] Based on the tridiagonal effective thermal conductivity stiffness matrix K ， Rearrange the above expression into matrix form:

[0078] ;

[0079] Using forward difference Decomposed into discrete formats based on time step:

[0080] ;

[0081] On each discrete grid, Multiply by grid step size Δz The diagonal sensible heat capacity mass matrix was then obtained. C ;

[0082] At any moment At that time, based on the diagonal sensible heat capacity mass matrix C and the effective thermal conductivity matrix of the tridiagonal K , the formula Rewritten as:

[0083] .

[0084] Based on the temperature values ​​of each grid cell, the grid cells that transition from a frozen state to a molten state are determined. Specifically, this includes: setting a phase transition temperature range; the phase transition temperature range lies between the solidus temperature and the liquidus temperature of the target site; determining the liquid water content fraction of each grid cell based on its temperature value, including: if the temperature value is greater than the liquidus temperature, the liquid water content fraction of the grid cell is set to 1; if the temperature value is less than the solidus temperature, the liquid water content fraction of the grid cell is set to 0; if the temperature value is within the phase transition temperature range, the liquid water content fraction of the grid cell changes linearly between 0 and 1 with the temperature value; comparing the liquid water content fraction of the current time step with that of the previous time step, and determining the grid cells whose liquid water content fraction changes from 0 to greater than 0 as those that transition from a frozen state to a molten state within the current time step.

[0085] Specifically, the "liquid water content fraction" is a dimensionless parameter characterizing the proportion of unfrozen liquid water in the total pore water mass within a grid cell, with a value range of [0,1]. In the thermodynamics of multiphase media in frozen soil, this parameter directly reflects the current phase composition and latent heat release / absorption state of the soil. In this embodiment, this parameter serves as a core intermediate variable for judging the phase change process: when its value is 0, it indicates that the grid cell is in a completely frozen state; when its value is 1, it indicates that it has completely thawed; when its value is between 0 and 1, it indicates that it is in a partially thawed state, and its value can reflect the degree of thawing.

[0086] Optionally, when the grid node region is a common frozen soil region, the overlying load is converted into excess pore water pressure, and the consolidation equation is solved to obtain the settlement of the common frozen soil region. Specifically, this includes: calculating the total weight of all soil layers above the thawing grid unit as the overlying self-weight stress of the grid unit; superimposing the overlying self-weight stress onto the pore water pressure of the current grid unit to obtain the excess pore water pressure; calculating the consolidation coefficient based on the permeability coefficient and volume compressibility coefficient of the grid unit, and solving the one-dimensional consolidation equation using the finite difference method based on the consolidation coefficient and excess pore water pressure to obtain the dissipation of excess pore water pressure of the grid unit within the current time step; multiplying the dissipation of excess pore water pressure of the grid unit by the volume compressibility coefficient to obtain the volumetric strain of the grid unit, which is the consolidation settlement of the grid unit within the current time step.

[0087] Optionally, if the grid cell belongs to an ice-rich region, a volume collapse mode is adopted, and the thickness of the melted ice within the grid cell is determined as the settlement of the grid cell. Specifically, this includes: the liquid water content fraction of the grid cell in the monitoring area type of an ice-rich region; if the liquid water content fraction is greater than 0, it is determined that the ice within the grid cell has undergone phase change melting; the ice melting thickness of the grid cell in the current time step is calculated, and the ice melting thickness is used as the collapse settlement of the grid cell.

[0088] Specifically, the thaw consolidation model is a mechanical response model established for permafrost regions with complete soil skeleton support capacity. Its essence is to incorporate the effective stress redistribution process caused by thawing into the classical consolidation theory framework: the overlying self-weight stress is applied as an external load to the top surface of the thawed layer, which is transferred through pore water to form excess pore water pressure, and then gradually dissipates through the synergistic effect of infiltration drainage and soil compression, ultimately manifesting as compressible deformation. The volume collapse model is a structural instability model established for ice-rich regions. Its core is to acknowledge that under high porosity conditions, the soil skeleton cannot maintain its original configuration after the ice melts, resulting in irreversible local volume shrinkage. This shrinkage can be directly quantified by the thickness of the melted ice without undergoing a drainage consolidation process. The two models correspond to different physical mechanisms, and their selection depends entirely on the previously determined region type, without relying on real-time judgment of intermediate variables such as temperature and water content, ensuring logical closed loop and engineering interpretability.

[0089] Specifically, the overlying self-weight stress is a vertical effective stress component in soil mechanics that characterizes the weight of the overlying soil layer on the underlying unit. Its physical meaning is the resultant force of gravity per unit area. In this embodiment, this stress serves as the initial mechanical excitation source driving the discharge of pore water and the compression of the soil skeleton. Its value is obtained by accumulating the dry density, volume, and gravitational acceleration of each grid unit above, without introducing additional parameters or model assumptions. The calculation method for the overlying self-weight stress can be: based on the thickness, dry density, and weight water content of each grid unit in the one-dimensional layered model, the saturated unit weight of each layer is calculated, and its gravity contribution is accumulated layer by layer along the depth direction. Excess pore water pressure is the additional pore water pressure exceeding the hydrostatic pressure as defined in the Terzaghi effective stress principle. Its generation mechanism is that after an external load (here, the overlying self-weight stress) is suddenly applied to saturated soil, due to the limited permeability of the soil, water cannot be discharged in time, resulting in pore water bearing part of the load. In this embodiment, this pressure constitutes the driving force for the subsequent consolidation process, directly determining the drainage rate and compression response intensity. This superposition operation can be achieved by algebraically adding the overlying self-weight stress value directly to the current pore water pressure value to form the initial excess pore water pressure distribution.

[0090] Specifically, the permeability coefficient is a physical parameter characterizing the soil's ability to allow water to flow through it, measured in m / s. A higher value indicates faster drainage. The volumetric compressibility coefficient is a parameter reflecting the rate of change of volumetric strain in soil under a unit effective stress increment, measured in MPa. -1 The higher the value, the more easily the soil is compressed; both the coefficient of compaction and the consolidation coefficient together determine the "consolidation coefficient," which is... ,in Permeability coefficient, The volume compressibility factor is 1. The specific weight of water; this relationship is a well-known definition of the consolidation coefficient in the art, without introducing a new formula derivation; in this embodiment, the consolidation coefficient is used to quantify the consolidation response rate of the mesh element in the current state; the one-dimensional consolidation equation refers to a partial differential equation with time as the independent variable and depth as the single spatial coordinate, and its standard form is... ,in The excess pore water pressure, The equation describes the evolution of excess pore water pressure over time and space. In this embodiment, the equation is solved by an explicit finite difference scheme, which has the advantages of high computational efficiency and easy embedding into a coupled framework. The solution process of the finite difference method can be as follows: the current grid cell and its upper and lower adjacent cells form a three-point template, the central difference approximates the second-order spatial derivative, the forward difference approximates the first-order time derivative, an explicit iterative scheme is constructed and the excess pore water pressure is updated.

[0091] Specifically, the dissipation of excess pore water pressure refers to the decrease in excess pore water pressure of the grid cell within the current time step, and its value is equal to the difference between the initial excess pore water pressure and the excess pore water pressure at the next moment; this value directly reflects the degree of water expulsion; volumetric strain is the volume shrinkage rate of a unit volume of soil during consolidation, and its physical meaning is the compressive deformation response of the soil skeleton under the action of effective stress growth; in this embodiment, this strain is expressed through a linear constitutive relation. Obtain, among which This refers to the dissipation of excess pore water pressure. The volumetric compressibility coefficient is derived from the fundamental assumptions of classical Terzaghi consolidation theory and does not introduce nonlinear corrections. Consolidation settlement is the vertical displacement of the grid cell due to drainage compression within the current time step, and its value is equal to the volumetric strain multiplied by the current vertical thickness of the grid cell. This thickness has been updated by the Lagrange method in the previous time step and is a dynamically evolving parameter. The latest updated volumetric compressibility coefficient of the current grid cell is directly multiplied by the excess pore water pressure dissipation calculated within the current time step.

[0092] Specifically, the ice melt thickness refers to the equivalent vertical thickness of the ice layer reduced by the grid cell due to phase change within the current time step. Its physical meaning is: the equivalent vertical scale obtained by converting the mass of melted ice within the cell to volume using the pure ice density and then dividing by its horizontal cross-sectional area. This parameter directly characterizes the geometric spatial loss caused by ice phase loss. In this embodiment, this thickness is equivalently mapped to the instantaneous settlement displacement of the grid cell in the vertical direction, without introducing consolidation time effects or drainage path assumptions, thus truly reflecting the catastrophic characteristic of ice-rich layers "melting and collapsing." The method for calculating the ice melt thickness can be: based on the difference in liquid water content fraction between the current time step and the previous time step, multiplied by the initial ice-containing thickness of the grid cell, the ice melt thickness within the current step is obtained.

[0093] This application uses the liquid water content fraction as a criterion for the melting state of ice-rich layers, and leverages the logic mechanism that a value greater than 0 triggers a collapse response to accurately capture the structural instability behavior of high-ice-bearing strata. A schematic diagram of the physical model for pure ice layer melting collapse and ordinary permafrost consolidation can be found in [reference needed]. Figure 3 Based on this, the thickness of the ice melt is directly defined as the collapse settlement, skipping the complex process of excess pore water pressure dissipation and stress redistribution in the traditional consolidation model, so that the settlement calculation and phase variables are strictly conserved. Finally, the collapse settlement and the consolidation settlement in the ordinary permafrost area are added together to form a complete physical expression of the total surface settlement, thus uniformly characterizing the two dominant mechanisms in the melting and settlement process of ice-rich permafrost under the thermo-mechanical coupling framework—progressive consolidation and sudden collapse.

[0094] Additionally, after the collapse settlement of the grid cell is determined, the process also includes: resetting the pore water pressure of the grid cell to zero to prevent the pore water pressure of the grid cell in the current time step from entering the calculation of the melting and consolidation mode in the next time step. See the overall execution flowchart. Figure 2 .

[0095] The pore water pressure of a grid cell refers to the hydrostatic pressure exerted by the pore water within that cell, reflecting the mechanical effect of undissipated water in the soil pores on the soil skeleton. Resetting to zero is a forced assignment operation, indicating that after completing the settlement calculation corresponding to the volume collapse mode, the pore water pressure value stored in the current grid cell is actively cleared, preventing it from participating in any subsequent mechanical response calculations driven by excess pore water pressure. "Preventing the pore water pressure value stored at the current moment from entering the calculation of the melting consolidation mode in the next time step" reflects the mechanism isolation logic—in ice-rich areas, the high initial porosity leads to structural collapse after ice melts. The water discharge process is instantaneous and irreversible, failing to meet the slow seepage and pressure diffusion prerequisites upon which consolidation theory relies. Therefore, the input path of the pore water pressure to subsequent consolidation equations is blocked. The reset operation can be executed immediately after the collapse settlement calculation is completed, updating the pore water pressure variable of the corresponding grid cell to 0 through a direct assignment statement.

[0096] S103: Sum the consolidation settlement and collapse settlement of all grid cells to obtain the total ground settlement in the cold region.

[0097] Optionally, the total ground subsidence in cold regions is calculated using an iterative calculation method. During the iterative calculation, the process of determining the time step for each iteration specifically includes: traversing all grid cells of the one-dimensional layered model to obtain the thermal diffusivity and consolidation coefficient of each grid cell; selecting the grid cell with the maximum thermal diffusivity and, in conjunction with the vertical dimensions of the grid cell, calculating the thermal critical time step that satisfies the stability of the explicit differential solution of the heat conduction equation; selecting the grid cell with the maximum consolidation equation and, in conjunction with the vertical dimensions of the grid cell, calculating the mechanical critical time step that satisfies the stability of the explicit differential solution of the consolidation equation; comparing the thermal critical time step with the mechanical critical time step, and selecting the minimum value as the simulation time step for the next iteration.

[0098] The "thermal diffusivity" is a physical quantity that characterizes a material's ability to diffuse heat under a unit temperature gradient. It is defined as the ratio of thermal conductivity to volumetric heat capacity. ,in For effective thermal conductivity, is the heat capacity per unit volume; this parameter reflects the rate of heat propagation in the soil and directly determines the upper limit of the time step in explicit heat conduction numerical simulations; the "thermal diffusivity coefficient" is used to quantify the instantaneous response of each grid cell to temperature disturbances, and its maximum value corresponds to the region with the most intense thermal response, constituting the dominant constraint for the numerical stability of the heat conduction subsystem; the vertical dimension of the grid cell, i.e., the grid step size, refers to the thickness of a single grid cell obtained by dividing along the depth direction, and is the basic scale for spatial discretization of the explicit difference scheme, its value directly affecting the discretization error and stability condition of the heat diffusivity term; the "thermal critical time step" is the theoretically maximum allowable time step to ensure the numerical stability of the explicit difference scheme of the heat conduction equation, and its calculation is based on the Courant-Friedrichs-Lewy (CFL) type stability criterion, expressed as ,in For grid step size, This is the maximum thermal diffusivity across the entire profile; this time step ensures that the temperature field update process will not cause false oscillations or divergences due to excessively large time step sizes.

[0099] The consolidation coefficient is a key parameter describing the coupling relationship between the pore water pressure dissipation rate and volumetric deformation in saturated soil under load. It is defined as follows: ,in Permeability coefficient, The volume compressibility factor is 1. is the unit weight of water; this parameter reflects the time-scale characteristics of the soil drainage consolidation process; in this embodiment, the consolidation coefficient is only activated after the grid element undergoes phase change melting, used to characterize the mechanical response hysteresis of the element after it changes from a frozen state to a saturated softened state; its maximum value appears in the region where meltwater is abundant and the skeleton is most significantly weakened, constituting the dominant constraint for the numerical stability of the consolidation subsystem; the "mechanical critical time step" is the theoretical maximum allowable time step to ensure the numerical stability of the explicit difference scheme of the one-dimensional consolidation equation, and its expression is as follows: ， in This is the maximum consolidation coefficient among all currently molten mesh elements; this time step ensures that the coupled evolution process of excess pore water pressure dissipation and settlement deformation remains numerically convergent.

[0100] Optionally, based on the consolidation settlement and collapse settlement of each grid cell, the vertical coordinates and physical property parameters of each grid cell are updated in real time using the Lagrangian method. Specifically, this includes using the Lagrangian method to reduce the vertical thickness of the corresponding grid cell based on the consolidation settlement and collapse settlement, and redetermining the depth coordinates of each grid cell; updating the porosity of each grid cell based on the volumetric strain generated by the settlement of each grid cell; wherein the porosity decreases with the increase of settlement; based on the updated porosity, recalculating the permeability and thermal conductivity of each grid cell using the empirical formula Kozeny-Carman equation or the geometric mean method, and substituting the updated physical property parameters into the thermo-mechanical coupling calculation of the next time step.

[0101] Specifically, the Lagrange method is a numerical modeling method that tracks material points. In this scheme, it is manifested as follows: each grid cell is considered a physical entity carrying its own mass and properties. When it settles, it not only changes its geometric coordinates but also simultaneously causes changes in physical properties such as a decrease in porosity, permeability, and thermal conductivity. This update process constitutes a key link in the thermo-mechanical bidirectional coupling. On the one hand, the coordinate change reconstructs the heat conduction path and boundary, affecting the subsequent temperature field distribution; on the other hand, the change in physical parameters is fed back into the calculation of the thermal diffusivity and consolidation coefficient in the next time step, giving the entire simulation process deformation-driven nonlinear evolution capabilities. Entering the calculation of the total surface settlement in the next time step means that the method continues to advance in a time-iterative manner until it covers the design prediction cycle. The entire process does not rely on empirical correction coefficients, and all outputs are direct derivations of physical equations and material constitutive relations. See [link to relevant documentation]. Figure 4 and Figure 5 These are, respectively, a comparison and verification diagram of settlement prediction results and measured data, and a comparison diagram of prediction results from different strata models.

[0102] Wherein, volumetric strain refers to the volume change per unit initial volume, and its value is equal to the ratio of the total settlement to the initial thickness of the corresponding grid cell; void ratio is defined as the ratio of pore volume to solid particle volume in the soil; in this embodiment, volumetric strain serves as the direct driving force for void ratio updates, and its mechanism is as follows: when the soil undergoes compressive settlement, solid particles are compacted, pore volume decreases, thereby leading to a decrease in void ratio; this relationship conforms to the basic physical constraints of soil mechanics, and in this embodiment, it is manifested as a monotonically decreasing function relationship; specifically, the updated void ratio... It can be represented as ,in The volumetric strain caused by settlement within the current time step. The porosity is the porosity at the previous time step; this expression is only used to illustrate the qualitative relationship between porosity and volumetric strain, and does not constitute a limitation on the specific mathematical form.

[0103] The process of updating the vertical coordinates using the Lagrange method is as follows: First, the consolidation settlement and collapse settlement of each grid cell are linearly superimposed to obtain the total settlement of the grid cell; then, the initial thickness of the grid cell is scaled proportionally, and the depth coordinates of its upper and lower interfaces are adjusted sequentially to ensure that the interfaces of adjacent grid cells are continuous and do not overlap.

[0104] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this invention.

Claims

1. A method for predicting settlement due to melting of ice-rich permafrost, characterized in that, Includes the following steps: A one-dimensional stratified model of the ground in the cold region was constructed, and the one-dimensional stratified model was discretized into several grid cells along the depth direction. Each grid cell was then divided into ordinary permafrost region and ice-rich region. Based on a one-dimensional layered model, the temperature value of each grid cell is determined; based on the temperature value of each grid cell, the grid cells that have changed from a frozen state to a thawing state are determined; based on the region type of the grid cells that have changed from a frozen state to a thawing state, the settlement calculation mode corresponding to the region type of the grid cell is determined. The settlement calculation mode includes: if the grid cell belongs to a normal frozen soil area, the thaw consolidation mode is adopted, the overlying self-weight stress of the grid cell is determined according to the soil gravity of all grid cells above the grid cell, and converted into the excess pore water pressure of the grid cell and substituted into the consolidation equation for solution to obtain the consolidation settlement of the grid cell; if the grid cell belongs to an ice-rich area, the volume collapse mode is adopted, and the thickness of the melted ice in the grid cell is determined as the collapse settlement of the grid cell. The total ground settlement in the cold region is obtained by summing the consolidation settlement and collapse settlement of all grid cells.

2. The method for predicting the thawing settlement of ice-rich permafrost as described in claim 1, characterized in that, The one-dimensional layered model is constructed based on geological survey data of cold regions; the geological survey data includes: dry density, weight water content, thermal conductivity, volume compressibility and permeability of the soil layer; the dry density and weight water content are used to calculate the initial void ratio and the overlying self-weight stress.

3. The method for predicting settlement of thawing permafrost as described in claim 1, characterized in that, The division of each grid cell into ordinary permafrost regions and ice-rich regions specifically includes: Obtain the initial porosity of each grid cell; When the initial porosity of the grid cells is greater than the preset pure ice threshold, the grid node region is determined as an ice-rich region. When the initial void ratio of the grid cells is less than or equal to the preset pure ice threshold, the grid node region is defined as a normal permafrost region.

4. The method for predicting settlement of ice-rich permafrost thawing as described in claim 1, characterized in that, The determination of the temperature value of each grid cell based on the one-dimensional hierarchical model specifically includes: The porosity, liquid water content, and specific heat capacity and thermal conductivity of soil particles, ice, liquid water and air of each grid cell are obtained, and the effective thermal conductivity and sensible heat capacity of each grid cell are calculated; wherein, the sensible heat capacity includes: the sensible heat capacity of the soil skeleton and the latent heat contribution within the phase change temperature range. A vectorized explicit finite difference equation is constructed. Using the temperature vector of the full profile at the current moment, combined with the effective thermal conductivity and sensible heat capacity, the temperature value of each grid cell of the full profile at the next moment is calculated in one go.

5. The method for predicting settlement of thawing permafrost as described in claim 1, characterized in that, The step of determining the grid cells that transition from the frozen state to the thawed state based on the temperature value of each grid cell specifically includes: A phase transition temperature range is set; the phase transition temperature range is between the solidus temperature and the liquidus temperature of the target site; Based on the temperature value of each grid cell, determine the liquid water content fraction of each grid cell, including: If the temperature value is greater than the liquidus temperature, then the liquid water content fraction of the grid cell is determined to be 1; If the temperature value is less than the solidus temperature, the liquid water content fraction of the grid cell is determined to be 0. If the temperature value is within the phase transition temperature range, the liquid water content fraction of the grid cell changes linearly with the temperature value between 0 and 1. Compare the liquid water content fraction of the current time step with that of the previous time step, and determine the grid cells whose liquid water content fraction changes from 0 to greater than 0 as the grid cells that have changed from frozen to thawed in the current time step.

6. The method for predicting settlement of ice-rich permafrost as described in claim 1, characterized in that, If the grid cell belongs to a common frozen soil region, a thawing consolidation mode is adopted. The overlying self-weight stress of the grid cell is determined based on the soil gravity of all grid cells above it, and then converted into the excess pore water pressure of the grid cell. This excess pore water pressure is substituted into the consolidation equation for solution, yielding the consolidation settlement of the grid cell. Specifically, this includes: Calculate the total weight of the soil layer in all grid cells above the grid cell that has melted, and use it as the overlying self-weight stress of the grid cell; The superposition of the overlying self-weight stress onto the pore water pressure of the current grid cell yields the excess pore water pressure. The consolidation coefficient is calculated based on the permeability and volume compressibility coefficient of the grid cells. Based on the consolidation coefficient and excess pore water pressure, the one-dimensional consolidation equation is solved using the finite difference method to obtain the amount of excess pore water pressure dissipated in the grid cells within the current time step. Multiplying the dissipation of excess pore water pressure in a grid cell by its volumetric compressibility coefficient yields the volumetric strain of the grid cell, which is the consolidation settlement of the grid cell within the current time step.

7. The method for predicting settlement of thawing permafrost as described in claim 1, characterized in that, If the grid cell belongs to an ice-rich region, a volumetric collapse mode is adopted, and the thickness of the melted ice within the grid cell is determined as the collapse settlement of the grid cell, specifically including: The monitoring area type is the liquid water content fraction of grid cells in ice-rich areas; If the liquid water content is greater than 0, it is determined that the ice in the grid cell has undergone phase change and melting. Calculate the ice melt thickness of the grid cell within the current time step, and use the ice melt thickness as the collapse settlement of that grid cell.

8. The method for predicting settlement of ice-rich permafrost as described in claim 1, characterized in that, After the collapse settlement of the grid cell is determined, the method further includes: resetting the pore water pressure of the grid cell to zero to prevent the pore water pressure of the grid cell in the current time step from entering the calculation of the melting and consolidation mode in the next time step.

9. The method for predicting settlement of thawing permafrost as described in claim 1, characterized in that, The calculation method for the total ground subsidence in the cold region is iterative calculation; the process of determining the time step used in each iteration during the iterative calculation specifically includes: Traverse all grid cells of the one-dimensional layered model to obtain the thermal diffusivity and consolidation coefficient of each grid cell; By selecting the mesh element with the maximum thermal diffusivity and combining it with the vertical dimension of the mesh element, the thermal critical time step that satisfies the stability of the explicit differential solution of the heat conduction equation is calculated. By selecting the mesh element of the maximum consolidation equation and combining the vertical dimension of the mesh element, the mechanical critical time step that satisfies the stability of the explicit differential solution of the consolidation equation is calculated. The thermal critical time step is compared with the mechanical critical time step, and the minimum value is selected as the simulation time step for the next iteration.

10. The method for predicting the thawing settlement of ice-rich permafrost as described in claim 1, characterized in that, Also includes: Based on the consolidation settlement and collapse settlement of each grid cell, the vertical coordinates and physical property parameters of each grid cell are updated in real time using the Lagrange method, specifically including: Using the Lagrange method, based on the consolidation settlement and collapse settlement of the corresponding grid cells, the vertical thickness of the corresponding grid cells is reduced, and the depth coordinates of each grid cell are redefined. The porosity of each grid cell is updated based on the volumetric strain generated by the settlement of each grid cell; wherein the porosity decreases as the settlement increases. Based on the updated porosity, the permeability and thermal conductivity of each grid cell are recalculated using the empirical formula Kozeny-Carman equation or the geometric mean method, and the updated physical property parameters are substituted into the thermo-mechanical coupling calculation for the next time step.

Images