Thanikella Vijayasri

Department of Civil Engineering, Indian Institute of Science, Bangalore, Karnataka 560012, India

Abstract: The seismic response analysis of a tailing dam is studied using a fully coupled effective stress approach in conjunction with an advanced multi yield surface plastic constitutive model for tailing material. Strain controlled static and cyclic triaxial tests were carried out to obtain the constitutive model for the tailing material. The tailing materials were collected from the Rampura Agucha tailing dam (Rajasthan State, India). A 2D nonlinear finite element (FE) model was then developed using different boundary conditions from the tailing embankment constructed using the downstream and upstream method of rising using OpenSees software. In first case, the model boundary was fixed in both the X and Y directions, and in the second case, viscous dashpots were introduced for both side and horizontal boundaries. The model was validated with experimental results on tailing material. Analyses were carried out considering five different earthquake motions, which were applied at the base. Comparisons of the different boundary conditions in terms of displacement flow vectors, pore pressure and stress-strain curves during shaking are presented. From the analysis, it was observed that the viscous boundary condition replicates the actual field conditions more accurately than the fixed boundary condition. In addition, it was found that the tailing embankment constructed by the downstream and upstream method of rising is not susceptible to liquefaction and lateral spreading for earthquake motions, even for a magnitude > 5.5.

Keywords: tailing embankment; cyclic triaxial test,; seismic analysis; liquefaction potential; transmitting and fixed boundary conditions

1 Introduction

Tailing dams are structures established to accommodate the waste tailings generated during mining projects. Three different methods are used for the construction of tailing dams: the upstream method,downstream method, and centreline method. The most suitable method is selected based on the availability of land, characteristics of the settled tailings, and the cost of construction. Most mine tailings are considered toxic in nature and are contained in a closed environment. In the case of failure of the tailing dam, the loss of control on the flooding of the liquid mine tailings pollutes the ecosystem in the downstream (D/S) areas, which consequently results in health risks and numerous social and economic crises in the vicinity. In fact, tailing dams are highly susceptible to instability problems when compared to water storage dams, which is attributable to causes like multi-stage raising, locally derived fill sand,lack of regulations, etc. (Ricoet al., 2008a). Over the past decades, the world has seen many disastrous accidents related to the large tailings and public attention has been drawn to the safety of tailing dams (Ricoet al., 2008b).The primary causes for the failure of tailing dams can be categorized as follows: too many elements like piping,slope instability, seismic liquefaction, excessive rainfall,foundation failure insufficient freeboard, seepage and rise in the phreatic surface, etc. Ricoet al. (2008b)studied the reasons for the failure of the dams based on a database of 147 cases of tailing dam failures around the globe, which he established. It was found that 85% of the failures were seen in active tailing dams.Meteorological factors such as uneven rainfall and snowfall are among the key causes of failure. Note also that the 66% of the failed tailing dams were constructed using the upstream (U/S) method of embankment rising.Hence, a comprehensive static and seismic analysis is necessary before the construction of any future tailing embankments.

In the past, studies on the seismic response of earth fill dams and embankments have been followed by simplified numerical methods such as the pseudo-static method of analysis and aimed to calculate the critical factors of safety that could cause instability of a dam slope. Later, sliding block analysis and Newmark′s block method are used to evaluate permanent displacements as a function of the critical acceleration. In addition,simplified numerical methods such as ground response analysis based on the equivalent viscoelastic approach has been reported (Schnabelet al., 1972; Idriss and Boulanger, 2010). Sitharam and Hegde (2016) carried out comprehensive stability analysis on a rock fill tailing dam and found that blasting activities did not affect its safety. More advanced shear beam analyses were developed to investigate the dynamic response of dams (Gazetas, 1987). Fully coupled effective stress formulation (Biot, 1941; Zienkiewiczet al., 1999) has also gained tremendous attention in the area of seismic investigations of earthen dams. Furthermore, reliabilitybased methods that consider the spatial variability of the mine tailings and the soil (e.g., Liet al., 2014; Sitharam and Hegde, 2019) have been investigated.

Subsequent finite element (FE) studies concentrated on the analysis of earth dams (Dakoulas, 2012;Rampelloet al., 2009) and elasto-plastic behavior of the soil materials (Prevost, 1978; Yanget al., 2003) using constitutive models of various levels of complexity. The FE method is considered to be one of the most powerful tools to analyze the response of dams and embankments during earthquakes. During the past few decades,nonlinear constitutive models have been implemented into numerical procedures to improve prediction of the hysteresis loop behavior of dams during seismic loading conditions (Elgamalet al., 2002; Eliaet al., 2011)

