Deep rock masses are subjected to high in situ stress and also characterized by a significant presence of discontinuities such as joints and faults. Many studies have focused on the blast-induced fracture in intact rock under in situ stress or in shallow discontinuous rock masses without in situ stress. The blasting fracture behavior of deep rock masses under the combined actions of in situ stress and discontinuities remains insufficiently understood. In this study, the split Hopkinson pressure bar tests under confining pressure are conducted to validate the Riedel-Hiermaier-Thoma (RHT) model parameters of granite. A series of dynamic FEM simulations with the RHT model are performed to investigate the explosion stress wave propagation and crack distribution of single-hole blasting under various in situ stress and discontinuity conditions. Based on numerical modeling, the mechanism of blast-induced fracture in a rock mass with discontinuity subjected to in situ stress is exposited, and the blasthole arrangement for utilizing the discontinuity to facilitate rock fragmentation by blasting is discussed. The results show that blasting in discontinuous rock masses mainly generates three types of cracks, namely radial cracks emanating from the blasthole, spalling cracks resulting from explosion stress wave reflection on the discontinuity surface, and wing cracks derived from the discontinuity ends. With increasing in situ stress, the radial cracks and wing cracks are significantly reduced; particularly, wing cracks completely disappear when the in situ stress reaches 30 MPa. The spalling cracks are the primary mechanism for the discontinuity to influence the blast-induced rock fracture under in situ stress. The area of the cracks generated between the blasthole and the discontinuity expands with the increase of the discontinuity length. For maximizing the facilitation of the discontinuity on the blast-induced rock facture, the blasthole is recommended to be arranged on the midvertical line of the discontinuity. The blasthole-discontinuity distance should be carefully designed to avoid excessive rock fragmentation near the blasthole and ensure the formation of the effective spalling cracks. Under high in situ stress, a smaller blasthole-discontinuity distance is advisable.

Progressive rock damage may develop before visible surface deformation, limiting the effectiveness of deformation-based monitoring alone. This study evaluates electro-mechanical impedance (EMI) sensing for tracking cyclic damage in sandstone, marble, and granite under uniaxial compression. Surface-bonded piezoelectric patches acquired conductance spectra, while digital image correlation (DIC) measured full-field deformation on the opposite specimen surface. A normalized cumulative energy dissipation-based damage indicator derived from cyclic stress-strain curves served as an independent mechanical reference. With increasing loading cycles, the principal conductance resonance generally shifted to lower frequencies and changed in amplitude, reflecting variations in the effective stiffness, damping, and coupling of the PZT-adhesive-rock system. Root mean square deviation (RMSD), correlation coefficient deviation (CCD), and mean absolute percentage deviation (MAPD) increased with cumulative damage and showed strong correlations with the energy dissipation-based damage indicator. Lithology influenced both mechanical and EMI responses. Porous sandstone exhibited progressive compaction and energy dissipation, whereas granite showed limited inelastic deformation followed by abrupt brittle failure. Under the present opposite-side configuration, EMI-based damage indicators changed progressively before clearly connected DIC strain-localization bands appeared in some granite specimens. EMI and DIC therefore provide complementary information. EMI is sensitive to local mechanical-impedance changes around the bonded sensor, whereas DIC resolves surface strain localization. The results support EMI as a complementary local sensing method for engineering-geological damage monitoring, although validation under confining stress, environmental variation, fractured rock masses, multi-sensor configurations, and field conditions remains necessary.

Rainfall infiltration significantly affects slope stability.Accurately describing the dynamic behavior of the rainfall infiltration boundary is crucial for analyzing seepage and stability of slopes.Although the unidirectional con-version from flow-controlled to pressure-controlled boundaries has been achieved,several challenges remain.These include the unclear mechanisms in the governing equation and the lack of bidirectional conversion of boundary condi-tions,especially regarding treatment of the rainfall infiltration boundary in fracture-pore media.To this end,based on Darcy's seepage theory,this study reformulates the pressure governing equation,and establishes a rainfall infiltra-tion governing equation that can achieve bidirectional conversion between flow-controlled and pressure-controlled boundaries.Then,rainfall infiltration boundary treatment methods are proposed for discrete fracture-pore media and dual-permeability media.Three challenges are addressed:(1)a position head gradient is introduced to eliminate the influence of semi-permeable layer thickness on the depth of surface water accumulation;(2)a bidirectional conver-sion mechanism with two judgment conditions is developed to achieve adaptive control of the rainfall infiltration boundary;and(3)a rainfall infiltration boundary treating method is proposed to reveal the influence of rainfall inten-sity on the water exchange between the dual media.The results show that the proposed method can effectively allevi-ate the errors induced by the conventional rigid,one-way conversion between rainfall infiltration boundaries.Further-more,fractures are shown to play a significant role in the slope rainfall infiltration process.Under the low-intensity rainfall condition,the matrix dominates water transport,while the fractures become the primary seepage channels under the high-intensity rainfall.These findings provide theoretical support for studying slope fracture seepage and failure mechanism.

