#### Abstract

The buildup of tailings in China has expanded dramatically with economic development and industrial demand, and the safety of tailings reservoirs has become increasingly serious. Due to the difficulty in finding a new reservoir site, the expansion approach of building a new tailings dam downstream of the original reservoir area was investigated. The stability of the tailings reservoir after expansion was calculated using the traditional dynamic and static stability solution method and taking into account the unpredictability of dam construction materials and tailings material parameters in the reservoir area. The results reveal that throughout the tailings accumulation process in the new reservoir, the tailings will build a back pressure slope at the original reservoir’s initial dam, which can considerably improve the original reservoir’s dynamic and static stability. The Monte Carlo method clearly outperforms older methods for tailing pond stability analysis. The results of this paper’s calculations will give a theoretical foundation and practical reference for the later management and maintenance of such tailings reservoirs, as well as fresh ideas and insights for comparable projects due to limited site selection.

#### 1. Introduction

A tailings reservoir is a massive geotechnical structure that mining companies use to store tailings throughout the ore extraction and separation process. As the greatest geotechnical structure, the primary function of the tailings dam is to block tailings. There are approximately 24,605 tailings ponds in the world [1], and tailings dam breakdown events may create irrevocable catastrophic consequences [2–8] and long-term environmental pollution [9–12].

The stability of the tailings pond is crucial for preventing accidents. Numerous scientists have researched the seismic stability of tailings ponds [13–16]. Fu-Sheng L I et al. [15] estimated the infiltration line of a tailings dam slope, examined the dam body condition and stability safety factor after heightening and expansion, and demonstrated that the elevation and expansion scheme of a tailings rockfill dam is safe and dependable. Wei et al. [17] analyzed the successful cases of increasing and extending the service life of an existing tailings pond using the upper bend drainage system. Ozcan et al. [18] analyzed a copper-zinc tailings dam elevated by 7 m in the Black Sea region of Turkey using static and pseudostatic analysis methods. The results indicate that the tailings dam has maintained its static and dynamic stability after being raised by 7 meters. Many scholars have used the Monte Carlo approach to calculate the stability of tailings ponds due to its less problem-limiting nature, high calculation accuracy, simple operation, and high adaptability [19–21]. Li et al. [22] examined the probability and sensitivity of spatial fluctuation of sand layer using the Monte Carlo test principle. It was believed that different sampling procedures created huge fluctuations in the chance of instability and the dependability index but did not appreciably alter the safety factor.

To sum up, there is no relevant literature utilizing the Monte Carlo approach to analyze the stability of the expansion mode of the new tailings dam downstream of the original reservoir area that forms twin reservoirs. This study employs the Monte Carlo approach to solve the dynamic and static stability of the tailings reservoir in the double reservoir area and compares it with the standard stability solution method to determine the “good formula” for the stability solution of the double reservoir area project. The research results give recommendations for the project’s future operation, maintenance, and safety warnings, as well as a theoretical foundation and numerical reference for similar initiatives in the future.

This paper’s outline is as follows: in the second section, the engineering situation of the tailings pond, the distribution characteristics of the tailings, the physical parameters of the tailings material, and the calculation model are described. The third section elaborates on the static and dynamic analysis theory and Monte Carlo theory utilized in this study. The fourth section presents the static and dynamic stability calculation findings of tailings based on the traditional approach and Monte Carlo theory. The fifth section of this paper describes the research findings. The sixth section is the conclusion.

#### 2. Enlargement Scheme of Tailings Pond

##### 2.1. Engineering Overview

The working area is located in the western slope of Ailao Mountains in Yuxi City, Yunnan Province, China. The geographical location is shown in Figure 1 [23]. The initial dam of the expanded tailings pond is located at about 640.0 m downstream of the original tailings dam site in Banmao ravine, which is the same as the final accumulation elevation of the built tailings pond design, and finally forms the whole. The proposed initial dam height is 64.0 m, the dam crest width is 5.0 m, the dam bottom width is 255.0 m, the dam crest length is 180.0 m, and the dam type is rockfill dam; the total storage capacity is 1757.3 million m^{3}, and the effective storage capacity is 1521.9 million m^{3}, which belongs to the second-class storage. The final accumulation elevation is 1260.0 m, the accumulation slope is 1 : 4.5, and the projected service life is 19.4 years.

