Centrifuge and Numerical Modeling of Liquefied Flow and Nonliquefied Slide Failures of Tailings Dams

The gold tailings used in this study was obtained from a gold mine situated east of Johannesburg, South Africa. The material was collected from the daywall of the tailings dam. Based on the particle size distribution (PSD), shown in Fig. 2, the gold tailings used in this study is classified as a silty sand. The tailings has a fines content ($D<63\mathrm{\mu m}$) of 26.5%, an average grain size $(D_{50})$ of 120 $\mu \mathrm{m}$, and a maximum particle size of 500 μm. The saturated coefficient of permeability, $k_s=1.4\times {10}^{-6}\mathrm{m}/\mathrm{s}$, was estimated using Hazen’s (1930) equation:where $D_{10}$ is the Hazen’s effective grain size (mm). The gold tailings has a specific gravity, $G_s$, of 2.65. The index properties of the gold tailings are summarised in Table 2.

The model slopes were constructed in a container with interior dimensions of $800\times 160\times 360\mathrm{mm}$ (Fig. 1), intending to model plane-strain conditions. The container was constructed from 50 mm thick aluminum panels and was fitted with a 50 mm thick glass panel to allow observations.

Twin sand filters comprised of a coarse sand filter, followed by a fine sand filter, were placed immediately upstream of the model slopes. This was done to prevent potential internal erosion in the dam slope (Ng et al. 2022) during the fluid level rise. The filters allowed for the viscous fluid to seep slowly and evenly before reaching the tailings slope. The permeability of the coarse sand filter was on the order of ${10}^{-4}\mathrm{m}/\mathrm{s}$ and that of the fine sand filter ${10}^{-5}\mathrm{m}/\mathrm{s}$.

The model slopes were prepared using moist tamping. This method was chosen because capillary effects between soil particles allow for the preparation of loose specimens (Castro 1969; Sladen et al. 1985; Kramer and Seed 1988; Konrad 1990; Ishihara 1993; Cai 2001). The moist tamping method also allows a contractive fabric with strain-softening tendencies, which is more susceptible to liquefaction, to be created (Chang et al. 2011). Both model slopes were formed by moist tamping the gold tailings in six 50 mm horizontal layers at an initial moisture content of 6.5%. The density was controlled by placing a predetermined amount of soil in each layer and compacting the material by hand in 50 mm layers. A temporary steel plate shutter was placed at the upstream end of the model at the boundary between the tailings and sand filters. Once a layer of tailings was compacted, the twin sand filters were placed at the upstream end of the models, also with a thickness of 50 mm each. Once six layers had been compacted, i.e., a thickness of 300 mm, the slope was cut by hand to an angle of 35° to the horizontal. The model slope in Test GT-1 was prepared at an initial void ratio, $e_0$, of 0.950 and in Test GT-2 of 0.845.

A groundwater system similar to Ng et al. (2022) was adopted and modified for this study. The system consisted of a fluid storage tank, a solenoid valve, a pressure transducer, a perforated steel plate covered by geotextile, an aluminum partition at the toe of the model slope, and an outlet. A perforated plate was installed upstream of the model slope to create a reservoir that allowed for the control of the hydraulic head upstream of the slope (Fig. 1). The perforated plate was covered with geotextile to prevent sand from filling the reservoir. A fluid storage tank, holding the viscous fluid, was mounted near the centrifuge rotation axis and was connected to the reservoir upstream of the slope via 8 mm plastic tubing. The fluid level in the reservoir was controlled using a solenoid valve. A pressure transducer installed at the inlet of the model container monitored the hydraulic head in the reservoir. The fluid level at the downstream end of the model slope was controlled using an aluminum partition plate with a height equal to that of the toe of the slope. As the fluid seeped through the model slope, the fluid level rose until it overflowed the partition, thus maintaining a downstream level constantly level with the toe, preventing inundation of the lower parts of the slope.

The intention with the centrifuge models was to investigate flowslide behavior during slope instability. It was therefore desirable to create instability-type events in the centrifuge. In order to successfully model the flowslide behavior of the slope and the subsequent dissipation of excess pore pressures, it was necessary to increase the viscosity of the pore fluid linearly with the scale factor (Askarinejad et al. 2014; Take et al. 2004; Zhang and Askarinejad 2021). The viscosity therefore had to be 60 times that of water to appropriately scale the dissipation of excess pore pressures. For both centrifuge tests, a mixture of glycerine and water was used. To determine the concentration glycerine required, a series of glycerine-water mixtures were prepared with 15%, 30%, 45%, 75%, 85%, and 92% glycerine. The viscosity of each mixture was measured at 25°C using a rheometer. In addition, the viscosities of the mixtures were measured at varying shear strain rates (${\mathrm{s}}^{-1}$). Fig. 3 shows the viscosity of the mixtures at 25°C and shear strain rates of 10 and $251{\mathrm{s}}^{-1}$. From the viscosity relationship shown in Fig. 3, a glycerine-water mixture with 86% glycerine and 14% water was chosen to obtain a fluid with a viscosity 60 times higher than water. The temperature within the centrifuge facility was set at 25°C using a dedicated ventilation system, which allowed for the viscosity of the fluid to remain unchanged. More details about the selection of a viscous fluid can be found in Crous et al. (2022).

To monitor the response of the pore fluid pressure, tensiometers developed at the University of Pretoria were constructed (Jacobsz 2018). The tensiometers are able to measure both positive and negative pore fluid pressures. The tensiometers were installed at various locations throughout the model slopes to give a comprehensive view of the slope response throughout the tests (Fig. 1). Some of the tensiometers dried prior to the centrifuge tests, so that some readings could not be used. The tensiometers that dried are indicated by the white circles in Fig. 1.

Two linear variable differential transformers (LVDTs) were installed at the crests of the model slopes to monitor crest settlement during the tests. To capture the rapid slope failures during the centrifuge tests, a continuous recording of the slope responses was required. Two video cameras were used to monitor and record the response of the model slopes during the centrifuge tests. One camera was placed in front of the window and a second camera was installed to give an oblique view of the slope face. Both video cameras had a display resolution of $\mathrm{1,280}\times 720$ pixels with a frame rate of 30 frames per second.

