Dynamic response of Zelazny Most tailings dam to mining induced extreme seismic event occurred in 2016

The paper deals with the stability analysis of tailings dam subjected to dynamic loading induced by mining shocks which occurred in neighbouring copper mine. The main goal of the paper was to model the dynamic response of the dam during two extreme paraseismic events which occurred in 2016 based on accelerograms recorded at the dam toe. Dynamic response of the tailings dam was calculated using finite element method and the implicit time-integration method implemented in commercial codes. The boundary condition corresponding to dynamic loading was determined by deconvolution procedure. The error analysis showed that most precise signal reproduction is achieved while using target signal with peak value reduced by 40% as a test signal. Both acceleration and displacement time-series were successfully reproduced. Moreover, the stability analysis was conducted for five independent signals with design peak horizontal acceleration and showed that no permanent displacements should occur . The temporary horizontal displacement of the dam crest should not exceed 13 mm, assuming equivalent linear material model.


Introduction
The paper deals with the stability analysis of tailings dam subjected to dynamic loading induced by mining shocks which occurred in neighbouring copper mine. The dam is a part of huge Tailings Storage Facility (TSF) Zelazny Most being the only place for deposition of post-flotation tailings for whole the Polish copper production, Fig. 1. A comprehensive description of geotechnical aspects of the development of one of the world's largest copper tailings depositories is given in [1].
Two relatively strong mining induced seismic events which occurred in March and October 2016 were the motive of presented analyses. The seismic energies of that events were equal to 1.4•10 7 J and 1.0•10 8 J, respectively and the hypocentres were localized at similar depth equal to 915 m and 963 m. In the paper the dynamic response of the tailings dam is assessed using finite element method and the implicit time-integration method. The analyses were carried out for the VIIIW cross-section of the Western Dam ( Fig. 1) where on March 8 th 2016 peak horizontal ground acceleration (PHA) equal to 683 mm/s 2 occurred which was the highest ever recorded acceleration at the TSF Zelazny Most. The epicentral distance to the seismic station located at the dam toe for analysed seismic events was equal to 1008 m and 1235 m, respectively.

Geotechnical characterisation
The tailings dam analyzed is being constructed by upstream method. The typical cross-section of the dam body is presented in Fig. 2. It consists of starter dam (2) and embankments/shells (1) which are formed as construction fill from coarsest tailings deposited near the dam crest. In order to improve the overall stability, in some parts of the facility the loading berm (5) is constructed. The tailings-water mixture is discharged from the dam crest creating a beach. In the process of the gravitational flow the segregation of the tailings occurs during which coarser tailings (3) settle close to the embankments whereas the finer ones flow towards the pond (4). In order to reflect different geotechnical properties the fine-grained tailings were divided into three parts, from which two are 80 m wide strips.
Due to the process of water infliltration large mass of tailings is fully saturated. In the analysis the position of phreatic line was assumed based on piezometric measurements, the elevation of circumferential and starter dam drains as well as the location of water line in the decant pond 200 m distant from the dam crest. At the moment of the analyzed extreme seismic shock the Western Dam was around 50 m high, the crest elevation was equal to 180 m asl and the tailings elevation was 2.5 m lower.

Fig. 2. VIIIW cross-section analysed.
The Zelazny Most facility is founded on a very complex geological structure consisting of Pliocene (Pl) and Quaternary (Q) clays and sands affected in the past by the three glaciations. In the analysed cross-section there exists 75 m deep layer of quaternary clays which includes relatively small lenses of quaternary fluvioglacial sands and clays. Pliocene clays are expected to be beneath this formation with around 17 m deep layer of tertiary sands. Moreover, near the surface there is 2 m deep layer of quaternary clays which had been weakened due to solifluction and cryoturbation processes occurring in interglacial periods.
The stiffness of these soils, even at great depths, is far lower than the stiffness of rock. The above has been included in the discrete model assumed for the stability analysis. In order to construct stable and efficient discrete model the geometry of the dam and the geological structure of the foundation were simplified, Fig. 2.

