Dynamic response and stability analysis of tailings storage facilities under seismic waves

Abstract

Earthquakes are one of the primary causes of instability in tailings storage facilities (TSFs), highlighting the necessity of studying dynamic responses and stability analysis under seismic wave action. This paper investigates the dynamic characteristics of a TSF proposed for decommissioning under seismic loading using the finite element time-history method. The study comprehensively analyzes dynamic deformation, seismic-induced permanent deformation, and pore water pressure changes, while identifying potential liquefaction zones and slope safety factors after seismic events. Results show that stress levels in the tailings dam remain low under applied dynamic loads, indicating a significant safety margin. Liquefaction zones are mainly concentrated in the shallow surface layers of the tailings beach, particularly in the fine sand and silt layers, with no liquefaction observed in the dam slope. The dam body’s deformation accumulates over the duration of seismic shaking, reaching a final value of 11.95 cm. The minimum slope safety factors during seismic events meet regulatory requirements, confirming dam slope stability. The findings provide a scientific basis for the design of TSF decommissioning projects.

Keywords

tailings dam time-history analysis dynamic characteristics seismic waves liquefaction evaluation

Introduction

Tailings storage facilities (TSFs) are sites used to store tailings discharged after ore processing in metal and non-metal mines.¹ The number and height of TSFs have been rapidly increasing to accommodate the growing demand for mineral resources and the resulting rise in tailings production. Among various hazards, earthquakes are one of the primary causes of tailings dam instability.^2,3 Due to the fine particle size and high specific gravity of tailings dam materials, these structures may appear stable under normal conditions but are highly sensitive to external disturbances.⁴ This sensitivity is particularly pronounced during high-intensity seismic events, which may trigger liquefaction and destructive deformation. For instance, the 2011 Tōhoku Earthquake in Japan (magnitude 9.0 on the Richter scale) caused structural instability in several tailings storage facilities, leading to tailings leaks that contaminated nearby farmland and water resources.^5–7

The dynamic response of tailings dams under seismic wave action is highly complex. It involves changes in acceleration, displacement, and stress, all of which directly impact their stability. Therefore, research on the dynamic response and stability of TSFs under seismic loading is essential.⁸ Scholars worldwide have extensively studied this issue using methods such as numerical simulation, physical modeling, and field monitoring.^9–11 Among these, numerical simulation has gained widespread application in the stability analysis of tailings dams due to its cost-effectiveness, efficiency, and high reproducibility. Numerical simulation methods primarily include the traditional pseudo-static method and the time-history analysis method.^12,13 The latter offers distinct advantages by accounting for seismic wave characteristics and integrating the dynamic properties of tailings materials (e.g., dynamic strength, pore pressure, shear modulus, and damping ratio), ensuring more representative simulation results.^14,15 Yin analyzed the permanent deformation and dynamic stability of a high tailings dam using the time-history analysis method, the global deformation analysis method, and the limit equilibrium method, with Xiaodae tailings pond as an example.¹⁶ Focusing on the Dashawo Tailings Reservoir heightening and expansion project, Wang adopted the time-course power analysis method to systematically study and comprehensively evaluate the dam’s dynamic response and stability before and after the expansion.¹⁴

Therefore, this study focuses on a proposed decommissioned tailings storage facility as the research object and employs the finite element time-history method to investigate the dynamic characteristics of the dam under seismic wave action. A comprehensive analysis is conducted on the dynamic deformation characteristics of the dam body, seismic-induced permanent deformation, and dynamic pore water pressure. Additionally, the study identifies potential seismic liquefaction zones and evaluates post-seismic slope safety factors, offering a scientific basis for designing TSF decommissioning projects.

Engineering overview

The tailings storage facility of a certain iron mine, which was designed in 1992 and put into operation in 1994, has now ceased operation. It is a valley-type tailings storage facility, surrounded by mountains on three sides. The surrounding mountainous terrain, characterized by dense vegetation and simple geological structures, is stable. No major faults, landslides, or debris flows have been identified in the area.

Geotechnical surveys and engineering assessments indicate that the dam crest has reached its highest elevation of 336 meters. This corresponds to a maximum dam height of 86 meters. The reservoir’s total storage capacity is 4.62 million cubic meters. With an effective storage coefficient of 0.75, this equates to an effective capacity of 3.465 million cubic meters. This capacity classifies the reservoir as a third-class facility.