Tests GT-1 and GT-2 followed the same general procedure. The slight differences between the two tests are discussed in detail subsequently. Both model slopes were accelerated to 60 g. As the model slopes were prepared at a loose state, the slopes experienced settlement during acceleration. Once 60 g was achieved, settlement was monitored until settlement stopped, upon which the wetting process was initiated. The solenoid valve was opened, which allowed the upstream reservoir to fill and the fluid to seep into the model slopes. Once the fluid level in the reservoir was close to the crest of the slopes, the solenoid valve was closed. The hydraulic head in the reservoir was continuously monitored and the solenoid valve was opened and closed accordingly to maintain a constant fluid level in the reservoir. The fluid seeped from the reservoir in the downstream direction, progressively raising the fluid level. The wetting process continued until the model slopes experienced a flow failure. Once the failure events had occurred, the centrifuge was stopped. The pore fluid pressures and crest settlements were monitored throughout the tests. For Test GT-1, the instrumentation sampling rate was 1 Hz. In retrospect, this was somewhat low, therefore the sampling rate was increased to 10 Hz for Test GT-2.

A bi-axisymmetric model triaxial specimen was used to investigate the simulated soil behavior during a CSD stress path. Fig. 4(a) shows the two-dimensional mesh, dimensions, and boundary conditions used for the element-scale numerical CSD stress path simulations. Fifteen noded triangular elements were used for the model. The dimensions of the model were $50\mathrm{mm}\times 25\mathrm{mm}$. The bottom and left boundaries that were normally fixed at the top and right boundaries were set as free. Line loads were applied to the top and right boundaries in order to control the confining and shear stress.

Fig. 4(b) shows the two-dimensional mesh, dimensions, and boundary conditions used for the numerical back-analysis of Test GT-2. Fifteen noded triangular elements were used for this model. The dimensions are at the model scale. The boundary conditions and dimensions were chosen to be as consistent with the centrifuge tests. The sand filters upstream of the model were not considered in the numerical analysis. The slope angle was set at 35° to replicate the initial geometry of the model slope. The initial water level was set at the bottom boundary of the model. The surfaces of the slope crest, slope face, and at the toe of the slope were set as seepage boundaries. The bottom and downstream vertical boundaries were set as impervious, except the groundwater flow boundary at the downstream end, which was set as open to model the boundary conditions of the centrifuge test.

The critical state framework provides a strong framework to formulate constitutive models that can adequately predict soil behavior associated with static liquefaction (Jefferies and Been 2016). Jefferies (1993) introduced such a constitutive model with NorSand. Despite containing “sand” in the name, NorSand is applicable to all soils where contact forces, instead of bonding, control particle to particle interactions (Shuttle and Jefferies 2010). The NorSand model is state dependent, and the response of a soil can be predicted considering the initial void ratio. The NorSand model parameters ($\mathrm{\Gamma }$, $\lambda$, $M_{tc}$,$N$, $H_0$, $H_y$, ${\chi }_{tc}$, $G_{ref}$, $m$, and $\upsilon$) were determined for the gold tailings from a series of triaxial tests.

The critical state line (CSL) in the state diagram ($e-\ln p^{\prime }$ space) was determined from a series of drained and undrained triaxial tests (Fig. 5). CU represents isotropically consolidated undrained triaxial compression tests and the number following denotes the initial confining pressure. CD represents isotropically consolidated drained triaxial tests and the number following denotes the initial confining pressure. The CSL in the state diagram can be expressed aswhere $\mathrm{\Gamma }$ = void ratio of the CSL at $p^{\prime }=1\mathrm{kPa}$; and $\lambda$ = slope of CSL in the state diagram. For the gold tailings used in this study, $\mathrm{\Gamma }=0.915$ and $\lambda =0.039$ provided the best fit.

The volumetric coupling coefficient, $N$, does not vary significantly and is often assumed to be 0.3 for simplicity (Shuttle and Jefferies 2010; Jefferies and Been 2016). To calibrate the state dilatancy parameter, ${\chi }_{tc}$, drained tests need to be conducted on dense specimens to get a dilative response. Despite preparing some triaxial specimens as dense as possible, it was not possible to achieve a dilative response using moist tamping. Thus, ${\chi }_{tc}$ could not be calibrated for the gold tailings used and it was assumed to be 3.2, following Reid and Smith (2021).

Fig. 6 shows the shear behavior of the gold tailings for isotropically consolidated undrained triaxial compression tests, along with the calibrated NorSand model predictions. Figs. 6(a and b) compare the measured and predicted stress-strain behavior and effective stress paths for the undrained triaxial tests, respectively. As bender elements were not available for the triaxial tests, the elastic shear modulus, $G_{ref}$, and the elastic exponent, $m$, were adopted, following Theron et al. (2004), as 60 MPa and 0.65, respectively. However, as noted by Shuttle and Jefferies (2016) and Reid and Smith (2021), there is an issue with the implementation of the elastic bulk modulus, $K_{ref}$, in NorSand when calibrating the model for silts. NorSand estimates $K_{ref}$ from the geophysically measured $G_{ref}$ usingwhich requires a constant and representative Poisson’s ratio, $\upsilon$. As discussed by Shuttle and Jefferies (2016), when a geophysically based $K_{ref}$ is used, the predicted response of the NorSand model produces a too-stiff initial response, which subsequently develops into a sharp transient loss of strength, followed by strong dilation to the correct critical state. Thus, the predicted value does not compare with the measured. This behavior is evident when comparing the measured and predicted values for specimen CU200 (Fig. 6). To circumvent this issue, Shuttle and Jefferies (2016) and Reid and Smith (2021) recommend that the elasticity values be reduced by a factor of up to four for some tailings when calibrating the NorSand model. Thus, a $G_{ref}$ of 15 MPa was adopted for this study to produce the best fit between the measured and predicted values. Table 3 summarizes the model parameters used.

From Fig. 6(a), it is evident that the predicted stress-strain behavior of the calibrated NorSand model is in good agreement with the measured values for loose gold tailings. It can be seen that the NorSand model can also adequately predict the effective stress paths of undrained triaxial tests [Fig. 6(b)]. As the mean effective stress, $p^{\prime }$, decreases, the deviatoric stress, $q$, increases until a peak is reached. After the peak, the deviatoric stress starts to decrease along with the mean effective stress. An important characteristic of the NorSand model is that it is able to capture the strain-softening behaviour of loose gold tailings during undrained shearing. After reaching a peak deviatoric stress, the shear strength decreases where large strains develop.

