Numerical simulation and safety distance analysis of slope instability of ionic rare earth tailings in different rainy seasons

Abstract

A significant number of ionic rare earth tailings slopes were generated during the process of rare earth mining and the factor of safety of these slopes continued to decrease under the action of rainfall. This study focussed on a specific tailings slope located in Dayu County, Ganzhou City. Through the utilization of MatDEM numerical simulation software, the variations in infiltration distribution and destabilization evolution of the tailings slope between different seasons were comparatively examined. The results indicated that during the rainy season (May to June), increased rainfall led to increased rate of water infiltration, consequently reducing the stability of the slope and exacerbating landslide damage. Furthermore, the study provided analyses of safety distances under different seasons of rainfall, taking into account the geographical location of the tailings slope in Dayu County, Ganzhou City, and the results of the analyses provided theoretical foundation into ensuring the safety of buildings and residents downstream of the slope.

Keywords:

• MatDEM
• rainfall
• safety distance
• tailings landslides

1. Introduction

China has a significant amount of rare earth resources (Wübbeke 2013), which accounts for approximately 23% of the world’s total reserves. Ionic rare earth resources are primarily found in southern mountainous areas such as Ganzhou City, Jiangxi Province and Longyan City, Fujian Province (Central People’s Government of the People’s Republic of China 2018). In pursuit of higher profit returns, some private companies have accumulated a large amount of tailings sand during the mining process, further increasing the potential risk of local tailings landslides (Nguyen et al. 2019; Ongpaporn et al. 2022). These companies use technical processes such as heap leaching and pond leaching, resulting in a relatively low output ratio of rare earths, with one tonne of rare earths generating more than 2,000 tonnes of tailing sand (Central People’s Government of the People’s Republic of China 2012). Because the heavy metals of tailing sand are excessive, tailing sand can only be stockpiled centrally (Yin et al. 2016). In recent years, with the government’s emphasis on protecting arable land, the pile of tailing sand is getting higher and higher (Central People’s Government of the People’s Republic of China 2017). In the case of continuous rainfall, the slopes are prone to sudden environmental pollution and landslides (Ran et al. 2018; Senthilkumar et al. 2018). Dayu County has several tailings, including one located in Zuo Ba town, as shown in Figure 1 below. Downstream of the tailings, there are residential houses, farmland, roads, equipment, and factories. During the rainy season, there is a significant risk of tailings landslides that could endanger the safety of the residents.

Figure 1. Location of the tailings in Zuo Ba town.

There are many factors that affect the stability of slopes, such as the type of rainfall, strength reduction factor, cracks and duration of rainfall (Komolvilas et al. 2021; Yuan et al. 2021; Kafle et al. 2022; Yamusa et al. 2022). Wu et al. (2022) analyzed the relationship between the degree of saturation and rainfall intensity based on different types of rainfall intensity and discovered that the higher the rainfall intensity, the faster the change in the degree of saturation. Luo et al. (2021) found a negative correlation between rainfall duration and safety factor by changing the duration of artificial rainfall. This relationship may be attributed to the depth of infiltration, where greater depth corresponds to lower safety factor for the slope (Nguyen et al. 2017). Cracksas a key factor influencing water infiltration (Chen et al. 2018), have a strong adsorption effect on water (Wang et al. 2022; Zhang et al. 2023), and water stored in cracks reduces the shear strength of slopes (Liu et al. 2005; Miao et al. 2022). Moreover, landslide susceptibility mapping plays a vital role in analyzing landslide stability. For instance, Miao et al. (2023) proposed a high-precision boosting-C5.0 machine learning model, which can serve as a foundation for geological hazard risk management and control.

2. Methods and materials

In this paper, the analysis related to tailings landslides was mainly done through numerical simulation, as shown in Figure 2. Before simulating landslides, the macromechanical parameters of discrete units were set according to the actual physical properties of tailing sand, and the infiltration process under different rainfall conditions was calculated using MatDEM based on Darcy’s law (Whitaker 1986). During the process, saturation increased with increasing water content and exceeded the saturation threshold for water content, which was around 30% (Gu et al. 2022). As saturation continued to increase, the connections between discrete units in the cracks began to break, forming unstable areas along the cracks. This ultimately led to the occurrence of tailings landslides.

