Applied Science, Faculty of
Civil Engineering, Department of
DSpace
UBCV
Lino Ramírez, Erick Alexander
Master of Applied Science - MASc
University of British Columbia
The release of tailings and impounded water after dam breaches has had catastrophic consequences for the population and the environment. These events have gained the attention of the geotechnical community because they occur without warning, causing massive damage. Flow liquefaction of tailings materials has been the cause of instability of recent tailings dam failures. These disasters were initiated by the triggering event that led the system to experience flow liquefaction and large deformations. To predict the evolution of flow-type instabilities, designers usually carry out dam-break analyses that rely on fluid mechanics methods since conventional slope stability methods such as limit equilibrium and finite elements only provide information until the onset of failure. Such decoupling of the geomechanics-based analysis for stability and fluid mechanics-based analysis for dam-break is not desirable, given the recent advances in critical state soil mechanics, strain-softening material models, and large deformation analysis methods. Hence, there is a need to develop analyses that can simulate the failure initiation and the post-failure stage observed during tailings dam disasters. In the present study, the triggering and the subsequent motion have been analyzed in a unified framework for the case of a saturated upstream tailings dam. Two different triggers are considered: the increment of pore water pressure and the sudden loss of strength in the foundation. A critical state bounding surface plasticity model is used for the tailings. This model can reproduce the contractive brittle behavior of loose saturated tailings. A two-phase single-point formulation of the material point method (MPM) has been considered to simulate the whole deformation process. This method has recently increased its popularity as a tool to deal with large deformation problems in geotechnical engineering. The influence of the initial void ratio on system performance is evaluated. The current study shows the capabilities of MPM along with an advanced constitutive model in simulating the complex stress-deformation mechanism observed in tailings dam failures.
Modeling of flow liquefaction and large deformations in tailings damsusing material point methodbyErick Alexander Lino RamírezB.Sc. Universidad Nacional de Ingeniería, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Civil Engineering)The University of British Columbia(Vancouver)October 2021© Erick Alexander Lino Ramírez, 2021The following individuals certify that they have read, and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, the thesis entitled:Modeling of flow liquefaction and large deformations in tailings dams using materialpoint methodsubmitted by Erick Alexander Lino Ramírez in partial fulfillment of the requirements for thedegree of Master of Applied Science in Civil Engineering.Examining Committee:Dr. Mahdi Taiebat, Civil Engineering, UBCSupervisorDr. Arcesio Lizcano, SRK Consulting CanadaAdditional ExamineriiAbstractThe release of tailings and impounded water after dam breaches has had catastrophic consequencesfor the population and the environment. These events have gained the attention of the geotechni-cal community because they occur without warning, causing massive damage. Flow liquefactionof tailings materials has been the cause of instability of recent tailings dam failures. These disas-ters were initiated by the triggering event that led the system to experience flow liquefaction andlarge deformations. To predict the evolution of flow-type instabilities, designers usually carry outdam-break analyses that rely on fluid mechanics methods since conventional slope stability meth-ods such as limit equilibrium and finite elements only provide information until the onset of fail-ure. Such decoupling of the geomechanics-based analysis for stability and fluid mechanics-basedanalysis for dam-break is not desirable, given the recent advances in critical state soil mechanics,strain-softening material models, and large deformation analysis methods. Hence, there is a need todevelop analyses that can simulate the failure initiation and the post-failure stage observed duringtailings dam disasters.In the present study, the triggering and the subsequent motion have been analyzed in a unifiedframework for the case of a saturated upstream tailings dam. Two different triggers are considered:the increment of pore water pressure and the sudden loss of strength in the foundation. A criti-cal state bounding surface plasticity model is used for the tailings. This model can reproduce thecontractive brittle behavior of loose saturated tailings. A two-phase single-point formulation of thematerial point method (MPM) has been considered to simulate the whole deformation process. Thismethod has recently increased its popularity as a tool to deal with large deformation problems ingeotechnical engineering. The influence of the initial void ratio on system performance is evalu-ated. The current study shows the capabilities of MPM along with an advanced constitutive modelin simulating the complex stress-deformation mechanism observed in tailings dam failures.iiiLay SummaryRecent tailings dam failures have impacted significantly the environment, human life, and property.This has prompted the geotechnical community to improve this situation, thus providing numericaltools to simulate such problem that typically involve large deformations. In this context, the maingoal of this research is to simulate the failure of a tailings dam using an advanced model and anumerical platform that can capture large deformations. A hypothetical tailings dam was analyzed.The simulation results provided valuable information to understand the deformation process. Suchinformation may help practitioners to evaluate the consequences of tailings dam disasters.ivPrefaceThe material presented in this thesis is an outcome of a research project on modeling and simulationof the spontaneous failure of tailings dams in the Theoretical and Applied Geomechanics (TAG)research group at UBC, in collaboration with SRK Consulting Canada Inc. (SRK).The implementation of the constitutive model and the Anura3D software used in this studyare freely distributed and were not developed by the author. The integration between these twocomponents to perform numerical simulations in a unified framework was done by the author andrepresent the main contribution of this research. The identification and design of the numericalmodels, the execution of numerical simulations, and the analysis and interpretation of results arethe original, unpublished, and independent work of the author, Erick Alexander Lino Ramírez, donein consultation with his research supervisor at UBC, Professor Mahdi Taiebat, and the industryadvisor, Dr. Arcesio Lizcano.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Failure of tailings dams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.4 Thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Numerical Methods for Large Deformations . . . . . . . . . . . . . . . . . . . . . . 82.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Mesh-based methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2.1 Arbitrary Lagrangian-Eulerian method . . . . . . . . . . . . . . . . . . . . 102.2.2 Coupled Eulerian-Lagrangian method . . . . . . . . . . . . . . . . . . . . 102.3 Particle-based methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3.1 Smoothed Particle Hydrodynamics (SPH) . . . . . . . . . . . . . . . . . . 102.3.2 The Material Point Method (MPM) . . . . . . . . . . . . . . . . . . . . . 11vi2.4 Why MPM? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 Modeling Flow Liquefaction with MPM . . . . . . . . . . . . . . . . . . . . . . . . . 173.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.2 Flow liquefaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.2.1 Critical state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.2.2 The state parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.2.3 Undrained softening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.3 Constitutive models and MPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.4 SANISAND . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.4.1 General description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.4.2 Model overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.4.3 Verification in Anura3D . . . . . . . . . . . . . . . . . . . . . . . . . . . 274 Description of the numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.1 Model setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.1.1 Geometry description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.1.2 Zoning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.1.3 Background mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.1.4 Sequence of modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.1.5 Definition of triggers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.3 Method of solving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.3.1 Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.3.2 Time increment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344.4 Material parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344.4.1 Upper and lower tailings . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.4.2 Raising lifts, foundation and buttress . . . . . . . . . . . . . . . . . . . . . 394.5 Initial and boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.5.1 Stage 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.5.2 Stage 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.5.3 Stage 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 445.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 445.2 Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45vii5.2.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 455.2.2 Deformation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 455.2.3 Time evolution of material state . . . . . . . . . . . . . . . . . . . . . . . 465.3 Case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 535.3.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 535.3.2 Deformation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 535.3.3 Time evolution of material state . . . . . . . . . . . . . . . . . . . . . . . 545.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 546 Assessment of flow liquefaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 596.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 596.2 Effect of initial void ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 606.2.1 Analyses with trigger 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 606.2.2 Analyses with trigger 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.3 Effect of trigger intensity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.3.1 Analyses with trigger 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.3.2 Analyses with trigger 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 626.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 627 Summary, Conclusions and Future Research . . . . . . . . . . . . . . . . . . . . . . 767.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 767.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 777.3 Recommendations for future research . . . . . . . . . . . . . . . . . . . . . . . . 77Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79A Use of SANISAND for different values of state parameter . . . . . . . . . . . . . . . 86A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86A.2 Selection of initial state parameter . . . . . . . . . . . . . . . . . . . . . . . . . . 86A.3 Element tests simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86A.4 Summary of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90viiiList of TablesTable 4.1 SANISAND model parameters for tailings . . . . . . . . . . . . . . . . . . . . 36Table 4.2 Physical parameters for tailings . . . . . . . . . . . . . . . . . . . . . . . . . . 39Table 4.3 Mechanical parameters for Raising Lift, Foundation and Buttress units . . . . . 39Table 4.4 Basic parameters for Raising Lifts, Foundation and Buttress . . . . . . . . . . . 39ixList of FiguresFigure 1.1 Cumulative released volume and life loss due to tailings dam failures. FromSantamarina et al. (2019) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1Figure 1.2 Aerial photo of Stava before (left) and after (right) the tailings dam failure.From WISE Uranium Project . . . . . . . . . . . . . . . . . . . . . . . . . . . 2Figure 1.3 Sullivan failure. Ground view of the movement two days after the event. FromDavies et al. (1998) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3Figure 1.4 Aerial view of the Merriespruit dam failure. From Fourie et al. (2001) . . . . . 4Figure 1.5 Fundão tailings dam before (left) and after (right) the failure. From Morgensternet al. (2016) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5Figure 1.6 Feijão tailings dam 11s after the first observed deformation. From Robertsonet al. (2019) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5Figure 2.1 Spatial discretization. (a) Lagrangian (b) Eulerian. From Yerro (2015) . . . . . 9Figure 2.2 MPM procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12Figure 3.1 Critical void ratio from direct shear tests (Casagrande, 1976) . . . . . . . . . . 19Figure 3.2 Critical state (CS) for the Stava tailings. Modified from Carrera et al. (2011) . . 19Figure 3.3 Definition of state parameter. Modified from Powrie (2018) . . . . . . . . . . 20Figure 3.4 Undrained softening. From Gens Solé (2019) . . . . . . . . . . . . . . . . . . 20Figure 3.5 Characteristic behavior of saturated sand under undrained monotonic loading.From Vaid and Chern (1985) . . . . . . . . . . . . . . . . . . . . . . . . . . . 22Figure 3.6 Schematic diagram of the SANISAND model (Yang et al. 2019) . . . . . . . . 24Figure 3.7 Undrained triaxial simulations on Toyoura Sand. (a) e0 = 0.833 (b) e0 = 0.907 . 28Figure 4.1 Geometry of the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 4.2 Background mesh: (a) Triangular elements (b) Initial material points distribution 31Figure 4.3 Simulation stages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32xFigure 4.4 Trigger application - Description . . . . . . . . . . . . . . . . . . . . . . . . . 33Figure 4.5 Constitutive models in model domain . . . . . . . . . . . . . . . . . . . . . . 34Figure 4.6 Initial void ratios assumed in the model: (a) Case 1 (b) Case 2 . . . . . . . . . 35Figure 4.7 SANISAND calibration for polymetallic tailings . . . . . . . . . . . . . . . . 37Figure 4.8 Critical state line for various tailings (continuous line) and sands (dashed lines).Modified from Jefferies and Been (2016). . . . . . . . . . . . . . . . . . . . . 38Figure 4.9 Initial conditions - Stage 1: (a) mean effective stress (b) pore-water pressure . . 40Figure 4.10 Solid boundary conditions - Stage 1: (a) Essential (b) Natural . . . . . . . . . 41Figure 4.11 Fluid boundary conditions - Stage 1: (a) Essential (b) Natural . . . . . . . . . 41Figure 4.12 Initial conditions - Stage 2: (a) mean effective stress (b) pore-water pressure . . 42Figure 4.13 Initial conditions - Stage 3: (a) mean effective stress (b) pore-water pressure . . 43Figure 4.14 Fluid boundary conditions - Stage 3: (a) Essential (b) Natural . . . . . . . . . 43Figure 5.1 Control points used in analysis of results . . . . . . . . . . . . . . . . . . . . . 44Figure 5.2 Displacement and velocities of control points P1, P2, P3 and P4 for case 1 . . . 45Figure 5.3 Variation of tailings dam deformed shape for case 1 . . . . . . . . . . . . . . . 48Figure 5.4 Evolution of deviatoric strains for case 1 . . . . . . . . . . . . . . . . . . . . . 49Figure 5.5 Evolution of mean effective stress for case 1 . . . . . . . . . . . . . . . . . . . 50Figure 5.6 Evolution of soil parameters at control point P2 for case 1 (a) Variation of voidratio, p and displacements (b) Variation of η , and stress development in thep−q and e− p plot during the deformation process . . . . . . . . . . . . . . . 51Figure 5.7 Evolution of soil parameters at control point P4 for case 1 (a) Variation of voidratio, p and displacements (b) Variation of η , and stress development in thep−q and e− p plot during the deformation process . . . . . . . . . . . . . . . 52Figure 5.8 Displacement and velocities of control points P1, P2, P3 and P4 for case 2 . . . 53Figure 5.9 Variation of tailings dam deformed shape for case 2 . . . . . . . . . . . . . . . 55Figure 5.10 Evolution of deviatoric strains for case 2 . . . . . . . . . . . . . . . . . . . . . 56Figure 5.11 Evolution of mean effective stress for case 2 . . . . . . . . . . . . . . . . . . . 57Figure 5.12 Evolution of soil parameters at control point P2 for case 2 (a) Variation of voidratio, p and displacements (b) Variation of η , and stress development in thep−q and e− p plot during the deformation process . . . . . . . . . . . . . . . 58Figure 6.1 Definition of the runout distance and selected control (moving) point for vul-nerability assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Figure 6.2 Comparison of the evolution of deviatoric strains with different initial void ra-tios: (a) e0 = 0.82 (b) e0 = 0.86 when trigger 1 is considered . . . . . . . . . . 64xiFigure 6.3 Comparison of the evolution of the tailings dam deformed shape for differentinitial void ratios: (a) e0 = 0.82 (b) e0 = 0.86 when trigger 1 is considered . . . 65Figure 6.4 Comparisons of the evolution of displacement, state parameter and mean pres-sure for different initial void ratios at control point P when trigger 1 is considered 66Figure 6.5 Comparison of the evolution of deviatoric strains with different initial void ra-tios: (a) e0 = 0.86 (b) e0 = 0.87 when trigger 2 is considered . . . . . . . . . . 67Figure 6.6 Comparison of the evolution of the tailings dam deformed shape for differentinitial void ratios: (a) e0 = 0.86 (b) e0 = 0.87 when trigger 2 is considered . . . 68Figure 6.7 Comparison of the evolution of deviatoric strains with different trigger intensi-ties: (a) pw = 200 kPa (b) pw = 250 kPa when e0 = 0.80 is considered . . . . . 69Figure 6.8 Comparison of the evolution of the tailings dam deformed shape with differenttrigger intensities: (a) pw = 200 kPa (b) pw = 250 kPa when e0 = 0.80 is considered 70Figure 6.9 Comparisons of the evolution of displacement, state parameter and mean pres-sure for different trigger intensities at control point P when e0 = 0.80 is consid-ered . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 6.10 Comparison of the evolution of deviatoric strains with different foundation thick-ness: (a) h = 1.0 m (b) h = 2.0 m when e0 = 0.86 is considered . . . . . . . . . 72Figure 6.11 Comparison of the evolution of the tailings dam deformed shape with differentfoundation thickness: (a) h = 1.0 m (b) h = 2.0 m when e0 = 0.86 is considered 73Figure 6.12 Relationship between runout and initial void ratio for different values of porewater pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74Figure 6.13 Relationship between runout and initial void ratio for different thicknesses inthe foundation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure A.1 Distribution of initial state parameter for different values of initial void ratio:(a) e0 = 0.78 (b) e0 = 0.80 (c) e0 = 0.82 (d) e0 = 0.84 (e) e0 = 0.87 . . . . . . . 87Figure A.2 Element test response for: (a) ψ = -0.05 (b) ψ = +0.05 . . . . . . . . . . . . . 88Figure A.3 Element test response for: (a) ψ = +0.10 (b) ψ = +0.15 . . . . . . . . . . . . . 89Figure A.4 Relationship between shear strength ratio and the initial state parameter of thelower tailings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90xiiAcknowledgmentsAbove all I would like to thank my wife Yuri Sulca for her deep love and her care for my devel-opment. Thank you for being so supportive and optimistic. But most of all, thank you for beingmy best friend during the last 13 years. Thank you for being part of this ongoing adventure, yoursupport and encouragement was in the end what made this thesis possible.I would like to express my sincere gratitude to my research supervisor Dr. Mahdi Taiebat forhis patience, guidance, and unfailing support throughout my study. I have benefited greatly not onlyfrom his vast knowledge in computational geomechanics but also from his kindness, considerationand sincere care in my struggling period. I am proud of, and grateful for, my time working with Dr.Taiebat.I also want to thank Dr. Arcesio Lizcano for our discussions and his encouragements with aperfect blend of insight and humor. His energetic, creative and open mind has benefited me notably.I also thank Luis Angel Aviles from Universitat Politècnica de Catalunya (UPC) for introducingme to the MPM world and for several hours of discussion in different chapters. He always had apositive attitude that inspired me.I am also grateful for the research partnership with SRK Consulting, and in particular to late Dr.Maritz Rykaart. He made all in his power to support financially this investigation. His memory andlegacy will be with me always. I also want to thank my SRK friends Osvaldo Ledesma and AndresBarrero for their comments and feedback on my research progress.I am indebted also to members of the Theoretical and Applied Geomechanics (TAG) researchgroup. I am greatly thankful to Dr. Ming Yang who helped me find my footing as I started thisprocess, Mr. Andrés Reyes for encouraging me in the hard times, Ms. Francesca Palmieri, Mr.Sounik Banerjee and Ms. Keith Pérez.My mother Lucy deserves my endless gratitude for here resilience, for teaching me the valuesof hard work, honesty and empathy. Thank you mom for working hard so that I could have theopportunity to follow my dreams. Thank you aunt Miriam for helping my mother to take care ofme. I also want to thank my sister Fiorella for her emotional support and her kindness.xiiiFinally, I would like to thank my late father Humberto for foster in me the passion towardengineering as a career.Financial support provided by Mitacs is gratefully acknowledged.xivChapter 1Introduction1.1 Failure of tailings damsOver the last 100 years the world has witnessed an increasing number of catastrophic failures ofmine tailings facilities. Such structures are built to ensure that the tailings are safely stored toprotect the environment. However, when tailings dams leak and spill tailings into the environment,the consequences are detrimental for the planet, profits, and people. Santamarina et al. (2019)showed that due to tailings dam failures the amount of released volume and the loss of human lifehave increased dramatically (see Figure 1.1). In this context, the failure of tailings dams represent aglobal issue that requires a proper understanding. A review of tailings dam failures is presented.SCIENCE sciencemag.orgGRAPHIC: C. BICKEL/SCIENCEgest to regulators and the public that the problem lies with the impounded materi-als themselves. However, in the absence of internal collapse or induced shear (for example, as a result of a seismic event), liquefaction and outflow of ponded ash and tailings occur after the dam has failed. Thus, liquefaction does not cause the fail-ure, but rather the disaster that follows.FAILURE MECHANISMSForensic analyses of ash pond and tailings dam failures are hindered by the massive destruction that washes away evidence from the failed zone. Nevertheless, forensic investigations of notable dam failures have identified several mechanisms, including overtopping due to water mismanagement (Merriespruit, South Africa, 1994) (6), shear failure of foundation soils (Mount Polley, Canada, 2014) (7), and shear of compressible low-permeability tailings placed near the perimeter of the impoundment (Samarco, Brazil, 2015) (8). These postfailure investi-gations have often found the convergence of more than one weakness in the design, construction, and/or operation of the dam.The Brumadinho dam was closely ob-served, and abundant laboratory tests, in situ tests, and monitoring data were avail-able. Furthermore, it failed 3 years after closure, contrary to the expectation that geotechnical structures become more stable over time. Indeed, an August 2018 audit concluded that the dam was stable (9). That an apparently closely monitored and inac-tive dam failed catastrophically less than half a year after being audited suggests that there may be gaps in the scientific under-standing of waste-storage facilities and of time-delayed triggers of failure. Filling these knowledge gaps will require changes in how coal ash and mine tailings are character-ized, improved test protocols, and enhanced physics-informed data interpretation.Tailings dams and ash ponds experience various time-dependent processes, such as gradual strains associated with pore-water pressure dissipation or creep, internal ero-sion and piping, and the migration of fine grains that can cause the gradual clogging of internal drainage pathways and the rise of fluid pressure within the impoundment. These processes occur throughout the entire life of the structure, and their effects can take years or even decades to become evident. By contrast, most laboratory experiments con-ducted on tailings and coal ash samples last from a few hours to a couple of weeks. This discrepancy in time scales points to the need for continuous monitoring of these structures to properly assess their evolution in time.The large mudflows that follow dam fail-ures imply the presence of loose, water-satu-rated sediments that want to contract upon shear: Water cannot drain fast enough, and grains become temporarily suspended, forming a dense fluid (liquefaction). How-ever, laboratory studies that attempt to emulate the hydraulic deposition used to deliver ash and tailings to the pond often fail to predict liquefaction under monotonic loading until impoundment depths exceed 15 to 20 m. I other words, these studies suggest that the hydraulic deposition of ash and tailings does not necessarily produce a contractive, liquefiable sediment.Thus, more-complex mechanisms are required to explain the time-delayed trig-ger and the mudflow that often follows ash ponds and tailings dam failures. These may involve coupled hydro-chemo-mechanical processes. For example, recent data show that fly ash (fine particles from burnt fuel) can experience early diagenetic cementation after placement (10). This process locks in high porosity during burial and results in a brittle material response that may favor liq-uefaction either by internal fabric collapse or during the shear that follows dam failure.Cemented layers have also been identi-fied in tailings dams (11), suggesting that mine tailings may be prone to similar ef-fects. Shear-wave velocity measurements sense even minor changes in small-strain stiffness and could help to detect incipient diagenetic cementation in the field (10).MANAGEMENT AND OPERATIONAs filling of the storage pond progresses, the retaining dams for tailings or coal combustion residuals are raised by build-ing new dykes on top of the preceding one. Each new dyke may be partially shifted toward the pond (upstream construction), aligned with the preceding dyke (center-line construction), or shifted away from the pond (downstream construction). In centerline and upstream constructions, part of the new dyke thus rests on im-pou ded material. The upstream met od (see the second figure) is the most econom-ical, but it requires appropriate construc-tion practices to create w de beaches made of the coarser tailings close to the dykes to improve drainage and stability, while the finer tailings are deposited in the pond further away from the dyke. Much of the discussion that followed the Brumadinho failure has focused on the upstream con-struction method used at the site (1).Analysis of case histories points to sub-standard management and operational practices, often driven by the pursuit of greater financial returns, which lead mine managers to compromise on the safety of tailings dams (6). In his 1980 classical es-say, the geotechnical engineer Ralph Peck noted the high probability of failure of earth dams (12). He remarked that failures were more often associated with disregard for proper engineering practices, rather than lack of knowledge. Similarly, the International Commission on Large Dams 1King Abdullah University of Science and Technology, Thuwal 23955-6900, Saudi Arabia. 2University of the Witwatersrand, Johannesburg 2000, South Africa. 3Geosyntec Atlanta, GA 30319, USA. Email: carlos.santamarina@kaust.edu.sa10 MAY 2019 • VOL 364 ISSUE 6440 5271910 20101990197019501930050100Released volume (103 m3)150200250300Philippines 1992Brazil 2015 Brazil 2015Canada 2014Philippines 1992Canada 20141910 201019901970195019300500100015002500200030003500Bulgaria 1966Italy 1985Brazil2019China 2008Mexico1937Human life lossBulgaria 1966Italy 1985Brazil2019China 2008Mexico1937Massive damage was caused by the failure of a mine tailings dam near the town of Brumadinho, Brazil. More than 232 lives were lost. Periodic failures of impounded industrial waste have catastrophic social, economic, and environmental consequences worldwide.Catastrophic dam failuresOver the past century, tailings dam and ash pond failures and the resulting fast-moving mudflows have led to a cumulative loss of almost 3000 lives. Data from (3).Published by AAASon June 24, 2020 http://science.sciencemag.org/Downloaded from Figure 1.1: Cumulative rel ased volume and life loss du to tailings dam fa ures. From San-tamarina et al. (2019)1Stava tailings dam failureThis event occurred on 19 July 1985 at 12:23 hours at Stava, near Tesero, in the Trento Region ofItaly. The failure resulted in 268 deaths and considerable damage to property. The upstream methodwas used as method of construction of the tailings facility. In total, 200,000 m3 of tailings flowed 4.2km downstream. The area affected was roughly 43.5 hectares. The resulting slurry wave travelledto Stava at a speed of 30 km/h during the first 26 seconds, later it reached an average of about 90km/h below Stava village.Chandler and Tosatti (1995) indicated that the mobility of the debris was consequence of theliquefaction of tailings once the support provided by the dam was removed. They concluded that thetailings dam was constructed with a low factor of safety and that the failure was probably triggeredby a blocked decant pipe located within the tailings. This trigger caused the water level to buildup inside the tailings facility. The precarious stability of the system could not withstand the rise inwater-levels, and it probably collapsed as a consequence.Figure 1.2: Aerial photo of Stava before (left) and after (right) the tailings dam failure. FromWISE Uranium Project2Sullivan mine tailings slideThis event happened on August 23, 1991, at about 14:00 hours near Kimberly in southeastern BritishColumbia, Canada. The failure took place in the southeastern portion of the Active Iron Pond dyke.This facility was built with conventional disposal of tailings (Davies et al., 1998). After the event,the toe area was displaced about 15 to 45 m downstream of its original position. However, there wasevidence that the affected area was up to about 100 m of the original toe. A full flowslide was notdeveloped, although 75,000 m3 of tailings were involved in the failure (Jefferies and Been, 2016).Davies et al. (1998) indicated that the tailings foundation of the Active Iron Pond dyke lique-fied. The authors also mentioned that potential triggering events for this failure were the additionalloading on the crest due to the construction of the dyke extension, and/or the dynamic loads im-posed near the toe due to the traffic of construction equipment. These triggers caused a toe failurefollowed by the liquefaction of tailings and subsequent sliding.i', , :I1· I I· , I I I I I I If I,' I, Figure 5. Aerial View of Slump Area Figure 6. Ground View of Slump Area 51st Canadian Geotechnical Conference 130 51..,m• Conference Canadienne De Geotechnique � I I I ! I Figure 7. Appearance of Area Following Remediation· 51st Canadian Geotechnical Conference 131 51klm• Conference Canadienne De Geotechnique I I I' ' Figure 1.3: Sullivan failure. Ground view of the movement two days after the event. FromDavies et al. (1998)Merriespruit tailings dam failureThis event occurred on 22 February 1994 at Merriespruit, a suburb in the gold mining town of Vir-ginia in the Free State Province of South Africa. The failure took place after a torrential downpourof about 50 mm. Approximately 600,000 m3 of tailings were released over a distance of 3 km,affecting the village of Merriespruit and claiming the lives of 17 people.3Fourie et al. (2001) argued that flow liquefaction was the cause of the Merriespruit flow failure.The authors established that the trigger mechanism was the storage of an excessive volume of wateron the top of the impoundment. This caused a breach in the outer confining slope by erosion.Laboratory tests and measured insitu void ratios on tailings adjacent to the failure showed that thetailings material had void ratios higher than the critical value and were likely in a metastable state.Fourie et al. (2001) concluded that the tailings that did flow in the Merriespruit dam were in a evenlooser state. Therefore, flowed when confinement was suddenly removed.Following the failure, a number of investigations werelaunched. Three different teams of geotechnical engineersundertook these investigations, one representing the state,one the owner, and one the operator of the dam. A compre-hensive evaluation of the findings of these three teams isgiven by Wagener et al. (1998), and only points relevant tothe discussion of static liquefaction as an explanation for thefailure are repeated here.Operational history of the damConstruction of the tailings dam commenced in 1978,with a low “starter wall” at the low end of the dam perime-ter. The dam was operated in a manner similar to that of thevast majority of gold tailings dams in South Africa, usingthe upstream daywall paddock construction method (ICOLD1982). A blanket drain was provided below sections of thestarter wall with a view to depressing the phreatic surface.Instead of a single ringmain around the dam, as is currentpractice in South Africa, three delivery pipelines were usedduring the last few years of operation of the dam, with deliv-ery stations at about 500 m intervals. Water was extractedfrom the top of the dam using permanent, embedded pen-stocks. No water storage dam was provided and the pen-stocks thus decanted into a sump, from where watergravitated back to the metallurgical plant on demand. Be-cause of the inadequate storage capacity of the sump, this ar-rangement often led to storage of water on top of the tailingsdam, which as described later was a major contributor to thefailure.© 2001 NRC Canada708 Can. Geotech. J. Vol. 38, 2001Fig. 1. Aerial view of the Merriespruit tailings dam failure showing the path of the destructive mudflow that occurred.Fig. 2. Partially destroyed house that was about 350 m from the breach location.I:\cgj\Cgj38\Cgj04\T00-112.vpTuesday, July 31, 2001 12:10:10 PMColor profile: Generic - CMYK US Negative ProofingComposite Default screenFigure 1.4: Aerial view of the Merriespruit dam failure. From Fourie et al. (2001)Fundão tailings dam failureThe failure of the Fundão tailings dam took place on November 5, 2015 at the state of MinaisGerais in Brazil, and resulted in the release of about 32 million m3 of tailings representing 61%of the impoundment contents. This event killed 19 people and caused long-term environmentaldamage to several water channels in the Doce River basin since a massive amount of iron minewaste was released.Morgenstern et al. (2016) mentioned that failure was initiated at the dam’s left abutment andoccurred due to flow liquefaction of the tailings i.e., the mass of soil lost strength and flowed like aliquid. This transformation from solid to liquid was abrupt and occurred in just seconds. The PanelReview (Morgenstern et al., 2016) concluded that the triggering event that initiated the failure wasthe lateral extrusion mechanism. This trig er was developed in the area where slimes rather t ansands were deposit d. In this way, the left abutme t of the Fundão tailings am was constructedabo slimes-rich deposits, thus facilitating the lateral extrusion m chanism.4Fundão Tailings Dam Review Panel Report on the Immediate Causes of the Failure of the Fundão Dam August 25, 2016 Page 18Figure 2-18 Failure initiation sequence By the time the events on Figure 2-18 had occurred, the growing failure would have become apparent to the observers on the crest at locations 4, 5, and 6 as it progressed back behind the crest and into the reservoir. Only then did the central and right sides of the dam begin to disintegrate. A conveyor crossing the Fundão stream channel about 1300 m downstream from the offset crest stopped functioning at 3:49PM, four minutes after failure is reported to have begun at 3:45PM. From this, it is ascertained that the flowslide was moving at about 11 m/s by the time it reached the conveyor. It is calculated 3 of tailngs dam failure statistics. In e ad been Figure Figure 2-19 Fundão damsite and reservoir (a) before, (b) after failureFigure 1.5: Fundão tailings dam before (left) and after (right) the failure. From Morgensternet al. (2016)Feijão tailings dam failureThis event occurred on January 25, 2019 at 12:38 p.m. at Córrego do Feijão Iron Ore Mine (DamI), located 9 km north-east of Brumadinho in the state of Minas Gerais, Brazil. The failure of thisfacility resulted in the release of 9.7 million m3 of tailings, representing roughly 75% of the storedvolume. More than 250 people passed away as a result of the collapse.Robertson et al. (2019) indicated that the material in the dam showed a sudden and significantloss of strength and rapidly become a liquid that flowed downstream. Based on these observations,the authors argued that the Feijaão failure was the result of flow liquefaction within the tailings inthe dam. The remaining component of the failure was the triggering event. Robertson et al. (2019)concluded that the Feijão tailings dam failed due to a critical combination of ongoing internal strainsdue to creep, and a strength reduction due to loss of suction in the unsaturated zone.R port of the Expert Panel on the Technical Causes of the Failure of F ijão Dam I 7CAM1 clearly shows that the failure extends across much of the face of the dam and occurs within the slope of the dam extending from the c est to just above the toe (Figure 9). The slope failure appears to be relativ ly shallow in depth. As the failed mass progressed downstream, the failure regressed into the stored tailings as a series of rapid slips, leaving near-vertical back scarps, as shown in Figures 10 through 12. The capture of further progression of the failure by CAM1 became obscured by the dust and mist that ros above the failing mass. Figure 10: Regression of the Failure 11 s After the First Observed Deformation, as Captured by CAM1 Figure 11: Regression of the Failure 18 s After the First Observed Deformation, as Captured by CAM1 Figure 1.6: Feijão tailings dam 11s after the first observed deformation. From Robertsonet al. (2019)51.2 MotivationRecent catastrophic tailings dam disasters have significantly impacted the environment, human life,and property. These events have attracted the attention of the geotechnical community to determinethe causes that led to the development of such failures, thus providing learning opportunities forpreventing these disasters from happening in the future. In this context, the flow liquefaction behav-ior observed in tailings represents a major concern in mining activities since this phenomenon hasbeen identified as a major contributor to the instability of tailings dams. Examples of tailings damflow liquefaction failures were described in Section 1.1.Predicting the post-failure large deformation of these facilities would allow for quantifying therisk and expected damage. However, a proper numerical simulation of tailings dam failures shouldhave two ingredients: an advance constitutive model that can capture the tailings flow liquefactionbehavior and a numerical platform that can handle the large deformations that are typically involvedduring these disasters. Alonso, (2021) highlighted the added advantage of establishing a joint con-sideration of triggering and motion in a unified framework since it allows linking the pre-failurebehavior with the subsequent “dynamic” stage. Some researchers have recently started to performsuch unified analyses; e.g., Yerro et al., (2016) and Cuomo et al., (2019); however, these analysesare limited to simple material models, and also still not common in the current engineering practice.1.3 ObjectivesThis research is intended to simulate the failure of a tailings dam by capturing the flow liquefac-tion and subsequent large deformations observed during the deformation process. To thoroughlydescribe this event, a material point method (MPM) simulation along with an anisotropic boundingsurface plasticity model is presented in this study. The integration between these two componentswas assessed through a base case study. In conjunction with the overarching goal, the objectives ofthe current study are to:• Integrate the constitutive model SANISAND into the material point method software Anura3D.• Design a stable geometry configuration and triggering events that led the system to experienceflow liquefaction and large deformations.• Evaluate the performance of the designed tailings dam against flow liquefaction by consider-ing different initial values of void ratio and trigger intensities.61.4 Thesis structureThis thesis is organized into seven chapters. The chapters are presented in a sequential manner:Chapter 2 presents a overview of the numerical methods considered in geomechanics to capturelarge deformations. This Chapter describes the MPM method and explain the reasons why thisplatform was considered.Chapter 3 describes the flow liquefaction phenomenon and provides information about the con-stitutive models used in MPM highlighting that the SANISAND model was selected for this study.This Chapter presents an overview of the selected model.Chapter 4 describes the boundary value problem considered in the simulations. Two cases arepresented. This chapter depicts the geometry, material parameters, sequence of modeling, definitionof triggers and the initial and boundary conditions of the problem. It also shows the numericalparameters considered in the analyses.Chapter 5 presents the simulation results for the two cases described in Chapter 4. This Chapterillustrates the deformation process in detail. The kinematics of the motion and stress path duringthe large deformations process are discussed.Chapter 6 presents the assessment of flow liquefaction of the tailings dam described in Chapter4. The role of the initial void ratio and trigger intensity are analyzed and the response of the tailingsdam is evaluated.Finally, the summary and conclusions of this study as well as the recommendations for futureresearch are presented in Chapter 7.7Chapter 2Numerical Methods for LargeDeformations2.1 IntroductionThe purpose of this chapter is to outline the existing numerical methods to simulate large defor-mations problems in geotechnical engineering. These are divided on the basis of the discretizationprocedure into mesh-based and particle-based methods. Each numerical scheme is mentioned and abrief description is provided. Among the different methods, the Material Point Method (MPM) wasselected in this thesis as the numerical modeling scheme to deal with large deformations. Section2.4 shows the main reasons why MPM was preferred.2.2 Mesh-based methodsMesh-based methods are the most popular for modeling deformation problems, and have been usedfor over fifty years. In these methods, the relationship between the deforming continuum and themesh determines the ability of the method to deal with large deformations (Yerro, 2015). Thisrelationship is determined by how the motion is described. In continuum mechanics, there are twotype of views to describe a motion: a Lagrangian description, in which the mesh deforms with thematerial, and a Eulerian description, in which the mesh is fixed and the material deforms in timewith respect to the grid.In the Lagrangian formulation (see Figure 2.1a), time and material coordinates are independent.This enables easy modeling of material free surfaces or multiple materials. Furthermore, it allowsthe modeling of history-dependent constitutive equations. However, an excessive deformation of8the soil may lead eventually to a high distortion of the computational mesh, thus affecting thecomputation due to ill-defined interpolations (Armero and Love, 2003).In the Eulerian formulation (see Figure 2.1b), spatial coordinates are used to track the flow of thematerial in time through the fixed mesh. This approach is used mostly in fluid mechanics and largedeformations can be easily captured. The Eulerian description results in significant difficulties forconventional problems in soil mechanics. The modeling of free boundaries and interfaces betweenthe materials limits its use (Armero and Love, 2003).Chapter 2. Introduction to Comput tional MechanicsFigure 2.1: Spatial discretisation scheme of lagrangian and eulerian mesh-based methods.The Finite Element Method (FEM) originated in the engineering literature and it isthe most used technique in geotechnics for finding approximate solutions to boundaryvalue problems (Zienkiewicz, 1977). The material geometry is replaced by an equivalentlagrangian finite element mesh. The shape of the elements can be very diverse and evencan have curved edges depending on the integration order. State variables are determinedat the nodes of the mesh by interpolation functions. This method is well suited to solvecouped hydromechanical problems, however, despite being a valuable tool in the field ofmodelling, the standard lagrangian FEM has serious limitations in problems involvinglarge deformations of the material. This is due to distortions suffered by the mesh andthat in some cases they lead to large inaccuracies (Armero and Love, 2003).In the field of computational fluid dynamics, the Finite Volume Method (FVM) is used(Eymard et al., 2000). Typically, it uses an eulerian description to describe the continuum.Each element of the computational grid represents a control volume in which the balanceequations are solved in it, hence conservation principles are imposed element by element.The determination of input and output flows through the contour of each cell is the mainpurpose of the calculation.One attempt to modelling large deformations maintaining well defined the boundariesarises from the arbitrary Lagrangian-Eulerian (ALE) methods (Donea et al., 1977, 2004).These make use of the two classical kinematical descriptions of motion, lagrangian andeluerian, by taking the advantages of both. Here, the nodes of the computational meshmay be moved with the continuum, or be held fixed. Therefore, to keep the computa-tional mesh as regular as possible and to avoid mesh entanglement during the calculation,remeshing algorithms are used. However, this requires greater computational time.10(a)Chap er 2. Introduction to Computatio al Mecha icsFigure 2.1: Spatial discretisation scheme of lagrangian and eulerian mesh-based methods.The Finite Element Method (FEM) originated in the engineering literature and it isthe most used technique in geotechnics for finding approximate solutions to boundaryvalue problems (Zienkiewicz, 1977). The material geometry is replaced by an equivalentlagrangian finite element mesh. The shape of the elements can be very diverse and evencan have curved edges depending on the integration order. State variables are determinedat the nodes of the mesh by interpolation functions. This method is well suited to solvecouped hydromechanical problems, however, despite being a valuable tool in the field ofmodelling, the standard lagrangian FEM has serious limitations in problems involvinglarge deformations of the material. This is due to distortions suffered by the mesh andthat in some cases they lead to large inaccuracies (Armero and Love, 2003).In the field of computational fluid dynamics, the Finite Volume Method (FVM) is used(Eymard et al., 2000). Typically, it uses an eulerian description to describe the continuum.Each element of the computational grid represents a control volume in which the balanceequations are solved in it, hence conservation principles are imposed element by element.The determination of input and output flows through the contour of each cell is the mainpurpose of the calculation.One attempt to modelling large deformations maintaining well defined the boundariesarises from the arbitrary Lagrangian-Eulerian (ALE) methods (Donea et al., 1977, 2004).These make use of the two classical kinematical descriptions of motion, lagrangian andeluerian, by taking the advantages of both. Here, the nodes of the computational meshmay be moved with the continuum, or be held fixed. Therefore, to keep the computa-tional mesh as regular as possible and to avoid mesh entanglement during the calculation,remeshing algorithms are used. However, this requires greater computational time.10(b)Figure 2.1: Spatial discretization. (a) Lagrangian (b) Eulerian. From Yerro (2015)The conventional mesh-based methods are the Lagrangian finite element method (FEM) andfinite difference method (FDM). Lagrangian FEM is well suited to solve coupled hydromechani-c l problems. In this technique, the finite interpolation unctions ar built upon the mesh, whichensures compatibility of the interpolation. However, this procedure is not always dv ntageous asone may experience mesh distortion, thus ending the computation altogether or resulting in drasticdeterioration of accuracy. To overcome this issue, adaptive remeshing FEM has become a neces-sity. However, when it comes to large deformations, this approach experiments difficu ties not o lyin the remeshing part, but also mapping the state variables from the old mesh to the new mesh.For this reason, frequent remeshing is not desirable (Li and Liu, 2002). In this context, alternativemesh-based approaches were developed to simulate large deformations problems. A brief overviewis presented.92.2.1 Arbitrary Lagrangian-Eulerian methodThe Arbitrary Lagrangian-Eulerian (ALE) method is one attempt to modeling large deformations,using the best features of the Lagrangian and Eulerian formulations. The purpose of the ALEformulation is to make the mesh independent of the material so that the mesh distortion can beminimized (Li and Liu, 2002). In this method, 3 steps can be considered. First, the material nodesare moved to new positions. Then, the mesh is modified by remeshing algorithms so that the elementquality is maintained until the end of the calculation. Finally, the current solution is transferredfrom the old to the rezoned mesh. In this method, remeshing algorithms plays a relevant role for thequality of results (Konkol, 2014). Application of the ALE method in geomechanics were developedby Di et al., (2007), and Aubram et al. (2015). For more details about ALE formulation the readeris referred to Donea et al., (2004).2.2.2 Coupled Eulerian-Lagrangian methodThe Coupled Eulerian-Lagrangian (CEL) method also considers the coupling between the La-grangian and Eulerian formulations. In this case, there are two kinds of meshes. The material(usually the soil) is discretized with Eulerian mesh, whereas the other part of the geometry is dis-cretized with Lagrangian mesh. The interaction between these two meshes is modeled using acontact algorithm defined by the user. In geotechnical engineering, the CEL method has been usedin penetration problems. See Qiu and Grabe (2011), and Qiu et al. (2011). For more details aboutthe CEL formulation the reader is referred to Noh, (1963).2.3 Particle-based methodsParticle methods are also called meshless or mesh-free methods (Yerro 2015). They discretize acontinuum body with a collection of particles and do not use a rigid mesh connectivity comparedwith the conventional mesh-based methods such as Lagrangian FEM.Li and Liu (2002) described the particle methods as the future generation of computationalmethods that will benefit applications in several fields of engineering. This section provides a de-scription of the most popular particle-based methods which have been applied in soil mechanics.2.3.1 Smoothed Particle Hydrodynamics (SPH)SPH is the earliest meshfree technique. It was developed by Lucy, (1977) and Gingold and Mon-aghan, (1977). SPH was used originally for astrophysical problems (Monaghan, 2005). However,since its development, SPH has been widely used for modeling large deformation problems in ge-omechanics. Soga et al., (2016) listed some of the geotechnical applications where SPH was used10successfully such as the modeling of river levee failures (Mori, 2008), and the simulation of seepagefailure and ground erosion (Maeda and Sakai, 2010). Bui et al., (2011) extended SPH to evaluatethe stability of a slope, and to a simulate discontinuous slope failure.In the SPH method, the motion of a continuum is modeled using a set of moving particles (orpoints) while the definition of a mesh is not required. These particles are assigned a constant mass,and carry all field variable information like density, pressure, velocity, etc., and moves with thematerial velocity. Each particle has a spatial distance, known as the smoothing length, over whichtheir properties are smoothed by a kernel function. The term ’kernel function’ refers to a weightfunction and defines how much of each field variable contributes to the field variable at a point (Buiet al., 2007).Soga et al., (2016) described SPH as a well adapted method to solve flow problems but haslimitations in the application of boundary conditions since it requires the use of ghost nodes. Fur-thermore, SPH is computationally expensive as a result of the search for the neighbouring particles.2.3.2 The Material Point Method (MPM)The Material Point Method (MPM) was developed to capture large deformations in history-dependentmaterials (Sulsky et al. 1994, 1995). It was presented as an extension of the particle-in-cell (PIC)method that was typically used to model fluid flow. In MPM, the continuum is represented byLagrangian points (called material points or particles) that move through an Eulerian fixed mesh.The material points (MPs) carry all physical properties of the continuum such as mass, momen-tum, material parameters, velocities, strains, stresses, external loads as well as state variables of theconstitutive model adopted, whereas the Eulerian mesh carry no permanent information.At the beginning of each time step, the information is transferred from the particles to the nodesof the computational mesh. Figure 2.2a shows the particles and nodes represented by red-coloredcircles (p) and empty squares (i), respectively. The equation of motion is solved for the nodal accel-eration values from which the nodal velocities are computed. The velocity vectors are representedby the blue arrows in Figure 2.2b. These nodal values are used to compute the strains and stressesat the material points as well as to update its position by mapping back the nodal solution to thematerial points (Figure 2.2c). At the end of the time step, the mesh carry no permanent informa-tion. Thus, it can be reset to the initial configuration (Figure 2.2d). Through this approach, MPMcombines the advantages of Eulerian (i.e., fixed finite element grid) and Lagrangian (i.e., movingmaterial points) formulations to simulate large deformations while avoiding the problem of meshdistortion.11ip(a) Quantities are mapped from MPs(p) to mesh nodes (i)pi(b) Solve equilibrium equations at thebackground mesh which deformspi(c) Map data from nodes to MPs.Then, update MPs positionpi(d) The background mesh can be re-defined to its original positionFigure 2.2: MPM procedureFor the sake of clarity, it is important to highlight that in MPM, the nodes are attached to thebackground grid and do not carry permanent information, while the material points represent thecontinuum and carry all physical properties. The background mesh is only used to solve the equationof motion with the information that was initially mapped from the material points. After solvingthe equation of motion, the nodes map back its information to the material points. At this point,the mesh carries no permanent information and can be reconstructed for computational convenience(Sulsky et al., 1995).12Historical developmentsSince the original, or standard, formulation of MPM (Sulsky et al., 1994), there have been manyimprovements and modifications of the method in order to overcome some numerical difficulties.These are related to convergence issues and non-physical noise caused by the interpolation betweennodes and material points, which is inherent to the method (Fern et al., 2019).Bardenhagen and Kober (2004) presented the Generalized Interpolation Material Point Method(GIMP). In this approach, the interpolation functions are not linear as opposed to standard MPM,for which they are linear. Also, particle characteristics functions are introduced to the method.The Dual Domain Material Point Method (DDMPM) is an approach to smooth the gradient of thenodal shape functions (Zhang et al., 2011). Unlike GIMP, no characteristic functions are used. Thisscheme is mostly used in the field of numerical mathematics.Beuth (2007, 2012) proposed a MPM mixed integration to mitigate the noise produced by gridcrossing. In this approach, the internal forces are computed using Gauss integration, which is basedon a fixed number of integration points. The fastest Gauss integration is done by using just oneintegration point. The reader is referred to Beuth (2012) and Fern et al. (2019) for more detailsabout MPM mixed integration.Integration schemesMost MPM implementations consider the dynamics in the formulation and use an explicit inte-gration scheme to discretize the time (Sulsky and Schreyer 1996, Nairn 2003, Zabala García 2010,Al-Kafaji 2013, Yerro 2015 and Ceccato et al. 2021). However, when problems to be solved in-volve very slow process, an implicit time integration scheme has been used by several researchers(Guilkey and Weiss 2001 and Beuth et al. 2008) to reduce computational cost. Fern et al. (2019)indicated that most of the dynamic problems can only be treated properly by an explicit integrationscheme as an implicit approach tends to smooth the solution in areas which might be of interest.Multi-phase formulations in MPMMany problems in geotechnical engineering involve the modeling of saturated soils where the in-teraction between the solid skeleton and the fluid is important. The application of MPM to suchmultiphase problems is recent. In the context of MPM, two different numerical approaches havebeen developed to model the behavior of saturated soils: the two-phase single-point and the two-phase double-point formulations. In the former, each material point represents a portion of saturatedsoil, carrying the information of both phases. The latter considers that the solid skeleton and thefluid phase are represented separately by two sets of material points, thus modeling the saturated13soil as a superposition of two continuous media.The first approach (two-phase single-point formulation) was presented by Zabala García (2010).He considered a u-p approximation for the governing equations where the solid displacement (u)and the pore pressure (p) are the primary variables. This approach was successfully used in theMPM simulation of the Aznalcóllar dam (Zabala and Alonso, 2011). However, the main drawbacksof this formulation are that it neglects the acceleration of fluid, suffer from pressure instability, andmust rely on Rayleigh damping to consider the energy dissipation that is proportional to the velocity.Al-Kafaji (2013) presented another version of the two-phase single-point approach based on thev-w formulation, where v and w denote the velocities of the solid and fluid phases, respectively. Inthis case, all acceleration terms are considered in the governing equations. This approach was suc-cessfully used in a variety of geotechnical situations such as: the modeling of the cone penetration inpartially drained conditions (Ceccato et al., 2016), the simulation of slope failures (Soga et al. 2016and Yerro et al. 2016), and to model the flow liquefaction behavior observed during a flume test(Cuomo et al. 2019). The reader is referred to Jassim et al. (2013) for a compact description of thetwo-phase single-point formulation based on the v-w form governing equations of saturated soil.The two-phase double-point formulation was initially presented by Bandara 2013, and laterextended by Abe et al. 2014 and Martinelli 2016. One important advantage of the double-pointformulation compared to the single-point is that the liquid material points embodies either liquidwithin the pores or free liquid, thus allowing users to simulate the free water in dams. The double-point formulation was considered to simulate a levee failure problem (Bandara and Soga, 2015), asubmerged column collapse (Martinelli and Rohe, 2015), a dike failure (Martinelli et al., 2017), anda levee overtopping problem (Liang et al., 2019).Finally, a recent publication (Ceccato et al., 2021) has developed a computational frameworkbased on MPM to study the coupled seepage-deformation process in unsaturated soils. The contin-uous development of multi-phase formulations in MPM illustrates its popularity in the geotechnicalcommunity.MPM softwareOne of the most popular MPM software available for researchers and practitioners is the open sourcecode Anura3D (Anura3D, 2021). Anura3D is a tool developed by the MPM Research Communityto simulate large deformations problems in geotechnical engineering. This includes the modelingof saturated soils through the two-phase single point formulation.In order to use Anura3D properly, pre-processing and post-processing software are required.The program GiD (Melendo et al., 2021) is the pre-processing software and allows users to definethe geometry, material properties, number of material points per element, boundary conditions as14well as generate unstructured and structured meshes for the problem. GiD is developed by CIMNE(International Centre for Numerical Methods in Engineering) and it is required to prepare the inputfiles for Anura3D. The program used for post-processing is Paraview (Ahrens et al., 2005) whichis an open-source multi-platform data analysis and visualization application that was developedto analyze extremely large datasets. Paraview collects the information provided by Anura3D tofacilitate the interpretation of results.The main features of the Anura3D software to simulate the failure a tailings dams are listed asfollows (Anura3D, 2021):• Allow 2D and 3D implementation of the material point method to capture large deformations.• Simulate the behavior of saturated soils considering the two-phase single-point formulationbased on the v-w approximation for governing equations (see Al-Kafaji 2013, Ceccato 2015and Yerro 2015).• Allow the implementation of user-defined constitutive models according to the ABAQUS usermaterial file standard (UMAT). Such implementation must be written in FORTRAN.• Allow the definition of essential and natural boundary conditions for the solid and fluidphases. For the solid phase, the essential boundary condition is the prescribed velocity whilethe natural boundary condition is the prescribed traction. In case of the fluid phase, the pre-scribed fluid velocity represents the essential boundary condition while the prescribed liquidpressure is the natural boundary condition.• Solve the discretized governing equations using an explicit integration scheme.2.4 Why MPM?MPM (through the software Anura3D) has been selected for modeling the failure of a tailingsdam. The main reasons for using MPM in large deformations problems were summarized by Ban-dara (2013), Yerro (2015) and Soga et al. (2016) as follows:• MPM uses the Lagrangian-Eulerian approach to simulate large deformation problems withoutmesh distortion.• MPM allows users to incorporate advanced history-dependent soil constitutive models sincea Lagrangian description of the motion is considered.15• MPM takes lesser computational time compared to other meshfree methods since particlesinside mesh elements can be easily identified.• MPM allows straightforward application of boundary conditions due to the presence of thebackground mesh. Therefore, the boundary conditions can be directly applied to grid nodesas in the FEM. This is the major advantage of MPM compared to other methods such as SPH.• MPM can be seen as an extension of FEM for large deformation problems. Therefore, theMPM implementation is intuitive for users of FEM.Finally, it is important to note that the number of publications related to MPM is growing fastdemonstrating the great interest of the geotechnical community on this modeling scheme. Duringthe last decade, MPM has been used in a variety of geotechnical situations such as:• The modeling of the dynamic evolution of landslides (Andersen and Andersen, 2010 andLlano-Serna et al., 2016).• The modeling of a river levee failure (Soga et al., 2016).• The run-out analysis of slopes and landslides (Yerro et al., 2016 and Yerro et al., 2019).• The simulation of liquefaction and retrogressive slope failures (Cuomo et al., 2019).• The modeling of earthquake-triggered landslide runout (Alsardi et al., 2021).16Chapter 3Modeling Flow Liquefaction with MPM3.1 IntroductionChapters 1 and 2 described recent failures of tailings dams and the numerical platforms that are ableto simulate large deformations. One important ingredient to provide a realistic simulation of tailingsdam failures is the constitutive model. The features that one is looking for in a model dependson the problem. For the case of tailings dam failures, recent studies (Morgenstern et al. 2016,Robertson et al. 2019) have noted that the tailings experienced large strains (likely reaching thecritical state) and flow liquefaction. To adequately simulate the tailings response, the constitutivemodel chosen must account for a critical state soil mechanics (CSSM) framework and having thecapability of modeling flow liquefaction. In this thesis, the SANISAND constitutive model (Dafaliasand Manzari, 2004) was chosen for the simulations. This chapter presents a description of the flowliquefaction phenomena and the SANISAND constitutive model used to capture such behavior.3.2 Flow liquefactionThe flow liquefaction behavior has been observed to occur in loose saturated clean sands (Cas-tro 1969, Ishihara 1993, Chu 1995, Leong et al. 2000), silty sands (Yamamuro and Lade 1998) andtailings (Schnaid et al. 2013) under undrained conditions. This phenomenon is characterized by asudden increase of pore water pressure with a significant loss of shear strength and occurs withoutwarning, thus the consequences of this phenomena often results in catastrophic failures. Given theimportance and frequent devastating consequences, a large amount of research has been performedon this topic.Critical state soil mechanics (CSSM) offers a clear framework to understand flow liquefaction17and to formulate constitutive models that can reproduce such behavior. This section describes rele-vant concepts of soil mechanics to understand the mechanism of flow liquefaction.3.2.1 Critical stateThe critical state soil mechanics (CSSM) is a general theoretical framework that allows to explainthe stress-strain behavior of soils. This theory is based on the concept of a critical density or criticalvoid ratio, first introduced by Casagrande (1936), who observed that, in shear box tests, loose sandscontracted and dense sands dilated until reaching approximately the same void ratio after largedisplacements (Figure 3.1). This final void ratio was named the critical void ratio and was associatedwith no further volume change once the sand had reached this final state.The concept proposed by Casagrande (1936) was further developed by Taylor (1948) who foundthe critical void ratio to be stress dependant. Roscoe et al. (1958) extended the concept to clayeysoils while Schofield and Wroth (1968) presented a more general approach where the conditionof a soil undergoing shear strains under constant stress and with no volume change is termed thecritical state of the soil. This concept is a key element of the CSSM framework. Furthermore,Schofield and Wroth (1968) introduced the concept of the critical state line or critical state locus(CSL) by establishing a relationship between critical void ratio and mean effective stress. The CSLis typically approximated as a straight line in the semi-log space. However, the CSL is often bilinearor curved as described in Carrera et al. (2011) for the Stava tailings (Figure 3.2). Gens Solé (2019)highlighted that early applications of the critical state or critical void ratio concepts were made toinvestigate the flow liquefaction of the Fort Peck dam (Casagrande 1976), thus establishing that theCSSM framework was used to describe the flow liquefaction phenomena from the start.3.2.2 The state parameterThe state parameter represents a important contribution for characterising the behavior of soils. Itwas first introduced by Been and Jefferies (1985) in order to provide a physical parameter that cancharacterize the mechanical behavior of granular materials. The state parameter Ψ is defined inthe e vs p′ diagram (Figure 3.3) as the difference between the current void ratio (e0) and the voidratio at the critical state (ec) at the same effective stress p′. Been and Jefferies (1985) highlightedthat the state parameter combines the effects of density and stress level in a unique parameter foreach material, thus providing a rational engineering approach to the construction of structures usingundensified hydraulic sand fills.The state parameter provides a method to relate shear parameters (dilatancy, friction angle, peakstress ratio). For instance, Been and Jefferies (1986) have presented data of peak angles of shearing18Ol�E.CT SHEA� TESTS DRY OR. DRAINED -ZERO PORE. PRes.svR,s QL � De ().f 6'n (c) y 1,0 JO K<j/cm2 (/09) IIJ � ... V) ci ct IL.I X W\ { PEAK 5T�ENGiH -r,...I(DENSE SANO --- --- ----ULTIMATE STRENGTH t'u LOOSE SAND (a) z' M' NO VOLU PISPLACEMENT d Loose. SAND RE"SPONSE - - -CRITICALVOID RATIO DENSE. SAND DILATIVE RESPONSE. (b) DISPLACEMENT iFigure 3.1: Critical void ratio from direct shear tests (Casagrande, 1976)1·00·80·60·4e0·1 1 10 100 1000 10000p: kPa(b)Compressive andgenerally strain hardeningInstability liquefaction zone1·2 1D-NCL sandPeak pointsCS sandInstability pointCompressive,strain softeningFigure 3.2: Critical state (CS) for the Stava tailings. Modified from Carrera et al. (2011)resistance as a function of state parameterΨ. Because of this capabilities, some constitutive modelsconsider the dependence of dilatancy and peak stress ratio on the state parameter since it unifies theconstitutive description of sands at various levels of density and confining pressure with the sameset of model constants (Dafalias and Manzari, 2004).19292 Soil mechanics: concepts and applications© 2010 Taylor & Francis Group, LLCBeen and Jefferies (1985) proposed the use of a state parameter ψ = v − vc to describe thestate of a sand, where vc is the speci c volume at the critical state (Figure 5.52). This is more satisfactory than the density index ID = (e − emin)/(emax − emin), because ψ takes into account the effects of the con ning pressure p′. Muir-Wood (1990) suggests that the in uence of mineralogy, angularity and particle size could be accounted for by normalising ψ by divid-ing it by (vmax − vmin).Jefferies (1993) went on to develop a critical state model for sand called Nor-Sand, which separates the volumetric state of the soil from the degree of overconsolidation and uses a hardening law based on ψ rather than on the speci c volume v to de ne the size of the cur-rent yield surface. A conceptually simpler example of a critical-state based model for sand, with essentially similar features to Nor-Sand, was proposed by Muir-Wood et al. (1994) and is summarised in Section 5.21.5.21 NON-LINEAR SOIL MODELSSoils only behave truly elastically at very small strains (Simpson et al., 1979). Furthermore, the stiffness of soil decreases rapidly with increasing strain following a change in the direc-tion of the stress path (e.g. loading followed by unloading): Jardine et al. (1984); Atkinson et al. (1990). If realistic predictions of the deformations associated with complex geotechni-cal constructions are required, it is necessary to take these factors into account in the model used to describe the stress-strain behaviour of the soil.One possibility is to relate the shear stiffness of the soil to the current stress and the strain following the last reversal in the direction of the stress path, by means of an empirical curve tted to triaxial stress-strain data (Duncan and Chang, 1970; Jardine et al., 1986). Usually the shape of the curve is hyperbolic, expressed in the formε( )σ − σ = ε+a b1 311(5.61)where σ1 and σ3 are the major and minor principal (total or effective) stresses, ε1 is the major principal strain (i.e. the axial strain in a triaxial compression test), and a and b are curve- t parameters determined empirically (Duncan and Chang, 1970, following Kondner, 1963; Kondner and Zelasko, 1963).A drawback of this approach is that it does not allow for the development of a peak strength, as the nal strength is approached asymptotically with increasing strain. This can Ψ = e0– ecp′ ln p′OCurrent meaneffectivestressCriticalstate lineTypicalstate pointee0ecFigure 5.52 Definition of the state parameter ψ for sand.Figure 3.3: Definition of state parameter. Modified from Powrie (2018)3.2.3 U rained softeningThe mechanism observed during flow liquefaction is illustrated by the undrained stress path andstrain curve of Figure 3.4. The key feature in the stress path is that the deviatoric stress reachesa peak and the material experience strain softening afterwards, thus reducing the shear strength ofsoils. Vaid and Chern (1985) not d that the stress-strain behavior of saturated san s under undrainedmonotonic loading was found to correspond to ne of the thre respons s illustrated in Figure 3.5.Keynot Lecture – Antonio Gens ECSMGE-2019 – Proceedings 6 IGS a) b) Figure 7. Undrained softening: a) effective stress path; b) stress-strain curve The association of flow liquefaction with undrained softening and pore pressure increase was already clearly identified in the stress-controlled tests carried out by Castro (1969) (Figure 8). It should be noted that, because of the stress-controlled nature of the tests, the strain rate after peak is several orders of magnitude larger than the strain rate before the peak, raising some questions about potential rate effects. The same undrained softening mechanism was also recognized by Bishop et al. (1969) and Bishop (1973) as the ultimate cause underlying the Aberfan disaster. To quantify the magnitude of softening, Bishop (1973) proposed the definition of a brittleness index that, for an undrained case, reads: ( ) ( )( )u f u rBu fc cIc−= (1) where (cu)f is the peak undrained shear strength and (cu)r is the residual (critical state) undrained shear strength. The peak point and the critical state have in fact very different characteristics as discussed in the next two subsections. Figure 8. Results of a stress controlled test on loose Banding sand (modified from Bishop 1973, data from Castro 1969): deviator stress-axial strain curve and pore pressure-axial strain curve 2.3 Critical/steady state As indicated above, the final state of the soil after undrained softening corresponds closely to the definition of critical state (Jefferies and Been 2006). In terms of effective stress, the mobilised friction angle corresponds to that of critical state (or constant volume). In the context of flow liquefaction, the final state has also been called steady state. It has been defined by Poulos (1981) as follows: ‘the steady state of deformation of any mass of particles is that state in which the mass is PeakCSPeakCSKeynote Lecture – Antonio Gens ECSMGE-2019 – Proceedings 6 IGS a) b) Figure 7. Undrained softening: a) effective stress path; b) stress-strain curve The association of flow liquefaction with undrained softening and pore pressure increase was already clearly identified in the stress-co trolled tests ca ried out by Castro (1969) (Figure 8). It should be noted that, because of the stress-controlled nature of the tests, the strain rate after peak is s ver l orders of magnitude larger than th strain rat before the peak, raising some questions about potential rat effects. The a undrained s ftening mechanism was lso recognized by Bishop et al. (1969) and Bishop (1973) as the ultimat cause underlying the Aberfan disaster. To quantify the magnitude of soft ning, Bishop (1973) proposed the definitionof a br ttleness index that, for an undrain d case,reads: ( ) ( )( )u f u rBu fc cIc−= (1) where (cu)f is the peak undrained shear strength and (cu)r is the residual (critical state) undrained shear strength. The peak point and the critical state have in fact very different characteristics as discussed in the next two subsections. Figure 8. Results of a stress controlled test on loose Banding sand (modified from Bishop 1973, data from Castro 1969): deviator stress-axial strain curve and pore pressure-axial strain curve 2.3 Critical/steady state As indicated above, the final state of the soil after undrained softening corresponds clos ly to the definition of critical state (Jefferies and Been 2006). In terms of effective stress, the mobilis d friction angle corresponds to that of critical state (or constant v lume). In the context of flow liquefactio , the final state as also been called steady state. I has been defined by Poulos (1981) as follows: ‘the steady state of eformation of any mass of particle is that state in which he mass is PeakCSPeakCSFigure 3.4: Undrai ed softenin . From Gens Solé (2019)20Type 1 (flow behavior)Type 1 is sometimes defined as static (flow) liquefaction (Yamamuro and Lade, 1998) or true liq-uefaction (Chern, 1985). It corresponds to the behavior described in Figure 3.4. In this case, loosesands are susceptible to lose almost all of their shear strength. Flow liquefaction happens when theinitial state is sufficiently located above the CSL in the e− p plane (Lashkari, 2016). If liquefactionoccurs (curve 1 in Figure 3.5), the sand eventually reaches the state state. The effective stress paththen terminates at m while steady state deformation occurs without limit (Vaid and Chern, 1985).Type 2 (limited flow behavior)Type 2 is sometimes called limited liquefaction (Vaid and Chern, 1985). This type represents thestrain softening response. Medium-loose sands display a initial remarkable contraction along witha peak in shear strength. Limited liquefaction happens when the sand state is initially placed veryclose to the CSL in the e− p plane (Lashkari, 2016). If limited liquefaction occurs, the arrestof strain softening response corresponds to a minimum in deviator stress (point n, Figure 3.5). Theeffective stress path at that point takes a sudden turnaround. This state at point n has been consideredas phase transformation (PT) by Ishihara et al. (1975).Type 3 (non-flow behavior)Type 3 represents the strain hardening response. Sand which responds to undrained loading in thismanner is called dilative. A state of PT can also be defned for such sand. Vaid and Chern (1985)noted that for mildly dilative sands, the effective stress path may show a recognizable turnaround instress path similar to the case of limited flow or limited liquefaction.. However, such a turnaroundmay not be discernible for highly dilative states. Sands in dense state exhibit a tendency towarddilation and strain hardening up to the critical state.On the other hand, it is important to mention that flow liquefaction has been mainly studiedunder undrained conditions. However, this behavior can be triggered under completely drainedconditions (Olson et al., 2000). In this context, Chu et al. (2003) highlighted that the instabilitybehavior of sands can be triggered under drained or undrained conditions defining instability as abehavior in which large plastic strains are generated rapidly due to the inability of a soil element tosustain a given stress or load. Leroueil et al. (2009) described a clear example of flow liquefactionunder drained conditions in the context of slopes subjected to pore water pressure increase. In thiscase, the total stress and shear stresses remain essentially constant but the mean effective stressdecrease which represents a drained test with constant q value and decreasing p′.21'O Figure 3.5: Characteristic behavior of saturated sand under undrained monotonic loading.From Vaid and Chern (1985)3.3 Constitutive models and MPMIn the MPM framework, the use of constitutive models that can capture the flow liquefaction be-havior described in Section 3.2 is limited. In fact, it is rare to find the use of advanced constitutivemodels in MPM. Yerro et al. (2015) and Yerro et al. (2016) used a Mohr-Coulomb Strain Softeningmodel to simulate progressive failures and the strength loss after peak strength conditions, Mastet al. (2015) suggested a non-associative variant of the Matsuoka-Nakai constitutive framework tosimulate column collapses, Fern and Soga (2016) used a Nor-Sand constitutive model to evaluatethe role of the constitutive models in the collapse of granular columns using the one-phase single-point MPM, and more recently Cuomo et al. (2019) selected a hypoplastic model to capture theflow liquefaction behavior observed during a well-documented flume test. In engineering practice,Robertson et al. (2019) used the Mohr-Coulomb model to perform a runout analysis of the recentlarge-scale tailings dam failure occurred at Brumadinho (Feijão Dam). It should be noted that, inMPM, the use of constitutive models that consider the critical state soil mechanics (CSSM) frame-work and capture flow liquefaction is limited. However, these features are required for a propersimulation of tailings dam failures.22This thesis considers the SANISAND model in the simulations since it is well-suited to capturethe features described before, and also provides additional advantages that make its use attractive inthe MPM framework. They are described as follows:• is based on a multiaxial formulation facilitating its use in 2D and 3D applications since it isnot limited to only triaxial stress paths.• is able to reproduce with one single set of parameters the monotonic behavior of sands atdifferent levels of pressures and densities.• includes the concept of the bounding surface plasticity, thus allowing the modeling of strainsoftening.3.4 SANISAND3.4.1 General descriptionThe SANISAND constitutive model was originally presented by Manzari and Dafalias (1997) con-sidering the two-surface plasticity formulation and the critical state soil mechanics (CSSM) frame-work. These publications formed the basis for a SANISAND class of models (simple anisotropicsand constitutive models) whose capabilities to reproduce the behavior of sands in monotonic load-ing have been recognized through different publications (Dafalias and Manzari 2004, Taiebat andDafalias 2008, and Taiebat et al. 2010).The SANISAND class was proposed in an effort to simulate the response of sands under mono-tonic and cyclic loading using the same constitutive parameters at any pressure and density (Man-zari and Dafalias 1997). Various extensions to the original formulation have been developed overthe years to capture different physical response observed in laboratory tests. In this context, thisstudy considers the version with fabric change effects (Dafalias and Manzari 2004) referred with thegeneric name of SANISAND hereafter. Recent publications developed by the Theoretical and Ap-plied Geomechanics (TAG) research group at UBC have added new ingredients to SANISAND. Thisis the case of the semifluidized concept that was introduced by Barrero et al. (2020) (SANISAND-Sf). The concepts of memory surface (Corti et al. 2016) and semifluidized state were incorporatedby Yang et al. (2020) (SANISAND-MSf).Yang et al. (2019) presented an schematic diagram of the SANISAND model in the stress space(Figure 3.6). The model comprises of four major cone-type surfaces. The narrow Lode angleindependent yield surface, and three Lode angle dependent bounding, critical and dilatancy surfaces23as illustrated in Figure 3.6. The yield surface has a small opening and is allowed to rotate in thestress space exhibiting isotropic and kinematic hardening. Furthermore, no plastic deformation willoccur under constant stress ratio as the stress rate remains inside the yield surface. The boundingsurface guarantees a smooth transition from elastic to the elastoplastic domain and controls thekinematic hardening of the yield surface. The dilatancy surface represents the phase transformationline described in Ishihara et al. (1975), within this surface the model response is contractive whileoutside the response is dilative. Finally, the opening of the bounding and dilatancy surface dependson the state parameter such that, at the critical state, the surfaces collapse onto the fixed critical stateenvelope.which the response is contractive and dilative, respectively. Theopening sizes of these latter surfaces are determined through the stateparameter, with the critical state envelope acting as an attractor. For astate parameter of zero, both surfaces coincide with the critical stateenvelope.3.3. ComparisonThe selected versions of SANICLAY and SANISAND models similarlyadopt the classical Lode angle dependent critical state theory, i.e., un-ique existence of critical state in the space of mean effective stress, p,and void ratio, e. The bounding surface formulation in both modelsrelies on the stress reversals as an indication for initiation of newloading process enabling both of these models to mimic the change ofstiffness in unloading and reloading events as observed in experiments.However, distance dependence of plastic modulus in their boundingsurface formulation differs in the sense that SANICLAY uses measure-ment of Euclidean distance in the stress space, while SANISAND usesmeasurement in the stress ratio space. This is a result of a basic dif-ference between clay and sand behaviors based on their simplifiedidealization; i.e. clays as assemblies of soft particles and sands as as-semblies of hard particles. In clays, plasticity can take place for loadingin any loading direction, i.e. either with changes of stress ratio or withconstant stress ratio, while in sands plasticity usually takes place forchanges of stress ratio, as long as particle crushing mechanism is notinvolved. Another substantial difference between these two models isthe flow rule that determines plastic strain increments. In SANICLAY,the flow rule is explicitly defined by means of plastic potential whosegradient provides the direction of plastic strain increment and therebythe ratio between its volumetric and deviatoric counterparts. In SANI-SAND, however, there is no explicit definition of plastic potential; in-stead the dilatancy surface is used for defining the plastic volumetricstrain increment, and a deviatoric non-associated flow rule is in-troduced for the plastic deviatoric strain increment in the stress ratiospace. Specific features are present in SANICLAY and SANISANDmodels with some success for reproducing strain softening and cyclicliquefaction, respectively, as observed in the experiments. In SANIC-LAY, a scalar-valued damage parameter is used evolving with magni-tude of plastic deviatoric strains. The damage parameter directly affectsplastic modulus, and consequently plastic volumetric and deviatoricstrain increments. In SANISAND, fabric dilatancy tensor is used evol-ving with the plastic volumetric strain increment with respect to theloading direction. The fabric dilatancy tensor directly affects the con-traction tendency in reverse loading that follows dilation regime.These two models have shown reasonable performance inreproducing some aspects of response in both monotonic and cyclicresponse of clays and sands, mainly in triaxial loading conditions [e.g.,[31,21,36,22,24]]. The validation studies are however still limited, andneed to be expanded to a broader range of loading conditions. With thisbackground in mind, in the present study these models are calibratedand used for evaluation of their performance in a number of verycomplex multidirectional cyclic shearing scenarios with reference to thetests described in the previous section.4. Calibration of the modelsThe SANICLAY and SANISAND models are calibrated using experi-mental data of Gulf of Mexico clay [37,18] and Monterey sand [38,17],respectively. The calibration process in each case is briefly discussed inthis section, with emphases on the parameters responsible for cyclicloading simulations. The parameters are calibrated essentially against acyclic 1-D linear test (DSS) in each case, and then the calibrated modelsare used in the next section for simulations of various multidirectionalcyclic shear tests in the corresponding database.4.1. Gulf of Mexico clayThe Gulf of Mexico clay database used in this study is from a testingprogram at the University of Texas, Austin [18]. Detailed description ofsampling and geotechnical characterization of the tested undisturbedsamples is available in [38,17]. The deposits in the area from whichsamples were taken consist of uniform surficial highly plastic Holoceneclays. Samples were retrieved from 30 different locations with varyingdegrees of seabed sloping. As the shallower samples were too soft, thetesting program was performed on samples below 7m depth. Theplasticity index for these samples was in the range of 40–55, consistingof 51–72% of clay particles.In addition to multidirectional cyclic shear tests, which are of directinterest in the present study, Rutherford [18] also conducted constantrate of strain consolidation tests of the samples, and using Casagrande'sgraphical method estimated a pre-consolidation stress of around 72 kPa.Murali [37] also conducted a series of undrained monotonic triaxialtests on these samples starting from either isotropic or K0 consolidationstates. Calibration of the model parameters is following the processdescribed by Seidalinov and Taiebat [24], unless otherwise specified.The calibrated model parameters are presented in Table 4. The slopes ofthe normal consolidation and the rebound lines in the e– pln space,denoted by λ and κ, respectively, are calibrated from the constant rateof strain consolidation tests. The experimental data is available ine– σln a space and the procedure described by Taiebat et al. [22] is fol-lowed to transform the data into the e– pln space. The critical statestress ratios Mc and Me, corresponding to the triaxial compression andextension, respectively, are calibrated from undrained triaxial com-pression and extension tests on CK0 specimens sheared from stressstates corresponding to a range of OCR values. The stress ratio N at theapex of the bounding surface is calibrated from the stress path recordedfor undrained triaxial compression test on CK0 normally consolidatedTable 4Model parameters of SANICLAY for Gulf of Mexico clay.Model parameters Symbol ValueElasticity κ 0.06ν 0.2CSL λ 0.247Mc 1.1Me 0.81Bounding surface N 1.1h0 80ad 4Rotational hardening C 3x 1.75Fig. 4. Schematic diagram of the SANISAND model in the principal stress space.M. Yang et al. Soil Dynamics and Earthquake Engineering 124 (2019) 230–258235Figure 3.6: Schematic diagram of the SANISAND model (Yang et al. 2019)3.4.2 Model overviewA compact version of SANISAND mathematical formulation is shown in this section based on thesummary presented in Taiebat et al. (2010). Both stress and strain quatities are assumed positive incompression. All stress components should be considered as effective stress. Finally, in terms ofnotation, tensor quantities are denoted by bold-faced symbols.24ElasticityThe nonlinear elastic response of SANISAND is assumed to be hypoelastic and defined bye˙ =s˙2Gε˙ev =−p˙K(3.1)where s is the deviatoric stress tensor; p the mean effective pressure; e the deviator strain tensor; εvthe volumetric strain; superscript e denoting elastic; superposed dot denoting the rate; and G, K thehypoelastic shear and bulk moduli, respectively, given byG = G0 patm(2.97− e)2(1+ e)(ppatm)0.5 K =2(1+ v)3(1−2v)G (3.2)where G0 is a dimensionless model parameter, v is the Poisson’s ratio, e is the void ratio and patm isthe atmospheric pressure.Yield surfaceThe yield surface is defined byf = [(s− pα ) : (s− pα )]1/2−√2/3mp = 0 (3.3)where m is a dimensionless material constant representing the opening size of the yield surface andα is the deviatoric back-stress ratio. This surface can be visualized as a cone in Figure 3.6. Thegradient of the yield surface in stress space is computed according to∂ f∂σ= n− 12(n : r)I n =r−α√2/3m(3.4)where n is the deviatoric unit tensor (n : n = 1) that is normal to f = 0, r is the stress-ratio tensordefined as r = s/p and I is the identity second-order tensor. The tensor n of Equation 3.3 is used tocalculate the Lode angle θ as:cos3θ =−√6tr(n3) (3.5)The critical stress ratio for a given Lode angle θ is denoted by M and has been obtained byinterpolation between values of Mc at θ = 0 (triaxial compression) and Me at θ = π/3 (triaxialextension) according toMCS = g(θ ,c)Mc g(θ ,c) =2c(1+ c)− (1− c)cos3θ c =MeMc(3.6)25where Mc and Me are model parameters.Bounding, dilatancy and critical surfacesSANISAND uses three concentric and homologous surfaces: the bounding, dilatancy, and critical,that are considered in the deviatoric stress ratio π−plane. The image points α bθ , α dθ , and α cθ onthese surfaces are defined by central mapping from the origin along the direction n asα bθ =√23αbθn αbθ = g(θ ,c)Mexp(−nbψ)−m (3.7)α dθ =√23αdθ n αdθ = g(θ ,c)Mexp(ndψ)−m (3.8)α cθ =√23αcθn αcθ = g(θ ,c)M−m (3.9)where ψ = e− ec is the state parameter; ec is the analytical expression of the critical state linedefined as:ec = ec0−λc( pcpatm )ε (3.10)with constants e0, λ and ε; and nb and nd are material constants.Plastic flowSANISAND includes a non-associative flow rule allowing realistic evaluations of plastic strain in-crements according toε˙ p = ⟨L⟩R (3.11)R = Bn+C(n2− 13 I)− 13 DI (3.12)B = 1+321− ccgcos3θ C = 3√321− ccg (3.13)where R is the direction of the ε˙ p; L is the loading index; and ⟨⟩ the Macauley brackets. Thedilatancy coefficient D in Equation 3.12 is defined byD = Ad(α dθ −α ) : n (3.14)Ad = A0(1+ ⟨z : n⟩) (3.15)26where A0 is a material constant. The so-called fabric-dilatancy internal variable z is introduced toaccount for the effect of fabric changes during the reverse loading followed by a plastic dilatantphase.Evolution lawsSANISAND has two tensorial internal variables, namely, the back-stress ratio tensor α and thefabric-dilatancy tensor z. The evolution law for the back-stress ratio α is a function of the distancebetween bounding and current back-stress ratio in the form ofα˙ = ⟨L⟩2/3h(α bθ −α ) (3.16)where h is the hardening coefficient. A constant value of h yield reasonable simulations. However,better results are obtained if h is made function of the current state variable ash =b0(α −α in) : n (3.17)b0 = G0h0(1− che)(p/patm)−1/2 (3.18)where α in is the initial value of α at initiation of a new loading process and is updated to the newvalue when the denominator of Equation 3.17 becomes negative. The parameters h0 and ch arematerial constants.Finally the evolution law for the fabric-dilatancy z is introduced asz˙ =−cz⟨−LD⟩(zmaxn+ z) (3.19)with cz and zmax as the material constants that control the maximum value of z and its pace ofevolution.3.4.3 Verification in Anura3DThe SANISAND implementation considered in this thesis was the freely distributed version of theuser-defined soil model subroutine UMAT written originally by M. Martinelli, C. Miriano and C.Tamagnini, and then updated by D. Mašín. This implementation is available at the soilmodels.infowebsite (Mašín, 2017). A contribution of this thesis was to include this implementation in thesoftware Anura3D to simulate the response of contractive tailings. To verify the correctness ofsuch implementation in Anura3D, a comparison between the undrained triaxial test simulationsperformed in the programs Anura3D and Incremental Driver is presented in Figure 3.7.27      q (%)q (kPa),QFUHPHQWDO'ULYHU$QXUD'    p (kPa)q (kPa),QFUHPHQWDO'ULYHU$QXUD'      q (%)u (kPa),QFUHPHQWDO'ULYHU$QXUD'(a)      q (%)q (kPa),QFUHPHQWDO'ULYHU$QXUD'    p (kPa)q (kPa),QFUHPHQWDO'ULYHU$QXUD'      q (%)u (kPa),QFUHPHQWDO'ULYHU$QXUD'(b)Figure 3.7: Undrained triaxial simulations on Toyoura Sand. (a) e0 = 0.833 (b) e0 = 0.90728Incremental Driver is a program written by Niemunis (Niemunis, 2017) that allows the testing,at one stress point, of the Abaqus implementations of different constitutive models (the so-calledUMAT’s) (Wichtmann et al., 2019). The main task of Incremental Driver is to force UMAT’s tofollow a prescribed loading path formulated in strain or stress components. For the verificationprocess, the undrained triaxial test path was considered.To complete the analysis, undrained triaxial tests were simulated using the FEM option of theAnura3D code. In this case, it is not possible to make simulations using only one stress point. Forthis reason, a element test model with two stress points was prepared in the software consideringappropriate initial and boundary conditions. Figure 3.7 shows the results of only one stress pointsince they were exactly the same.The model parameters considered in both simulations correspond to Toyoura sand and wereadopted from Dafalias and Manzari, (2004). Figure 3.7 shows the results of the verification pro-cess. It should be noted that there is good agreement between the simulations performed with twodifferent programs. Furthermore, Figure 3.7 illustrates the capabilities of SANISAND to capturethe strain softening behavior associated with the flow liquefaction phenomena.29Chapter 4Description of the numerical model4.1 Model setup4.1.1 Geometry descriptionThe geometry and layering of the hypothetical tailings dam analyzed in this study were motivatedby a historical tailings dam that experienced flow liquefaction (Davies et al. 1998). Such featuresare described in Figure 4.1. The distance between the toe of the buttress and the lower raising liftsis roughly 40 m. The thickness of the foundation is 1 m while its extension is 24 m. The thicknessof the lower and upper tailings are 7 and 8 m, respectively. The extension of the model is 130 m.The coordinate system is placed with its X and Y origins at the lower left corner of the model.              [ P \ P %XWWUHVV5DLVLQJ/LIWV8SSHU7DLOLQJV)RXQGDWLRQ /RZHU7DLOLQJV 5DLVLQJ/LIWVFigure 4.1: Geometry of the model4.1.2 ZoningFigure 4.1 describes all material zones considered in this study. The synthetic geometry includes abuttress located at the toe of the raising lifts. This buttress has been built to provide initial stability30to the system. Moreover, a foundation layer is considered to represent a weak foundation thatmay cause instability to the system. This mode of failure has been observed in recent tailings damdisasters such as Mount Polley (Morgenstern et al. 2015) and Cadia (Morgenstern et al. 2019).The tailings material was subdivided in two layers: upper and lower tailings which distinguishtwo type of tailings due to different deposition periods. The upper tailings is considered as an activeand relatively dense tailings whilst the lower tailings represents an old and liquefiable tailings.Furthermore, upper and lower raising lifts were considered to represent the dykes that are used tostorage the tailings.4.1.3 Background meshFor a 2D simulation in Anura3D, only triangular elements are available for the background mesh.This is a plain strain simulation. The mesh has been refined in the region where the flow lique-faction behavior is expected to get more accurate results and optimise the computational cost. Anunstructured mesh is considered with an element size of 1.5 m in the area where the refinement isnot considered while a size of 1.0 m is considered in the lower tailings. Initially, 3 material pointsper element are used within each element. Figure 4.2 shows the background mesh and the initialdistribution of material points considered for the calculations. The reader should note that the meshhas been extended more than the tailings structure itself. The MPM needs empty elements wherethe material points are expected to be as a consequence of the deformation process. The backgroundmesh is composed of 3,079 triangular elements and 1,652 nodes.              (a)              (b)Figure 4.2: Background mesh: (a) Triangular elements (b) Initial material points distribution314.1.4 Sequence of modelingThree simulation stages have been considered in the calculation procedure. First, the initial effectivestress and pore-water pressures were computed from a preliminary quasi-static simulation using theMohr-Coulomb constitutive model in the tailings units (Stage 1). Then, the model was switched toSANISAND and the equilibrium was verified with another quasi-static calculation (Stage 2). Thissequence was used to mitigate the numerical difficulties observed in the SANISAND constitutivemodel at initial low-stress levels. Finally, the trigger event is numerically applied and the deforma-tion process starts (Stage 3). Figure 4.3 describes each of them.Bring the system to equilibrium using Mohr-CoulombAssign stress considering weight, gravity and K0=0.5Change of constitutive model from Mohr-Coulomb to SANISAND7ULJJHUDSSOLFDWLRQBring the system to equilibrium using SANISANDSTAGE STAGE 1'HIRUPDWLRQSKDVHZLWK6$1,6$1'STAGE 3Figure 4.3: Simulation stages4.1.5 Definition of triggersA triggering event is defined as the event that causes large permanent deformations and destabi-lizes the tailings dam. Such event can be a moment in time or a continual change in the loadingconditions. In this study two triggers were used: the increment of pore water pressure at the baseof the model (trigger 1) and the sudden loss of shear strength in the foundation layer (trigger 2).Both triggers were applied at the beginning of stage 3 and are intended to cause the flow lique-faction of the lower tailings described in Figure 4.1, thus leading to the overall instability of thesystem. The first trigger represents an artificial recharge of the system (water surcharge) which isone of the most common reasons to cause slope instabilities. This trigger has been successfullyused in the MPM simulation of the Selborne experiment (Soga et al. 2016) and the simulation offlowslide flume tests (Cuomo et al. 2019). In both cases, the pore water pressure was applied inone boundary of the model. In this thesis, the applied pore water pressure changes with time, itstarts from 150 kPa and rises linearly during 1 sec, up to 200 kPa. After this time, the pressureon the boundary is maintained constant during all the calculation. The second trigger represents asystem with a soft foundation due to strain accumulation during construction or simply a soft layer32is located in the foundation. From a numerical point of view this trigger is considered by decreasingthe shear strength of the foundation unit in one single step, hence the shear strength loss is assumedto be instantaneous. This type of trigger was successfully considered in the MPM simulation of theOso Landslide (Yerro et al. 2019). In this study, the friction angle of the lower tailings is suddenlyreduced from 36° to 8°.The author recognizes the simplifications adopted in the trigger application. However, regard-less of how the trigger event is numerically imposed to the system, the main goal of this study is toset off the flow liquefaction phenomena in the lower tailings, thus establishing a joint considerationof the pre- and post-failure behavior of tailings structures.7ULJJHU LQFUHPHQWRISRUHZDWHUSUHVVXUH SZ P7ULJJHU /RVVRIVKHDUVWUHQJWKFigure 4.4: Trigger application - Description4.2 Governing equationsThe governing equations of the two-phase single-point MPM were considered in all simulationsstages described in Section 4.1.4. This approach was implemented in the software Anura3D andbriefly described in Section 2.3.2. The MPM formulation is fully dynamic being the momentumbalance the principal governing equation.4.3 Method of solvingAnura3D uses the Euler-Cromer (explicit) integration scheme to solve the discretized form of thegoverning equations. To complement this procedure, numerical parameters that affect the solutionand computational time of the simulations need to be defined.4.3.1 DampingIn order to accelerate the convergence to quasi-static equilibrium or reduce numerical instabilitiesduring calculations, a damping force has been included in the momentum balance equation. In thetwo-phase single-point formulation, the fluid and solid phases are damped separately. This force is33proportional to the corresponding unbalanced force by means a proportional factor α , and acts in theopposite direction of the velocity. A local damping coefficient of 0.75 was applied to stages 1 and 2(Section 4.1.4) to reach equilibrium as fast as possible without taking into account dynamic effects.However, for stage 3 a damping coefficient of 0.05 is considered to reduce numerical instabilities.4.3.2 Time incrementThe explicit integration scheme is conditionally stable since the time increment (or integration time)must be smaller than certain critical time which defines the biggest time increment that can be usedfor a stable calculation.In case of the two-phase single-point formulation, the definition of a critical time step is onerous.Anura3D defines this value considering the simplified analysis performed by Mieremet, 2015. Inthis criteria, the critical time step depends on the size of the mesh, the solid and fluid stiffness, theporosity, the density of the material and the permeability. The bulk modulus of water is reduced bya factor of 100 in order to increase the critical time step to gain computational efficiency with theexplicit integration scheme. Finally, a time increment of 1E-4 s was used in the simulations.4.4 Material parametersThe layers defined in Section 4.1.2 requires a set of parameters depending on the assigned con-stitutive model. Figure 4.5 illustrates the constitutive models assigned to each layer in the modeldomain. SANISAND was considered for the upper and lower tailings (Section 4.4.1) while a MohrCoulomb model with perfect plasticity and associated flow rule was used for the raising lifts, foun-dation and buttress (Section 4.4.2). Mohr Coulomb was used to simulate a rigid behavior that willnot experience softening whilst SANISAND was assigned to the tailings to simulate the contractivebrittle behavior that was observed in recent dam failures. The simulation considers the two-phasesingle point formulation for each layer. For this reason, physical properties related to the solid andfluid phases, and model parameters related to the assigned constitutive model must be defined.              0RKU&RXORPE6$1,6$1'Figure 4.5: Constitutive models in model domain344.4.1 Upper and lower tailingsThe SANISAND model parameters for the upper and lower tailings were obtained considering thecalibration process described in Taiebat et al. (2010) and the laboratory data provided by SRK Con-sulting which corresponds to tailings in Mexico . The data consists of three undrained triaxial testsperformed on polymetallic tailings at different levels of pressures and densities, thus covering a widerange of applicability. Figure 4.7 shows the calibration results considering the set of parameters de-scribed in Table 4.1. As a reference, Table 4.1 includes the SANISAND parameters for Toyoura,Nevada, Ottawa and Karlsruhe sands published by Dafalias and Manzari 2004, Taiebat et al. 2010,Ramirez et al. 2018 and Wichtmann et al. 2019, respectively. Figure 4.7 suggests that the upper andlower tailings may have a strain softening response for a certain range of pressures and initial voidratios.Two cases are considered in the simulations, one for each trigger defined in Section 4.1.5. Figure4.6 illustrates the initial void ratio values assumed for the upper and lower tailings for each scenario.Case 1 uses trigger 1 as triggering mechanism, while case 2 uses trigger 2. The initial void ratios aredescribed in Table 4.2, and were assumed to consider that the upper tailings is denser than the lowertailings. This assumption provides stability to the system. The initial void ratios were assigned onlywhere the SANISAND constitutive model was assumed.              6$1,6$1' e0=0.75 6$1,6$1' e0=0.82 (a)              6$1,6$1' e0=0.75 6$1,6$1' e0=0.87 (b)Figure 4.6: Initial void ratios assumed in the model: (a) Case 1 (b) Case 235It should be noted that the calibration process provides the critical state parameters that definesthe critical state line (CSL) used for the lower and upper tailings. Figure 4.8 illustrates the criticalstate line (CSL) for several tailings material including the one considered in this thesis. It is impor-tant to highlight that the CSL for tailings is steeper than the CSL for clean sands at relatively lowlevel of pressures. Figure 4.8 is a modification from Jefferies and Been (2016), the main changesare the inclusion of the CSL used in this study, the CSL obtained by Schnaid et al. (2013) for goldtailings, and the CSL obtained for Ottawa, Toyoura and Nevada sands which were described inTable 4.1.Finally, the basic properties described in Table 4.2 for the upper and lower tailings have beenassumed taking into account conventional values typically used in geotechnical engineering. It isimportant to mention that considering the parameters defined in Table 4.2, the Darcy permeabilityassumed for the tailings is 1E-2 m/s.Table 4.1: SANISAND model parameters for tailingsParameter Symbol ToyouraSandNevadaSandOttawaSandKarlsruheSandThisthesisElasticity G0 125 150 125 150 125ν 0.05 0.05 0.05 0.05 0.05Critical State Mc 1.25 1.14 1.26 1.34 1.50c 0.712 0.78 0.735 0.70 0.712ec 0.934 0.83 0.78 1.103 1.191λc 0.019 0.027 0.0287 0.122 0.545ξ 0.70 0.45 0.70 0.205 0.14Yield Surface m 0.01 0.02 0.02 0.05 0.02Dilatancy nd 3.5 1.05 2.5 2.0 1.1A0 0.704 0.81 0.626 0.90 0.2Kinematic nb 1.1 2.56 2.3 1.2 1.1Hardening h0 7.05 9.7 5.0 10.5 7.05ch 0.968 1.02 0.968 0.75 1.1Fabric dilatancy zmax 4.0 5.0 11.0 20.0 4.0cz 600 800 500 10000 60036      q (%)q (kPa)H S N3DH S N3DH S N3D6LPXODWLRQ    p (kPa)q (kPa)H S N3DH S N3DH S N3D6LPXODWLRQ      q (%)u (kPa)H S N3DH S N3DH S N3D6LPXODWLRQFigure 4.7: SANISAND calibration for polymetallic tailings370HDQ (IIHFWLYH 6WUHVV  N3D  9RLG5DWLR7DLOLQJVEHDFK$UJHQWLQDPL[HGEHDFK6NRXULHV2[LGHWDLOLQJV1HYDGD&RSSHUWDLOLQJV&RQJDGU\WDLOLQJV6XGEXU\ QLFNHO VDPSOH%7DLOLQJVVDQG2LOVDQGWDLOV6XGEXU\ QLFNHO VDPSOH$7RURPRFKR6FKQDLGHWDO2WWDZDVDQG)7R\RXUDVDQG1HYDGDVDQG7KLVVWXG\Figure 4.8: Critical state line for various tailings (continuous line) and sands (dashed lines).Modified from Jefferies and Been (2016).38Table 4.2: Physical parameters for tailingsParameters Symbol Upper and lowertailingsVoid ratio e (-) 0.75 | 0.82 and 0.86Solid density ρs (kg/m3) 2700Fluid density ρ f (kg/m3) 1000Intrinsic permeability fluid k f (m2) 1.0E-9Bulk modulus fluid K f (kPa) 2.15E4Dynamic viscosity fluid µ f (kPa/s) 1.0E-64.4.2 Raising lifts, foundation and buttressThe model parameters for the raising Lifts, foundation and buttress were described in Table 4.3 andcorrespond to the assigned Mohr-Coulomb constitutive model. For these materials, conventionalgeotechnical properties were selected according to the current engineering practice. The physicalparameters are defined in Table 4.4.Table 4.3: Mechanical parameters for Raising Lift, Foundation and Buttress unitsParameter Symbol Raising Lifts Foundation ButtressYoung Modulus E (kPa) 10000 10000 10000Poisson ratio ν (-) 0.3 0.3 0.3Friction angle φ (°) 38 36 40Cohesion c (kPa) 1 1 1Dilation angle ψ (°) 0 0 0Table 4.4: Basic parameters for Raising Lifts, Foundation and ButtressParameters Symbol Raising lifts, foundationand buttressVoid ratio e (-) 0.50Solid density ρs (kg/m3) 2700Fluid density ρ f (kg/m3) 1000Intrinsic permeability fluid k f (m2) 1.0E-9Bulk modulus fluid K f (kPa) 2.15E4Dynamic viscosity fluid µ f (kPa/s) 1.0E-6394.5 Initial and boundary conditionsIn order to solve the governing equations, initial and boundary conditions must be defined for eachstage of simulation outlined in Section 4.1.4.4.5.1 Stage 1The initial stresses are assigned considering the porosity, density, gravity and the K0 coefficient. Inthis thesis, a K0 value of 0.50 is considered. Figure 4.9 illustrates the stage 1 initial conditions.Figure 4.10 shows the essential and natural boundary conditions for the solid phase. Theyare established by defining the prescribed solid velocities at the borders. The essential boundaryconditions are simulated by rollers on the left and right sides of the model to constrain the horizontaldisplacements. At the top of the model rollers are considered to constrain the vertical displacementwhilst hinges at the bottom are used to prevent both vertical and horizontal displacements. Thenatural boundary conditions establishes zero traction at the surface of the model.Figure 4.11 illustrates the essential and natural boundary conditions for the fluid phase of stage1. These conditions are established by defining the prescribed water velocity at the borders. Theessential boundary conditions define the left and right sides of the model as impermeable bound-aries. The top and bottom sides are considered impermeable meaning that the fluid cannot movethrough these borders. The natural boundary conditions define zero water pressure at the surface ofthe model.       p′ N3D (a)         p w N3D (b)Figure 4.9: Initial conditions - Stage 1: (a) mean effective stress (b) pore-water pressure40vx=0vx=0vy=0vx= vy=0(VVHQWLDO(a)=01DWXUDO(b)Figure 4.10: Solid boundary conditions - Stage 1: (a) Essential (b) Naturalwx=0wx=0wy=0wy=0(VVHQWLDO(a)pw=01DWXUDO(b)Figure 4.11: Fluid boundary conditions - Stage 1: (a) Essential (b) Natural414.5.2 Stage 2The final stress state after the quasi-static simulation performed in stage 1 represent the initialstresses for stage 2. Figure 4.12 shows the initial state for stage 2. The solid and fluid boundaryconditions for this stage are the same that were described in Figures 4.10 and 4.11.       p′ N3D (a)       p w N3D (b)Figure 4.12: Initial conditions - Stage 2: (a) mean effective stress (b) pore-water pressure4.5.3 Stage 3At the beginning of this stage, the trigger event outlined in Section 4.1.5 is applied to the system. Thefinal stress state after the quasi-static simulation performed in stage 2 represent the initial stressesfor stage 3. Figure 4.13 shows the initial state for stage 3.The boundary conditions considered in this stage depend on the trigger application. Two dif-ferent triggers were described in Section 4.1.5: the increment of pore pressure at the base of themodel (trigger 1) and the sudden loss of strength in the foundation (trigger 2). In case trigger 1 isconsidered, the solid boundary conditions are the same as described in Figure 4.10. However, theboundary conditions for the fluid phase are different and are shown in Figure 4.14. If trigger 2 isused, the boundary conditions are the same as described in Figures 4.10 and 4.11.42       p′ N3D (a)       p w N3D (b)Figure 4.13: Initial conditions - Stage 3: (a) mean effective stress (b) pore-water pressurewx=0wx=0wy=0wy=0(VVHQWLDO(a)pw=0pw= variable with time according to Section 4.1.51DWXUDOP(b)Figure 4.14: Fluid boundary conditions - Stage 3: (a) Essential (b) Natural43Chapter 5Simulation results5.1 IntroductionThis Chapter shows the results for the two cases described in Chapter 4. The aim of this Chapter isto give a detailed failure description of the designed problem. These results allow to describe thefailure initiation as well as post-failure dynamics of the tailings dam. In the analysis, the response ofthe lower tailings is of interest since a higher initial void ratio was assumed, and a flow liquefactionresponse may take place. Figure 5.1 shows selected locations within the lower tailings that areconsidered to evaluate the kinematics, the deformation process and the stress path developmentduring the tailings dam motion. The results presented in this Chapter include contour plots whichshow the deformed shape, the evolution of deviatoric shear strain and mean pressures at differenttimes. Moreover, the evolution with time of displacements, velocities, void ratio and mean pressureare illustrated. It was observed that the description of the failure mechanism is the same for the twotriggers considered. For this reason, only the results for trigger 1 are described. However, the resultsfor the two cases are presented.      3 3 33x (m)y (m)Figure 5.1: Control points used in analysis of results445.2 Case 15.2.1 KinematicsThe main advantage of using MPM to simulate the failure of tailings dam is that the post-failurebehavior can be analyzed in conjunction with the pre-failure stage. The calculated displacementsand velocities of control points P1, P2, P3 and P4 are shown in Figure 5.2. When analyzing theseplots it can be noted that the TSF remains stable until t = 2 s. Then, the failure mechanism starts todevelop and the control points P1 and P2 increase quickly their velocity (accelerate) until the peakvelocity is attained at t = 6 s. Afterwards, the velocity decreases its value towards a new equilibriumstate. At t = 13 s the velocity of the control points tends to stabilize (velocity equal to zero). Notethat the time at which each control point starts to accelerate depends on its position within the lowertailings. The control point P3 starts to move at t = 4 s while the control point P4 remains withoutmovement until t = 6 s. Therefore, the TSF deformation process in the lower tailings starts betweenthe coordinates x = 40 and 60 m, and propagates to the right side where the control points P3 andP4 are located.        7LPH V 'LVSODFHPHQW P 3333        7LPH V 9HORFLW\ PV 3333Figure 5.2: Displacement and velocities of control points P1, P2, P3 and P4 for case 15.2.2 Deformation processTo describe the deformation process it is important to clarify the concept of failure. In the MPMliterature, several researchers such as Soga et al. (2016), Alonso (2021) and Ceccato et al. (2021)45use the term failure in the context of the classical or traditional geotechnical analyses (limit equi-librium) where the concept of factor of safety provides a numerical scale to quantify the safety ofgeotechnical structures. In this scenario, the term failure is used to describe the moment where thedriving forces are higher than the resistance forces i.e., when the system is not in equilibrium andstarts to develop displacements. This definition has been considered in this thesis. After the failureinitiation the post-failure dynamics proceeds. Yerro et al. (2016) and Cuomo et al. (2019) referredto the onset of failure as the onset of slope instability. Soga et al. (2016) refer to the deformationprocess, which represents the description of the static and subsequent dynamic stage after failure,as the process of instability. This study considers the term deformation process.Figure 5.3 shows the tailings dam deformed shape at different times when trigger 1 is applied.This plot illustrates the initial and final configuration when equilibrium has been re-established i.e.,when the velocity of the material points considered in the simulation is equal to zero. The horizontaland vertical axes represent the dimensions in m of the model. To present a detailed description ofthe deformation process Figures 5.4 and 5.5 illustrate the evolution of deviatoric strains and meaneffective stress at five different times. Figure 5.4 shows that at t = 4 s there is a complete developmentof shear surface within the tailings dam. As the tailings material flows, the thickness and extensionof the shear zone is increased. Figure 5.5 shows that at t = 4 s the mean effective stress in the shearsurface has been reduced due to the application of trigger 1. As the failure propagates the meanpressures are reduced in the mobilized tailings, facilitating the deformation process of the system.Finally, it is important to note that Figures 5.3, 5.4 and 5.5 describe the performance of the tailingsdam during the pre- failure, the onset of failure, and subsequent motion when trigger 1 is applied.5.2.3 Time evolution of material stateTo describe the time evolution of relevant soil parameters, it is important to recall that the onsetof flow liquefaction is used in this thesis to describe the initiation of strain softening which repre-sents the peak deviatoric stress in the p−q space under undrained loading (Vaid and Chern, 1985).Andrade et al. (2013) treat flow liquefaction as an instability. This thesis uses the term flow lique-faction.The results showed in this section represent the quantities obtained at selected control points. Infact, the control points described in Figure 5.1 are moving (tracking) points where the soil parame-ters are computed and saved at the end of each time step according to the MPM procedure describedin Section 2.3.2. This procedure implies the continuous interaction between the material points andthe nodes of the background mesh which is inherent to MPM.Figures 5.6 and 5.7 illustrates the behavior of tailings during the deformation process. For thispurpose, the control points P2 and P4 defined in Figure 5.1 have been selected to show the tailings46response, thus providing a detailed description of the pre-failure and subsequent motion when acritical state-based constitutive model is considered. A similar approach was used to describe thematerial behavior at large deformations by Soga et al. (2016) and Cuomo et al. (2019).To facilitate the description of the deformation process points A, B, C, D and E have beenincluded in Figures 5.6 and 5.7. These points represent relevant moments in time to describe thetailings response. Point A is located at the beginning of the simulation at t = 0 s. Point B representsthe onset of flow liquefaction. Point C represents the control point reaching the critical state. Atpoint D the control point stop moving. Finally, point E depicts the final condition at t = 16 s.For the control point P2, Figure 5.6a shows that from A to B there is a reduction of the meaneffective stress due to the pore water pressure applied on the base of the model. This reductionresults in a slight movement of the control point P2 and no change in the void ratio. At B, themean effective stress starts to reduce dramatically (due to the increment of pore pressure) and thedisplacements suddenly increase. For these reasons, point B can be considered as the onset offlow liquefaction. After B, the control point reaches the critical state at C. From B to C, a flowliquefaction behavior as described in Section 3.2 is observed. This is supported by the stress pathshowed in the (p− q) diagram and the sudden increment of displacement. Furthermore, Figure5.5 shows that the reduction of mean effective stress is not only at the control point P2 but also inmany parts of the lower tailings. At C, the material reaches the critical state. Figure 5.6a shows thatthe void ratio (e) and mean pressure (p) of the tailings material match the critical void ratio (ec)and the critical mean pressure (pc), respectively. Furthermore, Figure 5.6b shows that the tailingsstress ratio (η) match the M value computed at the critical state (MCS). It should be noted that theevolution of the void ratio (e), mean pressure (p) and stress ratio (η) oscillate around the criticalvalues. It is also necessary to specify that the values ec, pc and MCS were computed considering theequations described in Section 3.4.2 for the SANISAND constitutive model. From state C to D, thematerial remains at the critical state, thus experiencing large deformations. Finally, at state D thematerial point is about to stop and leave the critical state. Afterwards, the pore pressure dissipatesand the mean pressure starts to increase until reaches state E.In reference to control point P4, Figure 5.7a shows that from A to B the material is beingdensified under gravity loading, thus increasing the mean effective stress and reducing the initialvoid ratio. There is no displacement from A to B. Similar to previous case, B represents the onsetof flow liquefaction. After B, the mean pressure drops sharply reaching the critical state at C.From B to C, the control point P4 experimented flow liquefaction. From C to D the material is atthe critical state with the parameters e, p and η oscillating around the critical values. This stageis associated with large deformations. At D, the pore water pressure is dissipated and the meanpressure is increased until reaching E. From D to E, there is no displacement.47      W V      W V      W V      W V      W VFigure 5.3: Variation of tailings dam deformed shape for case 148q  00 1.00.5      W V      W V      W V      W V      W VFigure 5.4: Evolution of deviatoric strains for case 149p N3D 01100      W V      W V      W V      W V      W VFigure 5.5: Evolution of mean effective stress for case 150        7LPH V  9RLGUDWLR  HHF        7LPH V SSSF        7LPH V 'LVSODFHPHQW P 3BACDEABCDEA BCD E(a)        7LPH V 0  0&6/    p>N3D@q>N3D@p>N3D@ e>@&6/BAC DEABCDEA,BCDE(b)Figure 5.6: Evolution of soil parameters at control point P2 for case 1 (a) Variation of voidratio, p and displacements (b) Variation of η , and stress development in the p−qand e− p plot during the deformation process51        7LPH V  9RLGUDWLR  HHF        7LPH V SSSF        7LPH V 'LVSODFHPHQW P 3A B CD EABCDEDAEB C(a)        7LPH V 0  0&6/    p>N3D@q>N3D@p>N3D@ e>@&6/ABC DEABECDABECD(b)Figure 5.7: Evolution of soil parameters at control point P4 for case 1 (a) Variation of voidratio, p and displacements (b) Variation of η , and stress development in the p−qand e− p plot during the deformation process525.3 Case 25.3.1 KinematicsThe computed displacements and velocities of control points P1, P2, P3 and P4 are shown in Figure5.8. These results show that the control points P1, P2 and P3 started to move at t = 2 s while P4started its movement at t = 4 s. Afterwards, the control points increase their velocity reaching thepeak velocity at t = 4 s (P1 and P2) and t = 7 s (P3 and P4). Control points P3 and P4 present higherdisplacements than control points P1 and P2. After reaching the peak value, the control pointsreduce their velocity. At t = 16 s, the velocity of the control points is almost zero.        7LPH V 'LVSODFHPHQW P 3333        7LPH V 9HORFLW\ PV 3333Figure 5.8: Displacement and velocities of control points P1, P2, P3 and P4 for case 25.3.2 Deformation processFigure 5.9 shows the tailings dam deformed shape at different times when trigger 2 is applied.It should be noted that displacements presented in Figure 5.9 are higher than the displacementsobserved in previous case. Figure 5.10 illustrates the deviatoric strains. At t = 2 s, there is aconcentration of deviatoric shear strains in the foundation due to trigger application. Also, the shearsurface in the lower tailings started to develop. At t = 4, the development of the shear surfacehas been completed and the system started to move. Figure 5.11 shows the evolution of the meaneffective stress. It should be noted that at t = 2 the mean pressures in the lower tailings have alreadybeen reduced. This reduction of mean pressures propagates downstream as the material flows.535.3.3 Time evolution of material stateFigure 5.12 shows the evolution of material state for control point P2 when trigger 2 is applied.From A to B there is no displacement and the mean pressure increases. Point B can be consideredas the onset of flow liquefaction. After B, the material reaches the critical state at C. From B to C,the material experimented flow liquefaction . From C to D the material remains at the critical state.It should be noted that in this simulation, the oscillations around the critical values are higher thancase 1. From C to D large deformations are observed. At D, the mean pressure increases due topore water pressure dissipation until reaching point E.5.4 SummaryThis Chapter presented the results of the problem described in Chapter 4. The explanation of thekinematics, the deformation process and the evolution of relevant parameters provided a detaileddescription of the tailings dam failure. These results have shown the capabilities of an advancedconstitutive model such as SANISAND along with the MPM platform to simulate the designedproblem. The simulations captured the onset of failure, the subsequent propagation of the tailingsmaterial and its final configuration.54      W V      W V      W V      W V      W VFigure 5.9: Variation of tailings dam deformed shape for case 255q  00 1.00.5      W V      W V      W V      W V      W VFigure 5.10: Evolution of deviatoric strains for case 256p N3D 01100      W V      W V      W V      W V      W VFigure 5.11: Evolution of mean effective stress for case 257        7LPH V   9RLGUDWLR  HHF        7LPH V SSSF        7LPH V 'LVSODFHPHQW P 3A B CDEBCADECBDAE(a),M()MCSLCSLABC DEABCEDABCDE        7LPH V      p>N3D@q>N3D@p>N3D@  e>@(b)Figure 5.12: Evolution of soil parameters at control point P2 for case 2 (a) Variation of voidratio, p and displacements (b) Variation of η , and stress development in the p−qand e− p plot during the deformation process58Chapter 6Assessment of flow liquefaction6.1 IntroductionIn the current engineering practice is crucial to understand the potential consequences of a tailingsdam failure due to flow liquefaction. In this context, the numerical simulations presented in Chapter5 represent an interesting tool to address this issue since it allows to describe the failure of tailingsdams using the SANISAND model and the MPM framework. This Chapter assessed the perfor-mance of the tailings dam described in Chapter 4 against flow liquefaction. Two approaches wereconsidered. First, the role of the initial void ratio of the lower tailings is analyzed. Furthermore,different trigger intensities are considered to evaluate the response. In both approaches, the triggersdefined in Section 4.1.5 were considered. The response of the system was quantified numericallyby the runout distance. This parameter was defined as the distance between the toe of the initialbuttress and the toe of the buttress after failure once the equilibrium has been re-established (seeFigure 6.1). The results are presented for each approach. Comparisons between two different val-ues of initial void ratio and two different trigger intensities are described. Moreover, time series areused to describe the influence of the initial void ratio and trigger intensity at the selected controlpoint (see Figure 6.1). A total of 58 numerical simulations were performed. The results were com-piled in a unique plot to assess the performance against flow liquefaction. Finally, it is importantto mention that the results presented in this Chapter may have practical application to evaluate theconsequences of tailings dam failures.593       5XQRXWFigure 6.1: Definition of the runout distance and selected control (moving) point for vulnera-bility assessment6.2 Effect of initial void ratioThe model described in Chapter 4 was evaluated for varying initial void ratios of the lower tailings.These analyses were performed for each trigger defined in Section 4.1.5. Appendix A shows contourplots of the initial distribution of the state parameter when different values of initial void ratio areconsidered. This Appendix also shows the element test response of the tailings material for selectedvalues of state parameter and mean pressures.6.2.1 Analyses with trigger 1The influence of the initial void ratio (e0) when trigger 1 is considered is demonstrated in Figure6.2 for e0 = 0.82 and 0.86. In the analysis, the deviatoric strain (εq) has been chosen to describethe failure development. Figure 6.3 offers a complementary picture of the deformation process byillustrating the deformed shape at different times. The evolution of displacements together with thevariation of the state parameter (ψ) for the control point P described in Figure 6.1 are depicted inFigure 6.4. This plot also shows the (e− p) diagram for different values of void ratios.Notably, Figure 6.2 shows that the initial failure occurs at the lower tailings and propagatestowards the right side of the model and inward to the tailings facility along an approximate circularsurface. This surface represents the boundary between the mobilized material and the temporal staticzone. Figure 6.3 indicates that during t = 0-4 s, there is almost no deformation for e0 = 0.82, whilethere is an important settlement in the raising lift for e0 = 0.86. During t = 4-12 s, the deformationprocess proceeds and more displacements are observed when higher void ratios are considered. At60t = 12 s, the geometry is almost flat for e0 = 0.86 and the mobilized tailings affected the raising liftslocated downstream. It is clear that the higher the initial void ratio, the longer the runout distance.Figure 6.4 shows that the displacements of control point P started earlier when the initial voidratio is higher. Furthermore, it is clear that the large displacements observed are associated withthe fact that control point P is approaching the critical state i.e., when the state parameter (ψ) iszero. The higher the initial void ratio, the faster the material approaches the critical state. It shouldbe noted that lower values of mean pressures (p′) are obtained when the material is approachingthe critical state for higher initial void ratios. This implies lower values of shear strength, thusexplaining the sudden movement and larger displacements when higher void ratios are used in theanalysis. 6.4 shows that at the end of the simulation, when the control point P becomes steady again,the mean pressure is gradually recovered. Finally, the results presented in Figure 6.4 show strongfluctuations in the stresses when higher values of initial void ratio are used. This is associated withthe significant shearing that is expected for more contractive materials, the MPM algorithm and theadvanced model used in simulations.6.2.2 Analyses with trigger 2This Section shows the influence of the initial void ratio when trigger 2 is used in the analysis.Figures 6.5 and 6.6 describe the deformation process. It should be noted that there is a significantdifference in the response of the system when the initial void ratio of the lower tailings is e0 =0.86 and 0.87. Figure 6.5 demonstrates that trigger 2 is not sufficient to cause displacements in thetailings structure when e0 = 0.86 is assumed. In this scenario, there is strain accumulation in thefoundation due to the loose of shear strength, however, a failure mechanism was not observed. Onthe contrary, when the initial void ratio e0 is 0.87 the shear surface in the lower tailings has alreadybeen developed at t = 2 s. These results illustrates how the initial void ratio of the lower tailingsplays a decisive role in the stability of the tailings dam analyzed.6.3 Effect of trigger intensity6.3.1 Analyses with trigger 1To evaluate the effect of trigger intensity different values of maximum pore water pressure havebeen considered. The model presented in Chapter 5 assumed a maximum pore water pressure of200 kPa applied at the base of the model (see Figure 4.4). In this Section eight scenarios have beenconsidered which correspond to maximum pore water pressures of 150, 170, 180, 190, 195, 200,250 and 300 kPa.61Figure 6.7 compares the deformation process for the same initial void ratio e0 = 0.80 whentrigger intensities of 200 and 250 kPa are considered. Both cases show that the failure starts at thelower tailings and further propagates horizontally towards the downstream direction. The elapsedtime for the failure initialization is different for each case. When pw = 200 kPa the slip surface isalmost fully developed at t = 7 s, while for the pw = 250 kPa case the failure surface was developedat t = 0-4 s. Figure 6.8 shows that higher displacements are obtained when higher pore waterpressures are applied. During t = 0-7 s there is no deformation when pw = 200 kPa is used, whilethere is an important settlement for the same period of time when the trigger intensity is pw = 250kPa. These differences illustrate the effect of trigger intensity on the tailings dam performance.Figure 6.9 illustrates the displacement of control point P throughout the time for different triggerintensities. It is clear that the higher the intensity, the larger the displacement. It should be notedthat the movement of control point P starts as soon as there is a reduction in the state parameter. Thisreduction is associated with a significant drop of the mean effective stress. Lower mean pressuresare reached when higher trigger intensities are considered. The results presented in Figure 6.9 showstrong oscillations in the mean effective stresses when higher trigger intensities are considered.6.3.2 Analyses with trigger 2To evaluate the effect of trigger intensity two additional thickness of foundation were considered.The model described in Chapter 4 assumed a thickness foundation of 1.0 m. This Section also con-siders thicknesses of 2 and 3 m for the foundation. Similar to the simulation presented in Chapter 5,the motion was triggered by a sudden loss of shear strength in the foundation for all cases analyzed.Figure 6.10 compares the evolution of deviatoric strains for two different foundation thickness: 1and 2 m. It should be noted that for the case with 1-m thickness foundation, there is only strainaccumulation in the foundation and does not propagate to the tailings. However, for the case with2-m thickness foundation, there is a shear surface in the lower tailings that propagates downstreamas the tailings flows. Figure 6.11 compares the deformed shape for each trigger intensity. It can beobserved that, under the same initial void ratio e0 = 0.86, the system with 2-m thickness foundationpresented large deformations while the tailings dam with 1-m foundation was not affected by theapplied trigger. These results illustrates how the thickness foundation plays an important role in thestability of the tailings dam analyzed as long as its shear strength is suddenly reduced.6.4 SummaryThe results presented in this Chapter are summarized considering the runout distance. Figures 6.12and 6.13 show the correlation between runout and the initial void ratio of the lower tailings for each62trigger analyzed.Figures 6.12 shows that the higher the intensity of the triggering mechanism, the larger therunout distance for the same void ratio. Furthermore, the higher the void ratio, the larger the runoutfor the same trigger intensity. It should be noted that regardless of the trigger considered in theanalysis, the largest runout values computed during simulations correspond to initial void ratios of0.86 and 0.87. It should be mentioned that values of void ratio higher than e0 = 0.87 are not possibleto consider in the analysis since for such values the geometry described in Chapter 4 is not stableunder gravity loading and a trigger mechanism is not necessary to cause instability. Figures 6.13illustrates that the thickness foundation plays an important role on the tailings dam performance.It should be noted that when e0 = 0.85 the only simulation that presents a considerable runout isthe one that assumes 3-m thickness foundation. It is important to mention that when high values ofinitial void ratio are considered, large displacements are expected and the boundary of the modelmay affect the runout calculation.Finally, Figures 6.12 and 6.13 allow to define a zone where the system will experience largedeformations regardless of the trigger application. This zone can be considered between e0 = 0.86and 0.87. Furthermore, the zone between e0 = 0.80 and 0.85 represent the area where the applicationof trigger plays a crucial role in the stability. These results indicate that the initial void ratio shouldbe acknowledged as an important factor when it comes to analyzing the stability of tailings damsespecially when flow liquefaction behavior is expected.63q  00 1.00.5      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.2: Comparison of the evolution of deviatoric strains with different initial void ratios:(a) e0 = 0.82 (b) e0 = 0.86 when trigger 1 is considered64      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.3: Comparison of the evolution of the tailings dam deformed shape for different ini-tial void ratios: (a) e0 = 0.82 (b) e0 = 0.86 when trigger 1 is considered65        7LPH  VHF 7RWDO'LVSODFHPHQW  H H H H         7LPH  VHF 6WDWH3DUDPHWHU  H H H H &6/        7LPH  VHF 0HDQ(IIHFWLYH6WUHVV N3D H H H H Figure 6.4: Comparisons of the evolution of displacement, state parameter and mean pressurefor different initial void ratios at control point P when trigger 1 is considered66q  00 1.00.5      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.5: Comparison of the evolution of deviatoric strains with different initial void ratios:(a) e0 = 0.86 (b) e0 = 0.87 when trigger 2 is considered67      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.6: Comparison of the evolution of the tailings dam deformed shape for different ini-tial void ratios: (a) e0 = 0.86 (b) e0 = 0.87 when trigger 2 is considered68q  00 1.00.5      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.7: Comparison of the evolution of deviatoric strains with different trigger intensities:(a) pw = 200 kPa (b) pw = 250 kPa when e0 = 0.80 is considered69      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.8: Comparison of the evolution of the tailings dam deformed shape with differenttrigger intensities: (a) pw = 200 kPa (b) pw = 250 kPa when e0 = 0.80 is considered70        7LPH  VHF 7RWDO'LVSODFHPHQW  SZ  N3DSZ  N3DSZ  N3D        7LPH  VHF 6WDWH3DUDPHWHU  SZ  N3DSZ  N3DSZ  N3D&6/        7LPH  VHF 0HDQ(IIHFWLYH6WUHVV N3D SZ  N3DSZ  N3DSZ  N3DFigure 6.9: Comparisons of the evolution of displacement, state parameter and mean pressurefor different trigger intensities at control point P when e0 = 0.80 is considered71q  00 1.00.5      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.10: Comparison of the evolution of deviatoric strains with different foundation thick-ness: (a) h = 1.0 m (b) h = 2.0 m when e0 = 0.86 is considered72      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V      W V(a)      W V(b)Figure 6.11: Comparison of the evolution of the tailings dam deformed shape with differentfoundation thickness: (a) h = 1.0 m (b) h = 2.0 m when e0 = 0.86 is considered73       ,QLWLDO 9RLG 5DWLR   5XQRXW P SZ   N3DSZ   N3DSZ   N3DSZ    N3DSZ    N3DSZ   N3DSZ   N3DSZ   N3DFigure 6.12: Relationship between runout and initial void ratio for different values of porewater pressure74       ,QLWLDO 9RLG 5DWLR   5XQRXW P K  PK  PK  PFigure 6.13: Relationship between runout and initial void ratio for different thicknesses in thefoundation75Chapter 7Summary, Conclusions and FutureResearch7.1 SummaryRecent tailings dam failures have shown the need for performing large deformation simulationsto estimate the expected damage that may help risk assessment. To evaluate tailings dams per-formance, the current engineering practice uses either limit equilibrium or finite element methods.However, these approaches only provide information to the proximity of the failure and cannot gobeyond. Another aspect that should bear in mind is that the primary factor for causing the insta-bility in tailings dams is the flow liquefaction phenomena. Therefore, a proper simulation of thesecomplex evolution process should consider two components: a numerical platform that is able to re-produce large deformations and an advance constitutive model that can capture the flow liquefactionobserved in tailings.The current study focused on conducting numerical simulations that considered the triggeringand motion in a unified framework, thus capturing the failure initiation and subsequent runout. Theanalyses utilized the MPM code, Anura3D, in conjunction with an anisotropic bounding surfaceplasticity model, SANISAND. The reliability of the model implementation in the MPM code wasfirst assessed through verification process. The results obtained from element tests simulationsperformed using the program Incremental Driver were compared with numerical simulations ofundrained triaxial tests conducted in the software Anura3D. The comparison showed that the resultsusing the SANISAND implementation in Anura3D are in good agreement with the element testsresults obtained from Incremental Driver. Moreover, the verification results illustrated the predictivecapabilities of SANISAND to capture the flow liquefaction behavior observed in tailings.76Chapter 2 presented a overview of the numerical methods in geomechanics to deal with largedeformations. In this Chapter the MPM method was described and reasons for its use were pro-vided. Chapter 3 described the flow liquefaction phenomenon and the constitutive models used inMPM. The features that one is looking for in a model to simulate the failure of tailings dam werehighlighted and the SANISAND model was chosen for the simulations. Chapter 4 presented thehypothetical tailings dam assumed in the analysis. This geometry was motivated by a historicaltailings dam where flow liquefaction took place. The model parameters used in the simulationswere calibrated for a particular tailings. Chapter 5 presented the simulation results for the two casesdescribed in Chapter 4. The results showed that numerical simulations were able to capture the de-formation process of tailings dam failures. The results also suggested that the combination of MPMand SANISAND allow to simulate relevant aspects of the response such as the flow liquefactionphenomenon and the tailings behavior at the critical state. Chapter 6 described the assessment offlow liquefaction for the hypothetical tailings dam described in Chapter 4. Two approaches wereconsidered for the evaluation. For each case, the triggers defined in Chapter 4 were considered. Theresults suggested that the higher the trigger intensity, the larger the runout for the same initial voidratio. Furthermore, it should be noted that the higher the initial void ratio, the larger the runout forthe same trigger intensity. Finally, the results showed that there is a value of initial void ratio whereregardless of the triggering event the system will be experience large deformations. This shouldbe acknowledged as an relevant factor when it comes to analyzing tailings dam subjected to flowliquefaction.7.2 ConclusionsTo the author’s knowledge, this thesis is the first application of the SANISAND model in the studyof large deformation problems with a MPM code. This represents the main contribution of thisstudy. The numerical simulations illustrated the capabilities of MPM along with the SANISANDmodel to simulate the failure of tailings dams. The simulation results provided valuable informationabout the onset of flow liquefaction and subsequent large deformations. This is thought to be one ofthe key aspects achieved in this thesis. The results from this type of simulations may help engineersto evaluate the consequences of tailings dam disasters. Therefore, the numerical tool used in thisdocument may have practical application.7.3 Recommendations for future researchThe oscillations in stresses obtained in this study are mostly related to the performance of theSANISAND model in the MPM framework especially when higher initial void ratios and higher77trigger intensities are considered. The oscillations can be attributed to the high deviatoric strainsinvolved during the process or to the MPM algorithm. To reduce the stress oscillations more inves-tigation is desired. In this context, SANISAND could be used in another MPM platforms such thegeneralized interpolation material point (GIMP) method or the dual domain material point (DDMP)method. Such algorithms do not consider linear interpolation functions. Therefore, a better responseof the SANISAND model may be expected.This thesis can represent the first step to continue using the SANISAND model in the MPMframework. One application that might provide valuable information for practitioners and re-searchers is the simulation of the cone penetration tests (CPTu) in loose sands. This study canalso encourage to use the SANICLAY model in the MPM framework. Applications for this modelmay be the simulation of brittle clays where large deformations are expected.78BibliographyK. Abe, K. Soga, and S. Bandara. Material point method for coupled hydromechanical problems.Journal of Geotechnical and Geoenvironmental Engineering, 140(3):04013033, 2014.J. Ahrens, B. Geveci, and C. Law. Paraview: An end-user tool for large data visualization. Thevisualization handbook, 717(8), 2005.I. K. Al-Kafaji. Formulation of a dynamic material point method (MPM) for geomechanicalproblems. PhD Thesis, University of Stuttgart, 2013.E. E. Alonso. Triggering and motion of landslides. Géotechnique, 71(1):3–59, 2021.A. Alsardi, J. Copana, and A. Yerro. Modelling earthquake-triggered landslide runout with thematerial point method. Proceedings of the Institution of Civil Engineers-GeotechnicalEngineering, pages 1–14, 2021.S. Andersen and L. Andersen. Modelling of landslides with the material-point method.Computational Geosciences, 14(1):137–147, 2010.J. E. Andrade, A. M. Ramos, and A. Lizcano. Criterion for flow liquefaction instability. ActaGeotechnica, 8(5):525–535, 2013.Anura3D. Anura3D - MPM Research Community. Available from http://www.mpm-dredge.eu/,2021.F. Armero and E. Love. An arbitrary lagrangian–eulerian finite element method for finite strainplasticity. International journal for numerical methods in engineering, 57(4):471–508, 2003.D. Aubram, F. Rackwitz, P. Wriggers, and S. A. Savidis. An ale method for penetration into sandutilizing optimization-based mesh motion. Computers and Geotechnics, 65:241–249, 2015.S. Bandara and K. Soga. Coupling of soil deformation and pore fluid flow using material pointmethod. Computers and geotechnics, 63:199–214, 2015.S. S. Bandara. Material point method to simulate large deformation problems in fluid-saturatedgranular medium. PhD thesis, University of Cambridge Cambridge, UK, 2013.79S. G. Bardenhagen and E. M. Kober. The generalized interpolation material point method.Computer Modeling in Engineering and Sciences, 5(6):477–496, 2004.A. R. Barrero, M. Taiebat, and Y. F. Dafalias. Modeling cyclic shearing of sands in thesemifluidized state. International Journal for Numerical and Analytical Methods inGeomechanics, 44(3):371–388, 2020.K. Been and M. Jefferies. Discussion: a state parameter for sands. Géotechnique, 36(1):123–132,1986.K. Been and M. G. Jefferies. A state parameter for sands. Géotechnique, 35(2):99–112, 1985.L. Beuth. Formulation and application of a quasi-static material point method. PhD Thesis,University of Stuttgart, 2012.L. Beuth, T. Benz, P. Vermeer, C. Coetzee, P. Bonnier, and P. Van Den Berg. Formulation andvalidation of a quasi-static material point method. In Proceedings of the 10th internationalsymposium on numerical methods in Geomechanics (NUMOG), Rhodes, Greece, pages189–195, 2007.L. Beuth, T. Benz, P. A. Vermeer, and Z. Wie˛ckowski. Large deformation analysis using aquasi-static material point method. 2008.H. H. Bui, K. Sako, and R. Fukagawa. Numerical simulation of soil–water interaction usingsmoothed particle hydrodynamics (sph) method. Journal of Terramechanics, 44(5):339–346,2007.H. H. Bui, R. Fukagawa, K. Sako, and J. C. Wells. Slope stability analysis and discontinuous slopefailure simulation by elasto-plastic smoothed particle hydrodynamics (sph). Géotechnique, 61(7):565–574, 2011.A. Carrera, M. Coop, and R. Lancellotta. Influence of grading on the mechanical behaviour ofstava tailings. Géotechnique, 61(11):935–946, 2011.A. Casagrande. Characteristics of cohesionless soils affecting the stability of slopes and earth fills.J. Boston Society of Civil Engineers, 23(1):13–32, 1936.A. Casagrande. Liquefaction and cyclic deformation of sands-a critical review. Harvard SoilMechanics Series, Harvard University, Cambridge, Massachusetts., (88), 1976.G. Castro. Liquefaction of sands. PhD thesis, Harvard University, Cambridge, MA, 1969.F. Ceccato. Study of large deformation geomechanical problems with the material point method.PhD Thesis, Universita Degli Studi Di Padova, 2015.F. Ceccato, L. Beuth, and P. Simonini. Analysis of piezocone penetration under different drainageconditions with the two-phase material point method. Journal of Geotechnical andGeoenvironmental Engineering, 142(12):04016066, 2016.80F. Ceccato, A. Yerro, V. Girardi, and P. Simonini. Two-phase dynamic mpm formulation forunsaturated soil. Computers and Geotechnics, 129:103876, 2021.R. J. Chandler and G. Tosatti. The stava tailings dams failure, Italy, July 1985. Proceedings of theinstitution of civil engineers-geotechnical engineering, 113(2):67–79, 1995.J.-C. Chern. Undrained response of saturated sands with emphasis on liquefaction and cyclicmobility. PhD thesis, University of British Columbia, 1985.J. Chu. An experimental examination of the critical state and other similar concepts for granularsoils. Canadian Geotechnical Journal, 32(6):1065–1075, 1995.J. Chu, S. Leroueil, and W. Leong. Unstable behaviour of sand and its implication for slopeinstability. Canadian Geotechnical Journal, 40(5):873–885, 2003.R. Corti, A. Diambra, D. M. Wood, D. E. Escribano, and D. F. Nash. Memory surface hardeningmodel for granular soils under repeated loading conditions. Journal of Engineering Mechanics,142(12):04016102, 2016.S. Cuomo, P. Ghasemi, M. Martinelli, and M. Calvello. Simulation of liquefaction andretrogressive slope failure in loose coarse-grained material. International Journal ofGeomechanics, 19(10):04019116, 2019.Y. F. Dafalias and M. T. Manzari. Simple plasticity sand model accounting for fabric changeeffects. Journal of Engineering mechanics, 130(6):622–634, 2004.M. Davies, B. Chin, and B. Dawson. Static liquefaction slump of mine tailings - A case history.pages 123–132. 51st Canadian geotechnical conference, 1998.Y. Di, J. Yang, and T. Sato. An operator-split ale model for large deformation analysis ofgeomaterials. International Journal for Numerical and Analytical Methods in Geomechanics, 31(12):1375–1399, 2007.J. Donea, A. Huerta, J.-P. Ponthot, and A. Rodríguez-Ferran. Arbitrary l agrangian–e ulerianmethods. Encyclopedia of computational mechanics, 2004.E. J. Fern and K. Soga. The role of constitutive models in mpm simulations of granular columncollapses. Acta Geotechnica, 11(3):659–678, 2016.J. Fern, A. Rohe, K. Soga, and E. Alonso. The material point method for geotechnicalengineering: a practical guide. CRC Press, 2019.A. Fourie, G. Blight, and G. Papageorgiou. Static liquefaction as a possible explanation for themerriespruit tailings dam failure. Canadian Geotechnical Journal, 38(4):707–719, 2001.A. Gens Solé. Hydraulic fills with special focus on liquefaction. In Proceedings of the XVIIECSMGE-2019: Geotechnical Engineering foundation of the future, pages 1–31. InternationalSociety for Soil Mechanics and Geotechnical Engineering (ISSMGE), 2019.81R. A. Gingold and J. J. Monaghan. Smoothed particle hydrodynamics: theory and application tonon-spherical stars. Monthly notices of the royal astronomical society, 181(3):375–389, 1977.J. Guilkey and J. Weiss. An implicit time integration strategy for use with the material pointmethod. In Proceedings from the First MIT Conference on Computational Fluid and SolidMechanics, volume 2, pages 216–219, 2001.K. Ishihara. Liquefaction and flow failure during earthquakes. Geotechnique, 43(3):351–451,1993.K. Ishihara, F. Tatsuoka, and S. Yasuda. Undrained deformation and liquefaction of sand undercyclic stresses. Soils and foundations, 15(1):29–44, 1975.I. Jassim, D. Stolle, and P. Vermeer. Two-phase dynamic analysis by material point method.International journal for numerical and analytical methods in geomechanics, 37(15):2502–2522, 2013.M. Jefferies and K. Been. Soil liquefaction: a critical state approach. CRC press, 2016.J. Konkol. Numerical solutions for large deformation problems in geotechnical engineering. PhDInterdisciplinary Journal, 2014:49–55, 2014.A. Lashkari. Prediction of flow liquefaction instability of clean and silty sands. Acta Geotechnica,11(5):987–1014, 2016.W. Leong, J. Chu, and C. Teh. Liquefaction and instability of a granular fill material. GeotechnicalTesting Journal, 23(2):178–192, 2000.S. Leroueil, J. Chu, and D. Wanatowski. Slope instability due to pore water pressure increase. InThe first Italian workshop on landslides, Naples, pages 81–90, 2009.S. Li and W. K. Liu. Meshfree and particle methods and their applications. Appl. Mech. Rev., 55(1):1–34, 2002.D. Liang, X. Zhao, and A. Rohe. Simulation of seepage and overtopping induced dam failures. In2nd International Conference on the Material Point Method for Modelling Soil-Water-StructureInteraction, pages 278–83, 2019.M. A. Llano-Serna, M. M. Farias, and D. M. Pedroso. An assessment of the material point methodfor modelling large scale run-out processes in landslides. Landslides, 13(5):1057–1066, 2016.L. B. Lucy. A numerical approach to the testing of the fission hypothesis. The astronomicaljournal, 82:1013–1024, 1977.K. Maeda and H. Sakai. Seepage failure and erosion of ground with air bubble dynamics. InGeoenvironmental engineering and geotechnics: progress in modeling and applications, pages261–266. 2010.82M. T. Manzari and Y. F. Dafalias. A critical state two-surface plasticity model for sands.Geotechnique, 47(2):255–272, 1997.M. Martinelli. Soil-water interaction with material point method. double-point formulation. Reporton EU-FP7 research project MPM-Dredge PIAP-GA-2012, 324522, 2016.M. Martinelli and A. Rohe. Modelling fluidisation and sedimentation using material point method.In 1st Pan-American Congress on Computational Mechanics, pages 1–12, 2015.M. Martinelli, A. Rohe, and K. Soga. Modeling dike failure using the material point method.Procedia Engineering, 175:341–348, 2017.D. Mašín. SANISAND Abaqus UMAT and Plaxis implementations.https://soilmodels.com/download/plaxis-umat-sanisand/, 2017. Accessed: 2021-07-26.C. M. Mast, P. Arduino, P. Mackenzie-Helnwein, and G. R. Miller. Simulating granular columncollapse using the material point method. Acta Geotechnica, 10(1):101–116, 2015.A. Melendo, A. Coll, M. Pasenau, E. Escolano, and A. Monros. GID – The personal pre and postprocessor. https://www.gidhome.com/, 2021. Accessed: 2021-10-13.M. Mieremet. Numerical stability for velocity-based 2-phase formulation for geotechnicaldynamic analysis. Reports of the Delft Institute of Applied Mathematics, 15-03, 2015.J. J. Monaghan. Smoothed particle hydrodynamics. Reports on progress in physics, 68(8):1703,2005.N. R. Morgenstern, S. G. Vick, and D. V. Zyl. Report on Mount Polley tailings facility breach.2015.N. R. Morgenstern, S. G. Vick, C. B. Viotti, and B. D. Watts. Report on the immediate causes ofthe failure of the Fundão dam. Fundão Tailings Dam Review Panel, 2016.N. R. Morgenstern, M. Jefferies, D. V. Zyl, and J. Wates. Report on NTSF Embankment Failure.Independent Technical Review Board, 2019.H. Mori. SPH method to simulate river levee failures. PhD thesis, University of Cambridge, 2008.J. A. Nairn. Material point method calculations with explicit cracks. Computer Modeling inEngineering and Sciences, 4(6):649–664, 2003.A. Niemunis. Incremental driver, user’s manual. University of Karlsruhe KIT, Germany, 2008.A. Niemunis. Download Package of Incremental Driver.https://soilmodels.com/download/incrementaldriver-zip/, 2017. Accessed: 2021-07-26.W. F. Noh. Cel: A time-dependent, two-space-dimensional, coupled eulerian-lagrange code.Technical report, Lawrence Radiation Lab., Univ. of California, Livermore, 1963.83S. M. Olson, T. D. Stark, W. H. Walton, and G. Castro. 1907 static liquefaction flow failure of thenorth dike of wachusett dam. Journal of Geotechnical and Geoenvironmental engineering, 126(12):1184–1193, 2000.W. Powrie. Soil mechanics: concepts and applications. CRC Press, 2018.G. Qiu and J. Grabe. Explicit modeling of cone and strip footing penetration under drained andundrained conditions using a visco-hypoplastic model. geotechnik, 34(3):205–217, 2011.G. Qiu, S. Henke, and J. Grabe. Application of a coupled eulerian–lagrangian approach ongeomechanical problems involving large deformations. Computers and Geotechnics, 38(1):30–39, 2011.J. Ramirez, A. R. Barrero, L. Chen, S. Dashti, A. Ghofrani, M. Taiebat, and P. Arduino. Siteresponse in a layered liquefiable deposit: evaluation of different numerical tools andmethodologies with centrifuge experimental results. Journal of Geotechnical andGeoenvironmental Engineering, 144(10):04018073, 2018.P. K. Robertson, L. de Melo, D. J. Williams, and G. W. Wilson. Report of the Expert Panel on theTechnical Causes of the Failure of Feijao dam I. 2019.K. H. Roscoe, A. Schofield, and a. P. Wroth. On the yielding of soils. Geotechnique, 8(1):22–53,1958.J. C. Santamarina, L. A. Torres-Cruz, and R. C. Bachus. Why coal ash and tailings dam disastersoccur. Science, 364(6440):526–528, 2019.F. Schnaid, J. Bedin, A. Viana da Fonseca, and L. de Moura Costa Filho. Stiffness and strengthgoverning the static liquefaction of tailings. Journal of Geotechnical and GeoenvironmentalEngineering, 139(12):2136–2144, 2013.A. N. Schofield and P. Wroth. Critical state soil mechanics, volume 310. McGraw-hill London,1968.K. Soga, E. Alonso, A. Yerro, K. Kumar, and S. Bandara. Trends in large-deformation analysis oflandslide mass movements with particular emphasis on the material point method.Géotechnique, 66(3):248–273, 2016.D. Sulsky and H. L. Schreyer. Axisymmetric form of the material point method with applicationsto upsetting and taylor impact problems. Computer Methods in Applied Mechanics andEngineering, 139(1-4):409–429, 1996.D. Sulsky, Z. Chen, and H. L. Schreyer. A particle method for history-dependent materials.Computer methods in applied mechanics and engineering, 118(1-2):179–196, 1994.D. Sulsky, S.-J. Zhou, and H. L. Schreyer. Application of a particle-in-cell method to solidmechanics. Computer physics communications, 87(1-2):236–252, 1995.84M. Taiebat and Y. F. Dafalias. Sanisand: Simple anisotropic sand plasticity model. InternationalJournal for Numerical and Analytical Methods in Geomechanics, 32(8):915–948, 2008.M. Taiebat, B. Jeremic´, Y. F. Dafalias, A. M. Kaynia, and Z. Cheng. Propagation of seismic wavesthrough liquefied soils. Soil Dynamics and Earthquake Engineering, 30(4):236–257, 2010.D. W. Taylor. Fundamentals of soil mechanics, volume 66. LWW, 1948.Y. Vaid and J. Chern. Cyclic and monotonic undrained response of saturated sands. In Advances inthe art of testing soils under cyclic conditions, pages 120–147. ASCE, 1985.T. Wichtmann, W. Fuentes, and T. Triantafyllidis. Inspection of three sophisticated constitutivemodels based on monotonic and cyclic tests on fine sand: Hypoplasticity vs. sanisand vs. isa.Soil Dynamics and Earthquake Engineering, 124:172–183, 2019.J. A. Yamamuro and P. V. Lade. Steady-state concepts and static liquefaction of silty sands.Journal of geotechnical and geoenvironmental engineering, 124(9):868–877, 1998.M. Yang, G. Seidalinov, and M. Taiebat. Multidirectional cyclic shearing of clays and sands:Evaluation of two bounding surface plasticity models. Soil Dynamics and EarthquakeEngineering, 124:230–258, 2019.M. Yang, M. Taiebat, and Y. F. Dafalias. Sanisand-msf: a sand plasticity model with memorysurface and semifluidised state. Géotechnique, pages 1–20, 2020.A. Yerro. MPM modelling of landslides in brittle and unsaturated soils . PhD Thesis, UniversitatPolitecnica de Catalunya, 2015.A. Yerro, E. Alonso, and N. Pinyol. The material point method for unsaturated soils.Géotechnique, 65(3):201–217, 2015.A. Yerro, E. Alonso, and N. Pinyol. Run-out of landslides in brittle soils. Computers andGeotechnics, 80:427–439, 2016.A. Yerro, K. Soga, and J. Bray. Runout evaluation of oso landslide with the material point method.Canadian Geotechnical Journal, 56(9):1304–1317, 2019.F. Zabala and E. Alonso. Progressive failure of aznalcóllar dam using the material point method.Géotechnique, 61(9):795–808, 2011.F. Zabala García. Modelación de problemas geotécnicos hidromecánicos utilizando el método delpunto material. PhD Thesis, Universitat Politecnica de Catalunya, 2010.D. Z. Zhang, X. Ma, and P. T. Giguere. Material point method enhanced by modified gradient ofshape function. Journal of Computational Physics, 230(16):6379–6398, 2011.85Appendix AUse of SANISAND for different values ofstate parameterA.1 IntroductionThis Appendix describes the element test response for different initial state parameters when theSANISAND constitutive model is used in the simulations. To ensure the use of the same initialstate parameter, different combinations of void ratios and mean pressures are used. The softwareIncremental Driver is considered in the simulations (Niemunis [2008]).A.2 Selection of initial state parameterThe range of initial state parameters depends on the initial void ratio used in the analyses. FigureA.1 shows the state parameter initial distribution for different values of initial void ratio assumed inthe lower tailings. It should be noted that state parameter varies between -0.10 and +0.20. There-fore, state parameters equal to -0.05, 0.05, 0.10 and 0.15 are considered in the Incremental Driversimulations.A.3 Element tests simulationsFigures A.2 and A.3 illustrate the element test response for different initial state parameters whenundrained triaxial conditions are considered. It should be noted that the SANISAND model is ableto capture the flow liquefaction response observed in tailings. Furthermore, the higher the stateparameter, the lower deviatoric stress values are reached for the same level of mean pressure.86       (a)       (b)       (c)       (d)       (e)Figure A.1: Distribution of initial state parameter for different values of initial void ratio: (a)e0 = 0.78 (b) e0 = 0.80 (c) e0 = 0.82 (d) e0 = 0.84 (e) e0 = 0.8787      q (%)q (kPa)H  H H  H       p (kPa)q (kPa)H  H H  H p (kPa) Void Ratio (-)&6/ ,QLWLDO)LQDO(a)      q (%)q (kPa)H  H H  H       p (kPa)q (kPa)H  H H  H p (kPa) Void Ratio (-)&6/ ,QLWLDO)LQDO(b)Figure A.2: Element test response for: (a) ψ = -0.05 (b) ψ = +0.0588      q (%)q (kPa)H H  H H        p (kPa)q (kPa)H H  H H  p (kPa) Void Ratio (-)&6/ ,QLWLDO)LQDO(a)      q (%)q (kPa)H  H H  H       p (kPa)q (kPa)H  H H  H p (kPa) Void Ratio (-)&6/ ,QLWLDO)LQDO(b)Figure A.3: Element test response for: (a) ψ = +0.10 (b) ψ = +0.1589A.4 Summary of resultsFigure A.4 summarizes the results considering the undrained shear strength ratio ( q2p0 ) obtainedfrom the simulations. The results suggest that the undrained strength ratio is only independent ofthe initial pressure when the initial condition of the soil represents the critical state (ψ = 0). Foranother initial state parameters, the strength ratio is pressure-dependent. It also should be notedthat lower strength ratio values are obtained when higher state parameters and lower pressures areconsidered.      ,QLWLDO 6WDWH 3DUDPHWHU    q2p0 ( )S   N3DS   N3DS   N3DS   N3DFigure A.4: Relationship between shear strength ratio and the initial state parameter of thelower tailings90
Thesis/Dissertation
eng
Civil Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Taiebat, Mahdi
Modeling of flow liquefaction and large deformations in tailings dams using material point method
Text