##### 2.2. Physical and Mechanical Properties of Each Component Material in Tailings Pond

In this paper, field density measurement, consolidation test, undrained static triaxial shear test, and dynamic triaxial test are carried out. The test materials are mainly composed of tailings dam accumulation material, subdam, tailings silt, and tailings silty clay. The dam is built from the local soil in Yunnan. Tailings silt and silty clay are taken from the tailings site. Tailings dam accumulation material is a kind of mixed material. According to the results of field compaction test at dam site, the rockfill material is tested indoor, and the average value of field dam slope material is graded.

Tailings dam body and dam material have nonlinear characteristics. Its deformation not only changes with the load size, but also is related to the stress path of loading. The stress-strain relationship is obviously nonlinear. In order to better describe the nonlinear characteristics of tailings dam and dam material, the nonlinear model is used to calculate the material characteristics. In nonlinear models, the Duncan-Chang Hyperbolic model is commonly utilized because of its clear concept [24], extensive expertise in parameter determination, and greater consistency with engineering practice. Therefore, in this paper, the Duncan–Chang model is used to calculate the nonlinear stress-strain relationship of dam materials. The parameters of Duncan–Chang model determined by three-axis experiments are arranged as shown in Table 1.

The maximum dynamic shear modulus is the main indicator of material dynamic parameters. At present, there are primarily three techniques for computing the dynamic shear modulus: field shear wave test calculation [25], calculation according to the undrained strength value [26], and empirical calculation considering physical properties [27]. In this paper, the empirical method is used to calculate the dynamic shear modulus of tailings dam materials. In the calculation process, the dynamic shear strain is 10^{−6}, and the relationship between the maximum dynamic shear modulus and confining pressure of silty tail soil and silty tail clay is shown in Figure 2.

The dynamic shear modulus ratio and damping ratio are crucial factors for evaluating the seismic safety of an engineering site and analyzing the seismic response of a soil layer. In this paper, the dynamic shear modulus ratio and damping ratio are used as two indices to measure the dynamic parameters of tailings dam materials. The dynamic shear modulus ratio and damping ratio of initial dam, subdam, tailings silt, and tailings silt clay are calculated. Through the calculation and analysis of the accumulation height of each component in the tailings pond, the confining pressure of the initial dam rockfill is set to 200 kPa and 800 kPa, and the consolidation ratio Kc is 2; the confining pressures of subdam, tailing silt, and silty clay are set to 100 kPa, 150 kPa, and 200 kPa, respectively, and the consolidation ratio Kc is 1. The calculation results are shown in Figure 3.

**(a)**

**(b)**

**(c)**

**(d)**

##### 2.3. Calculation Model of Tailings Pond Expansion

For the first case in China, a new tailings dam is constructed downstream of the old reservoir area, thereby upgrading the existing tailings reservoir from a third-level reservoir to a second-level reservoir. First, this extension is controlled by the local topography and geomorphology; second, the site selection in the area downstream of the original reservoir can focus on mitigating the reservoir’s environmental impact. Simultaneously, the development of additional tailings in the area downstream of the original reservoir can utilize the same infrastructure and significantly reduce economic expenses.

Figure 4 is the calculation model of three working conditions; that is, the tailings beach in the new reservoir area accumulates to the initial dam top of the original reservoir area (condition 1), the tailings beach in the new reservoir area accumulates to the central point of the atomic dam in the original reservoir area (condition 2), and the beach in the new reservoir area accumulates to the top of the atomic dam in the original reservoir area (condition 3). Under the three working conditions, the infiltration line height is inferred by field hole measurement.