Earth dams using Jiangxi lateritic soil have been widely constructed in China. The water retention behavior of lateritic soil varied when subjected to reservoir fluctuations, influencing the safety and stability of earth dam. For this reason, the water retention behavior of Jiangxi lateritic soil and its microscopic mechanism under drying–wetting (D–W) cycles were investigated. The numbers of D–W cycles N = 1, 3 and 5 were considered. At the macro-scale, the filter paper method was adopted for the soil–water retention curve (SWRC) measurements. At the micro-scale, the mercury intrusion porosimetry and scanning electron microscope tests were conducted for microstructural observation. The results showed that under a given N value, the double S-shaped SWRC of lateritic soil was exhibited, with two air entry values (AEV1 and AEV2) identified. This was accompanied by the double peak hysteresis of SWRC at these two AEVS. The bi-modal pore size distribution (PSD) curve was observed, which was separated into four ranges. It was found that the double S-shaped SWRC was divided into four zones, which were dependent on the corresponding four ranges of PSD curve, respectively. Meanwhile, the d1 (the dominant macro-pores size), dinf (pore diameter at the inflection point) and d2 (the dominant micro-pores size) in PSD curve were consistent with AEV1, ψinf (the suction at the inflection point) and AEV2 in SWRC by Young–Laplace equation, respectively. With increasing N, the decrease in water retention capacity of four zones of SWRC was explained by the corresponding increase in pore size of four ranges.

Accurately inferring the joint distributions of geotechnical parameters is essential for geometric structure modelling and reliability assessment of geotechnical structures. In engineering practice, only sparse data can be acquired for a specific site, posing a challenge in modelling the joint distributions of geotechnical parameters. To address this challenge, this study proposes a Bootstrap-enhanced Bayesian Updating with Structural reliability (BBUS) method for the probabilistic characterisation of correlated geotechnical parameters with sparse data. The proposed method employs the parametric bootstrap technique to reconstruct the likelihood function, enabling the efficient generation of robust posterior samples. These posterior samples are subsequently used to infer the posterior predictive distributions of correlated geotechnical parameters through random sampling. The effectiveness of the proposed method is validated through two case studies. The results demonstrate that the proposed method can accurately model the joint distributions of correlated geotechnical parameters with sparse data, yielding posterior predictive distributions that closely match the observed data. Notably, the proposed method is versatile enough to be extended to correlated non-normal parameters and high-dimensional parameter scenarios. These findings also confirm the potential of the proposed BBUS method for the probabilistic characterisation of correlated non-normal geotechnical parameters with sparse data.

Due to the structurally control effect by bedding planes, the deformation and failure of layered rocks during tunnel excavation involve complex mechanical processes, including fracturing, bending and block spalling. When simulating the large deformation and its discontinuous behavior of this type of rock mass, the calculations based on the implicit numerical method have certain limitations in terms of convergence. To accurately reveal the progressive failure mechanisms under excavation, this study developed a numerical model for simulating quasi-static excavation process through Python and Fortran, which is suitable for the explicit cohesive zone model and finite discrete element method (CZM-FDEM). The influence of cohesive element stiffness on in-situ stress equilibrium is analyzed, and numerical verification is conducted through typical opening problems. The results indicate that the embedding of cohesive elements reduces the overall stiffness of the model and affects the accuracy of in-situ stress equilibrium. When the ratio of cohesive element stiffness Knto the elastic modulus of solid elements Es exceeds 108 and 560, the deviations in stress equilibrium fall within 8 % and 2 %, respectively. Quasi-static simulation of tunnel excavation is achieved through the dynamic relaxation method combined with explicit solution strategy. The results show high consistency with implicit methods and the Kirsch analytical solution. Regarding excavation responses of surrounding rocks, the effects of stratification features and tunnel section shape on the failure mechanisms of surrounding rock are systematically analyzed. Due to the obstruction to stress transmission and deformation continuity, the bedding planes significantly control the failure mode. As the layer thickness increases, the anisotropy of the surrounding rock gradually weakens, which results in the excavation damaged zone (EDZ) distribution approaching that of homogeneous rock. Meanwhile, the overall failure trend in surrounding rocks approximately parallel to the bedding planes with different bedding inclinations. The arched tunnel is more prone to stress concentration at the corner and arch shoulder. In this case, the side wall forms typical V-shaped failure. The explicit CZM-FDEM approach proposed in this study can provide a reliable means for stability assessment of layered tunnel engineering.