The original structure was a permeable rock-fill dam, measuring 27 meters in height and 135 meters in crest length, with an outer slope ratio of approximately 1:1.75 and a crest width of 3.5 meters. The dam’s expansion utilized upstream tailings embankment methods, with sub-dam slope ratios ranging from 1:4.0 to 1:4.3. The extended dam section has an average outer slope ratio of 1:4.5 and a final crest length of 301 meters. The layout and design of the tailings pond are illustrated in Figure 1.

Figure 1.

Plan layout of the tailings storage facility.

Dynamic finite element calculation theory

Dynamic constitutive model

The Hardin-Drnevich model is employed to calculate the dynamic shear modulus and damping ratio, which serve as the foundation for dynamic analysis in finite element calculations.¹⁷

The model includes the following relationships:

a. Dynamic Modulus:

G = \frac{G_{\max}}{1 + γ / γ_{r}}

(1)

b. Damping Ratio:

λ = λ_{\max} \frac{γ / γ_{r}}{1 + γ / γ_{r}}

(2)

c. Maximum Shear Modulus:

G_{\max} = K_{2} P_{a} {({σ^{'}}_{m} / P_{a})}^{n}

(3)

where

γ_{r}

is the reference shear strain,

γ_{r} = τ_{\max} / G_{\max}

;

σ_{m}^{'}

is the mean effective stress;

P_{a}

is the atmospheric pressure; and

K_{2}

and

n

are experimental parameters.

These equations provide the relationships needed to model the stress-strain behavior of materials under dynamic loading conditions. The Hardin-Drnevich model is particularly suitable for representing soil behavior in seismic analyses, as it captures the nonlinear stress-strain response and energy dissipation through damping.

Dynamic control equation and solution steps

The dynamic control equation is expressed as:

[M] {\ddot{δ} (t)} + [C] {\dot{δ} (t)} + [K] {δ (t)} = {F (t)}

(4)

where

δ (t)

\dot{δ} (t)

, and

\ddot{δ} (t)

represent the node displacement, velocity and acceleration, respectively;

F (t)

is the nodal dynamic load, determined by the seismic acceleration;

[M]

is the mass matrix;

[K]

is the stiffness matrix;

[C]

is the damping matrix.

The damping matrix $[C]$ is assumed using Rayleigh damping:

[C] = λ ω [M] + \frac{λ}{ω} [K]

(5)

where

ω

is the fundamental frequency and

λ

is the damping ratio.

To solve the dynamic control equation, the Wilson-θ method is applied. The main calculation steps are as follows:

Step 1: Static Nonlinear Calculations

Perform static nonlinear calculations to determine the static stress distribution of each element before the earthquake.

Step 2: Seismic Process Segmentation

Divide the seismic process into multiple intervals based on the amplitude of seismic waves; assume that the dynamic parameters of each element, specifically the shear modulus $G$ and damping ratio $λ$ , remain constant within each interval; use the results from static analysis and the dynamic constitutive model to calculate the initial values of $G$ and $λ$ for each element at the start of the dynamic analysis.

Step 3: Time-Step Calculations

Use a time step $Δ t$ = 0.01∼0.02 seconds within each interval; Apply the Wilson-θ method to solve the dynamic control equation (4) and calculate the dynamic shear strain $γ$ process for each element in the interval; Take 65% of the maximum shear strain ( $0.65 \cdot γ_{\max}$ ) during the interval as the average shear strain $\bar{γ}$ for the element.

Step 4: Iterative Parameter Update

Based on the calculated average shear strain $\bar{γ}$ , update the shear modulus $G_{i}$ and damping ratio $λ_{i}$ for each element using the dynamic constitutive model; iterate within the current period until the updated shear modulus and damping ratio meet the required accuracy; use the final $G_{i}$ and $λ_{i}$ values as the starting shear modulus $G_{i + 1}$ and damping ratio $λ_{i + 1}$ for the next interval; calculate the fundamental frequency of the tailings dam for use in the next interval’s iterations.

Step 5: Dynamic Response Calculation

Repeat Steps 3 and 4 for each interval until the end of the seismic process is reached; obtain the required dynamic response quantities for the entire earthquake.

Calculation method of dynamic pore water pressure

The calculation formula for dynamic pore water pressure is as follows:

Isotropic Pressure Consolidation:

\frac{u_{d}}{σ_{3 c}^{'}} = \frac{2}{π} \arcsin {(\frac{N}{N_{L}})}^{\frac{1}{2 α}}

(6)

Anisotropic Pressure Consolidation:

\frac{u_{d}}{σ_{3 c}^{'}} = \frac{2}{π} (1 - m \frac{τ_{0}}{σ_{0}^{'}}) \arcsin {(\frac{N}{N_{L}})}^{\frac{1}{2 α}}

(7)

Differentiating equation (7) gives the incremental calculation formula for pore water pressure:

Δ u_{d} = \frac{σ_{0}^{'} (1 - m \frac{τ_{0}}{σ_{0}^{'}})}{α π N_{L} \sqrt{1 - {(\frac{N}{N_{L}})}^{\frac{1}{α}}}} {(\frac{N}{N_{L}})}^{\frac{1}{2 α} - 1} Δ N

(8)

where

τ_{0}

represents the initial shear stress;

m

is an empirical coefficient reflecting the decrease in pore pressure, typically ranging between 1.1 and 1.3;

Δ N

denotes the incremental equivalent number of vibration cycles within the specified time period; and

N_{L}

represents the total number of vibration cycles required to reach liquefaction.

$N_{L}$ can be obtained either from liquefaction test curves or calculated using the following semi-logarithmic formula.

N_{L} = 10^{\frac{b - \bar{τ} / {σ^{'}}_{0}}{a}}

(9)

where

a

and

b

are empirical constants associated with the soil’s anti-liquefaction properties;

\bar{τ}

represents the average shear stress, typically taken as 0.65 times the maximum shear stress within the period.

Calculation method for seismic permanent deformation

The assessment of seismic permanent deformation on slopes is conducted using the Newmark methodology,^18,19 as shown in Figure 2. In scenarios where the potential sliding plane forms a circular arc, the sliding moment caused by the mass under static conditions is calculated using the following formula:

W b = R \sum τ d s

(10)

Figure 2.

Newmark sliding block method.

Under dynamic conditions, the sliding moment of the mass must take into account the effect of seismic loads. The sliding moment in this case is expressed as:

M_{D} (t) = W b + k W h = R \sum S_{q} d s

(11)

where

k

is the seismic response coefficient, defined as

k = \frac{a}{g}

From equations (9) and (10), the following can be derived:

k = \frac{b}{h} (\frac{\sum S_{q} d s}{\sum τ d s} - 1) = \frac{b}{h} (\frac{\bar{s}}{\bar{τ}} - 1)

(12)

At a specific moment during an earthquake, when $h = d$ , the value of $k$ reaches its minimum. Upon reaching the critical sliding state, that is, when k = k_y, the resistance moment becomes equal to the sliding moment. At this point, the resistance moment is:

M_{R} (t) = W b + k_{y} W h

(13)

Throughout the earthquake event, the state of the sliding mass can be categorized based on the value of $k$ relative to $k_{y}$ :

The mass reaches a limit equilibrium state when $k = k_{y}$ ;

The mass remains stable when $k < k_{y}$ ;

The mass initiates sliding when $k > k_{y}$ .

Assuming the rotational angular acceleration of the sliding mass is $\ddot{θ}$ , the dynamic equilibrium equation governing the behavior of the sliding mass can be expressed as follows:

I \ddot{θ} = M_{D} (t) - M_{R} (t)

(14)

where

I

is the rotational inertia of the sliding mass, defined as

I = W d^{2} / g

;

d

is the distance from the center of gravity of the sliding body to the center of the circle; and

W

is the total weight of the sliding body.

From equation (13), the angular acceleration $\ddot{θ}$ is determined, and by performing a second integration over time, the rotational angular displacement of the sliding body can be calculated:

θ = \iint \ddot{θ} d t d t

(15)

Once the rotational angular displacement of the sliding mass is determined, multiplying it by the radius of the sliding arc gives the rotational arc length at any given point along the sliding surface. By applying the geometric principles of circular motion, the corresponding horizontal and vertical displacements for any point on the sliding arc can then be accurately calculated.

Calculation conditions and parameters

Dam body calculation section and finite element division

Based on the engineering survey report, the cross-sectional configuration and tailings zoning of the main dam at its final elevation of 336 meters were determined, as shown in Figure 3. The finite element mesh for the tailings dam cross-section was primarily generated using four-node quadrilateral elements, with triangular elements used locally for smooth transitions. The finite element computational mesh is shown in Figure 4, consisting of 3452 nodes and 3284 elements.

Figure 3.

Tailings dam calculation section diagram.

Figure 4.

Finite element calculation grid of tailings dam.

Static calculation parameters

The nonlinear stress-strain relationship of the dam construction material is characterized using the Duncan E-B model. The selection of static calculation parameters was based on comparable engineering projects, with the specific values presented in Table 1.^20–22

Table 1.

Duncan E-B model calculation parameters.

Parameter name	Density, g/cm³	C, kPa	φ,°	K	n	R _f	K _b	m
Initial dam	2.05	0	33	800	0.36	0.72	320	0.25
Silt fill layer ①	1.95	10	16	160	0.50	0.81	70	0.21
Tail fine sand layer ②	1.75	0.2	26.5	217	0.07	0.87	41	0.33
Tail silt layer ②₁	1.95	5	18	180	0.50	0.81	80	0.21
Tail silty clay layer ②₂	1.95	12	12	160	0.50	0.81	70	0.21
Tail fine sand layer ③	1.80	0.2	28	160	0.19	0.82	41	0.21
Tail silt sand layer ③₁	1.85	0.3	26.5	300	0.36	0.72	240	0.25
Tail silt layer ④	1.96	6	23	180	0.50	0.81	80	0.21
Tail silty clay layer ④₁	2.00	6	13	160	0.50	0.81	70	0.21

Seismic acceleration time history curve

At the project site, the peak ground acceleration is recorded as 0.15 g, corresponding to a seismic fortification intensity of VII degrees and a characteristic seismic response spectrum period of 0.45 s. Three groups of artificially fitted seismic waves are used in the calculation, with all amplitudes taken as 0.15 g. When calculating the dynamics of the tailings dam, both horizontal and vertical acceleration curves are input simultaneously, with the peak value of the vertical acceleration curve taken as 2/3 of the peak value in the horizontal direction. The calculation results show that the contour distribution of the different seismic waves is similar, with only a slight difference in the extremes of the physical quantities. Due to space limitations, this paper only presents the results of analyzing and calculating a group of seismic waves with larger extreme values of physical quantities.

Figures 5 and 6 display the synthesized seismic waves, the standard spectrum curve, and the seismic acceleration time history curves.

Figure 5.

Artificially synthesized seismic waves and standard spectrum curve.

Figure 6.

Artificially synthesized seismic acceleration time history curve. (a) Horizontal, (b) vertical.

Dynamic calculation parameters

Maximum dynamic shear modulus

For the test soil material, the dynamic shear modulus at a shear strain of 10⁻⁶ is considered as the maximum dynamic shear modulus

G_{\max}

. The relationship between the dynamic shear modulus

G_{\max}

and the mean effective stress

σ_{0}^{'}

is linear when plotted on a double logarithmic scale, following a power function form:

G_{\max} = K P_{a} {({σ^{'}}_{0} / P_{a})}^{n}

(16)

where

σ_{0}^{'}

represents the mean effective stress;

P_{a}

is the atmospheric pressure; and

K

and

n

are experimental parameters, as shown in Table 2.

Table 2.

Maximum dynamic shear modulus relationship with confining pressure.

Soil sample name	Confining pressure/kPa	$G_{\max}$ /MPa
Tail fine sand	100	26.8
	200	27.4
	300	31.7
Tail silt	100	17.6
	200	29.0
	300	38.7

Relationship between dynamic shear modulus ratio, damping ratio, and dynamic shear strain

According to the dynamic test report, the test results of the dynamic shear modulus ratio $G / G_{\max}$ and damping ratio $λ$ for tail silt and tail fine sand are shown in Figure 7.

Figure 7.

Relationship between dynamic shear modulus ratio, damping ratio, and dynamic shear strain. (a) Relationship curve between $G / G_{\max} - γ_{d}$ and $λ - γ_{d}$ for tail silt at $σ_{3}$ = 100 kPa, (b) relationship curve between $G / G_{\max} - γ_{d}$ and $λ - γ_{d}$ for tail silt at $σ_{3}$ = 200 kPa, (c) relationship curve between $G / G_{\max} - γ_{d}$ and $λ - γ_{d}$ for tail silt at $σ_{3}$ = 300 kPa, (d) relationship curve between $G / G_{\max} - γ_{d}$ and $λ - γ_{d}$ for tail fine sand at $σ_{3}$ = 100 kPa, (e) relationship curve between $G / G_{\max} - γ_{d}$ and $λ - γ_{d}$ for tail fine sand at $σ_{3}$ = 200 kPa, (f) relationship curve between $G / G_{\max} - γ_{d}$ and $λ - γ_{d}$ for tail fine sand at $σ_{3}$ = 300 kPa.

Finite element result analysis

Static calculation results

Figure 8 presents the isocontour map showing the static stress and stress levels within the tailings dam. The pattern of stress distribution and the stress level isocontours under static conditions appears logical and consistent. The dam body’s maximum vertical stress is recorded as 998.39 kPa, excluding the foundation. These stress isocontours, largely parallel to the dam’s slope, indicate that the dam body’s vertical stress predominantly results from the weight, or height, of the “soil column.” At the base of the upstream side of the dam, the highest shear stress is observed at 130.12 kPa, again excluding the foundation. With the maximum stress level in the dam body at 0.83, there is no indication of a risk of shear failure or plastic deformation. Notably, areas with extreme values are primarily located near the dam’s bottom, away from the external slope, suggesting that the slope stability of the dam is not compromised.

Figure 8.

Isocontour map of static stress and stress level in tailings dam (stress unit: kPa). (a) Vertical stress, (b) shear stress, (c) stress level.

Dynamic calculation results

Typically, the seismic acceleration time history curve reflects the surface conditions of the dam bedrock or the site itself. When directly applied to dynamic analyses of the dam grid, including the bedrock, this can lead to an overestimation of the dam’s seismic response due to the amplification effect of the bedrock. Additionally, it is observed that during seismic events, the deformation experienced by the dam body significantly exceeds that of the bedrock. Consequently, the seismic response analysis model is designed to focus exclusively on the area encompassing the tailings dam.

Dynamic displacement response

Figure 9 shows the isocontour map of dynamic displacement of the tailings dam under seismic dynamic loads.

Figure 9.

Isocontour map of dynamic displacement in tailings dam (unit: cm). (a) Horizontal, (b) vertical.

Figure 10 shows the locations of typical nodes at elevations 227 m (Point A), 312.1 m (Point B), and 336.3 m (Point C). Figure 11 presents the corresponding time history curves of dynamic displacement at these nodes under seismic dynamic load.

Figure 10.

Typical node locations of the dam cross-section.

Figure 11.

Time history curves of dynamic displacement at typical nodes. (a) Horizontal dynamic displacement at Point A, (b) vertical dynamic displacement at Point A, (c) horizontal dynamic displacement at Point B, (d) vertical dynamic displacement at Point B, (e) horizontal dynamic displacement at Point C, (f) vertical dynamic displacement at Point C.

In the context of seismic activity, the tailings dam exhibits a pattern in the distribution of dynamic displacement isocontours, with dynamic displacement values gradually increasing from the dam’s base to its apex. The analysis also reveals that horizontal displacement predominates over vertical displacement. Specifically, the peak values of horizontal and vertical dynamic displacements recorded for the entire structure are 6.40 cm and 1.06 cm, respectively, under the influence of seismic waves.

Acceleration response

Figure 12 shows the acceleration isocontour map of the dam body, representing the maximum absolute value of acceleration throughout the entire seismic duration.

Figure 12.

Acceleration isocontour map of tailings dam (unit: m/s²). (a) Horizontal, (b) vertical.

Depicted in Figure 13 are the acceleration time history curves for selected nodes at various elevations: 277.0 m (Point A), 312.1 m (Point B), and 336.3 m (Point C), under seismic dynamic loading. The precise locations of these nodes are shown in Figure 10.

Figure 13.

Time history curves of acceleration at typical nodes. (a) Horizontal acceleration at Point A, (b) vertical acceleration at Point A, (c) horizontal acceleration at Point B, (d) vertical acceleration at Point B, (e) horizontal acceleration at Point C, (f) vertical acceleration at Point C.

During seismic activity, the tailings dam exhibits a logical distribution of acceleration isocontours under various conditions, with the highest values predominantly near the dam’s crest. The dam experiences a maximum horizontal seismic acceleration of 4.06 m/s², with an amplification factor of 2.71. The maximum vertical acceleration is 2.64 m/s², with an amplification factor of 2.64.

Analysis of liquefaction potential

Figure 14 shows the distribution of dynamic pore pressure and liquefaction zones in the tailings dam under seismic action, with liquefaction areas highlighted in red.

Figure 14.

Isocontour map of dynamic pore pressure and liquefaction zones (pore pressure unit: kPa). (a) Dynamic pore pressure, (b) liquefaction zones.

During seismic wave activity, the peak dynamic pore pressure in the dam reaches 92.25 kPa. Liquefaction primarily occurs in the shallow layers of the beach area (consisting of fine sand and silt) within the tailings dam, with no liquefaction observed on the slope. This indicates that the overall structure of the tailings dam remains secure.^23–25

Permanent deformation

Figure 15 shows the distribution of cumulative permanent deformation in the tailings dam under seismic wave action.

Figure 15.

Cumulative curve of seismic permanent deformation of tailings dam.

During seismic events, the tailings dam experiences cumulative permanent deformation. When the acceleration of soil elements exceeds the yield acceleration, incremental residual deformations accumulate, reaching their peak by the end of the earthquake. Under seismic wave activity, the tailings dam undergoes permanent deformation of 11.95 cm.

Post-earthquake analysis of the dam body’s anti-slide stability

Figure 16 illustrates the change in the anti-slide stability safety factor of the tailings dam slope over the course of the earthquake, influenced by seismic wave action.

Figure 16.

Safety factor for anti-slide stability of tailings dam slope.

The sliding safety factor is calculated using the Bishop method, as expressed by the following formula:

F = \frac{\sum (σ_{n i} \tan φ + c) d l_{i}}{\sum τ_{i} d l_{i}}

(17)

where

σ_{n i}

represents the effective normal stress on soil slice

i

;

φ

denotes the friction angle of the soil;

c

is the cohesion of the soil;

τ_{i}

refers to the tangential reaction force acting on soil slice

i

; and

d l_{i}

represents the tangential length of the soil slice

i

During seismic wave activity, the minimum safety factor of the dam slope remains at 1.09. This value meets regulatory standards throughout the earthquake, indicating that the dam slope remains stable and resistant to sliding under seismic conditions.^26,27

Conclusions and recommendations

The nonlinear finite element static and dynamic analyses of the tailings dam project lead to the following conclusions:

Static conditions: Stress patterns and levels in the dam exhibit a generally rational distribution. The highest vertical stress reaches 998.39 kPa, with isocontours aligned parallel to the slope. The maximum shear stress, located near the upstream dam bottom, is 130.12 kPa. The stress level peak of 0.83 in the dam body suggests a low risk of shear failure, with high-stress regions primarily at the dam’s base and away from the outer slope, ensuring slope safety.

Seismic events: During seismic activity, the tailings dam exhibits a progressive increase in dynamic displacement from the base to the crest. Specifically, under an earthquake with a peak acceleration of 0.15 g, the dam experiences maximum horizontal and vertical dynamic displacements of 6.40 cm and 1.06 cm, respectively. The acceleration isocontours within the dam show a logical pattern, peaking near the crest. The highest horizontal and vertical accelerations are 4.06 m/s² and 2.64 m/s², respectively, with amplification factors of 2.71 and 2.64. The peak dynamic pore pressure reaches 92.25 kPa. Liquefaction primarily occurs in the shallow layers of the beach area (comprising fine sand and silt), but not on the slope, highlighting the dam’s overall safety. The cumulative permanent deformation due to seismic activity is 11.95 cm, and the lowest safety factor of the slope complies with regulatory standards, confirming the dam’s stability.

Despite seismic activity, the tailings dam remains stable. Liquefaction is confined to the shallow surface layers of the beach areas, specifically in fine sand and silt zones. To further enhance safety, seismic reinforcement in these specific areas is recommended.

Statements and declarations

Footnotes

Conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors acknowledge the General Project of the Science and Technology Plan of the Beijing Municipal Commission of Education (KM202310009005), the National Natural Science Foundation of China (No. 52309131, No. 52304202), the Guizhou Science and Technology Plan Project (Qianke Science Support [2023] General 122), the Yuxiu Innovation Project of NCUT (No. 2024NCUTYXCX209), and the Open Project of the Key Laboratory of Xinjiang Coal Resources Green Mining (No. KLXGY-KB2409).

ORCID iD

Songli Jin

References

Fotovvat

Sadrekarimi

. Instability of a gold mine tailings subjected to different stress paths. J Saf Sci Technol 2022; 148(5): 04022020.

Clarkson

Williams

. Critical review of tailings dam monitoring best practice. Int J Min Reclamat Environ 2020; 34(2): 119–148.

Tao

, et al. Experimental study on overtopping dam-break of a tailing reservoir under extreme conditions. Environ Sci Pollut Res 2024; 31(5): 6874–6890.

Lin

Wang

Liu

, et al. Regional distribution and causes of global mine tailings dam failures. Metals 2022; 12(6): 905.

Wang

Yang

Karen

, et al. Status and development for the prevention and management of tailings dam failure accidents. Chin J Eng 2018; 40(5): 526–539.

Koppe

. Lessons learned from the two major tailings dam accidents in Brazil. Mine Water Environ 2021; 40(1): 166–173.

Shahriari

Aydin

. Lessons learned from analysis of Los Frailes tailing dam failure. Adv Saf Manag Hum Factors 2018; 604(32): 309–317.

Guo

Liu

Chen

. Numerical simulation on seismic permanent displacement of clay anti-seepage tailing pond. J Saf Sci Technol 2015; 11(1): 5–10.

Qin

. Experimental study of a tailings impoundment dam failure due to overtopping. Mine Water Environ 2018; 37(2): 272–280.

10.

Ferrer

Sáez

Ledezma

. Numerical modeling of cracking pattern's influence on the dynamic response of thickened tailings disposals: a periodic approach. Earthq Eng Eng Vib 2018; 17(1): 179–190.

11.

Salam

Xiao

Khosravifar

, et al. Seismic stability of coal tailings dams with spatially variable and liquefiable coal tailings using pore pressure plasticity models. Comput Geotech 2021; 132: 104017.

12.

Aswathi

Jakka

. Seismic analysis of hybrid tailings dams: insights into stability and responses. Bull Eng Geol Environ 2024; 83(2): 56.

13.

Naeini

Akhtarpour

. Numerical analysis of seismic stability of a high centerline tailings dam. Soil Dynam Earthq Eng 2018; 107(9): 179–194.

14.

Wang

Yin

Wei

, et al. Study of the dynamic stability of tailings dam based on time-history analysis method. J China Univ Min Technol 2018; 47(2): 271–279.

15.

Liu

Wang

, et al. Feasibility and stability analysis of coal-gangue-based tailings dam on subsidence land. Fresenius Environ Bull 2018; 27(1): 37–45.

16.

Yin

Wang

Wei

, et al. Analysis of the permanent deformation and stability of high tailings dam under earthquake action. Rock Soil Mech 2018; 39(10): 3717–3726.

17.

Wang

Zhang

, et al. Stability analysis of iron ore tailings dams under action of seepage and earthquakes. Geol Explor 2018; 54(3): 614–622.

18.

Nimbalkar

Annapareddy

VSR

Pain

. A simplified approach to assess seismic stability of tailings dams. J Rock Mech Geotech Eng 2018; 10(6): 1082–1090.

19.

Jin

Song

Liang

, et al. Dynamic characteristics of tailings reservoir under seismic load. Environ Earth Sci 2018; 77(18): 654.

20.

Liu

Sun

, et al. Variations in permeability and mechanical properties of basaltic rocks induced by carbon mineralization. Sustainability 2022; 14(14): 15195.

21.

Liu

Sun

Liu

, et al. Damage properties of basaltic rocks under multilevel stress loading and microwave irradiation. Eng Fract Mech 2023; 289(289): 109391.

22.

Sun

Liu

. Uniaxial compressive strength determination of rocks using X-ray computed tomography and convolutional neural networks. Rock Mech Rock Eng 2021; 54(54): 4225–4237.

23.

Yuan

Xiong

Fang

, et al. Analysis of the working state and stability safety of No. 2 mud dam in Xiaomaliu Tailing Pond. Eng J Wuhan Univ 2022; 55(2): 123.

24.

Liu

Nagula

Jostad

, et al. Considerations for using critical state soil mechanics based constitutive models for capturing static liquefaction failure of tailings dams. Comput Geotech 2024; 167: 106089.

25.

Zeng

Liu

, et al. Enhanced method of slices with experimental evidence for predicting static liquefaction on tailings dams. Bull Eng Geol Environ 2024; 83(6): 227.

26.

Sun

Liu

You

, et al. Yielding performance of an intelligent terminal structure for cable bolts. J Intell 2024; 2(2): 9180034–9180114.

27.

Sun

Liu

, et al. Mechanical testing and numerical simulation of intelligent terminal structure of rockbolt used as a deformable support and for safety monitoring in rock engineering. Tunn Undergr Space Technol 2024; 148: 105769.