#### 3. Stability Analysis Theory and Material Parameters of Tailings Pond

##### 3.1. Static Stability Analysis

Finite element and limit equilibrium methods are the most commonly employed for calculating the static stability of tailings ponds [28, 29]. The most extensively utilized technique among these approaches is the Swedish arc technique. Here, the static stability of the tailings pond is studied using the Swedish arc method and the finite element strength reduction method, and the following equation is established:

Swedish arc method is as follows:

Finite element strength reduction method is as follows:

In the formula,

is strip gravity, kN; *c*_{i} is cohesion, kPa; *Ф*_{i} is internal friction angle, °; *b*_{i} is the width of soil strip, m; *n* is the number of blocks; *θ*_{i} is the bottom slope angle, °; *γ* is the weight of soil, kN/m^{3}; *h*_{i} is the average height of soil strip, m; *γ* is the weight of water, kN/m^{3}; is the average height below the soil infiltration line, m; *τ*_{r} is the reduced shear strength, kPa; *τ*_{f} is shear strength before reduction, kPa; *σ* is the total stress, kPa; and *F*_{s} is the reduction coefficient.

##### 3.2. Dynamic Stability Analysis

###### 3.2.1. Establishment of a Dynamic Model considering Confining Pressure

In dynamic calculations and analyses, dam body and foundation materials are typically considered viscoelastic. The dynamic shear modulus and damping ratio *D* are used to reflect the nonlinearity and hysteresis of the dynamic stress-strain relationship and are expressed as the relationship between dynamic shear modulus, damping ratio, and dynamic shear strain. In this paper, the Hardin–Drnevich model [30] is used to calculate the dynamic shear modulus and damping ratio, as follows:

Dynamic shear modulus is as follows:

Damping ratio is as follows:

Maximum shear modulus is as follows:

In the formula, *γ*_{r} is the parameter shear strain, ; *γ* is the total shear strain; *K*_{2}, *n* are test parameters; *Pa* is atmospheric pressure, MPa; and is the average effective stress, kPa.

###### 3.2.2. Dynamic Control Equation and Solution

Wilson method is used to solve the dynamic control equation [31]. In the calculation process, it is assumed that the dynamic shear modulus and damping ratio *D* of each unit in each period remain unchanged, and the time step is set as *∆t* = 0.01 ∼ 0.02 s. The 0.65 times of the maximum shear strain *γ*_{max} in this period is taken as the average shear strain *̄* of the unit in this period.

The dynamic control equation is

In the formula, *δ*(*t*), , and are node displacement, velocity, and acceleration; *F*(*t*) is the dynamic load of the node, determined by the seismic acceleration; [*M*] is the unit mass matrix; [*K*] is the element stiffness matrix; and [*C*] is the element damping matrix.

###### 3.2.3. Earthquake Permanent Deformation

The Newmark approach [32, 33], simplified analysis method [34], equivalent inertia force method [35], and probability analysis method [36, 37] are the most common seismic permanent deformation calculation methods. When Newmark method is used to calculate the permanent deformation of soil, the following fundamental assumptions are made: (1) the effect of earthquake on dam body can be equivalent to a quasistatic inertia force with a defined magnitude and direction; obvious potential sliding surface will be formed when dams are damaged in earthquakes. (2) When the seismic action exceeds its ultimate seismic capacity, rigid-plastic sliding occurs along the potential sliding surface; the strength of the sliding surface will not degrade significantly during the earthquake. (3) The quasistatic horizontal seismic acceleration applied when the potential landslide is in the critical limit equilibrium state or the antisliding safety factor is 1 is called the yield seismic acceleration *a*_{y} = *k*_{y} of the dam body, where the ratio of the yield seismic acceleration to the gravity acceleration is called the yield acceleration coefficient *k*_{y}. (4) When the total sliding force generated by the potential landslide under the action of seismic inertia force exceeds the total antisliding force on the sliding body, the potential landslide will slide instantaneously along the potential sliding surface. The sliding deformation generated by this instantaneous overload will stop when the acceleration is reversed, and the velocity of the landslide reduces to zero. The seismic slip is the accumulation of the sliding displacement caused by all instantaneous overload pulses during the earthquake duration.