Also shown in Fig. 6(b) is the CSL determined from a series of drained and undrained triaxial tests. The slope of the CSL under triaxial compression, $M_{tc}$, is 1.327, which corresponds to a critical state friction angle, ${Ø}_c^{\prime }$, of 33°. It should be noted that the numerical model used in this study is under plane strain conditions. While a constant value of $M_{tc}$ can be determined for triaxial compression conditions, as noted by Bishop (1966), the intermediate principal stress has an effect on the critical friction ratio, $M$. The stress state in plane strain varies from one plane strain state to another, because the stress state develops to accommodate the imposed strain condition. The stress state is often denoted by the Lode angle, $\theta$, which is a stress invariant that ranges from $\theta =+30°$ in triaxial compression to $\theta =-30°$ in triaxial extension, with all other stress states lying between these limits (Jefferies and Been 2016). Thus, the critical state friction ratio as a function of Lode angle, $M(\theta )$, is used as a general concept for the critical state friction ratio, where $M_{tc}$ is used as an input parameter into NorSand, to scale $M(\theta )$ (Woudstra 2021).

Based on the definition by Lade (1992), instability lines (ILs) can be found by connecting the peaks of the undrained stress paths within the origin of the $q-p^{\prime }$ plane. The average instability line (IL) for the gold tailings used in this study is shown in Fig. 6(b). As shown by Yang (2002), the slope of the IL, $M_{IL}$, is state dependent (i.e., the slope varies with the state parameter). Thus, the $M_{IL}$ throughout the model slopes will vary. However, the average $M_{IL}$ for the three undrained triaxial tests conducted shown is 0.800. This corresponds to an average mobilised friction angle at the peak, ${Ø}_{IL}^{\prime }$, of 21° applicable to undrained conditions. In order to induce instability, or trigger static liquefaction, the slope needs to be in a loose, saturated state and the soil must have a contractive tendency, which will cause strain-softening during undrained shearing (Take and Beddoe 2014). In addition, the minimum slope angle required for a slope to become unstable under undrained conditions, $\alpha$, can be determined using the following equation (Lade and Yamamuro 2011):where $Ø$ = the friction angle of the soil. For the gold tailings used, this translates to a minimum slope angle of 16.7°.

The implementation of NorSand in PLAXIS 2D requires that a uniform state parameter, ${\psi }_0$, be assigned to the material. To mitigate this issue, multiple layers with different ${\psi }_0$ values can be created and assigned to the model. However, for simplicity, a ${\psi }_0$ of 0.122 was adopted for the back-analysis of Test GT-2, which corresponds to the state parameter at an average depth of 125 mm in the model slope. It should be noted that $K_0$ consolidation will take place in situ. This can lead to shear densification within the slope. Thus, in order to study the impact of shear densification on the instability behavior of the numerical model, a parametric study was conducted. Two additional analyses were conducted using the same calibrated model and numerical procedure, but with ${\psi }_0$ values of 0.03 and 0.06. The ${\psi }_0$ used for the CSD stress path simulations were varied between –0.10 and 0.20, to determine whether a trend between the modified state parameter (the difference between the void ratio at the onset of instability and at the critical state, $\overline{\psi }$) and $M_{IL}$ can be established [after Chu et al. (2003)].

The unsaturated and saturated unit weights were set as 15.0 and $18.6\mathrm{kN}/{\mathrm{m}}^3$, respectively. These values correspond to the respective unit weights of the model slope in Test GT-2. The saturated permeability of the gold tailings was set as $2.3\times {10}^{-8}\mathrm{m}/\mathrm{s}$ (i.e., 60 times less than that of the tailings) for all analyses to satisfy the scaling law for seepage flow due to the viscous fluid used in the test. The initial pore fluid pressure measured by the tensiometers in the centrifuge tests was approximately $-15\mathrm{kPa}$. Thus, a suction of 15 kPa was applied to the model slope in the numerical analysis. For the initial stages of the numerical analyses, modelling the acceleration of the centrifuge, a Mohr-Coulomb model was used with a gold tailings friction angle, $Ø$, of 33°, as well as a $\upsilon$ of 0.3.

The numerical simulation of a CSD test simulates the stress path that soil elements followed during, and the numerical back-analysis of, Test GT-2. Three stages were used to simulate the CSD stress path. As will be discussed later, the initial stages of the back-analysis of Test GT-2 required the use of the Mohr-Coulomb model to obtain convergence. Thus, for Stage 1 of the CSD test simulation, initial conditions were set by applying gravity loading to the model using the Mohr-Coulomb model. During Stage 2, the model was anisotropically sheared under drained conditions using the Mohr-Coulomb model by linearly increasing ${\sigma }_1^{\prime }$ and ${\sigma }_3^{\prime }$ by 120 and 45 kPa, respectively. This resulted in a stress path similar to what a soil element experiences during centrifugal acceleration (to be discussed). After this stage, the material model was substituted for the calibrated NorSand model for Stage 3. During Stage 3, the model was subjected to a CSD stress path by linearly reducing ${\sigma }_1^{\prime }$ and ${\sigma }_3^{\prime }$ at the same rate, keeping the deviatoric stress constant whilst reducing the mean effective stress, until collapse was triggered.

Three stages were used for the numerical back-analysis of Test GT-2. During Stage 1, initial conditions were set by applying gravity loading to the model using the Mohr-Coulomb model. During Stage 2, the gravitational acceleration of the model slope was simulated by increasing the gravity level from 1 g to 60 g while still using the Mohr-Coulomb model in order to obtain convergence. The gravity level was progressively raised in order to mimic the centrifuge modelling process. After increasing the gravity level to 60 g, the material model was substituted for the calibrated NorSand model for Stage 3. During Stage 3, a hydraulic head of 280 mm was applied at the upstream side of the model [Fig. 4(b)] to allow the water to seep into the slope, modelled using a fully coupled transient flow-deformation analysis until slope failure was triggered. No additional stages were required, as the model predicted the drained instability of the slope during Stage 3. Stress points at the locations of T1, T3-T7, and T10 [Fig. 1(b)] were monitored during the numerical back-analysis, allowing for comparison with the measured pore fluid pressure response during Test GT-2.