Figure 2. The process of numerical analysis.

2.1. Mechanical analysis of discrete unit connection

In MatDEM, discrete units are connected to each other by springs and generate the normal force ($F_n$) and shear force ($F_s$) (Liu 2021). Based on the Mohr-Coulomb criterion, the maximum shear force ($F_{s\max }$) that the springs can withstand can be determined, and a damage criterion for discrete units can be established. When the shear force ($F_s$) exceeds the maximum shear force ($F_{s\max }$), or normal displacement ($X_n$) is less than 0 or greater than or equal to the fracture displacement ($X_b$), the tangential or normal springs will break, resulting in the disconnection of discrete units, as represented by the following equations. (1) $$F_{\mathit{smax}}=F_{so}-{\mathrm{\mu }}_p{F}_n$$(1) (2) $$F_s=\{\begin{array}{c}{K}_s{X}_s, F_s\le F_{\mathit{smax}}\left(\mathit{Connections} \mathit{Status}\right)\\ -{\mathrm{\mu }}_p{F}_n,F_s>F_{\mathit{smax}}\left(\mathit{Disconnection} \mathit{Status}\right)\end{array}$$(2) (3) $$F_n=\{\begin{array}{c}{K}_n{X}_n,X_n<X_b\left(\mathit{Connections} \mathit{Status}\right)\\ K_n{X}_n,X_n<0\left(\mathit{Disconnection} \mathit{Status}\right)\\ 0,X_n\ge X_b\left(\mathit{Disconnection} \mathit{Status}\right)\end{array}$$(3) where $F_{so}$ is the shear resistance between discrete units, and ${\mathrm{\mu }}_p$ is the friction coefficient between discrete units. $K_n$ is the normal stiffness, $K_s$ is the tangential stiffness, and $X_s$ is the tangential displacement.

2.2. Materials of the discrete units

In stacking models, the stacking and cementation of discrete units change their mechanical parameters, and the values of compressive strength, tensile strength, and modulus of elasticity of the materials are also affected by the phenomenon of "size effect" (Zhou et al. 2016). In MatDEM, the macroscopic mechanical parameters of discrete units are usually smaller than the actual values, e.g. Young’s modulus($\mathit{E}$) is usually 20% to 40% of the actual values, and compressive strengths (${\mathit{C}}_{\mathit{\mu }}$) and tensile strengths (${\mathit{T}}_{\mathit{u}}$) are 10% to 20% of the actual values (Liu et al. 2017).

In this simulation, the Young’s modulus of discrete units was 36% of the actual values, and the compressive and tensile strengths were 18% of the actual values. Since MatDEM is currently unable to adjust the Poisson’s ratio and the coefficient of internal friction of the material, the Poisson’s ratio $(\mathit{v})$ and the coefficient of internal friction $({\mathit{\mu }}_{\mathit{i}})$ of discrete units were kept the same as the actual values. The relevant mechanical parameters are shown in Table 1 below.

Table 1. Mechanical parameters of real materials and discrete units.

	Type	$\mathit{E}$ (GPa)	$\mathit{v}$	${\mathit{C}}_{\mathit{\mu }}$(MPa)	${\mathit{T}}_{\mathit{u}}$(MPa)	${\mathit{\mu }}_{\mathbf{i}}$
Actual material	Tailings	30.6	0.23	15.3	3.5	0.42
	Granite	65.2	0.27	202.8	20.6	0.72
Discrete units	Layer1	11.01	0.23	2.75	0.63	0.42
	Layer2	23.47	0.27	36.50	3.71	0.72

In addition, there are conversion equations between the macroscopic mechanical parameters and the micromechanical parameters, which can be used to obtain the five micromechanical parameters of discrete units, as shown in the following Table 2:

Table 2. The microscopic mechanical parameters of discrete units.

	${\mathit{K}}_{\mathit{n}}$(GN·m^-1)	${\mathit{K}}_{\mathit{s}}$(GN·m^-1)	${\mathit{X}}_{\mathit{b}}$(m)	${\mathit{F}}_{\mathit{so}}$(GN)	${\mathit{u}}_{\mathit{p}}$
Layer 1	182.3	63.8	1.75*10^-3	1.01	0.045
Layer 2	392.5	137.8	4.77*10^-3	2.78	0.076

The five conversion formulas: (4) $$K_n=\frac{\sqrt{2}Ed}{4(1-2v)}$$(4) (5) $$K_s=\frac{\sqrt{2}(1-5v)Ed}{4(1+v)(1-2v)}$$(5) (6) $$X_b=\frac{3K_n+K_s}{{6\sqrt{2}K}_n\left(K_n+K_s\right)}.T_u{d}^2$$(6) (7) $$F_{so}=\frac{1-\sqrt{2}{\mathrm{\mu }}_p}{6}·C_{\mathrm{\mu }}{d}^2$$(7) (8) $${\mathrm{\mu }}_p=\frac{-2\sqrt{2}+\sqrt{2}I}{2+2I}$$(8) (9) $$I={(\sqrt{1+{{\mathrm{\mu }}_i}^2}+{\mathrm{\mu }}_i)}^2$$(9) where $d$ is the diameter of discrete units and $I$ is the scaling factor.

2.3. Moisture diffusion model

Based on Darcy’s law, the water infiltration equation for unsaturated soil is shown below (Ross 1990). (10) $$\frac{\partial \theta }{\partial t}=\frac{\partial q}{\partial x}=\frac{\partial }{\partial x}\left[K\right(h\left)\frac{\partial h}{\partial x}\right]$$(10) where $\theta$ is the water content (%), $q$ is the soil water flux density (cm s⁻¹), $K$ is the unsaturated hydraulic conductivity (cm s⁻¹), $h$ is the water pressure head (cm), $x$ is the horizontal distance (cm) and $t$ is the time of infiltration (s).

In this numerical simulation, MatDEM implemented the infiltration process of water moving from the upper boundary to the lower boundary by selecting the path with the shortest distance for water infiltration. The equation governing the motion of water during vertical infiltration is as follows (Wang et al. 2018). (11) $$\{\begin{array}{c}\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left(D\left(\theta \right)\frac{\partial \theta }{\partial z}\right)+\frac{\partial K\left(\theta \right)}{\partial z}\\ \theta (z,0)={\theta }_0\left(z\right)\\ \theta (0,t)={\theta }_s\\ \frac{\partial }{\partial n}\theta (L,t)=0\end{array}$$(11) where $L$ is the boundary, ${\theta }_0$ is the initial moisture content (%), $D(\theta )$ is the water diffusivity, $z$ is the vertical distance (cm) and $n$ is the shape parameter of the soil pore distribution. In addition, according to the two most common B–C and V–G models (Kosugi 1994), the specific expressions for the nonlinear term $D(\theta )$ and the hydraulic conductivity of unsaturated soil $K(\theta )$ are obtained.

B–C model: (12) $$\{\begin{array}{c}K(\theta )=K_s(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}{)}^{m/n}\\ D(\theta )=K(\theta )\cdot \frac{dh}{d\theta }=K_s\cdot h_d\cdot (-\frac{1}{n})\cdot ({\theta }_s-{\theta }_r{)}^{(1-m)/n}\cdot (\theta -{\theta }_r{)}^{(m-1-n)/n}\end{array}$$(12)

V–G model: (13) $$\{\begin{array}{c}K(\theta )=K_s[\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}{]}^{1/2}1-1-(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}{)}^{1/m}{]}^m{]}^2\\ \begin{array}{cc}D(\theta )& =K(\theta )\cdot \frac{dh}{d\theta }=\frac{K_s(1\text{-m})(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}{)}^{1/2-1/m}}{\alpha m({\theta }_s-{\theta }_r)}\cdot [{(1-{(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r})}^{1/m})}^{-m}+{(1-{(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r})}^{1/m})}^m-2]\end{array}\end{array}$$(13) where ${\theta }_r$ is the residual moisture content (%), ${\theta }_s$ is the saturation moisture content (%), and $K_s$ is the saturated hydraulic conductivity (cm s⁻¹). $\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}={(\frac{1}{\alpha h})}^n,$ $m=2+3n$ (B–C model), $\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}={\left[\frac{1}{1+(\alpha h{)}^n}\right]}^m,$ $m\text{= 1}\mathrm{–}\text{1/}n$ (V–G model), $\alpha$ is the parameter of the shape of the characteristic soil moisture curve.