The present study has two parts. In the first part, cyclic triaxial experiments were conducted on tailing materials taken from the Rampura Agucha (RA) mine in the Bhilwara District of Rajasthan, India. Tests were carried out at varying cyclic axial-strain amplitudes (0.4%,1.0%, 1.5% and 2%), and cell pressures (70 kPa and 100 kPa) at a loading frequency of (0.3 Hz). In the second part, a robust nonlinear FE model of a tailing dam was developed, its seismic performance and liquefaction potential were assessed during various seismic motions relevant to its seismicity regions. Nonlinear backbone curves (Gsecandξ) developed experimentally were implemented in the 2D finite element program OpenSees through MATLAB subroutine codes. Emphasis was given to various response parameters in terms of pore pressure, hysteresis loop behavior, peak displacement,acceleration profiles, cyclic stress ratio and liquefaction potential with varying boundary conditions constructed by upstream and downstream methods of rising.

2 Details of tailing dam and materials

In this study, a tailing dam has been chosen for analysis purpose. It is located in the Rampura Agucha(RA) mine in the Bhilwara District of Rajasthan, India.The RA mine is the world′s largest zinc mine with a production capacity of 6.15 million MT per annum. The mine produces zinc as a primary product and lead as a secondary product. The tailing pond has an approximate plan area of 1200 m × 800 m shown in Fig. 1(a). The tailing dam was constructed incrementally in different phases using waste rocks generated from the mining projects. The typical cross section of the tailing dam constructed using the upstream and downstream method of rising is shown in Fig. 1(b) and 1(c).The total height of dam is about 57 m above the natural ground level.The dam was setup on low permeable soil underlain by rocky strata. The U/S slope is 2H:1V and the D/S slope is 36° with the horizontal (angle of repose of the waste muck). Tailing materials for laboratory testing have been collected at 3.5 m below the top of the crest (443.5 m RL). Index properties and engineering properties of the tailing material calculated from the experimental tests are shown in Table 1(b).

The RA mine tailing dam is composed of mine muck,impervious soil and tailing material. The main dam composed of waste mine muck with large particle size distribution is not susceptible to liquefaction. Distribution of different sizes and weights of the mine muck are given in Table 1(a). A 1.5 mm thick high-density polyethylene(HDPE) geomembrane has been used in earlier risings(after phase III) to the U/S slope. The average friction angle of the roughened geomembrane is 37o. Figure 1 shows the details of the U/S slope. The soil cover used along with the geomembrane on the upstream side is a cohesive soil which is only 0.5 mm thick and acts as a liquefaction resistant material.

Table 1(b) Index properties of tailing material

Fig. 1 (a) Plan view of the Rampur Agucha mine; (b) schematic diagram of tailing dam constructed by downstream and (c) upstream rising conditions along with geomembrane layer shown in Fig. 1(b)

3 Experimental investigations

To investigate the physical properties and engineering properties of the tailing samples, specific gravity, gradation curve, Atterberg limits, vibratory table, compaction, consolidated undrained (CU) triaxial shear tests, consolidation and permeability tests have been carried out. Samples were collected at 3.5 m below the top of the crest (443.5 m RL). Tests were conducted according to the Bureau of Indian Standards (1980a;1980b; 1981; 1985). The cyclic triaxial tests were carried out under strain controlled conditions to determine the dynamic properties and liquefaction potential of tailing samples based on ASTM (1992; 1996).

3.1 Material characterization

The tests described above were performed in the laboratory to determine their physical properties and engineering properties. From the results, it is seen that the specific gravity of tailing material was found to be 2.74,which is higher than average. Due to its higher specific gravity, dry density of tailing material is increased, which makes it suitable as a backfill material for retaining structures and earthen dams. From gradation analysis,the percentage of sand, silt, and clay fraction in the tailing material was found to be 51.12%, 48.88%, and 0,respectively. The maximum dry unit weight and optimum moisture content of the tailing material was 18 kN/m3and 22%, respectively. The index properties of tailing materials yielded from various tests are shown in Table 1(b).Other material properties such as mine muck, impervious soil and foundation soil, which were used in the tailing dam are shown in Figs. 1(b) and 1(c). These material properties are taken from the report supplied by M/s Hindustan Zinc Limited (Vedanta) Rampura-Agucha Mine, Rajasthan.