##### 3.3. Monte Carlo Method

The Monte Carlo approach is fundamentally a probability analysis technique. In the Monte Carlo sampling experiment, the material parameters are taken from the interval based on the assumption that the physical parameters of the material follow a particular probability distribution. The precision of the Monte Carlo sampling experiment is dependent on the sample size of the material characteristics and the number of extractions. When all requirements are qualified, the precision of the Monte Carlo sampling experiment can be ensured. In this article, 10,000 extractions are performed in order to calculate the stability of the new and original tailings using the Monte Carlo approach.

When the Monte Carlo method is used to analyze the stability of the tailings pond, the mean and standard deviation of the physical properties and mechanical parameters of the tailings material in Table 2 are input into the GeoStudio software. According to the normal distribution of the output parameters, the random values are substituted for *n* times in the solution equation of the dynamic and static stability. The *n* safety factors are *Z*_{1}, *Z*_{2}, *Z*_{3}, …, *Z*_{n}; if there are *m* safety factors greater than 1, the probability of dam failure can be expressed as follows:

At this time, the mean *μ*_{z} and standard deviation *σ*_{z} of the safety factor of the tailings dam are

##### 3.4. Seismic Acceleration Time-History Curve

According to the engineering geological survey study, the Simao area has experienced tectonic earthquakes in the past. According to records, Honghe, Ailaoshan, Jiujia-Anding, three NW-trending faults, have rarely seen earthquakes, which were mainly concentrated in the NW-trending secondary fault zone on the west side of the Amojiang fault zone, constituting Pu’er-Simao seismic tectonic belt. The earthquakes mainly occurred in the south of latitude 23°40′. The seismic activities that occurred from 1983 to 1984 also concentrated in the area from Pu’er to Simao. The epicenter extended northward to the vicinity of Mohei. There were 9 strong earthquakes with magnitude *M* ≥ 6, 4 moderate earthquakes with magnitude 6 > *M* ≥ 5, and 12 earthquakes with magnitude 5 > *M* ≥ 4. All of them were weak earthquakes with magnitude below 4. The epicenters of the above strong and weak earthquakes are distant from the engineering region and have minimal impact on its construction.

According to the feasibility study report of tailings pond expansion, the peak ground motion acceleration of the site is 0.10 g, the seismic fortification intensity is VII, and the characteristic period of the ground motion response spectrum is 0.45 s. In this calculation, two seismic waves are used, the time-history curves of El Centro waveform and Taft waveform are taken, and the maximum amplitude is reduced to 0.1 g and 0.15 g, as shown in Figure 5. In the dynamic calculation of tailings pond, the horizontal acceleration time history curve is input.

#### 4. Calculation Results and Analysis

##### 4.1. Static Stability

In order to analyze the stability of the tailings pond after adopting the new expansion method, three working conditions models of the tailings reservoir that the beach surface of the new reservoir area accumulates to different positions of the original reservoir area are designed. Since the new and the original two reservoir areas merge under working condition, there is no stability calculation result of the original reservoir area under . Critical sliding surface for slope stability calculation by limit equilibrium method is selected based on experience; that is, failure surface is generated from the outer slope angle or outer slope surface of the initial dam and is penetrated at the top of the subdam.