The slope was divided into two parts, the regular and irregular areas, based on the boundary conditions of the numerical simulation. The area to the left of the dividing line represents the regular area, where water infiltrates vertically. The corresponding equation is stated below. (14) $$\{\begin{array}{c}\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial x}\left(D\left(\theta \right)\frac{\partial \theta }{\partial x}\right)+\frac{\partial }{\partial z}\left(D\left(\theta \right)\frac{\partial \theta }{\partial z}\right)-\frac{\partial K\left(\theta \right)}{\partial z}\\ \theta (x,z,0)={\theta }_0,\left(x,z\right)\in [0,2]\times [0,8.5]\\ \theta (x,z,t)={\theta }_s,x=0,z=0\\ \frac{\partial }{\partial z}\theta (x,z,t)=0,(x,z)\in (0,2]\times [0,8.5]\end{array}$$(14)

The area to the right of the zoning line represented the irregular area, where water infiltrated into the slope’s interior. The corresponding equation is stated below. (15) $$\{\begin{array}{c}\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial x}\left(D\left(\theta \right)\frac{\partial \theta }{\partial x}\right)+\frac{\partial }{\partial z}\left(D\left(\theta \right)\frac{\partial \theta }{\partial z}\right)-\frac{\partial K\left(\theta \right)}{\partial z}\\ \theta (x,z,0)={\theta }_0,\left(x,z\right)\in \Omega ,\Omega \in [2,6.7]\times [0,8.5]\\ \theta (x,z,t)={\theta }_s,x=2,z=-\tan \hspace{0.17em}(\alpha )\cdot (x-1),x\in [2,6.7]\\ \frac{\partial }{\partial z}\theta (x,z,t)=0,z=-\tan \hspace{0.17em}({\theta }_L)\cdot (x-1),x\in [2,6.7],{\theta }_L=60^{\circ }\end{array}$$(15)