3.2 Test program and methodology

To investigate the effective parameters, dynamic properties and number of cycles for liquefaction evaluation of the tailing samples, consolidated undrained static and cyclic triaxial tests were performed in the laboratory. Sample preparation was carried out on the tailing material as per the procedure given in the Bureau of Indian Standards (1981). Cylindrical specimens of 50 mm × 100 mm were prepared targeting a maximum relative density of 65%. The dry deposition method was used for the sample preparation. The sample was prepared in three layers, and porous stone and filter paper were placed on top and bottom of the samples.The assembly was sealed with an O-ring. The saturation process was continued until B (pore pressure parameter)reached beyond 0.98. After the saturation process,consolidation was carried out and volume change was recorded. Finally, shearing was carried out under the desired effective cell pressures of 70 kPa, 100 kPa and 150 kPa.

To determine the dynamic properties and liquefaction potential of the tailing samples, cyclic triaxial tests were performed under strain controlled conditions on 50 mm ×100 mm samples. Shearing has been carried out at cell pressures of 70 kPa, 100 kPa and 150 kPa, strain amplitudes of 0.5%, 1.0%, and 1.5% and at a loading frequency of 0.3 Hz. From the gradation curve, it is observed that the tailing material is well graded sand with fines. Generally, fine particles require sufficient time for pore water pressure equalization, taking into consideration samples tested at a frequency of 0.3 Hz.

4 Experimental results and discussions

4.1 Static triaxial tests

When the sample is subjected to shearing under the desired confining stress, pore pressure may not dissipate under a given drainage path and time. In such situations,static triaxial consolidated undrained tests play a main role. A series of static CU tests have been carried out by considering the effect of various cell pressures (70 kPa,100 kPa and 150 kPa). All tests were kept at a relative density of 65% and at a strain rate of 0.6 mm/min was maintained. When soils are tested at 85% relative density, they show dilative behavior where liquefaction did not occur. Therefore, the relative density of 65%(medium sand) was used in this study. The results for a selection of typical tests showing different behaviors are illustrated in Fig. 2. Figures 2(a) and 2(b) present the deviator stress variations and excess pore pressures with axial strain for tailing samples considering the effect of confining pressures. From the results, strain hardening behavior has been observed on the tailing samples. The undrained behavior of a dilatant sample that reaches a steady state at large strains is shown in Fig. 2(a). When the shearing pore pressure attains a negative value, the samples behave in a dilative state under undrained tests as shown in Fig. 2(b). The peak deviator stress of the tailing samples varies from 180 kPa to 270 kPa with axial strain of 28% to 30% and cell pressures from 70 kpa to 150 kPa, respectively. The percentage increase of peak deviator stress is about 17%–50% as compared with the cell pressure of 150 kPa.

Fig. 2 Effect of confining pressures on (a) deviator stress versus axial strain; (b) pore pressure versus axial strain; (c) effective stress path of tailing samples tested at 65% relative density

Figure 2(c) shows the effective stress path, from the figure, it can be seen that for higher levels of confining pressures, the mean effective stress is found to be increasing. The effective friction angle 36° and cohesion 13 kPa for tailing samples was found from undrained triaxial experiments. Tailing samples dilated under compressive shearing and progressed towards stable at critical states, as indicated in Fig. 2(c). Thep–qdiagrams plotted are based on Massachusetts Institute of Technology stress field.

4.2 Cyclic triaxial test results

The cyclic behavior and liquefaction potential of tailing samples are investigated by taking the effect of varying axial-strain amplitudes (0.5%, 1.0% and 1.5%) and cell pressures (70 kPa, 100 kPa and 150 kPa)at a loading frequency of 0.3 Hz into account. The loading frequency was selected depending on the type of the material to allow sufficient time for pore-water pressure equalization. The number of cycles to produce initiation of initial liquefaction increases as the loading frequencies decrease. Harmonic loading is applied using a sine wave with frequencies of 0.3 Hz. Typical plots for the variations of excess pore pressure ratio, mean effective stress and deviator stress with number of cycles are shown in Figs. 3 and 4.

4.2.1 Effect of axial strain (%) on liquefaction potential of tailing material