Seismic loading
The region of the Zelazny Most facility is characterized by very low natural seismicity [2,3], however due to mining activity it is exposed to mining-induced tremors with close epicentral distances. Seismic activity is monitored by local seismic network SEJS-NET [4]. Seismic stations, that record three components of earthquake vibrations, are located at the toe of the dam (P) and on its slope (K) installed in six cross-sections, Fig. 2. The X axis is directed downstream the dam.
In order to reduce low frequency noise and to achieve physically permissible velocity and displacement time series [5,6,7] the signals recorded by seismic station located in VIIIW cross-section on 8 th of March and 29 th November 2016 were modified using polynomial trend elimination technique. The frequencies above 10 Hz were cut as much higher than the fundamental frequency of the dam. The sign of the X signal was changed to be compatible with the axis direction of the numerical model. Moreover, in order to reduce the computational effort the signals were cut covering 0.05% to 99% Arias intensity [8].
Modified horizontal acceleration time-histories and its frequency spectrums are shown in Fig. 3 and Fig. 4. Due to the high amplitudes of waves with frequencies close to 10 Hz the parameters of the signals were noticeably modified. For example the PGAx of the signal measured at the toe of the dam on 8 th March has been decreased by 10.6% and the vertical one by 44.0% to be equal 610 mm/s2 and 225 mm/s2, respectively. The dominant frequencies of modified horizontal and vertical accelerograms are equal to 2.59 Hz (3.98 Hz) and 9.02 Hz, respectively. In both analysed seismic events the vertical peak ground acceleration was over 2 times lower, and the dominant frequency was three times higher than respective horizontal one.
Some parameters that characterise analysed horizontal acceleration time-histories are shown in Table  1. The highest value of PGAx was recorded on March 8th, however the highest value of Arias intensity IA equal to 501 mm/s was registered on November 29th. Strong motion duration tA did not exceed 5 s. The dominant frequency of signals recorded by both seismic stations were equal. Based on signals recorded during past shocks with energy higher than 110 7 J seismic hazard analysis was conducted for the mining activity plan in the period 2012 -2042. The forecast of peak horizontal ground acceleration value (PHA) for the analysed cross-section was estimated to be 0.13g [9]. This value with 10% of exceedance over a period of 30 years has been assumed as design acceleration. The analyses were performed for five independent accelerograms recorded in close distance to the cross-section analysed.
In order to determine the signal applied at the bottom boundary of the discrete model, and which corresponds to the measured or designed PHA related to free surface motions at the toe of the dam the deconvolution procedure has been implemented.

Discrete model
The dynamic analysis was performed using two dimensional finite element commercial codes. Full dynamic approach was conducted by PLAXIS 2D (2011), whereas simplified Newmark approach was calculated in GEOSTUDIO (2007) [10,11,12,13,14,15,16]. The same geometry of the model was implemented into both codes, the size of which was 1575.5x200 m (7 widths x 4 heights of the dam slope). The plane strain condition was assumed. In PLAXIS 2D viscous elements on the side boundaries of the model were used, whereas in GEOSTUDIO model boundaries are rigid causing some wave reflections. Despite significant model dimensions, seismic wave propagation along model base and its potential transformation are not considered, and seismic loading is uniformly applied along whole bottom boundary of the model. In both cases prescribed horizontal displacement time-series is used as the boundary condition, which is automatically calculated from predefined accelerogram. The base of the model is rigid reflecting all the dynamic energy backwards to the model.
In both cases stiffness dependency on stress and strains was incorporated [17,18], however with some differences depending on the model applied. The equivalent linear model (EQ) was used in GEOSTUDIO and Hardening Soil with small strain stiffness (HSs) in PLAXIS.
Each model was discretized taking into account the variable stiffness of the soil and the 10 Hz frequency wave length to be described by at least 5 nodes. In case of GEOSTUDIO the 4-nodes quad elements were used, whereas in PLAXIS 15-nodes triangular elements. The distances between the nodes in the direction of wave propagation were ranged from 3 to 8 m.
According to the Nyquist's criterion, time step equal to 0.05 s is sufficient to reproduce waves with frequency lower than 10 Hz. However the stability requirement based on Courant-Friedrich-Levy condition indicated that the time step should not be larger than 0.02 s. The dynamic calculations in GEOSTUDIO were conducted with time step equal to 0.0052 s (original signal was recorded with time step equal to 0.0026 s). The PLAXIS code has its intrinsic criterion which forced the calculation time step to be 0.002 s.