Figs. 7(a and b) show the measured pore fluid pressures normalized by the initial vertical effective stress, and the measured crest settlements normalized by the initial slope height during Test GT-1. Nine tensiometers (T1-T9) were installed in the slope [Fig. 1(a)], but T9 cavitated prior to the test. All functioning tensiometers measured a pore fluid pressure between $-10$ and $-14\mathrm{kPa}$ at the start of the test. When the slope was at the design centrifugal acceleration, the pore fluid pressures started to gradually decrease [Fig. 7(a)] due to the downward migration of moisture under high acceleration, as was also observed by Poulose et al. (2000). This trend was temporarily interrupted by the stoppages. As the water table neared the tensiometers, the pore fluid pressure started to increase. Once the phreatic surface reached the location of a tensiometer, the pore fluid pressure increased further, measuring positive pore pressures. The measured pore fluid pressure at T5, which was installed at the toe of the slope [Fig. 1(a)], plateaued after roughly 7 h [Fig. 7(a)] as the soil above T5 became completely saturated. Thus, as the water table continued to rise throughout the slope upstream of the toe, the head at T5 stayed constant.

As seen in Fig. 7(b), once the viscous fluid started to seep into the slope, rapid settlement was induced at the crest of the slope, indicating volumetric reduction of the slope due to the loss of matric suction as the slope was wetted (Jennings and Burland 1962; Chang et al. 1981; Sladen et al. 1985; Ng et al. 1998; Vaid and Sivathayalan 2000). After 24 h of wetting, the model slope experienced a flow failure. At this point, rapid fluctuations in the pore pressures were measured, along with sudden crest settlements. After the failure of the slope, excess pore pressures dissipated rapidly. The plateaus in the LVDT readings denote the end of travel.

The images of the failure sequence of the slope in Test GT-1 (Fig. 8) were captured by the video camera, looking downward onto the slope face. The camera recording the side view malfunctioned during the test.

The failure of the model slope was a retrogressive failure that initiated at the toe and progressed upstream until a flow slide failure was triggered. For discussion purposes, the reference time in this section is taken as the time when the first break in the slope occurred [i.e., Fig. 8(a)]. The failure process of the model slope can be divided into four main stages.

Due to large seepage forces, sloughing occurred at the toe and large tensile cracks formed above the toe of the slope before the failure event. A large soil mass near the toe broke away from the slope where the cracks formed, initiating the failure process. This break is defined as Stage 1. The initial break was slow moving, having an estimated velocity of $8\mathrm{mm}/\mathrm{s}$ based on the time of sliding (1.9 s). This break resulted in a loss of confinement at the toe of the slope. This loss of confinement resulted in the soil upstream in the model slope losing confining stress so that it was no longer supported by the material at the toe. This triggered a second soil mass to break away from the slope, defined as Stage 2 [Fig. 8(b)]. This break was more rapid than the initial break, with an estimated peak velocity of $23\mathrm{mm}/\mathrm{s}$. During Stages 1 and 2 no deformation was induced in the upstream section of the slope. However, Stage 3 was initiated after a more rapid third break that occurred 6.5 s after the initiation of the first break [Fig. 8(c)]. The third break had an estimated peak velocity of $58\mathrm{mm}/\mathrm{s}$. Just 0.9 s after the third break occurred, the remainder of the slope rapidly collapsed, which resulted in a flowslide failure [Stage 4, Fig. 8(d)]. This final collapse had an estimated peak velocity of $102\mathrm{mm}/\mathrm{s}$. Kennedy et al. (2020) observed a similar retrogressive landslide failure in a sensitive clay model slope. It is hypothesized that the Merriespruit tailings dam experienced a similar retrogressive failure, causing static liquefaction of the tailings and a large flow failure (Wagener et al. 1998; Fourie et al. 2001; Mànica et al. 2022). Cuomo et al. (2019) simulated a flume test that experienced a similar retrogressive slope failure that transitioned into a liquefaction failure, using the material point method.

Side views of the slope before and after the failure are shown in Figs. 9(a and b), respectively. The observed failure surface cutting through the locations of T4, T5, T6, and T8 is indicated. Only these tensiometers are indicated as only these sensors measured significant pore pressure fluctuations during the failure event (discussed subsequently). Only T5 is indicated in Fig. 9(b), as it was difficult to estimate the post failure locations of T4, T6, and T8.

Fig. 10 shows the measured excess pore pressure ratios, $r_u=\mathrm{\Delta }{u}_w/{\sigma }_v^{\prime }$, and normalised crest settlement immediately before, during, and after the failure of the model slope in Test GT-1, where $\mathrm{\Delta }{u}_w$ is the excess pore pressure and ${\sigma }_v^{\prime }$ is the initial vertical effective stress estimated from the slope profile before failure [Fig. 9(a)] at the location of each tensiometer [Fig. 1(a)]. Time zero in Fig. 10 is set at five seconds before the initiation of the failure event. T1, T2, T3, and T7 all measured negative pore pressures at the time of failure [Fig. 7(a)], indicating that the tensiometers were located in the unsaturated soil of the model slope. These tensiometers did not measure any significant changes in pore pressure during the failure event. Thus, only the measured responses of T4, T5, T6, and T8 are presented in Fig. 10.

As shown in Fig. 8(a), the first break occurred at the toe of the slope, where T4 and T5 were located. This is reflected by the measured changes in pore fluid pressure (Fig. 10). Only T4 and T5 initially measured a response, whilst T6 and T8 did not measure any changes. In addition, at the moment T4 and T5 measured a response, LVDT1 and LVDT2 had not yet measured any crest settlements, indicating that the soil near the crest did not experience any deformation, yet both T4 and T5 initially measured a decrease in pore fluid pressure (i.e., negative excess pore pressure). It is believed that the dilatant response was caused by tension in the tensiometer cables. As the slope started to move, the tensiometer cables would have been pulled, causing some dilation around the tensiometers, resulting in negative excess pore pressures. However, after the slight drop in pore pressure, both tensiometers measured positive excess pore fluid pressures. This fluctuation in the pore fluid pressure confirmed that the soil sheared under at least partially undrained conditions. The positive excess pore fluid pressures indicate the soil’s contractive tendency. T4 measured a peak excess pore fluid pressure of 13.2 kPa $(r_u=0.31)$ and T5 measured a peak of 49.8 kPa $(r_u=1.17)$. T4 measured a $r_u<1$. However, and unlike level ground, full liquefaction is not a prerequisite for slope failure. In contrast, T5 measured an $r_u=1$. After this point, the excess pore pressure dropped slightly before increasing again, this time reaching a peak $r_u=1.17$.