From Figs. 3(a) and 3(b) and Figs. 4(a) and 4(b), it is seen that as the number of cycles goes up, the porepressure ratio increases, leading to a decrease in mean effective stress and deviator stress of the tailing samples.From Figs. 3(c) and 3(d) and Figs. 4(c) and 4(d), it is seen that due to the increase in the excess pore water pressure, the hysteresis loop between axial strain and deviator stress became flatter. In addition, the effective stress path moves towards the left until a final steady state is attained. The liquefaction potential of the tailing sample has been studied under the application of a loading frequency is 0.3 Hz and relative density of 0.65%, respectively, subjected to axial strains of 0.5%–1.5%. The results are tabulated in Table 2 and Fig. 3.It is observed from the results that an increase in axial strain amplitude leads to the decrease in liquefaction resistance. It decreases about 15%–42%, 12%–20% and 9.3%–27.5% when subjected to confining pressures of 70 kPa, 100 kPa and 150 kPa, respectively.

4.2.2 Effect of confining pressure on the liquefaction potential of tailing material

Fig. 3 Variation of (a) deviator stress with time (b) pore pressure ratio with time (c) stress-strain loop and (d) stress path on tailing samples considering the effect of axial strain amplitudes (tested at maximum dry density 18 kN/m3)

The effect of cell pressure on the potential of liquefaction of tailing samples have been studied under a relative density of 65% and loading frequency of 0.3 Hz.It is observed from the results (Table 2 and Fig. 4) that the number of cycles for initiation of liquefaction varies from 19 to 11, 25 to 20 and 40 to 29 for tailing samples tested at confining pressures of 70 kPa, 100 kPa and 150 kPa, respectively, with axial strains of 0.5%–1.5%.The liquefaction resistance was higher at high confining pressure and low cyclic axial strain when compared to low confining pressure and higher cyclic axial strain.It increased about 31%–52%, 37.5%–100% and 81%–260% for axial strains of 0.5%–1.5% (Table 2).

5 Dynamic properties of compacted tailing material

The modulus reduction and damping ratio with number of cycles have been estimated considering the variation of axial-strain amplitudes (0.5%, 1.0% and 1.5%), cell pressures (70 kPa, 100 kPa and 150 kPa)and loading frequencies. In general,Gsecincreases andξdecreases with a rise in cell pressures and with reduction in axial strains. From Table 2, it is observed that the rise in confining pressure from 70 kPa to 150 kPa results in the rise of the modulus reduction ratio. It increases about 10%–18% for 0.5%, 2.4%–10.3% for 1.0% and 2.9%–14% for 1.5%. Similarly,ξdecreases as the confining pressure increases. It decreases about 11%–20% for 0.5%, 2.8%–11.5% for 1.0% and 3%–15% for 1.5%.

6 Dynamic FE formulation

The saturated soil system is modelled as a two-phase material based on the Biot (1941) theory for porous media. Theu–pformulation uses pore fluid (p) and solid phase displacement (u) as primary variables assuming that the acceleration of the pore fluid relative to the soil acceleration is negligible. The u-p formulation can be written as follows (Biot, 1941; Zienkiewiczet al., 1999):

Caughey and Penzien (1960) showed that the damping matrix could be expressed as a linear combination of mass and stiffness matrices. The Rayleigh damping matrix is given by

whereα,βare called mass and stiffness proportional damping parameters. The relation between the damping ratio and Rayleigh damping parameters is given by

ω1is the fundamental frequency of the embankment(0.95 Hz),ωnis the excitation frequency andξis the damping ratio. Here, the Newmark algorithm was used to generate the numerical solution in the time domain.

Fig. 4 Variation of (a) deviator stress (b) pore pressure ratio with time (c) stress-strain loop and (d) stress path of tailing samples considering the effect of confining pressures (tested at maximum dry density 18 kN/m3)

7 Finite-element model geometry and boundary conditions

The majority of the (HZL) RA mine tailing dam has a downstream method of rising. There is only a small section where land is not available to adopt the upstream method of raising. A nonlinear FEM of the tailing dam is formulated in the FE framework of OpenSees as shown in Fig. 8. For meshing of the embankment, four node quadu-pelements are used. Each node of this element has two translational and one pore pressure degrees of freedom. The approximate size of the element for the dam constructed with upstream and downstream rising and foundation soil was taken as 2.2 m. A total of 5,349 nodes, 5,155 elements were created for downstream rising and 3,240 nodes, 3,104 elements were created for upstream rising for the finite element modelling.Weathered rock about 10 m below the foundation soil have been considered for the analyses. The source to site distance from seismic zone II to V was 750 km to 900 km.When calculating the liquefaction potentialRand magnitude, most will be important parameters. But the earthquake motion has been scaled to seismic zone II from V and it is applied at the bedrock with a distance of 75 m from below ground level. Before describing the modelling in more detail, a brief discussion about the two types of boundary conditions is as follows.