Rocks release subtle geochemical warning signals before breaking. These signals, coming from naturally occurring nuclides (e.g., radon, helium, argon, and thoron), have often been reported before earthquakes, volcanic eruptions, landslides, and rock and ice avalanches. However, despite their high sensitivity to deformation, their detectability, as well as myriad promising observations over half a century, nuclide signals are still far from being applied to geohazard prediction or widely used for monitoring. Here, we first develop a decomposition and interpretation method for nuclide signals. By analyzing nuclide signal time series observed from a month-long laboratory rock failure experiment and year-long slope deformation in a field setting, we identify a universal paradigm unit of nuclide signal evolution. We find that this paradigm unit is characterized by two core characteristics: a transient pulse and equilibrium fluctuation which are intrinsically correlated to rupture area and crack aperture, respectively. Through analytical derivation and pore-scale simulations, we establish the constitutive equations that link these characteristic nuclide signals to key rupture structural parameters. Rooted in these constitutive relations, we further develop a diagnostic theory of rock rupture via nuclide signals. We apply the model to track rock failures at the laboratory and field scale. The proposed nuclide signal decomposition and rupturing model enable the unification of discrete signal units emitted by individual microrupturing events, with the integrated signal evolution observed during macroscopic failure. This integration may serve as a foundation for both the mesoscopic assessment of rock damage and the early warning of geohazards induced by rock ruptures.

Understanding fluid flow and heat transfer in fractured rock masses is essential for advancing geothermal energy recovery. Fluid flow and heat transfer in single fractures, considering mechanical effects such as normal and shear displacements, have been extensively studied, yet the behavior in deformable intersected fractures remains underexplored. This study demonstrates, for the first time, heat transfer characteristics in intersected fractures under varying mechanical conditions, focusing on the effect of shear displacement. We first develop a computational model to simulate the shear behavior of intersected fractures under varying mechanical boundary conditions, validated against laboratory shear tests conducted in the present study. We then simulate fluid flow and heat transfer processes in the intersected fractures after shear. The results show that increasing shear displacement shifts the fracture intersection, leading to a significant rise in the cumulative energy proportion at the outlet of main flow fracture-from 63.1% to 96.5% as shear displacement increases from 0.5 mm to 5 mm. Furthermore, overlooking the effect of Constant Normal Stiffness (CNS) conditions in geothermal simulations of deep fractured rock masses could lead to an overestimation of production energy. These findings enhance the understanding of heat transfer behavior in natural fractured rock masses, contributing to more efficient geothermal energy extraction.

Understanding contact characteristics in rock fractures is essential for predicting fluid flow and heat transfer in geothermal and other subsurface energy systems. This study systematically investigates how contact obstacles, characterised by contact area ratio (c) and spatial distribution mode (uniform versus normal), govern coupled thermo‑hydraulic processes in single fractures. A three‑dimensional parallel‑plate fracture model incorporating discrete contact elements is developed to isolate contact effects while eliminating geometric confounders such as roughness and aperture variability. Finite‑element simulations are performed for contact ratios ranging from c= 0% to 80% under different contact distributions. The results show that increasing contact area ratio markedly reduces hydraulic aperture and flow channel connectivity, leading to pronounced changes in flow regime and heat transfer behaviour. A critical threshold at c≈ 40% is identified, beyond which flow transitions from channelised to highly tortuous patterns, significantly delaying thermal breakthrough. Moreover, contact distribution exerts a strong control on flow path geometry: uniform contact distributions promote relatively straight preferential pathways and more efficient advective heat transport, whereas normal distributions induce higher tortuosity and enhanced conductive heat exchange with the fracture walls. Based on these mechanisms, a modified analytical solution for temperature evolution is proposed by introducing a hydraulic aperture formulation that explicitly accounts for both contact ratio and flow tortuosity. The analytical predictions show excellent agreement with numerical results. These findings provide mechanistic insights into contact‑controlled thermo‑hydraulic behaviour in fractured rock and offer general implications for the design and optimisation of geothermal and other geo‑energy extraction systems.