An aspect to consider is the rate of shearing. Looking at Fig. 10, the pore pressures increased by 76 kPa ($r_u=1.7$) at T5 within a span of three seconds. This indicates that the rate of shearing in the region of T5 was rapid. This rapid shearing resulted in the generation of significant positive excess pore pressures, and prevented the dissipation of the excess pore pressures. This led to a degree of localized liquefaction around the toe of the slope. Once the soil in the region experienced a degree of localized liquefaction, the tailings could not support the tailings upstream of the toe. As a result, the tailings upstream of the toe experienced a loss of confinement and a loss of support, which triggered the retrogressive failure (Fig. 8).

As noted previously, T6 and T8 initially did not measure a response at the initiation of the slope failure. However, the retrogressive failure of the slope (Fig. 8) was triggered due to a degree of localized liquefaction in the region of T5, and T6 and T8 only measured a response once the failure had progressed upstream and reached the crest of the slope (Fig. 10). At the time that T6 and T8 measured a response, both LVDT1 and LVDT2 measured changes in crest settlement as well. This indicates that, as soon as the crest started to deform, the tailings at T6 and T8 started to shear. Like T4 and T5, T8 initially measured negative excess pore pressures, which was soon followed by positive excess pore pressures, indicating the tailings had a contractive tendency. In contrast, T6 only measured negative excess pore pressures before returning to the original value. T6 measured a minimum excess pore pressure of $-33.1\mathrm{kPa}$ $(r_u=-0.39)$ and T8 measured a peak excess pore pressure of 28.0 kPa $(r_u=0.33)$. The negative $r_u$ measured by T6 could be due to the tensiometer being pulled from the soil during the failure event.

The measured positive excess pore pressures, specifically around the toe region, led to a significant reduction in effective stress, and thus in the shear strength of the slope. This reduction in shear strength due to the rapid undrained shearing triggered a degree of localized liquefaction at the toe of the slope. Eckersley (1990), Wang and Sassa (2000), Moriwaki et al. (2004), Take and Beddoe (2014), and Ng et al. (2022) also reported positive excess pore pressures during slope failures in flume and centrifuge tests. It was concluded that the excess pore pressures were generated due to the contractive tendency of the loose soils used in the tests. Furthermore, the magnitude of excess pore pressures generated depended on the rate of deformation as well as the soil state (i.e., stress state and void ratio). Based on the observed failure mechanism (Fig. 8), the rapid rate of shearing, the large measured positive excess pore pressure ratios at the toe of the slope (Fig. 10), and the long flow out distance of the failure, it can be concluded that the retrogressive failure was triggered by drained instability at the toe of the slope, which triggered the localized liquefaction at the toe.

Figs. 11(a and b) respectively show the measured pore fluid pressures normalized by the initial vertical effective stress and crest settlements normalized with the initial slope height during the rise of the water table in Test GT-2. Ten tensiometers (T1–T10) were installed in the slope [Fig. 1(b)], but T2, T8, and T9 dried prior to the test and their readings are therefore not shown. The initial pore pressure response during Test GT-2 was similar to that of Test GT-1. However, unlike the observations in Test GT-1, the measured pore fluid pressure continued to decrease around the toe of the slope (i.e., T5). In contrast to Test GT-1, the centrifuge was not stopped during this test. This resulted in the fluid level being raised more rapidly than in Test GT-1. It is thus possible that the degree of saturation of the soil in the toe region was low, resulting in negative pore pressures being measured in the toe region. Owing to the lower degree of saturation in the toe region, it was more difficult to induce localized liquefaction in the toe region. It is also possible that T5 dried out before the phreatic surface reached the tensiometer, which would result in a slow response.

As seen in Fig. 11(b), the crest of the slope settled gradually as the viscous fluid seeped into the slope, which indicates that volume reduction occurred due to a loss of matric suction. The rate of settlement accelerated when the pore pressures near the base of the slope began to rise (T1, T3, and T5). However, the initial rate of settlement was lower than that in Test GT-1, due to the model slope in Test GT-2 being denser that the slope in Test GT-1. The rate of settlement started to reduce throughout the test but increased at around 14 h into the test. At this point cracks started to form in the crest of the slope, visible on the side view (discussed subsequently). After 16.4 h of wetting, the model experienced a sudden slide failure, accompanied by sudden pore pressure fluctuations [Fig. 11(a)] and crest settlement [Fig. 11(b)]. However, the trigger and failure mechanism was different than Test GT-1 (discussed subsequently). The end of the LVDT measurements denote the end of travel.

Images of the slope failure sequence in Test GT-2 as captured by the video camera from the side are presented in Fig. 12. For ease of interpretation, the reference time is set at the time of initiation of the failure event. Fig. 12(a) shows the slope profile and the location of the tensiometers prior to the failure, as well as the observed failure surface. Fig. 12(b) shows the initiation of two tension cracks that formed beneath the crest of the slope, indicating the start of slope failure. At 0.95 s, the cracks had propagated downwards and the crack furthest upstream started to transition into a slip surface [Fig. 12(c)]. Furthermore, a crack had formed at the toe of the slope. At 1.11 s, the slip surface had progressed further downwards and the crack at the toe started to widen [Fig. 12(d)]. The soil mass above the slip surface started to move downstream, causing heave at the toe. At 1.29 s the crack at the crest widened [Fig. 12(e)]. The viscous fluid started to flow into the crack, further destabilizing the sliding soil mass. At 1.68 s the soil at the toe of the slope broke away from the larger sliding soil mass [Fig. 12(f)]. The heaved soil at the toe of the slope was pushed downstream by the larger sliding soil mass. At 5.85 s the fluid that had filled the crack at the crest started to overtop the slope and flowed downstream along the slope surface, ponding at the toe [Fig. 12(g)]. At 9.95 s the sliding soil mass had settled at the downstream side of the model container [Fig. 12(h)]. Some of the tailings that were located upstream of the crest fluidized and settled in the cavity at the upstream end of the slope. The estimated peak velocity of the failure was approximately $35\mathrm{mm}/\mathrm{s}$, which is significantly slower than the peak velocity during the failure in Test GT-1.