7.1 Fixed boundary condition

This is the simplest and most computationally efficient boundary condition. For a plane strain problem at the boundary, displacement in both theXandYdirections is zero. However, when the domain is fixed at the end, energy waves are reflected back at the boundary and the solution in the computational domain becomes polluted.

7.2 Viscous dash pot boundary condition

The most common way to restrict the infinite computational domain to a finite domain is with artificial boundaries. Lysmer and Kuhlemeyer (1969) introduced a viscous boundary. This boundary is a system of two series of viscous dampers oriented in the normal and tangential direction with a dashpot coefficient ofAρVP,AρVS, respectively,

According to the unbounded domain boundary conditions, the finite element method could be expressed as a combination of normal stressσ=ρVPw′ and shear stressτ=ρVSu′ .

8 Constitutive model

In the undrained condition, the behavior of saturated tailing material is simulated using fluid solid coupled Matrix. Solid phase of the tailing material, mine muck and foundation soil were modelled using a multi yield surface plasticity model (pressure dependent) and for impervious soil using the pressure independent multi yield material model. Here the nonlinearity of the tailing material has been modelled considering the yield surface as yield criteria (Prevost, 1978). Constitutive equations are presented in Appendix A. The Drucker-Prager type of yield surfaces are used. Input parameters in the simulation are as shown in Table 3. From the top of the embankment extending up to a depth of 18 (in RL 451.5–433.5) m,Material 1 is employed; Material 2 is employed for a depth from 18 m to 40 m (in RL 433.5–413.5), whereas for the lower half of the embankment (in RL 413.5–394.5), the property of Material 3 was applied as shown in Table 3. The average effective confining pressure was calculated for Material 1 (70 kPa), Material 2 (100 kPa)and Material 3 (150 kPa). This arrangement was made for the tailing material in both the upstream and downstream rising conditions. The experimental backbone curves of the tailing material were implemented in the 2D finite element program OpenSees. The form of the nonlinear hyperbolic equation used in OpenSees is shown below(Eqs. (6), (7) and (8)).

The four essential model input parameters are: the reference shear modulus (G); the maximum shear strain(γ); pressure depend coefficient (d), which is the fitted hyperbolic parameter, taken as 0.5); and the stress ratio.They are assigned based on calibration of the cyclic triaxial experimental data.

whereτis the shear stress,Gis the maximum shear modulus,ξis the damping ratio,γis the shear strain,γris the reference shear strain,cs1,cs2are the critical state parameters calculated from isotropic compression,pais the atmospheric pressure,Pris the reference pressure, andc,dandd1are the hyperbolic curve fitting parameters.

8.1 Tailing material – geomembrane connectivity

The geomembrane layer has been modelled as linear elastic truss elements. Uniaxial elastic material was used to assign the truss elements. To provide connectivity between the tailing element nodes and geomembrane nodes, zero-length contact elements were used. The frictional angle (δ) of 37 has been recorded between two materials. The tensile strength of the geomembrane was found to be 484 MPa.

9 Model validation

In this section, a comparison has been made between the numerical model and experimental results of the tailing sample. To develop the numerical model, a sensitivity analysis was carried out to accurately obtain the size of the element. Ten elements have been chosen for the study of the elemental model. The four node quadu-pelement was employed to investigate the seismic response of the saturated tailing sample. Figure 5 shows a comparison of the results including (a) deviator stress versus time (b) mean effective stress versus time, and (c)pore pressure ratio versus time at a cell pressure of 70 kPa, loading frequency of 0.3 Hz and at a relative densityof 65%. From Fig. 5, it is understood that the numerical model shows good agreement with the experimental results. Maximum deviation between experimental and numerical results is with an error of 15%. Similarly,Fig. 6 shows a comparison between the experimentally observed backbone curves with the numerical model.The observed backbone curves were recorded at different depths that approximately corresponded to the overburden pressure close to the experimental confining pressures. From Fig. 6(a) it is noted that the trend of the stress-strain behavior of the tailing material from the experiment match reasonably well with the numerical model. Maximum deviation between the experimental and numerical (stiffness and damping ratio) values were found to be 8%–20%, 25%–27%, for 70 kPa, 2%–15%,1.5%–15% for 100 kPa and 14%–18%, 9.6%–16% for 150 kPa.