Accurately modelling the interactions between continuous and discontinuous materials is essential for advancing engineering solutions across a wide range of fields. Owing to fundamental differences in the governing equations of discontinuous and continuous numerical models, distinct solution schemes have been developed, presenting challenges for their coupling. This paper proposes a unified scheme for modelling continuous and discontinuous materials, integrating them monolithically by adopting a variational framework for an implicit Discrete Element Method (DEM) and the Finite Element Method (FEM). The seamlessly coupled FEM-DEM model is formulated as a convex optimisation problem, which can be transformed into a standard second-order cone program (SOCP) and solved efficiently using modern optimisation algorithms. The proposed approach has been validated through a series of numerical examples. Additionally, a numerical model test on granular flow against an elastic barrier has been conducted, demonstrating the scheme's capability in modelling the impact dynamics of granular flows and the resulting deformation and stress distribution in structures, which has significant implications for engineering design involving granular material handling.

Sandstone-mudstone interbedded rock masses are widely distributed in reservoir areas of southwest China. As the reservoir water level fluctuates periodically, the mechanical properties of the interface between sandstone and mudstone tend to deteriorate, thereby undermining the stability of interbedded rock masses on bank slopes. This study examined the macroscopic shear performance and microscopic deterioration mechanism of the sandstone-mudstone interface subjected to cyclic drying-wetting treatments via direct shear tests, mineral composition analysis, and microstructure observation. The test findings indicated that both the peak shear strength and the friction angle of the interface decreased in a negative exponential manner as the number of drying-wetting cycles increased. The shear stiffness increased linearly as the normal stress rose during the drying-wetting treatments. Meanwhile, the microscopic investigation of the interface after shear tests showed that the increase in normal stress and the water-induced softening of minerals would enhance the interfacial wear effect. Significantly, most scratched mineral grains were scraped from the mudstone walls, which indicated that the water-caused degradation of the shear resistance of the interface between sandstone and mudstone was mainly caused by the shrinkage and expansion of the mudstone wall. Based on the evolution of shear stress, a shear damage constitutive model for the sandstone-mudstone interface was developed by drawing on the Saint-Venant model and statistical damage mechanics. The proposed model considered the important factors, including the water-induced damage and the applied normal stress, which was validated by test results and proved to have a satisfactory predictive resolution.

The study of the shear behavior of bonded rock-cement interface is important for understanding the strength and stability of grouted rock masses. This research aims to reveal the failure mechanism behind the shear property of bonded rock-cement interfaces. For the study, sandstone and granite joint blocks with specific morphology were fabricated by using a three-dimensional (3D) engraving technique. Bonded rock-cement joints with asperity inclination angles of 15°, 30°, and 45° were prepared. Shear tests were performed on these bonded rock-cement joints to investigate the shear response and failure modes considering the effect of applied normal stress and interface morphology. Meanwhile, the two-dimensional particle flow code (PFC2D) was utilized to model the entire shear process of bonded rock-cement interfaces. The macroscopic shear behavior and mesoscopic failure mechanism were comprehensively investigated by the laboratory test and numerical simulation. The results showed that the shear stress-displacement curves of bonded rock-cement joints exhibit two distinct peaks, and the shear stress evolution can be categorized into four stages including elastic growth, rapid stress drop, secondary stress growth, and progressive softening. Significantly, the number of acoustic emission events also exhibits two distinct peaks related to the double peak of the shear stress curves. The failure of bonded rock-cement interfaces is mainly induced by shear fractures, while the failure of rock and cement blocks is primarily caused by tensile fractures. The number of shear cracks in the bonded rock-cement interfaces reaches the peak when the shear stress reaches the primary peak; whereas as the shear stress continuously approaches the residual stage, the fracture of the bonded rock-cement joints is primarily characterized by tensile cracks in the blocks.