Figure 6 is the calculation results of the new and original reservoir areas based on the traditional method (the left is the Swedish arc method, and the right is the finite element strength reduction method) under three working conditions. The figure demonstrates that the safety factors of the new and original reservoir areas solved by the new tailings reservoir expansion method of double reservoir areas are far more than the standard design value, and the minimum value is the result of the original reservoir area under the Swedish arc method #1 condition. According to the calculation results of the finite element strength reduction method, it can be seen that there are two shear stress concentration areas in the original reservoir area under working condition, and the infiltration line can be obtained by comprehensive engineering speculation. The stress concentration area is near the intersection point of the infiltration line and the tail silt clay area of the original reservoir area. This is due to the saturation state of the tailings below the infiltration line, which leads to the decrease of the effective stress and then affects the stress distribution in the nearby area. Therefore, the static stability of the original reservoir area should be paid attention to in the construction project of the new reservoir area of the tailings reservoir in the double reservoir area. With the gradual increase of the beach surface in the new reservoir area, the stress concentration point in the original reservoir area disappears. This is because the gradual increase of the beach surface in the new reservoir area causes the back pressure slope protection at the initial dam and the subdam of the original reservoir area, thereby improving the stability of the original reservoir area. The two traditional methods calculate that the safety factor of the original reservoir area increases with the increase of the beach surface in the new reservoir area, which also proves this point. Consequently, compared to the way of extending the expansion at the original dam, the novel expansion method of the double reservoir area significantly enhances the storage capacity and is beneficial to enhancing the stability of the tailings reservoir. In the process of the gradual rise of the beach surface in the new reservoir area and the formation of the whole reservoir area, the intensive shear stress gradient of the reservoir area is gradually transferred from the original reservoir area to the new reservoir area and finally uniformly distributed at the dam foundation of the new reservoir area. Therefore, in the practical engineering of the double reservoir area, when the reservoir area meets operational condition , the monitoring center of the stability of the reservoir area should be transferred from the original reservoir area to the new reservoir area.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Figure 7 is the Monte Carlo calculation results of the new reservoir area under condition (MC-Swedish arc method on the left side and MC-strength reduction method on the right side). The calculation results of the new and original reservoir areas under other conditions are similar to those in Figure 7, so they are not listed. The findings of the MC-Swedish arc method and the MC-strength reduction approach are comparable to the normal distribution, whereas the standard calculation method yields a fixed value. This is because the traditional calculation method is to find the critical sliding surface in the calculation area and solve the stability, ignoring the discreteness and randomness of the soil, while the Monte Carlo method is to solve all the sliding surfaces in the calculation area. Therefore, the Monte Carlo method’s computation findings have greater practical relevance. When the minimum safety factor in the calculation area is obtained, the composition and distribution of the safety factor in the area can be mastered. Under three working conditions, Monte Carlo method provides comparable average value and safety factor as traditional calculation methods. However, in the original reservoir area, the safety factor calculated by Monte Carlo method increases abruptly, which is much larger than that of the traditional calculation method. This paper analyses that the safety factor calculated by Monte Carlo method is a closed interval, and the average value is calculated. The safety factor calculated by the Sweden arc method and finite element method is the minimum value of the critical sliding surface. With the progressive increase of the tailings beach in the new reservoir area, the back pressure slope is formed on the outer slope of the dam at the initial stage of the original reservoir area. Therefore, the safety factor of the other three conditions suddenly increases; MC-Swedish arc method and MC-strength reduction method calculate that the failure probabilities of the new and original reservoir areas under three working conditions are 0, which shows that the new and original reservoir areas under this condition have good static stability.

Table 3 shows the calculation results of the traditional method, MC-Swedish arc method, and MC-strength reduction method in the new and original reservoir areas under three working conditions. It can be seen from the table that several methods in the new reservoir area under three working conditions have good consistency, and the calculation error is minor. However, in the original reservoir area, the calculation results of Monte Carlo method are sharply increased and should be compared with other methods in practical engineering. The Monte Carlo method is combined with the Swedish arc method and the finite element strength reduction method to solve the static stability of the new original tailings pond after expansion. This calculation method compensates for the randomness and spatial variability of the tailings and dam material parameters that are overlooked in the traditional calculation method to solve the static stability of the tailings pond and replaces the fixed value with the normal distribution of the calculated material parameters to analyze and evaluate the slope stability. Therefore, under the same design calculation model, the calculation results of this method are obviously more realistic.