Table 3 Input parameters used in pressure-dependent and independent multi-yield surface model for different materials

Fig. 5 Comparison of experiment with element model results of undrained cyclic triaxial test including (a) deviator stress versus time (b) mean effective stress versus time,(c) pore pressure ratio versus time tested at confining pressure of 70 kPa

Fig. 6 Comparison of experiment with element model results of undrained cyclic triaxial test (a) backbone curves(b) Gsec vs. shear strain (c) damping ratio vs. shear strain

A synopsis of the results of the fully coupled effective stress is presented in the following sections for both upstream and downstream cases using varying boundary conditions in terms of displacement flow vectors, accelerations, pore pressures and hysteresis loops recorded along the tailing column (shown in Figs. 8(a), 8(b)).

10 Results of numerical analyses and discussion

The duration of strong motion can have an influence on earthquake damage. Many physical processes such as stiffness, strength and the build-up of pore pressure in loose and saturated sands, ashes and tailing material samples are sensitive to the number of stress reversals during an earthquake. Chung and Krinitzsky (1977)estimated the earthquake duration for soil (8–25 s) and rock sites (5–22 s) with varyingMwof 4.5–7.0. A set of five ground motions (Chamba, Chamoli, IndiaBurma and NortheastIndia) selected from the Himalayan frontal thrust with varying earthquake magnitudes of 4.5–6.5M and time duration varying from 8–25 s were selected for seismic analysis. Seismic motion data were taken from the COSMOS database (2010a, 2010b). As per seismic design code IS (1893) Part 1 (2002), the RA- mine tailing dam was located in Seismicity Region II, which has an expected seismic intensity of 0.10 g.To represent the regional seismicity, each of the seismic motion data points of amplitude are multiplied with a scaling factor of seismic intensity of zone II to obtain a seismicity intensity of 0.10 g. Scaled ground motions are summarized in Table 4. The acceleration response spectra of the input motions are shown in Fig. 7.The fundamental period of the tailing embankment is also overlain in the response spectra plot to show that, although each scaled motion has the same PGA,different motions exhibit different spectral acceleration demand for the embankment. To find the displacement and acceleration along the depth of the tailing material,the tailing column has been considered (shown in Fig. 8).

10.1 Acceleration profiles

The horizontal and vertical acceleration profiles along the entire depth of the tailing column obtained from FEM analysis for both upstream and downstream rising for fixed boundary conditions are presented in Fig. 9. The amplification in the foundation soil is not notable; however, the response is significantly amplified in the tailing portion for most of the motions for both the cases (dam constructed with downstream and upstream method of rising). From Fig. 9, it is also shown that the most critical acceleration is generated by the IndiaBurma input motion, whereas minute deformations are generated on the application of the Northeast earthquake for both the cases. From Fig. 9,the maximum peak ground horizontal acceleration at the base and surface for 0.10 g input bedrock motion is 0.14 g and 0.2 g for upstream rising and 0.2 g and 0.22 g for downstream rising conditions, indicating an amplification factor of 1.4, 2 for upstream rising and 2.0, 2.2 for downstream rising. The vertical acceleration shows even higher amplification, 0.2 g and 0.45 g for upstream rising and 0.11 g and 0.13 g for downstream rising conditions, indicating an amplification factor of 2, 4.5 for upstream rising and 1.1, 1.3 for downstream rising. It may be concluded that the chances of failure are greater for dams constructed by downstream rising conditions.

Figure 10 shows a comparison of peak horizontal and vertical acceleration profiles along the entire depth of the tailing column constructed by downstream rising for different boundary conditions. It is observed that the viscous boundary condition gives higher acceleration when compared with the fixed boundary condition. This may be because FBC boundary conditions in both theXandYdirections are fixed. Whereas for viscous boundary condition, some displacement was specified at the boundary in terms of tangential fluid velocity as well as normal fluid velocity. Due to the specified displacements,accelerations developed at the boundary. From the results, it can be concluded that viscous boundary conditions replicate the actual field situation gives accurate results, whereas the fixed boundary condition provides inaccurate results and is not acceptable. The increase of horizontal and vertical accelerations is about 2.3%–9.5%, 7.6%–15%, for Chamba 9.5%–13%,4%–30%, for Northeast and 33%–45%,15%–20%, for IndiaBurma earthquake motions when compared to the fixed boundary condition.