2.4. Rainfall analysis

Rainfall is one of the main causes of landslides disasters (Kuradusenge et al. 2020; Ahmed 2021), and rainfall is directly related to seasonal changes (Li et al. 2016). Warm and humid air currents from the Pacific Ocean pass through several provinces in southern China, resulting in a 2-month-long period of rain known as the rainy season. In addition, the time of occurrence of the rainy varies from region to region. Based on meteorological data from Weather Atlas (Weather Atlas 2002–2023), the average rainfall data for Ganzhou city were compiled. As depicted in Figure 3, the rainy season in Ganzhou occurred between May and June each year. Among these, the average rainfall during the non-rainy season was 73.4 mm, while the average rainfall during the rainy season reached 210 mm, which was approximately 2.86 times higher than the average rainfall during the non-rainy months. Prolonged rainfall led to the accumulation of rainwater on the slope surface and increased the saturation of tailing sand, thereby reducing the slope’s stability (Gasmo et al. 2000).

Figure 3. The average rainfall in Ganzhou city.

2.5. Slope stability analysis

Rainfall causes an increase in the self-weight of tailing sand, which disrupts the original mechanical equilibrium. Without considering the equilibrium of forces between the slices, the slope is divided into equally spaced slices according to Swedish slice method (Enoki et al. 1990), which can be used to find the integral ratio of the shear stress to the shear strength on each slice can be found., i.e. the safety factor ($F_s$). The formula for calculating the stabilisation coefficient ($F_s$) is as follows. (16) $$F_s=\frac{\sum _{i=1}^N (c_i{l}_i+w_i\mathrm{c}\mathrm{os}{\theta }_i\mathrm{t}\mathrm{an}{\phi }_i)}{\sum _{i=1}^N w_i\mathrm{s}\mathrm{in}{\theta }_i}$$(16) $c_i,$ ${\phi }_i$ are the cohesion and friction angle of the ith slice, respectively, and N is the number of divided slices. $l_i,$ ${\theta }_i$ and $w_i$ are the bottom length, bottom angle, and weight of the ith slice, respectively. Under continuous rainfall, the self-weight of the slices increases as the saturation of tailing sand increases, but the cohesion and angle of internal friction of tailing sand decreases (Zhang and Ding 2012). Thus the factor of safety decreases with increasing rainfall (Raj and Sengupta 2014), and a larger increase in rainfall leads to a faster decrease in the safety factor, resulting in greater slope instability.