Earthen dams made of Jiangxi lateritic soil are widely built in Jiangxi Province, China. Field observations showed uneven settlements and cracks in the earthen dams, which were attributed to the dynamic water-level fluctuation in the reservoir. Under this circumstance, the initiation and propagation of cracks of Jiangxi lateritic soil can be accelerated by the drying-wetting (D-W) cycles, threatening the safety and stability of earthen dams. For this reason, the dynamic characteristics of cracks in Jiangxi lateritic soil under D-W cycles and its microstructure mechanism were investigated in this study, for up to 5 cycles (N = 5). The microstructure of Jiangxi lateritic soil was measured using mercury intrusion porosimetry (MIP) and scanning electron microscopy (SEM), and its effect on the crack patterns was quantitatively analyzed through image processing technique. The results showed that: (1) The drying-induced desiccation cracks with increasing N was divided into three stages: the crack-generating stage (N ​= ​0–1), the crack-propagating stage (N ​= ​1–3) and the crack-stable stage (N ​= ​3–5). The initiation and propagation of cracks showed a strong correlation with microstructure damage (e.g., aggregate decomposition and pore expansion), which resulted from D-W cycles. With the penetration of large pores, the cracks were generated; (2) The wetting-induced healing behavior was categorized into two zones: the first zone corresponded to the healing of sub-cracks, while the second zone corresponded to that of primary cracks. With increasing N, the full-healing of primary cracks (N ​= ​2) was converted to the partial healing of primary cracks (N ​= ​3–4) and wetting-induced cracks (N ​= ​5); (3) The crack dynamic hysteresis (CDH) behavior consists of two stages, which were separated by a threshold water content (wth). With increasing N, the wth value decreased, indicating that more residual cracks, which were not healed in the wetting process, were accumulated. This study addressed the effect of D-W cycles on the cracks dynamic characteristics of Jiangxi lateritic soil, which can be helpful in the design of geotechnical engineering.

Grouting has been widely used in the reinforcement of jointed rock masses, and the interface between the rock and cement serves as a crucial binary interface in controlling the strength effect after grouting. In order to understand the shear failure behavior of the rock-cement interface, a series of unbonded sawtooth-shaped rock-cement interface specimens were fabricated and tested. Direct shear tests were conducted under a constant normal load to investigate the shear mechanical characteristics, considering the effect of the rock materials, interface morphology, and normal stress level. In addition, to investigate the mesoscopic failure behavior in detail, the particle flow code framework (PFC2D) was employed to simulate the shear behavior of unbonded rock-cement interfaces under combined compression and shear load action. It was found that the shear failure behavior was mainly concentrated at the cement block asperities, and the failure mode changed from wearing of the asperities on the cement block to cutting failure with increasing asperity inclination. The mesoscopic failure of the unbonded rock-cement interface was dominated by tensile failure of the cement block asperities. The progressive shear failure of the unbonded rock-cement interface showed as an extension of the damage from the surface to the interior of the asperities. Based on the discovered shear mechanism, a theoretical model for predicting the shear strength of unbonded rock-cement interfaces was derived by adopting the deduction method of the peak dilation angle. The model has clear physical significance because it considered both the interface morphology and shear properties, and was verified as having high prediction precision, providing theoretical guidance for studying the shear behavior of rock-cement interfaces.

Studying droplet behavior at intersections within three-dimensional fractured media is crucial for predicting fluid flow and solute transport. However, the relationship between droplet splitting at these intersections and macroscopic quantities like flow rate is still unclear. In this study, we developed and validated a theoretical model to describe droplet splitting in three-dimensional fractures using laboratory visualization experiments. Our findings show a strong alignment between the model's predictions and the experimental data, indicating the model's effectiveness in simulating real-world droplet behavior. We observed droplet behavior and transformation within a single three-dimensional fracture, noting that size is influenced by flow rate, aperture, contact angle, and inclination. Using these variables, we determined the capillary barrier size at intersections. We used our model to predict droplet splitting at intersections, defining three modes: Type I, full invasion of the horizontal fracture; Type II, partial invasion; and Type III, complete crossing of the fracture. Sensitivity analysis showed that droplet splitting shifts from Type I to Type II and stabilizes at Type III as the forward contact angle, flow rate, and inclination angle increase. Notably, there is a non-linear relationship between the splitting ratio and fracture aperture, especially during the Type Ito Type II transition. Our findings improve understanding and contribute to more accurate predictions of unsaturated flow in fractured rock formations, impacting groundwater remediation, oil and gas production, and geothermal energy extraction.