Fig. 7 Acceleration response spectra of the input motions

Table 4 Details of input ground motions

In addition, Figs. 11 (a) and 11(b) shows acceleration at the surface for different boundary conditions of the tailing column in the time domain for the Chamba and Northeast earthquake motions. Figures 11(c) and 11(d)shows the spectral accelerations at the crest and bedrock subjected to the Northeast and Chamba motions. For these two motions, the acceleration response spectra at the base are approximately the same (～ 0.45 g). Whereas at the crest, the spectral acceleration varies from 0.48 g and 0.5 g for fixed boundary condition (FBC) and 0.75 g(1.75 times higher) and 0.7 g for viscous boundary condition(BC) for the Chamba and Northeast, motions respectively.From the above explanation, it is clear that the viscous boundary condition replicates the actual field conditions.

Fig. 8 (a,b) FE mesh model of tailing dams constructed by downstream and upstream rising with applying fixed boundary conditions;(c,d) FE mesh model of tailing dams constructed by downstream and upstream rising conditions with applying viscous boundary conditions

10.2 Displacement flow vectors

Typical results showing the variations of horizontal and vertical displacement flow vectors for a tailing embankment rising by downstream conditions, under the excitation of the Chamba motion with varying boundary conditions, are shown in Fig. 12. From the results, it can be seen that from the foundation soil, more rapid horizontal deformations occur as the tailing material height is further incremented. It is also found that the viscous boundary condition gives higher displacement than the fixed boundary condition, which is the consequence of the phenomenon explained earlier. It may be concluded that viscous dampers with very large computational domain are expected to yield accurate results, which could be used as a benchmark to verify the results of other boundary conditions.

From the fixed boundary conditions, 107 mm of the maximum vertical displacement and 10 mm of the maximum horizontal displacement in the downstream condition have been obtained in the tailing material by applying the earthquake motion of the Chamba record.Similarly, for the same simulations of the viscous boundary condition, 126 mm of maximum vertical displacement and 58.5 mm of maximum horizontal displacement was observed. The increase of settlement of the dam is approximately 17% for vertical as compared with the fixed boundary condition. Smaller permanent vertical and horizontal settlements were observed during the Northeast earthquake motion in the analyses.

Fig. 9 Horizontal and vertical acceleration profiles (a, b) upstream rising, (c, d) downstream rising along the depth of tailing column for fixed boundary condition (tailing column shown in Fig. 8(a, b))

Fig. 10 Comparison of horizontal and vertical acceleration profiles along the depth of tailing column of dam constructed by downstream rising with varying boundary conditions

Fig. 11 (a, b) Acceleration time histories and (c, d) response spectra for Chamba and Northeast motions with varying boundary conditions

Fig. 12 (a, b) Dam constructed by downstream rising flow (display) vectors with varying boundary conditions (a) X- flow vectors(FBC), (b) X-flow vectors (Viscous BC); (c, d) downstream rising flow (display) vectors with varying boundary conditions(c) Y-disp (FBC), (d) Y-disp (viscous BC)

This is consistent with the results of the FE method,which has shown the strong effect of the input motion parameters, varying boundaries and magnitude of earthquakes on the predicted slope displacements. From this study, the seismic performance of the tailing dam constructed by downstream and upstream rising can be considered satisfactory from a life safety perspective during the most frequent occurring earthquake motions,because the maximum displacements are considerably lower.

10.3 Excess pore pressure profiles

The pore pressure profiles are recorded along the depth of the tailing column for downstream rising with varying boundary conditions for the Chamba,Northeast and India-Bangladesh earthquake motions as seen in Fig. 13. From the figure, it is observed that pore water pressure increases along the depth of the material column for both the cases (FB, viscous FBC).It is also seen that the viscous boundary condition gives a higher pore water pressure when compared with the fixed boundary condition. This may be a consequence of the phenomenon explained earlier. For fixed boundary conditions, the free surface of the material column shows excess pore pressure, which equals 10 kPa for both the Chamba and Northeast motions and at the base,it shows 30 kPa and 40 kPa. As for the India-Bangladesh motion, it shows 35 kPa for free surface and 80 kPa for the base. This is related to the various initial stress states reached at various locations in the dam before the shaking. Similarly, for viscous BC, the maximum value at the free surface is equal to 10 kPa and at the base, it is 125 kPa and 120 kPa for the Chamba and Northeast motions. For the India-Bangladesh motion, it is 35 kPa at the free surface and 100 kPa at the base. The increase is about 76%, 66% and 20% for the Chamba, Northeast and Bangladesh earthquake motions, respectively, when compared with the fixed boundary condition.