The failure observed in Test GT-2 is different from the observed failure in Test GT-1. The failure in Test GT-1 initiated at the toe of the slope and progressed upstream, which was confirmed by the delayed response of the instruments in the slope (Fig. 10). In Test GT-2 the model slope experienced a slip failure triggered by drained instability due to the rising fluid level, and the soil mass failed as one intact mass. The failure observed was similar to the slide-to-flow transformation described by Bolton et al. (2003) and Take et al. (2004). A drained slip failure was triggered by rising pore pressure. The increasing pore pressure results in a CSD stress path, where the mean effective stress is slowly reduced whilst the deviatoric stress is kept constant. As pore pressure increased in the slope, the safety factor against slope failure reduced towards unity, by which time a sudden slope failure occurred, providing little prior warning. The slip failure triggered at least partially undrained shearing within the slope, which generated positive excess pore pressures. The geometry of the slip failure and the slow moving soil mass then transitioned into a more rapid flowslide failure. However, the flowslide failure was still slower than the observed failure in Test GT-1.

Fig. 13 shows the $r_u$ and normalised crest settlements immediately before, during, and after the slope failure in Test GT-2. The vertical effective stress, ${\sigma }_v^{\prime }$, was estimated from the slope profile before failure [Fig. 12(a)] and the location of each tensiometer [Fig. 1(b)]. Time zero in Fig. 13 was set as five seconds before the initiation of the failure event. T2, T8, and T9 dried out before the failure event and T1, T3, and T5 measured negative pore pressures at the time of failure [Fig. 11(a)]. Thus, only the measured responses from T4, T6, T7, and T10 are presented in Fig. 13.

From Fig. 13 it can be seen that, as soon as crest settlements were measured, the tensiometers started to measure a change in pore pressures. Looking at the measured crest settlement, the collapse occurred fairly rapidly, as the crest settled more than 18 mm within 1 s. The slope actually experienced more settlement during this time, as the plateau of the measured crest settlements indicates the end of travel of the LVDTs. Similar to the measured pore pressure response in Test GT-1, the tensiometers initially measured a drop in pore pressure at the start of the failure event due to tension in the cables. The tensiometers were located close to the observed failure surface [Fig. 12(a)] and, as the soil mass moved, the cables of the tensiometers were pulled. However, positive excess pore pressures were measured soon after this initial dilation, indicating that the tailings had a contractive response during shearing. This fluctuation in measured pore pressures confirms that the failure of the sliding soil mass caused at least partial undrained shearing in the model slope.

Like what was observed in Test GT-1, T6 only measured negative excess pore pressures during the failure event. T6 measured a minimum $\mathrm{\Delta }{u}_w=-59.1\mathrm{kPa}$ $(r_u=-0.69)$ (Fig. 13). The negative $r_u$ measured could be due to the tensiometer being pulled from the soil during the failure event. T4 and T10 measured the highest peak excess pore pressure ratios: $\mathrm{\Delta }{u}_w=19.0\mathrm{kPa}$ $(r_u=0.44)$ and 54.2 kPa $(r_u=0.42)$, respectively. T7 measured a peak $\mathrm{\Delta }{u}_w=38.5\mathrm{kPa}$ $(r_u=0.30)$. These peak $r_u$ values are lower than the largest peak $r_u$ measured in Test GT-1, for two reasons. First, the slope in Test GT-2 was prepared in a denser state than in Test GT-1. Thus, it can be expected that the magnitude of positive excess pore pressures generated will be smaller. Second, the failure in Test GT-2 has a slower peak rate of shearing than that in Test GT-1. As the soil was sheared, positive excess pore pressures were generated within the tailings. However, as the excess pore pressures were generated, some excess pore pressures dissipated. Owing to the slower rate of shearing, more excess pore pressures could dissipate during the failure event when compared with Test GT-1. However, it should be noted that, even if partial drainage was present, which would lead to the dissipation of positive excess pore pressures, soils that exhibit quasi steady state behavior, as well as strain softening behavior, can still result in liquefaction events (Jefferies and Been 2016).

In order to validate the numerical model used for the numerical back-analysis of Test GT-2, a CSD stress path was simulated using the model shown in Fig. 14(a). Fig. 14 shows the simulated response of a CSD specimen with a ${\psi }_0=+0.10$ at the start of unloading. Fig. 14(a) shows the simulated stress path in the $p^{\prime }-q$ plane. The soil was sheared anisotropically under drained conditions until Point A, using the Mohr-Coulomb model. This simulated the stress path that a soil element within a model slope experiences during centrifugal acceleration. The average IL for the gold tailings used is also plotted in Fig. 14(a), which shows that Point A is above the IL. After reaching Point A, the material model was switched to NorSand and the unloading process was started, subjecting the model to a CSD stress path. The mean effective stress continued to decrease until instability was triggered at Point B, at which the soil collapsed to Point C. Instability was only triggered at $M_{IL}=1.96$, which is above the CSL in the $p^{\prime }-q$ plane. This is due to the lack of an internal cap in the current NorSand implementation in PLAXIS 2D, which would limit the hardening of the yield surface if implemented (Jefferies 1997). However, Fig. 14(b) shows the soil response during the CSD simulation in the $e-p^{\prime }$ space. During unloading (from Point A to Point B), there was a slight contraction after the initiation of instability at Point B, followed by rapid undrained shearing and a collapse to Point C on the CSL. Fig. 14(c) shows the CSD simulation in the $p^{\prime }-q-e$ space along with the CSL. It is evident that the sudden contraction at Point B triggered the rapid collapse of the soil towards the CSL. During the collapse of the soil, there was a rapid accumulation of axial strain, which resulted in the undrained shearing of the soil. This rapid collapse and undrained shearing triggered the generation of 34 kPa of positive excess pore pressure, which ultimately resulted in the sudden loss of shear resistance of the soil. This behavior was consistent with specimens with positive $\overline{\psi }$ values (i.e., contractive tendencies). For specimens with negative $\overline{\psi }$ values (i.e., dilative tendencies), slower failures were observed with very little change in pore pressures.