The existing quantitative evaluation methods of joint roughness are rich in achievements, among them most methods mainly focus on the influence of either the inclination or the height of joint asperity, and seldom consider both the two factors in the overall roughness quantification of rock joints. In this study, a batch of granite and sandstone joints were prepared for three-dimensional morphology analysis and were subjected to direct shear tests. The Structure-from-Motion (SfM) photogrammetry technique was adopted to digitally reconstruct the three-dimensional joint morphology, which provided data basis for roughness quantification. The non-stationary joint morphology was identified and was removed from the original morphology. Based on the stationary morphology feature, the influence of asperity amplitude on roughness estimation was investigated and a significant effect was revealed. A new statistical roughness parameter, the amplitude-weighted average asperity inclination θaw, was proposed that consider extra the contribution of the asperity height feature. In order to quantitatively estimate the JRC for a certain rock joint, a prediction model was suggested to assess the JRC of two-dimensional (2D) joint profile using the new roughness parameter θaw. The model was validated to be effective in predicting JRC of joint profiles from published studies. Subsequently, another model was established to predict the JRC of three-dimensional (3D) joint surface by incorporating a three-dimensional influence factor f3D into the 2D model. This 3D model was verified to have high prediction accuracy in quantifying the roughness of rough joints, through both test results and published data.

The mixing intensity plays a crucial role in determining the properties of cement-based materials. However, systematic research on this topic has been hindered by the lack of stable and reliable high-intensity mixing devices, leaving the quantitative influence on the hydration process and the underlying mechanism unclear. To address this gap, this study independently developed a high-intensity mixing device for cement-based materials preparation with a maximum no-load speed of 24,000. The findings revealed that high-intensity mixing (3000 rpm) remarkably enhanced both fluidity and compressive strength of cement paste. Compared with conventional mixing methods, the high-intensity mixed samples exhibited an improvement of 50.39 % in fluidity, along with 50.12 % and 21.44 % increases in 1-day and 28-day compressive strengths, respectively. Particle dispersion characterization revealed that high-intensity mixing effectively improved the dispersion and uniformity of cement particle, which was beneficial for a more intense hydration inside the samples. Thermogravimetric (TG) and 29Si nuclear magnetic resonance (NMR) analyses verified the hydration enhancement effect that high-intensity mixed samples achieved higher hydration degrees and longer average chain lengths. 1H NMR was further conducted to verify the reduced pore volume and decreased porosity in high-intensity mixed samples, thereby significantly enhancing the performance of cement-based materials. This study provides a novel and practical high-intensity mixing technology for the performance improvement of cement-based materials with broaden industrial application scenario.

Backward erosion piping (BEP) at the base of levees is one of the main causes of levee failure. BEP typically occurs in two-layer levees systems and involves multiple stages throughout its development. While random field methods have been widely applied to model soil spatial variability, existing approaches predominantly focus on single-layer soil structures. Analyzing the impact of only single-layer spatial variability on piping development, especially during the initial uplift stage, may not adequately capture the complexities of actual conditions. To address this, this paper proposes a general framework that can be used to generate both single-layer or two-layer random field models. The Karhunen-Loe`ve (K-L) series expansion method has been refined to generate two independent random fields, thereby facilitating the characterization of spatial variability in both the upper and lower soil layers. Monte Carlo Simulation (MCS) is employed to generate realizations of the random fields, which are then incorporated into numerical models to evaluate the effects of spatial variability on the whole process of BEP. The results show that single layer models tend to overestimate or underestimate the probability of piping initiation. The stage of backward erosion is influenced by the spatial variability of the upper and lower layered soil parameters. The single-layer random field model has a tendency to overestimate the probability of slope failure. In addition, the results of several stability analyses under the effects of extreme rainfall and upstream water level rise show that the probability of slope failure increases significantly during this process.

Reasonable division methods of landslide susceptibility indexes (LSIs) are crucial for producing landslide susceptibility levels (LSLs), including very low, low, moderate, high, and very high levels. However, few studies have systematically compared division methods such as natural break, equal interval, quantile, geometric interval, and K-means. Moreover, these methods start from LSIs but ignore the nonlinear correlation between known landslides and LSIs. To address this, the natural break-frequency ratio (FR) method is proposed, combining the natural break method for LSLs division with the FR method. First, the five conventional methods divide LSIs predicted by three machine learning models in An’yuan County, China. Then, the natural break-FR method is proposed to divide the same LSIs and compared with these methods. The natural break-FR, equal interval and K-means method yielded the largest sum of landslide ratio in very high and high susceptibility level, showing these methods can use high and very high susceptibility levels to predict as many landslides as possible. Finally, statistical perspectives of known landslide identification, division area proportion, and landslide ratio are applied to discuss how to select a suitable division method. Results show different division methods have comparative effects on final LSLs. The landslide ratios of equal interval, K-means, and natural break methods at high and very high susceptibility levels are greater than the former methods. The natural break-FR method performs best with MLP and SVM, but in the more precise RF model, the equal interval method outperforms it, followed by the natural break-FR method.