Fig. 13 Excess pore pressure profile along the depth of tailing column for downstream rising with varying boundary conditions

10.4 Stress-strain behavior

Figure 14 shows the hysteresis loop of two different Gauss elements at different depths of the tailing column with varying boundary conditions for downstream rising subjected to different earthquake motions. The Gauss elements are located at the surface and base inside the material column. From the figure, it can be seen that the stress-strain curve is representing a nonlinear material behavior for different motions and the shear strains become very large and plastic behavior is clearly visible in both the boundary condition approaches. It is also observed that as the overburden pressure increases, the shear stress also increases. Not surprisingly, the initial stiffness of the tailing material is found to escalate with the increase in the overburden pressure. Also, it is observed that the shear stress and shear strain values of the viscous boundary condition gives higher than the fixed boundary condition, which may be a consequence of the phenomenon explained earlier. The increase is about 18% and 28% for the Chamba and Northeast motions when compared to the fixed boundary conditions.

Fig. 14 Stress-strain hysteresis loop at various depths (a, b) surface, (c, d) at base of tailing column constructed by downstream rising with varying boundary conditions

Fig. 15 CSR and FSL profiles of tailing column (column shown in Fig. 8 (a, b)) constructed by (a, b) downstream and (c, d) upstream rising conditions

10.5 CSR and FSL profiles

Figure 15 shows the profiles of the cyclic stress ratio (CSR) and factor of safety against liquefaction(FSL) along the depth of the tailing column with varying boundary conditions constructed by the upstream and downstream method of rising. TheCSRis obtained by using the peak acceleration at each depth via a conventional method (Seed and Idriss, 1982). The water table is located at 441.5 m (in RL) shown in theFSLprofile. From the results, it can be observed that the tailing material constructed by upstream and downstream rising is not significantly affected by liquefaction throughout its depth for all the earthquake motions with magnitudes between 4.5–6.4 with a mean FSLvalue range of 2–2.5 for FBC upstream, 1.5–2.5 for viscous BC upstream,2–2.6 for FBC downstream, and 1.6–2.6 for viscous BC downstream.

11 Conclusions

The effective angle of friction and cohesion of the tailing samples are 36° and 13 kPa, respectively,which indicates good friction angle and low cohesion.Liquefaction resistance decreases as the cyclic axial strain amplitudes increase, and decreases about 15%–42%, 12%–20% and 9.3%–27.5% for 70 kPa, 100 kPa and 150 kPa, respectively. The tailing material shows higher liquefaction resistance and stiffness modulus.The modulus reduction ratio increases as the confining pressure increases, and increases about 10%–18%for 0.5%, 2.4%–10.3% for 1.0% and 2.9%–14% for 1.5%, respectively. Similarly, the increase in confining pressure leads to a decrease inξ, about 11%–20% for 0.5%, 2.8%–11.5% for 1.0% and 3%–15% for 1.5%,respectively, for the tailing samples. The numerical model shows good agreement with the experimental results on the tailing sample with a percentage deviation of 15%. The maximum deviation between experimental and analyticalGsecandξvalues are 8%–20%, 25%–27%,for 70 kPa, 2%–15%, 1.5%–15% for 100 kPa and 14%–18%, 9.6%–16% for 150 kPa.

A notable amount of amplification was observed in the acceleration response of most of the tailing embankment. Also, from the displacement flow vectors,it is found that the viscous boundary condition provides accurate results whereas for fixed boundary conditions,the results are not as accurate since the domain is suddenly fixed at the end and energy waves are reflected back at the boundary. It may be concluded that the viscous boundary condition replicates the actual field conditions.It was observed that that the stress-strain curve represents a nonlinear material behavior for different motions and the shear strains becomes very large and plastic behavior is clearly visible for both the boundary conditions.From the above results, it can be concluded that the seismic performance of the tailing dam constructed by downstream and upstream methods of rising can be considered satisfactory from a life safety perspective during the most demanding earthquake motions,because the maximum displacements are considerably lower. Through numerical simulation, it is found that the tailing dam is not susceptible to liquefaction from earthquake motions even for a magnitude > 5.5. Finally,based on the comparison of the tailing dams constructed by downstream and upstream methods of rising, the dam rising by the upstream building method is more stable.