Adopting the relationship between the slope on the instability line and the modified state parameter (after Chu et al. 2003), a trend between the simulated $M_{IL}$ versus the corresponding $\overline{\psi }$ can be determined (Fig. 15). As shown in Fig. 15, the change in $M_{IL}$ is relatively moderate for specimens with negative $\overline{\psi }$ (dilative specimens), which corresponds to the work performed by Chu et al. (2003). A more drastic change in $M_{IL}$ is observed for specimens with positive $\overline{\psi }$, which again is consistent with the common observation that the void ratio has a greater influence on the $M_{IL}$ for loose specimens than for dense specimens (Chu et al. 2003). One thing to note is that the $M_{IL}$ for all these simulations are higher than what we would expect from laboratory tests, which is due to the lack of an internal cap in the current NorSand implementation in PLAXIS 2D. Nevertheless, this indicates that the onset of instability during a CSD stress path can be predicted when using NorSand, albeit at a higher $M_{IL}$ values.

Based on the measured results and the video evidence of the failure event, it is hypothesized that the failure of the slope was a slide-to-flow failure (Bolton et al. 2003; Take et al. 2004). To confirm this hypothesized trigger and failure mechanism, a numerical back-analysis of Test GT-2 was conducted using the model shown in Fig. 4(b). In the back-analysis, a head of 280 mm was applied at the upstream end of the slope. The water was allowed to seep into the slope during the fully coupled, flow-deformation stage of the analysis, progressively raising the water table until the failure of the slope was triggered. Failure here is defined as when the solution failed to reach convergence. This is an indication that large plastic strains developed rapidly.

Fig. 16 shows the variations of normalized computed pore fluid pressures by ${\sigma }_v^{\prime }$ with time at the monitored stress points over the duration of the entire test. Comparing Fig. 16 with Fig. 11(a), one notable difference between the numerical and physical models is that the numerical model is not able to capture well the downward migration of pore fluid and the associated increase in matric suction that occurs in the unsaturated zone under high acceleration. A gradual increase in suction was noted during the course of the test on data from tensiometers T1, T3, T4, and T5 installed in unsaturated tailings [Fig. 11(a)]. This increase in suction would have resulted from the downward migration of moisture under high acceleration, along with the tailings, and consequently the tensiometers, drying out over the approximately 1 h test duration. The air-entry value of the tensiometers was 50 kPa, implying that excessive dry-out would not have been accompanied by cavitation, which makes the identification of excessive dry-out difficult. As the water table rose during the test, the measured pore pressures increased, but to values less than predicted numerically. It is believed that this was the result of excessive dry-out and slow tensiometer resaturation as the phreatic surface rose, resulting in undermeasurement of pore pressures. However, the magnitude of the positive pore pressures in the saturated zones were better matched.

Fig. 17 shows the computed stress paths of all the monitored stress points during the numerical back-analysis. It can be seen that all of the monitored stress points were loaded anisotropically during Stage 2 of the numerical back-analysis. During this stage, the deviatoric stress and mean effective stress increased linearly at all points, with the computed stress ratio $(q/p^{\prime })$ varying between 1.02 and 1.23, depending on the location within the slope. After the model reached the test acceleration level, wetting was initiated, and the fluid level was raised until the slope experienced instability. No change in the stress state was observed during the switch from the Mohr-Coulomb model to the NorSand model. T7 and T10 each experienced a slight decrease in deviatoric stress as the water table approached their location. It is possible that a certain amount of shear stress was mobilized during the acceleration stage. As the water table rose, the shear strength of the underlying material reduced somewhat first, resulting in the loss of some of the mobilized shear stress in the material overhead. Once the phreatic surface reached the location of these stress points, the effective stress started to reduce at a constant shear stress. The soil around the locations of T3, T4, T5, and T6 all experienced a slight decrease in mean effective stress once the phreatic surface reached their locations with an approximately constant deviatoric stress. Once the phreatic surface reached the location shown in Fig. 18, the slope collapsed. The initiation of instability is indicated by the triangular markers in Fig. 17. Once instability was triggered, the soil around T3, T4, T5, T6, T7, and T10 experienced a sudden reduction in mean effective stress, along with the rapid generation of shear strains and positive excess pore pressures (to be discussed), with a slight decrease in deviatoric stress. The soil around the location of T1 did not experience much change in the stress state, as it was located in the unsaturated zone above the phreatic surface (Fig. 18). At this point, the numerical back-analysis could not obtain convergence, due to the rapid generation of large shear strains, indicating a runaway slip failure. The stress state at the onset of instability denotes the state of the soil element at which it is no longer able to support the external shear stress imposed on it under undrained or partially undrained conditions (Ng et al. 2004). The soil elements experienced instability at different mobilized friction angles, ranging from 27.3° to 32.0°. At stress ratios above the IL, collapse is triggered due to a small change in the stress ratio, which leads to the triggering of instability within the slope.

Fig. 18 shows the computed shear strain distribution during the failure of the slope. It can be seen that the location of the phreatic surface, which triggered the failure of the slope, is in good agreement with the observed phreatic surface just before failure in Test GT-2 [Fig. 12(a)]. In addition, the failure that was triggered during the back-analysis was a slide-to-flow failure like the failure observed in Test GT-2. Here, a conventional slip failure was predicted. The overall shape, depth and location of the slip failure is also similar to the observed slip failure from Test GT-2. The slip failure has an approximately logarithmic spiral shape that starts just past the toe, progresses upstream and extends slightly upstream of the crest of the slope. When comparing the observed failure surface [Fig. 12(a)] with the computed failure surface (Fig. 18), it can also be seen that the locations of the tensiometers relative to the failure surface are in good agreement with the physical model.