2.6. Numerical model and monitoring sites

The upper limit of the number of units in Madam is approximately 6 million (Liu et al. 2017–2021), Moreover, a higher number of units or smaller particle sizes of the units impose greater performance requirements on the computer. The particle size of tailing sand exhibits a relatively wide gradation, with 80% of the particles ranging from 0.05 to 0.5 mm. However, such a wide gradation can hinder simulation efficiency, necessitating a reduction in the range of tailing sand particle size. The particle size of discrete units in the model ranges from 0.028 m to 0.033 m, as shown in Figure 4. In addition to the grain size, the slope angle is also a worthwhile aspect to consider in establishing the slope model. In general, the greater the slope angle, the more likely tailings landslides will occur. Slopes generally range from 30° to 60°, and in order to analyze the maximum damage caused by the tailings landslides to people and buildings, a slope with an angle of 60° was chosen for modelling. Regarding stratification, the slope was divided into tailing sand and bedrock layers. Tailing sand has a loose structure and the density of 1.48–2.56 g/cm³, and layer 2 is used to replace tailing sand in the numerical simulation. Since the radius of discrete units in layer 2 was different from that of tailing sand, the hydraulic parameter values of layer 2 were different from the actual hydraulic parameter values of tailing sand. The hydraulic parameter values of layer 2 were ${\theta }_r$=0.1, ${\theta }_s$=0.9, $K_s$=5.5, $\alpha$=0.71 and $n$=0.6. The bedrock is mainly granite, which has higher compressive and tensile strength, and the density of 2.60–2.82 g/cm³, and layer 1 was used to replace the bedrock in the numerical simulation. Layer 1 had very little effect on the slope, and it was set up mainly to have a supporting effect on the particles of the model and to prevent discrete units from flying out of the interior of the model. Therefore, the values of hydraulic parameters in layer 1 were the same as those in layer 2. In addition, in order to analyse the process of water infiltration, seven detection zones were set up, as shown in Figure 5. Four vertical monitoring zones with X values of 0.2–0.3 m and Z values of 8.2–7.8, 7.2–6.8, 6.2–5.8, and 5.2–4.8 m, respectively, were used to study the vertical infiltration of water. Three horizontal monitoring zone with Z values of 3.9–4.1 m and X values of 1.5–1.7, 2.7–2.9, and 3.9–4.1 m, respectively, were used to study the horizontal infiltration of water.

Figure 4. The particle size distribution.

Figure 5. The initial model diagram.

3. Results

Numerical simulations of tailings landslides under rainfall were conducted in this study, and the process took approximately 36 s. Twelve MatDEM files and eight data tables were generated for analysis, containing key information such as degree of saturation, displacement, velocity, and impact force. The data were extracted, classified, filtered, and analyzed to derive the evolution patterns of the displacement field, moisture field, and impact force.

3.1. Moisture field data analysis

Figure 6(a) shows that Monitoring Point 1 reached a saturation of 78% at 15.9 s and remained relatively stable thereafter. Monitoring Point 2 began to show an increase in saturation at 5.8 s, while Monitoring Point 3 started to exhibit an increase at 17.9 s. Monitoring Point 4 experienced minimal changes in saturation. In Figure 6(b), Monitoring Point 1 approached a saturation level of approximately 78% at 36 s. The saturation of monitoring point 2 started to increase at 7.7 s, and the degree of saturation increased to about 32% at 36 s. During different rainy seasons, there was a significant disparity in the saturation levels of Monitoring Point 2 at 36 s, with a difference in saturation of approximately 42%. The saturation of monitoring point 3 and point 4 remained essentially unchanged, indicating that there was no water infiltration into monitoring point 3 during the non-rainy season. This indicates that the depth of water infiltration was shallower than 2.5 m.