##### 4.2. Dynamic Stability

In order to systematically study the dynamic characteristics of a tailings reservoir in Yunnan, El Centro wave and Taft wave are selected to carry out the dynamic calculation of the height of the accumulation dam in three new reservoirs under six calculation conditions (confining pressure 800 kPa). Table 4 outlines the precise computation conditions.

###### 4.2.1. Dynamic Displacement Response

The extreme value of dynamic displacement (i.e., the maximum absolute value of dynamic displacement in the whole earthquake duration) is a decisive index for evaluating the dynamic stability of tailings pond. The horizontal displacement, vertical displacement, and total displacement of condition A are introduced as examples, as shown in Figure 8.

**(a)**

**(b)**

**(c)**

The same method is used to simulate other working conditions, and the horizontal displacement, vertical displacement, and total displacement of the six working conditions are counted in Table 5.

From the dynamic displacement nephogram, it can be seen that under the action of seismic waves, the horizontal displacement to the contact between tailings beach and subdam is larger, especially in the original reservoir area; the vertical displacement of the reservoir area is small, and the total displacement is highly consistent with the horizontal displacement. This is because the input acceleration time history curve is horizontal in the dynamic calculation of the tailings pond.

Table 5 demonstrates that the new and old tailings reservoirs have distinct degrees of displacement when subjected to the action of two types of seismic waves. The minimum value of total displacement under the action of El Centro wave is 2.40 cm when the beach surface of the new reservoir area accumulates to the top of the initial dam in the original reservoir area, and the maximum value is 4.87 cm when the beach surface of the new reservoir area accumulates to the top of the subdam. When the total displacement of the reservoir is compared at different points of the tailings pond, it is discovered that the total displacement of the reservoir gradually increases with the increase of the tailings beach under the action of the El Centro wave. The minimum value of total displacement under Taft wave action is 2.02 cm when the new reservoir beach accumulates to the midpoint of the original reservoir dam height, and the maximum value is 2.17 cm when the new reservoir beach accumulates to the original reservoir’s initial dam crest. When the variation law of the total displacement of the reservoir under the action of the two waveforms is compared, it is clear that the reservoir changes under the action of the El Centro wave. The El Centro waveform is closer to the eigenvalue period of the expansion tailings reservoir, yet its safety performance remains within the controllable range.

###### 4.2.2. Acceleration Response

The acceleration extremum and amplification factor of the dam under six working conditions are calculated. The acceleration extremum and amplification factor of condition A are introduced as examples, as shown in Figure 9.

**(a)**

**(b)**

The same method is used to simulate other working conditions. The horizontal acceleration and vertical acceleration under six working conditions are counted, and the corresponding magnification is calculated in Table 6.

The comprehensive analysis of horizontal acceleration and vertical acceleration extremum nephogram and magnification summary table shows that different waveforms have different effects on new and original tailings reservoirs, especially to prevent seismic waves similar to El Centro waveform, which has a greater impact; under the action of earthquake, the overall distribution of acceleration nephogram of tailings reservoir under different working conditions is reasonable, and the extreme value is mainly located at the initial dam of the new reservoir area, but due to the small value, its influence can be ignored; with the increase of the stacking height of the tailings beach in the new reservoir area, the seismic characteristic values of the new and original tailings reservoirs also change, resulting in different effects of the same waveform on different heights. Therefore, it is recommended that in the construction process of the new and original tailings, special attention should be paid to the seismic detection and early warning. With the increase of tailings beach in the new reservoir area, the extreme value of horizontal acceleration gradually increases, which is consistent with the conclusion in [38].

###### 4.2.3. Seismic Stability

The Newmark method introduced above is used to calculate the seismic stability under three working conditions, and the variation law of the minimum safety factor and the sliding surface where the minimum safety factor is located in the seismic process of the reservoir area is obtained so as to comprehensively understand the seismic stability of the double reservoir area in the construction process.