Materials
To model normalized shear modulus degradation curves G/Gmax depending on confining pressure (p') as well as plasticity index (PI) the formula given by Ishibashi and Zhang [19] was adopted in Quake/W for both cohesive and cohesionless soils. G/Gmax curves were determined mostly based on plasticity index, however for Tertiary clays and tailings they were fitted to the experimental TX shearing test results. Curves applied for selected materials at selected depth, e.g. PlCl_H -pliocene clay below the H depth, are presented in Fig. 5. Based on determined G/Gmax curves the reference value of shear strains 0.7 characterising stiffness reduction curve in HSs model was calculated. The dependency of maximum shear modulus Gmax on stress state was included. In EQ model Gmax depends on vertical effective stress, however in HSs the Gmax depends on minor effective stress.
Respective curves for Tertiary clays and tailings were approximated based on numerous series of shear wave velocity measurements for undisturbed soils samples, made in TX with piezoelectric bender elements at different level of confining stress. For cohesionless soils the empirical formula given in [20,21] was used. In order to model internal soil damping the Rayleigh damping model was applied in both codes however the values of the α and β parameters were established in different way. Constant damping ratio of 3% has been assumed for foundation soils and 5% for dam body and tailings [22]. Some additional hysteretic damping is included in the PLAXIS discrete model by using the HSs material model. The Poisson ratio for dynamic calculations was defined as 0.495 for all soils under phreatic line which is the closest allowable approach to the constant volume condition in computational model.
Shear strength of the majority of soils in stability analysis has been modelled by Coulomb-Mohr yield condition, excluding fine tailings under phreatic line for which linear variation of strength with mean effective stress was assumed (SHANSEP model). Angle of internal friction of the dam's materials is equal to 34 and of foundation soils is varying from 17 to 34.

Dynamic boundary condition
The kinematic boundary condition was defined using the deconvolution procedure described by Kramer in 1996 [23]. First, the transmittance functions were calculated based on the response of the model to assumed horizontal accelerogram. Then the boundary condition was calculated from recorded acceleration time histories, both vertical and horizontal ones. The main assumptions of the method described above are: linearity and constant parameters of the dynamic system which is not hold in the case of HSs model. In case of equivalent linear model there is also problem with the calculation of the equivalent shear stiffness. Thus during the study different types of test accelerograms were assumed, i.a. chirp signal which introduces another frequency in every subsequent time step. Application of such a signal has allowed to determine one transmittance function for multiple signals and thus has reduced the computational effort.

Calculation steps in simplified approach
Simplified dynamic displacement analysis of the dam is performed based on Newmark's approach using GEOSTUDIO software. In first step the initial static equilibrium stress state within the dam mass before earthquake is established. For the sake of simplicity the stage construction and overconsolidation was neglected. Subsequently, the simplified one-phase dynamic calculations are done using equivalent linear procedure. Due to these assumptions, only elastic dynamic strains are expected and the generation of pore water pressure is not analysed. After that, Newmark safety analyses using Slope/W were conducted. Stability calculations were performed for Coulomb-Mohr model and Stress History And Normalized Soil Engineering Properties (SHANSEP). Over one hundred various slip surfaces were examined covering all admissible failure mechanisms: from shallow to deep surfaces passing through the soft clays and boundary between saturated and unsaturated tailings. Finally, calculated permanent displacements are assessed in reference to the assumed critical value.

Calculation steps in full dynamic approach
In case of full dynamic approach overconsolidation of subsoils and stage construction processes of the dam were included. Simplified construction process was implemented which consists of continuous construction of 10 m high embankments and deposition of similar thickness of tailings during totally 730 days (2 years) as well as assuming the equivalent break for the period 2190 days (6 years) corresponding to the periods where the dams were not being raised. Both stages are being modeled in the Consolidation in a classical mode module. After that dynamic loading stage is started. Assessment of dynamic stability of the dam is performed in selected points (Fig. 7) based on the analysis of the results of following magnitudes: value of peak horizontal displacements (R1-R10), value of peak shear stress and deviatoric strains as well as the ratio of dynamic to static maximum shear stresses (P1-P10).

Results
The dynamic response of the tailings dam was successfully reproduced for both material models. Fit has been confirmed by comparison of calculated and measured spectrum of acceleration amplitudes at points corresponding to the location of both seismic stations, Fig. 8. Good agreement was achieved for equivalent linear model while using chirp signal as test signal. However, mean value of relative error of amplitude spectrum is 50% higher than for target signal. In the case of HSs material model quite good fit was achieved by incorporating recorded signal. In both cases appropriate selection of PHA turned out to be very important. The PHA of test signals was reduced by 40% with respect to the measured one.
Errors of the reproduced acceleration and displacement time-series at the seismic station located at the dam toe for both seismic events calculated in GEOSTUDIO were shown in Table 2. Relative error of calculated and measured peak horizontal acceleration pax for March seismic event was equal to 1.2%. Mean value of relative errors |Fa| calculated for every frequency was equal to 0.31. L 2 -norm ||Fa|| of absolute errors vector was equal to 0.023. Relative error of horizontal displacement pdx was not larger than 2.2%.
As it was mentioned above, the signals in PLAXIS were reproduced with larger error, e.g. relative error of pax was equal to 4.46%, and mean value of relative errors |Fa| calculated for every frequencies was equal to 4,6.  The stability analysis carried out for five independent accelerograms scaled to design peak acceleration did not result in any permanent displacements D calculated by simplified approach, Table 3. Maximum value of peak horizontal displacement pux of the dam crest caused by analyzed dynamic loading was equal to 12.9 mm.
Stability analysis with full dynamic approach has been carried out based on the assessment of the resulting variables calculated at selected points. In the article selected results are presented only, Table 3. Maximum value of permanent horizontal displacement due to dynamic loading applied was equal 9.3 mm. Maximum value of peak shear strains p calculated at selected location P1-P10 were in range 2.1•10 -5 to 9.1•10 -5 . The highest additional shear stress with reference to static shear stress maxd/s was calculated at point P3 and it was equal to 10.1%.

Conclusions
The dynamic response of the tailings dam was successfully reproduced. Both simplified and full stability analyses of the dam were performed and showed that the stability of the dam was preserved. For predicted design peak ground acceleration derived for the mining plan period from 2012 to 2042 no permanent displacements should occur. For equivalent linear material model the temporary horizontal displacement of the dam crest should not exceed 13 mm, assuming.
Some contributions to the deconvolution procedure as to test function selection have been proposed. According to that it is suggested to use the chirp signal with the peak value equal to 50% of peak value of the reproduced signal. However to gain precise reproduced signal it is recommended to use target signal with peak value reduced by 40%.