Figure 6. The variation of the degree of saturation at vertical monitoring points. (a)The rainy season; (b) the non-rainy season.

In addition to the monitoring points of the degree of saturation in the vertical direction, three horizontal monitoring points were also set up, namely monitoring point 5, monitoring point 6, and monitoring point 7. Figure 7(a) illustrates the changes in saturation levels for the horizontal monitoring points during the rainy season. At 36 s, the saturation of monitoring point 7 was 79%. The saturation of monitoring point 5 and monitoring point 6 increased from 0 s and exceeded 30% at 28.8 and 23.1 s, respectively. Due to the faster vertical infiltration of water compared to the horizontal infiltration in the cracks, the saturation of monitoring point 5 and point 6 was greater than that at monitoring point 7 before 11.4 s. However, the saturation of monitoring point 5 and point 6 was less than that of monitoring point 7 after 11.4 s, this was because the locations of monitoring point 5 and point 6 are closer to the interior of the slope than monitoring point 7. Figure 7(b) depicts the variation of saturation levels for the horizontal monitoring points during the non-rainy season. The saturation changed curves for monitoring point 5, point 6, and point 7 intersected at 20.1 s, which was delayed by 8.7 s, compared to the rainy season. The saturation of monitoring point 7 exceeded 30% at 25.4 s, which is delayed by 12.2 s, compared to the rainy season. The delay indicated that the rate of water infiltration under the non-rainy season was smaller than the rate of water infiltration under the rainy season.

Figure 7. The variation of the degree of saturation of horizontal monitoring points. (a) The rainy season; (b)the non-rainy season.

3.2. Analysis of water infiltration

Rainfall is a key condition for triggering landslides (Argyriou et al. 2022). With continuous rainfall, water infiltrated from discrete units with high saturation into discrete units with low saturation. Figure 8 illustrates the temporal distribution of saturation. During the rainy season, it could be observed that water infiltrated more rapidly, reaching a depth of approximately 2 m. During the non-rainy season, the amount of water infiltration was less compared to the rainy season and the water infiltrated at a relatively slower rate to a depth of about 1.5 m. Figure 8 also reveals that the presence of cracks enhanced the rate of water infiltration.

Figure 8. The time-range distribution of saturation. (a) The rainy season; (b) the non-rainy season.

3.3. Connection evolution law

The connections between discrete units in this simulation are shown in Figure 9 below. The green area indicates that discrete units are connected to each other, and the blue area indicates that discrete units are disconnected from each other. As water started to seep into the interior of the slope, the connections of discrete units gradually broke down in the cracks, and tailings landslides started to occur. When the saturation of discrete units in cracks to 30%, the connection began to break. Thus, the change of connection could correspond to the increase of saturation of discrete units in the cracks. In the rainy season, as shown in Figure 9(a), the connection of crack 1 was almost completely disconnected at 9 s, the connection of crack 2 was completely disconnected at 27 s, the connection of crack 3 was partially disconnected at 36 s, and the slope underwent an overall landslide. In the non-rainy season, as shown in Figure 9(b), the connection of crack 1 started to disconnect at 9 s and was almost completely disconnected at 18 s. The connections of crack 2 and crack 3 were partially disconnected at 36 s, and only partial landslide occurred on the slope.

Figure 9. The changes in connections between discrete units during different rainy seasons: (a) the rainy season; (b)the non-rainy season.

3.4. Safety distance analysis

In order to analyze the safety distance more accurately, the intervals of 7–9 m, 9–11 m, 11–13 m, 13–15 m, 15–17 m, 17–19 m, and 19–21 m in the x axis were set as the spacing areas. From the start of tailings landslides until their completion, the area and impact force of discrete units in the monitoring areas were calculated at intervals of 0.1 s. These calculations are presented in Figures 10 and 11. Based on the data from Figure 10, the maximum area of discrete units within each spacing area during various rainy seasons was determined, and the average height of the maximum area was calculated as the depth of inundation according to the area formula. The related research (Pollock and Wartman 2020) suggested that when the inundation depth ranged from 0.9 to 5.9 m, the potential harm caused by tailings landslides to individuals increases. Therefore, a critical value of 0.9 m was selected to determine the depth of inundation.