Figure 10 shows the minimum safety factor and critical slip surface ((a), (c), (e), (g), and (i) show El Centro wave in the earthquake of the new and original tailings pond under six working conditions; (b), (d), (f), (h), and (j) show Taft wave). It can be seen from the figure that the safety factors calculated by El Centro and Taft waves are similar; under the action of seismic waves, the safety factor of the original reservoir area of working conditions A and a will drop sharply, so the stability of the original reservoir area of working condition should be paid more attention in practical engineering. With the increase of the new reservoir beach, the back pressure slope will be formed in the original reservoir area, so the stability of the whole tailings reservoir can be restored to a higher level. Under the action of El Centro wave and Taft wave, the critical sliding surfaces under six working conditions have good consistency; the safety factor of El Centro wave is less than that of Taft wave, which further shows that the El Centro waveform is closer to the eigenvalue period of the expanded tailings pond.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

Table 7 lists the calculation results of dynamic stability of new and original tailings reservoirs under different working conditions by the Monte Carlo method. Based on Newmark method, the randomness of tailings and dam construction material parameters is considered. The simulation results are as follows.

It can be seen from Table 7 that when the Monte Carlo method is used to solve the dynamic stability of the new and original tailings reservoirs under different working conditions, the same situation appears as that when the static stability is solved, and there is a sudden increase in the calculated value of the safety factor at the original reservoir area of working conditions A, B, a, and b. The reasons for the above analysis have been explained in this paper, which is not repeated here. The calculation results of Monte Carlo method and Newmark method show that the new and original tailings reservoirs have high stability under El Centro wave and Taft wave under different working conditions. When the Monte Carlo method is used to calculate the dynamic and static stability of the double reservoir area, the calculation results of the new reservoir area are in good agreement with those of other methods, whereas the calculation value of the original reservoir area is higher. Attention should be paid to the comprehensive comparison with other methods in practical engineering application.

#### 5. Discussion

The limitation of this study is that the linkage effect of tailings and dam material parameters with dam instability damage is not considered. It has been shown that when a slope has a tendency to destabilize, parameters such as soil pore space and permeability coefficient will increase with time, and then the values of parameters such as soil cohesion and internal friction angle and elastic modulus will decrease, and eventually the soil will be destabilized due to a vicious cycle. Therefore, the subsequent study can solve the stability of tailing ponds from the perspective of destabilization damage and material parameter degradation linkage, which will be more in line with the engineering reality.

#### 6. Conclusions

In this study, we used the traditional dynamic and static stability solution method and Monte Carlo method to carry out the stability numerical calculation of the double reservoir area (three static conditions and six dynamic conditions) generated by the new expansion method. The results of the calculations indicate that the new expansion method is highly feasible and stable, and it will form a reverse pressure slope at the initial dam of the original reservoir area while significantly increasing the capacity of the tailings reservoir, which is more conducive to the stability of the tailings reservoir.

The traditional dynamic and static stability methods (Newmark method and Swedish arc method) disregard the randomness and spatial variability of tailings material when calculating the stability of a tailings pond and use the critical sliding surface in the calculation area as an index to evaluate the stability of the region, which makes it difficult to meet the engineering needs at this stage. The Monte Carlo method replaces the fixed value with the tailings material parameters obeying a specific distribution in the calculation process and takes into account the randomness of the calculation material parameters, and therefore it can exclude other influencing factors when calculating the failure probability of the tailings pond under different working conditions by sampling. Under the condition that the required number of samples and extraction times are met, Monte Carlo technique accuracy can also be guaranteed. Therefore, the Monte Carlo approach has greater reference importance and practical significance than traditional methods.

Monte Carlo approach predicts similar stability of the new reservoir area as other traditional methods, indicating its feasibility in solving the stability of the new reservoir area under the condition of double reservoir areas. However, this method may overestimate the stability of the original reservoir area, and therefore comparisons with other calculation methods are required.

#### Data Availability

All data supporting the findings in this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.