Figs. 19(a and b) respectively show the normalized pore fluid pressure response versus time around the time of failure at the monitored points in the physical and numerical models. The pore pressure responses at the time of failure are difficult to compare in detail because the numerical model can only yield values up to the point of nonconvergence. However, similar trends are seen up to that point, with pore pressures in both the physical and numerical models increasing. In addition, the computed pore pressures are generally higher than the measured pore pressures. As previously discussed, it is possible that the lower measured values are due to the dry-out and incomplete resaturation of the tensiometers. However, a qualitative comparison between the measured and computed pore pressures appears reasonable.

For comparison with the physical model, Fig. 20 shows the computed $r_u$ during the failure of the slope during the back-analysis of Test GT-2. During the numerical back-analysis, the predicted failure occurs instantaneously [Fig. 19(b)]. Thus, the computed $r_u$ is plotted against peak shear stain (e.g., deformation). It can be seen that for all the monitored stress points only a contractive response (i.e., positive excess pore pressures) was predicted during the numerical back-analysis. This is the major difference between the measured and computed $r_u$ values. The measured response of the slope first indicated a dilative response, which was followed by a contractive response (Fig. 13). As mentioned, it is possible that the shearing of the soil during the slope failure caused tension in the cables of the tensiometers. The tension can readily explain the initial dilation. However, shortly after the initial dilation, positive excess pore pressures were generated.

During Test GT-2, T4 and T10 measured the peak $r_u$ values, reaching an $r_u$ of 0.44 and 0.42, respectively (Fig. 13). The computed peak $r_u$ at T4 and T10 were 0.28 and 0.39, respectively (Fig. 20). The computed $r_u$ at T4 is somewhat less than the measured, possibly due to the location of T4. As T4 was situated above the failure surface, the soil around T4 experienced a lesser amount of shearing, which would cause smaller excess pore pressures to be generated. The computed $r_u$ at T10 is close to the measured value. Furthermore, the computed $r_u$ at T7 reached a peak value of 0.33, which is also close to the measured $r_u$ of 0.30.

T5 did not measure any significant changes in the pore fluid pressure during the slope failure in Test GT-2. The measured value of T5 at the time of failure indicated that the soil was still unsaturated. As the seepage in Test GT-2 happened much quicker than in Test GT-1, it is possible that the soil around the toe of the slope was partially saturated. Another reason could be that T5 dried before the fluid level reached T5, thus showing a slow response. The numerical back-analysis of Test GT-2 indicated that the soil at the toe of the slope (around T5) should have been fully saturated and have experienced the largest amount of shearing during the failure event (Fig. 18). This is confirmed by the computed $r_u$ at T5, which reached a value of 0.49, which is the largest computed $r_u$ in the slope during the failure event.

The computed $r_u$ at T6 also differs from the measured value. During Test GT-2, T6 measured only negative excess pore pressures (Fig. 13). The cable of T6 passed through the soil mass that failed in Test GT-2, which caused tension in the cable, which would have generated negative excess pore pressures. During the back-analysis of Test GT-2, T6 reached a peak computed $r_u$ of 0.45. The failure caused the undrained shearing of the soil in close proximity to the slip surface, which caused positive excess pore pressures to be generated. Like the measured $r_u$, negligible changes in the pore pressure during the failure of the slope were computed for T1 and T3.

As mentioned, $K_0$ consolidation can lead to shear densification within the slope. To study the impact of shear densification on the instability behavior of the model, two additional analyses were conducted using the same methodology as the numerical back analysis of Test GT-2, but with ${\psi }_0$ values of 0.03 and 0.06. During both analyses the slopes failed when the phreatic surface was approximately in the same location as that shown in Fig. 18, only slightly higher. The overall shape and depth of the slip surfaces during both analyses were also similar to those shown in Fig. 18. The pore pressure responses during failure for both analyses were also similar to those computed during the back analysis of Test GT-2 (Fig. 20). During both analyses positive excess pore pressures were generated during the failure events due to the triggering of undrained shearing, indicating a decrease in shear strength. However, the magnitudes of excess pore pressures generated were reduced.

For the analysis with ${\psi }_0=0.06$, the peak $r_u$ values at T4, T5, and T6 reduced by between 19% and 23%, compared with the analysis with ${\psi }_0=0.122$. The $r_u$ values at T7 and T10 reduced by 28% and 29%, respectively. For the analysis with ${\psi }_0=0.03$, the peak $r_u$ at T5 reduced by 25%, whereas at T4, T6, T7, and T8 the $r_u$ values reduced by between 48% and 55% when compared with the analysis with ${\psi }_0=0.122$.

These analyses indicate that shear densification will result in lower peak $r_u$ values. However, as the ${\psi }_0$ values were still positive (i.e., the soil has a contractive tendency), the positive excess pore pressures generated during the failure events still resulted in a reduction in shear strength, causing a slide-to-flow failure.

Considering the observed and computed failure mechanisms, as well as the measured and computed $r_u$ values during the failure of GT-2, it can be concluded that the rising water table triggered the drained failure of the slope, which in turn triggered the undrained shearing of the soil. This undrained shearing caused the generation of positive excess pore pressures, which reduced the effective stress. This reduction in effective stress caused a decrease in the shear strength of the soil, especially along the slip surface, which led to the failure of the slope. Based on the observed and predicted failure mechanisms (Figs. 12 and 18), it can be concluded that the slope in Test GT-2 experienced a slide-to-flow failure, similar to that described by Bolton et al. (2003) and Take et al. (2004). It should be pointed out that, although a drained slope stability analysis would predict this specific trigger and subsequent failure, tailings storage facilities with adequate drained factors of safety might have inadequate undrained factors of safety when subjected to a trigger capable of initiating an undrained response.

For interest’s sake, the yield strength ratio [$s_u(\mathrm{yield})/{\sigma }_v^{\prime }$] at which failure occurred was determined, after Olson and Stark 2003. Here $s_u(\mathrm{yield})$ is the peak shear strength available during undrained loading and ${\sigma }_v^{\prime }$ the pre-failure vertical effective stress. Although not strictly applicable, as the failure trigger was drained instability followed by undrained shearing post failure, an equivalent $s_u(\mathrm{yield})/{\sigma }_v^{\prime }$ of 0.54 was calculated for Test GT-2. This value is much higher than the typical range of values found by Olson and Stark (2003), as the failure was not the result of undrained instability.