Figure 10. Area of discrete units in the spacing areas: (a) The rainy season; (b) the nonrainy season.

Figure 11. Impact force of discrete units in the spacing areas: (a)The rainy season; (b) the nonrainy season.

In addition, the maximum force of impact of discrete units in the spacing areas of different rainy seasons could be obtained from Figure 11, which was useful for analysing the safety distances of buildings. Since the vast majority of buildings are nowadays reinforced concrete structures, the object of analysis in this paper was the reinforced concrete structural type of buildings. Previous research (Kang and Kim 2016) indicated that impact pressures exceeding 100 kPa resulted in extensive damage to buildings, while impact pressures ranging from 35 to 100 kPa led to moderate damage. Impact pressures below 35 kPa were associated with slight damage to structures. The results of the specific analysis are shown in Table 3.

Table 3. The safety distance analysis for landslides.

Seasons	Data contents	7–9 (m)	9–11 (m)	11–13 (m)	13–15 (m)	15–17 (m)	17–19 (m)
The rainy season	The average height/m	1.355	1.28	0.965	0.716	0.713	0.712
	Results	>0.9	>0.9	>0.9	<0.9	<0.9	<0.9
	Impact force/kPa	502.6	445.01	356.72	336.1	148.2	155.52
	Results	>100	>100	>100	>100	>100	>100
The non-rainy season	The average height/m	0.665	0.55	0.365	0.3	0.2	0.3
	Results	<0.9	<0.9	<0.9	<0.9	<0.9	<0.9
	Impact force/kPa	176.2	126.3	86.89	75.71	29.43	24.63
	Results	>100	>100	35 ∼ 100	35 ∼ 100	<35	<35

4. Discussion

In this paper, the processes of tailings landslides during different rainy seasons were analyzed using MatDEM. Compared with previous experiments and numerical simulations (Montrasio et al. 2016), the relationship between tailings landslides and rainfall amounts was investigated, and this relationship was verified through numerical simulations. Furthermore, the hazard of tailings landslides was further analyzed based on previous studies. While the study of tailings landslides was still in the preliminary stage, this study still provided a theoretical basis for local slope management.

Unfortunately, further studies on slopes with different gradients had not been carried out in this paper, resulting in limited data availability from modeling. In addition, the model only considered average monthly rainfall and did not consider the daily fluctuation of rainfall.

For the prevention of tailings landslides, an important direction would be to develop new techniques or materials for slope reinforcement (Ngo et al. 2019, 2023).

5. Conclusion

In this study, water infiltration and tailings landslides under different rainy seasons were analyzed using numerical simulations. Discrete units were constructed with specific elastic modulus and strength properties through a combination of conversion formulas for discrete units and numerical simulation experiments, enabling automated modeling of discrete units. Furthermore, additional analyses were conducted to determine the safety distances downstream of the tailings slope based on the numerical modeling results. Regarding the depth of landslide inundation, the inundation depth was less than 0.9 m within the 13–19 m during the rainy season, while in the non-rainy season, the inundation depth was less than 0.9 m across all spacing areas. In terms of building damage, extensive damage was observed in all spacing zones during the wet season. During the non-rainy season, buildings within the range of 7–11 m suffered extensive damage, those within the range of 11-15 m suffered moderate damage, and those within the range of 15–19 m suffered slight damage. The numerical simulations demonstrated differences in water infiltration and tailings landslides between the rainy season and the non-rainy season, with higher water infiltration rates and more severe damage observed during the rainy season. Therefore, special attention should be given to ensure the safety of tailings slopes during the rainy season.

Acknowledgments

The authors thank the anonymous reviewers for their helpful suggestions.

Data availability statement

Data is available on request from the authors.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was funded by China Postdoctoral Science Foundation, Grant No. 2022M722924; National Natural Science Foundation of China, Grant No. 52364013; “Youth Innovation Talent Introduction and Education Plan" of Shandong Colleges and Universities (Grant No. Lu Jiao Ke Han [2021] No. 51), and Jiangxi Province Thousand Talents Program Project, Grant No. jxsq2019102082.