On some uncertainties related to static liquefaction triggering assessments Open Access

The grave risk posed by static liquefaction to the safety of tailings storage facilities (TSFs) has been made clear by a series of recent failures where, in each case, the trigger mechanism was clearly shown to be of static origin (AECOM, 2009; Jefferies et al., 2019; Morgenstern et al., 2016; Robertson et al., 2019). For each of these failures, a robust panel investigation was carried out, where the contribution of the historical development and conditions within the TSF on the stress path leading to static liquefaction triggering (SLT) was examined. This process usually involved stress–deformation modelling, often utilising the constitutive model NorSand (Jefferies, 1993) based on the state parameter Ψ (Been and Jefferies, 1985). Publication of these modelling exercises, along with greater focus on SLT, has served to increase the application of such techniques in tailings engineering practice. However, despite the clear benefits of the techniques now available and their greater application, a significant number of uncertainties remain relating to various aspects of the SLT analysis process that warrant further investigation and discussion.

Of the uncertainties relevant in a static liquefaction assessment, those related to the distribution and variability of in situ Ψ appear to have received the most attention (Hicks and Onisiphorou, 2005; Robertson et al., 2019). Furthermore, tentative means to account for various uncertainties in a probabilistic framework have been made, for example by Stoutjesdijk et al. (1998) who investigated the effects of different in situ states and slope geometry on whether zones of an idealised slope were in a potentially unstable state, with similar methods employed by Mathijssen et al. (2015) in the design of dredge spoil storage areas. However, despite these advances, several additional uncertainties related to SLT assessments have not yet been studied in detail. For example, in all of the works related to static liquefaction cited above and in the recent work of forensic investigation panels, the instability stress ratio η_IL (= q/p′ at occurrence of instability) relevant to in situ conditions was developed based on triaxial compression (TX-C) tests (i.e. with a principal stress angle (α) of zero degrees from vertical), and with apparently differing views as to how the axisymmetric conditions of TX-C test should be related to in situ plane strain conditions. Furthermore, Stoutjesdijk et al. (1998) and Mathijssen et al. (2015) did not account for whether a potentially unstable condition in one zone was likely to propagate elsewhere owing to brittle stress redistribution – something that is likely a key consideration in the failure process of brittle soils (Been, 2016; Jefferies, 2018) and which can be readily observed in a video of the recent Brumadinho failure (YouTube, 2021).

The effects of the following uncertainties related to SLT assessments were examined in this work.

• The geostatic stress ratio K₀, including its effect on both the in situ stress ratio (η = q/p′) and the principal stress angle from vertical (α) below a slope.

• Whether plane strain conditions should be assumed to modify or scale η_IL compared with axisymmetric tests.

• The poorly understood potential for α rotating away from the vertical to affect the relevant value of η_IL under drained SLT.

These factors were first examined by means of a literature review and then implemented in a series of numerical models of an idealised slope to examine their effects. Given the urgent need for reliable static liquefaction assessments of TSFs in the mining industry, the authors are of the view that better quantification of uncertainties is crucial for engineering practice. This paper serves to address some of these in a preliminary manner for the purpose of illustration, with the aim of providing an improvement to SLT assessments and impetus for further works to reduce the uncertainties outlined.

The importance of the inherent uncertainty in K₀ (⁠$=σ′h/σ′v$⁠) in geotechnical practice as it relates to Ψ interpretation from cone penetration tests (CPTs) has been discussed in detail by Jefferies and Been (2015). However, to the authors’ knowledge, the effect of K₀ on in situ η, and hence the potential for SLT, has not received detailed examination. A schematic illustration of the effect of different values of K₀ on an element of soil under vertical effective stress of 100 kPa is provided in Figure 1. As shown in the figure, the value of K₀ has a large effect on how ‘close’ the element may be to the initiation of instability.

Most laboratory element tests carried out to measure K₀ on sandy soils and sandy or silty tailings – which tend to be reconstituted samples undergoing virgin consolidation – suggest K₀ values similar to the Jaky equation (⁠$1−sin⁡ϕ$⁠) (i.e. ∼0.5) are reasonable. Alternatively, based on tests involving repeated cycles of low-amplitude cyclic strains (Youd and Craven, 1975; Zhu and Clark, 1994), Jefferies and Been (2015) argue that, over time, a higher K₀ is likely to develop in situ. For example, it would not be unexpected for a TSF to be exposed to many low-amplitude cyclic events over the decades of its operational life owing to small seismic events and/or blasting at nearby mine works – although a similar observation would be true for essentially all soil deposits worldwide, at least in seismically active regions. There indeed is some evidence of higher K₀ values from self-boring pressure meter tests in hydraulically placed deposits (e.g. Ghafghazi and Shuttle, 2008). Alternatively, some extrusion-type deformation processes appear likely to reduce K₀ in situ (e.g. Morgenstern et al., 2016), again in a way that cannot be informed by typical element-test measurement of K₀. Finally, there are persuasive arguments that many historical measurements of K₀ in element tests suffer limitations (Talesnick et al., 2020), thus adding to the uncertainty of this parameter.

The different views on relevant K₀ values can be illustrated through a review of the assumptions adopted and/or values produced through modelling in three recent TSF failure investigations. These are summarised in Table 1, where the K₀ values assumed in CPT interpretation by each investigation panel are presented – clearly indicating that there is a lack of consensus even among failure investigation panels. Further, the NorSand constitutive model (Jefferies, 1993), adopted by all of the panels listed in Table 1 to varying degrees, generally provides estimates of K₀ > 0.7 for the inputs developed in these investigations (Jefferies et al., 2019; Morgenstern et al., 2016), which are higher than typical values seen in K₀ triaxial tests (or similar) on reconstituted sands or sand/silt tailings. Indeed, direct comparisons of measured K₀ values with predictions from NorSand as calibrated for the same sand have demonstrated the higher values that result from NorSand (Castonguay and Konrad, 2016). Clearly, the state of current knowledge is such that the value of K₀ must be recognised as a source of uncertainty in the context of estimating in situ η.

Prior to proceeding, it should be noted that while K₀ represents a useful ‘shorthand’ for in situ horizontal stress conditions, it is most useful in the context of material below level ground (BLG) where stresses can be assumed axisymmetric and α is vertical. Below a slope, where α has rotated (and thus the vertical effective stress is not the major principal stress), the values K_c (σ′₃/σ′₁) and η are more useful. However, to describe and compare different models in this paper, the term K₀ is used, referring therefore to the portion of each model BLG conditions not affected by the slope geometry.

An important difference between typical laboratory tests to infer η_IL in the laboratory (i.e. triaxial tests) and in situ conditions is that the triaxial device tests under axisymmetric conditions, whereas locations in situ below a slope are plane strain. A useful parameter to describe and compare intermediate principal stress conditions across a range of different situations is the intermediate principal stress ratio defined by Bishop (1966):

There is a reasonable quantity of plane strain and hollow-cylinder tortional shear (HCTS) data to allow the effects of b on various aspects of soil behaviour to be inferred. Of particular note here is the clear evidence of a dependence of the critical-state friction ratio (M) on b (e.g. Jefferies and Shuttle, 2011). Importantly, in the range of b likely relevant to below-slope plane strain conditions (0.2–0.4?), a reduction in M of 10–20% could be expected, where the term M_b can be used to refer to the relevant critical-state friction ratio under a particular value of b. This observed effect of b on M_b has been directly incorporated in constitutive models, for example by means of the expression proposed by Jefferies and Shuttle (2011), illustrated in Figure 2.

Turning to the focus of this paper – the initiation of instability – the most comprehensive and relevant dataset is that reported by Wanatowski (2005), where comparisons between the initiation of instability under undrained and constant shear drained (CSD) tests were carried out in TX-C and plane strain tests. One of the primary outcomes of that work is summarised in Figure 3, which shows the trends of η_IL against Ψ (at the occurrence of instability, referred to as the modified state parameter by Wanatowski (2005)) for TX-C and plane strain tests carried out by Wanatowski (2005), as synthesised by Wanatowski and Chu (2012). Importantly, scaling η_IL by the relevant value of M_b for a particular test type or b value (TX-C or plane strain) was found to produce a reasonably unified trend. This finding regarding the scaling of η_IL based on intermediate principal stress conditions was applied by the Cadia panel (Jefferies et al., 2019) when identifying what series of events would lead to η_IL being reached in tailings below the TSF slope; for example, at a given value of Ψ, η_IL is reduced by the same ratio as M_b suggested by the relationship shown in Figure 2. Alternatively, in preceding work by the Fundão panel (Morgenstern et al., 2016), it appears that the value of η_IL obtained from TX-C tests was used directly in the interrogation of a plane strain numerical model without any scaling.

While the data of Wanatowski and Chu (2012) (Figure 3) show a clear dependence of η_IL on b and the possibilities of unifying the datasets through scaling to M_b, an important additional factor is that, in the same overall study (Wanatowski, 2005; Wanatowski and Chu, 2007), a dependence of the critical-state line (CSL) elevation in the e–log p′ space on b was also demonstrated – with the plane strain CSL (CSL_PS) being at about 0.04 ‘higher’ elevation than the TX-C CSL (CSL_TX) in the compression plane. This observed dependence of CSL elevation has itself been incorporated in some numerical models of element tests with a plane strain condition (Jefferies et al., 2015). It is crucial to note that, for the data presented in Figure 3, the Ψ for each test type was calculated on the basis of the CSL relevant for that test type – that is plane strain tests against the CSL_PS and TX-C tests against CSL_TX. Alternatively, if Ψ for all of the tests were calculated based on the CSL_TX – implicit to many published numerical SLT analyses, as discussed subsequently – this would ‘shift’ the plane strain tests to the right by ∼0.04 Ψ on Figure 3, thus arguing against the scaling of η_IL by b on the basis of the reduction seen with M_b.

Returning to the work of the Fundão and Cadia panels, where NorSand was implemented in a boundary element problem to deduce events that led to a value of η_IL sufficient to trigger liquefaction, there appears to be no mention of implementing a dependence of CSL elevation in the compression plane on b. A review of the relevant model outputs and descriptions of the modelling in these works do not provide any indication this was carried out in either case. Alternatively, it appears clear that the Fundão panel did not scale η_IL with b, while the Cadia panel did. This summary therefore suggests an ongoing uncertainty and lack of consensus in engineering practice as to the implications of the available experimental data on how – or if – η_IL should be scaled with b.

A final uncertainty comes from the fact that other plane strain tests (e.g. Cornforth (1961), as interpreted by Jefferies and |Shuttle (2002, 2011)) do not appear to suggest a dependence of CSL elevation on b, while some discrete-element modelling (Huang et al., 2014) does suggest such a dependence, but of a much lower magnitude than the observations of Wanatowski and Chu (2007).

In summary, it appears difficult to conclude definitively whether η_IL should be scaled with b in numerical SLT analyses, or even what the current state of engineering practice is in this regard. As such, this is included as an uncertainty for examination in the work that follows.

It has long been recognised that soils deposited under a vertical gravity field can develop an anisotropic fabric that results in their engineering properties differing depending on the direction of shearing (Arthur and Menzies, 1972; Oda, 1981). However, despite decades of recognition of the effects of such anisotropy on the undrained response of sands (Shibuya et al., 2003; Sivathayalan and Vaid, 2002; Uthayakumar and Vaid, 1998), clays (Seah, 1990) and tailings (Reid et al. 2021), incorporation of this aspect does not seem to be common in current tailings engineering practice –as suggested by the work of recent TSF failure panels. However, given the widespread evidence for the effect of anisotropy on undrained strength and the initiation of instability, it seems difficult to argue against its inclusion in static liquefaction assessments (e.g. Reid, 2020). Indeed, the anisotropic nature of sandy soils is directly incorporated into common ‘screening level’ static liquefaction assessment tools (Sadrekarimi, 2014, 2016).

For plane strain below-slope conditions, the principal stress angle from vertical (α) can be calculated by:

where τ_xy is the shear stress acting on the horizontal plane (because of the slope geometry) and σ_v and σ_h are vertical and horizontal stresses, respectively. It is important to recognise that the denominator being a relationship between horizontal and vertical stress means that K₀ (itself a major uncertainty, as discussed previously) plays a role in the amount of α rotation under a given amount of τ_xy.

To provide context as to the likely effect of α rotation on η_IL, the data summarised by Reid (2020) sourced from relevant studies in the literature (Shibuya et al., 2003; Sivathayalan and Vaid, 2002; Uthayakumar and Vaid, 1998) are presented in Figure 4. This shows the variation in η_IL with α for three different HCTS programmes. The clear effect of α on η_IL is indicated by each test programme. Two important aspects of the data presented in Figure 4 are as follows.

• (a)

  Each dataset used to allow interrogation of the effect of α was carried out at the same value of b, and thus the trends seen are independent of potential effects of the intermediate principal stress ratio.

• (b)

  The HCTS device does not suffer the limitations of other devices such as the direct simple shear device, which – while finding increased application in static liquefaction studies (e.g. Reid and Fourie, 2019) – does not allow all principal stresses to be measured, thus giving an erroneous perception of maximum shear stress and η_IL (e.g. Jefferies and Been, 2015).

While the data in Figure 4 are persuasive, some limitations should be noted.

• The results are only for undrained tests, as no drained SLT test results at a range of α values using HCTS appear to have been published. However, given the general agreement between undrained and drained tests to infer η_IL at a given Ψ (Chu et al., 2012; Daouadji et al., 2010; Dong et al., 2016; Imam et al., 2002; Wanatowski and Chu, 2012), the tendency for α to also affect η_IL under drained SLT seems highly likely.

• The tests conducted by Uthayakumar and Vaid (1998) and Shibuya et al. (2003) were consolidated isotropically, then sheared undrained under different α values. This is unrealistic compared with the drained rotation of principal stress that would occur as a TSF slope is developed over a period of decades, and where the drained rotation of α itself is likely to affect the anisotropy of the soil. Alternatively, the tests carried out by Sivathayalan and Vaid (2002) were consolidated anisotropically under the same α as subsequent undrained shear – a much more representative test method when compared with in situ behaviour. As such, the authors suggest that the data of Sivathayalan and Vaid (2002) represent the best currently available dataset to infer the effect of α on η_IL. This dataset was therefore given precedence in the current study when developing the relationships to be used subsequently in this work.

Additional trendlines were added to Figure 4, based on the following assumptions.

• Assumption 1: α has no effect on η_IL (implicit in recent failure panel works for example).

• Assumption 2: η_IL reduces by 0.5% per degree that α rotates from vertical.

These two assumptions were adopted in the numerical modelling carried out in this work.

To frame the discussion, an idealised TSF was developed to allow interrogation of the effects of the various uncertainties discussed earlier in the paper. This was carried out using the finite-difference code FLAC v8.0 (Itasca, 2016). The geometry of the idealised TSF is outlined in Figure 5, with a height of 50 m, a slope of 1V:2.5H and initially being ‘dry’ in the sense of zero pore pressures throughout the model (although consideration of suction effects was generally excluded in the current study). The TSF sits on a 25 m stiff foundation, assigned strength and stiffness parameters such that the tailings behaviour was not influenced by foundation deformations (i.e. an extrusion-type contribution to the stress path is avoided) and to produce a K₀ BLG of 0.7, consistent with a stiff overconsolidated soil. Zone sizes in the model were 1 m square. A 20 m area at the toe of the TSF, representing a starter embankment, was assumed to be dilative (i.e. strength characterised by the indicated strength values regardless of phreatic conditions in all simulations).

The model was first developed with the tailings behaviour simulated using NorSand with the parameters developed by the Cadia failure investigation panel (Jefferies et al., 2019), as listed in Table 2. This constitutive model was adopted as it represents current industry state of practice for such assessments, as evidenced through applications in recent failure investigations in addition to Cadia (Morgenstern et al., 2016; Robertson et al., 2019). The version of NorSand used was the user-defined model version outlined by Shuttle and Jefferies (2018); this was adopted as it appears to be that used by the Cadia panel. Alternatively, the version of NorSand included in the FLAC software code (Itasca, 2020) appears to have subtle differences to that implemented by Shuttle and Jefferies (2018) and the relevant panels.

As the NorSand model is built up in layers, a ‘seed’ Ψ value is required for each layer upon initial placement. A seed value of +0.06 was adopted herein, consistent with the modelling of the Cadia panel.

The condition of the model after placement of all 50 layers is shown in Figure 6 in terms of η, Ψ and α. Below the slope, an increasing stress ratio and rotated principal stress angle occur by virtue of the slope geometry. The figure shows that increased η led to densification of the tailings in this area, with Ψ values as ‘low’ as +0.03 near the toe (compared with the seed value of +0.06). Similar below-slope reductions in Ψ were observed by the Cadia panel, who termed this behaviour shear densification. An element simulation is provided in Appendix 1 for further demonstration of this effect.

It is noted that the shear densification that occurs below the slope leads to higher values of η_IL in this area by virtue of the denser state of the material, as discussed later. Also noteworthy are the low values of η away from the slope, which occur due to the relatively high K₀ value produced (∼0.80 BLG). Alternatively, values of α up to 50° occur below the lower portion of the slope; that is, a magnitude of rotation such that it could be expected to affect η_IL based on Figure 4. The magnitude of rotation is affected significantly, as noted previously, by the K₀ value that develops in the model.

To interrogate the model and assess SLT, a series of NorSand element tests were carried out at various stresses and initial K₀ values. These models were run using the NorSand VBA program for TX-C element simulations (e.g. Shuttle and Jefferies, 2016) using the Cadia parameters outlined in Table 2. The results of these simulations are summarised as plots of Ψ against η_IL in Figure 7. Data from the panel's element testing are included. These simulations and experimental data show the expected general trend of a decreasing η_IL with increasing Ψ (e.g. Chu et al., 2003). The relationship adopted for use in this paper is also indicated in Figure 7 – the effect of other potential fits to the experimental and element simulation data presented are addressed in subsequent parametric analyses. It is noted that the data considered in Figure 7 are all for TX-C conditions; inclusion of the potential effects of α and b are discussed subsequently.

For the purpose of assessing static liquefaction failure, a phreatic surface approaching the slope was used as the relevant trigger mechanism – a stress path referred to as CSD. The CSD stress path was selected as it is amenable to simple quantification in terms of a phreatic surface level or location – which is beneficial for comparing uncertainty between various models – and it avoids inclusion of additional foundation effects and potential strain-softening uncertainties that might be required if using an extrusion-type triggering process (e.g. Morgenstern et al., 2016). Furthermore, it is a stress path commonly attributed as the trigger to brittle failures such as Aberfan (e.g. Jefferies and Been, 2015) and Wachusett dam (e.g. Olson et al., 2000). The phreatic surface changes applied to the model to produce a CSD stress path within the tailings are schematically illustrated in Figure 8. Pore pressures of zero were maintained in all areas above the relevant phreatic surface for a particular stage of the modelling process. A hydrostatic pore pressure was assumed below the phreatic surface.

The phreatic surface was shifted towards the toe of the slope in 1 m horizontal increments. When the potential for failure under a given phreatic level was to be checked, the following steps were carried out.

• (a)

  The model conditions were saved, after which all of the tailings were assigned the Mohr–Coulomb constitutive model with ϕ = 40° and c′ = 5 kPa. The Ψ values that had developed within each element from previous model steps were retained for use in subsequent calculations. The Mohr–Coulomb model was used for the SLT assessment process for two reasons: to match the parameters used in subsequent parametric analyses based on the Mohr–Coulomb model and to prevent drained changes in Ψ from occurring after some elements undergo SLT and shed shear stresses onto adjacent elements, as would occur – unrealistically – in the framework employed herein where drained conditions (i.e. fluid bulk modulus of zero) were used. Drained conditions were used for modelling the CSD path, given difficulties seen with the numerical implementation of NorSand in FLAC for undrained conditions (Shuttle and Jefferies, 2018) and general uncertainties related to the modelling of the unloading path as it approaches failure (e.g. Jefferies and Been, 2015). It is noted that this process of using NorSand to calculate in situ stresses, while employing the Mohr–Coulomb model to then further check stability of the slope, is consistent with that employed by the Cadia panel (Jefferies et al., 2019).

• (b)

  Each element below the phreatic surface was examined to see if η > η_IL for the relevant value of Ψ. The following four models were used.

• •

  Assuming either that η_IL should or should not be scaled to the same magnitude as that of M_b based on the relationship proposed by Jefferies and Shuttle (2011) and presented in Figure 2 (i.e. b scaling).

• •

  Adopting either assumption 1 or assumption 2 with respect to the potential effects of α on η_IL (i.e. α scaling).

• (c)

  For each element where η > η_IL, a liquefied strength ratio (⁠$sr/σ′vo$⁠) was calculated based on the Ψ of the element, which was used to calculate the relevant liquefied shear strength. That shear strength value was then applied as a cohesion input to the relevant element. For such elements the friction angle was also reduce to zero. The relationship between Ψ and the liquefied strength ratio adopted in the NorSand models is provided in Appendix 2.

• (d)

  The model was then re-run and additional checks were carried out to see if more zones now had η > η_IL (i.e. brittle stress redistribution).

• (e)

  If no additional zones underwent an initiation of instability (i.e. η > η_IL) and the model could be brought to equilibrium with all of the triggered zones reduced to a liquefied strength, the saved model from the commencement of the process (step (a)) was reactivated and the phreatic surface was again raised and the process repeated until the slope failed or until the phreatic surface reached the toe of the slope.

It is acknowledged that the Mohr–Coulomb model employed in the above steps means the models would not be expected to produce reliably all the aspects of tailings stress–strain behaviour during brittle strength loss. However, similar approaches were used successfully by the Cadia panel and the ‘manual’ switching to the Mohr–Coulomb model with liquefied strength inputs also finds usage in much state-of-practice seismic analysis (e.g. Chowdhury et al., 2019).

The pseudo-codes of the modelling process are provided in Appendix 3, with the full Flac input code used to carry out these models also provided as online supplementary material.

The outcome of the triggering assessment is summarised in Table 3 and Figure 9 in terms of the final phreatic surfaces achieved in the four models, along with the areas of the models where η > η_IL at the first occurrence of SLT under the failure phreatic surface and the maximum extent of triggered zones once brittle stress redistribution had completed (and ongoing deformations were occurring in the model). While the methods used here to infer when the slope would undergo failure do not produce a ‘failure surface’ as such, an approximate extent of the initial failed area can be inferred by demarcating the boundary between zones of significant displacement. This boundary is shown in Figure 9 for each model. It is noted that this only represents the initial area of the slope that is undergoing significant displacements. In a real slope, subsequent retrogressive failures would likely occur. However, attempting to model such high-strain aspects of the predicted failures – say through a rezoning scheme (e.g. Reid, 2013) – was beyond the scope of the current work.

What the results make clear is the significant variation in predicted failure of the slope based on the assumptions made regarding the effect of α on η_IL. The effect of b is smaller overall, having a more significant effect for models 1 and 3 (where α effects are excluded) as the relevant zones where SLT occurs in those models are further below the slope where higher b values are present (refer to Figure 6(d)). It is emphasised that the existing literature and the apparent state of engineering practice could seemingly provide support for any of these assumption combinations.

To provide reference data to interrogate the ‘realism’ of the NorSand model outcomes, three additional models were used where failure of the slope was identified using the shear strength reduction (SSR) approach in Flac to identify when a factor of safety (FoS) of unity would be obtained for the following scenarios.

• Assuming all tailings below the phreatic surface undergo SLT with a resulting liquefied strength ratio (⁠$sr/σ′vo$⁠) of 0.10. This represents the conservative assumption, increasingly advocated for standard engineering practice (Been, 2016; Jefferies, 2018), given the inherent difficulties in predicting SLT in contractive, brittle materials.

• Carrying out a more conventional SLT analysis procedure, with assumed yield stress ratios (⁠$su(yield)/σ′vo$⁠) of 0.20, 0.25 and 0.30. This represents a simplified version of common yield stress ratio type techniques (Olson and Stark, 2003; Sadrekarimi, 2016).

The $su(yield)/σ′vo$ values of 0.20, 0.25. and 0.30 used for comparison were selected as they (a) represent typical ranges for the back-analysis of static liquefaction flow slides, including those triggered by a CSD path (e.g. Olson et al., 2000), (b) are consistent with the undrained strength ratio of 0.20 inferred for near- and below-slope conditions by the Cadia panel and (c) agree with the range of η_IL values adopted in the parametric analyses (refer to Appendix 4 for derivations). While it is noted that application of a single value of $su(yield)/σ′vo$ does not itself capture the effects of brittle stress redistribution, the use of a single average value to characterise stress conditions within loose slopes at the initiation of failure (as done through back-analysis works to support triggering screening methods (e.g. Olson and Stark, 2003)) typically indicates average values in the range of those adopted here in the SSR analyses. Therefore, an expectation that a more advanced SLT analysis based on numerical techniques should provide results of a similar magnitude to yield stress ratio-type approaches seems justifiable.

The results of these SSR models are presented along with the NorSand model failure surfaces in Figure 10, with the phreatic surface distance at failure for each model annotated for comparison. The results indicate that the NorSand models generally ‘straddle’ the SSR models using a range of $su(yield)/σ′vo$ ratios seen in previous liquefaction case histories (e.g. Olson and Stark, 2003) – with the NorSand models including α indicating failure earlier than the SSR approach, while those not including α fail later. Indeed, the failures indicated by the models including α are such that they begin to approach the SSR model with the assumption that all saturated tailings trigger with a subsequent $sr/σ′vo$ of 0.10. This result might suggest that the magnitude of the effect of α on η_IL is lower in reality than that made with assumption 2 (i.e. 0.5% reduction per degree that α rotates from vertical), at least on the basis of the outcomes of the NorSand constitutive model as applied herein and compared with the SSR approach.

One area requiring additional focus in the context of the discussion of α effects is the in situ value of K₀; for example, as NorSand with the Cadia parameters produces relatively high K₀ values (∼0.80 BLG), this leads to significant principal stress rotation. As such, further parametric assessments under a wider range of K₀ values than can be produced with the NorSand model were carried out.

Based on the previous discussion around the effects of K₀, b and α, a further series of models was developed to allow a parametric investigation. In these simulations, the Mohr–Coulomb constitutive model was used for the tailings throughout the simulation process, with the same strength values of ϕ′ = 40° and c′ = 5 kPa as used in the SLT assessments carried out as a part of the previous NorSand modelling. To obtain a range of plausible K₀ values for the sensitivity analyses carried out, Poisson's ratio (ν) was varied as necessary, based on the following expression for elastic models under level-ground conditions:

It is also noted that specification of a particular ν will also dictate intermediate principal stress development, and thus b. This means that selection of a particular ν controls both K₀ and b, and the resulting interaction of these values will have a major impact on when particular elements of a given model will trigger. For example, for an elastic model, σ₂ can be calculated as:

The corresponding values of ν used in the analyses and the resulting BLG K₀ values used in the parametric analyses are outlined in Table 4. The foundation parameters were taken as the same as those used in the NorSand models, while elastic stiffness values for the tailings were ten times less than those for the foundation.

Given the relatively high strengths used in the Mohr–Coulomb model of the tailings, the behaviour of the model is elastic, and thus Equations 3 and 4 control the development of K₀ and intermediate principal stress, respectively. It is noted that Equation 3 applies to confined vertical loading (i.e. an oedometer test or BLG conditions). As noted previously, the K₀ of each model will generally increase closer to and below the slope, owing to principal stress rotation.

It is emphasised that while modification of Poisson's ratio to achieve a desired K₀ is simple and effective for the purposes of this study, this approach can hardly be considered leading practice for stress–deformation modelling, which is generally based on more complex constitutive models that can reproduce far more subtle and fundamental behaviours. However, while complex constitutive models are undoubtedly useful, many will result in a particular value of K₀ when provided with a set of material-specific inputs – as discussed previously in the context of the Cadia parameters in NorSand.

Other relevant parameters applied to the tailings material in the model, for the purposes of the parametric SLT assessment and to closely align with the NorSand modelling as far as practicable, are as follows.

• η_IL(TX-C) values of 0.8, 0.9 and 1.0 to cover the conceivable range for the Cadia tailings under typical Ψ values relevant to below-slope conditions (refer to Figure 6). It is noted that assigning a constant value of η_IL(TX-C) implies an assumption of a constant Ψ in the areas below the slope relevant to failure in these parametric models.

• Scaling (or not) of η_IL(TX-C) on the basis of both α and b, in the same manner as previously described for the NorSand models

• A liquefied strength ratio (⁠$sr/σ′vo$⁠) of 0.1, based on the average Ψ value indicated in the triggered region of the NorSand models (refer to Figures 6 and 9) and the Ψ−$sr/σ′vo$ relationship outlined in Appendix 2.

These inputs resulted in 48 different models, as summarised in Table 5. The same idealised slope and CSD path as previously outlined were adopted. Similar SLT assessment procedures were carried out as used for the NorSand models with slight simplifications – for example, switching constitutive models while tracking Ψ was not required (see Appendix 3 for the pseudo codes).

The significant differences in the model stress conditions that result from different Poisson's ratio inputs are highlighted in Figure 11 based on η_IL, α and b for each model after initial placement, prior to implementation of the CSD path. The previous NorSand model is included for comparison. Consistent with previous discussions, models with a low Poisson's ratio (and thus lower K₀) result in higher η. Furthermore, as indicated in Figure 11, the models with higher K₀ produced greater values of α at a given point below the slope.

The results of all the parametric models are summarised as a histogram in Figure 12 and a cumulative distribution in Figure 13, in both cases in terms of the phreatic surface distance from the toe at which failure was predicted. The distance for each model is listed in Table 5. A wide range of results is evident, with failure predicted when the phreatic surface had advanced to anywhere from 8 m to 55 m from the toe.

A first observation of the results in Figures 12 and 13 would be that, logically, none of the models fail until the phreatic surface has come to at least within 57 m from the toe because, until this point, the slope is stable even if a liquefied strength of $sr/σ′vo$ is 0.10 is assumed in all elements below the phreatic surface. However, the majority of the models predicted failure with the phreatic surface further from the toe than the $su(yield)/σ′vo$ = 0.20 SSR model. This is an interesting outcome because, in back-analyses of statically triggered flow liquefaction case histories, $su(yield)/σ′vo$ values are rarely – if ever – below 0.20 (e.g. Olson and Stark, 2003); indeed, this has led the value of 0.20 to be suggested as a useful screening value in itself (Morris, 2008) This suggests that some aspect of the parametric analyses carried out may be weighting the results to be too conservative compared with actual field-scale behaviour, despite the η_IL(TX-C) values clearly being plausible. The following hypotheses as to why this may be occurring are suggested.

• In the modelling framework, the assumption that triggering is instantaneous and assured when an element reaches a state such that η > η_IL may be slightly unrealistic. For example, it is possible for loose soils to be at a stress state above η_IL under certain stress paths; for example, where p′ is increasing under drained loading, at which point almost any subsequent perturbation causing a decreasing p′ could lead to instability. If such behaviour is typical in situ, it would make the current analysis framework slightly conservative.

• It could potentially be argued that the results of the current study provide further evidence of the higher K₀ values likely to develop within full-scale TSFs and other loose deposits, as models with low K₀ values generally fail with the phreatic surface at the greatest distances from the toe; in other words, excluding the K₀ = 0.5 models would shift the median result of the numerical SLT analyses to much closer to the SSR $su(yield)/σ′vo$ = 0.20 analysis. The implications of K₀ in this context are examined in more detail subsequently.

• The use of η_IL(TX-C) = 0.8 may be slightly too low for below-slope conditions because, while such a value is clearly plausible for loose soils prepared in element tests, shear densification under a slope may make such loose states less likely in areas where SLT will initiate. For example, if another set of models assuming η_IL(TX-C) = 1.1 were to be carried out – perhaps ‘just’ plausible on the basis of the η_IL–Ψ relationship in Figure 7 and the below-slope Ψ values produced in Figure 6 – this would serve to shift the median result to within the $su(yield)/σ′vo$ 0.20–0.25 range inferred from many flow liquefaction case histories.

The effect of K₀ is specifically examined in Figure 14, where the distance of the phreatic surface from the toe at failure is presented for all models in which η_IL(TX-C) = 1.0 is assumed (models 1–8; see Table 5). This value of η_IL(TX-C) was selected as it is likely the most realistic when considering the below-slope shear densification outlined in the NorSand model and the resultant range of Ψ values in the area most likely to trigger. Other than K₀, the potential effect of α on η_IL and the potential for b to scale η_IL, each of these models had identical inputs and methods. Included for comparison in Figure 14 are the SSR results assuming $su(yield)/σ′vo$ values of 0.20 and 0.25.

The results demonstrate the significant effect that different K₀ values can produce. For example, for the model sets where b scaling is not considered, the predicted phreatic surface distance from toe at failure can vary from 8 m to 23 m (when not considering α) and from 27 m to 45 m (when considering α). Furthermore, as in the NorSand models, the significant effect of considering α in such an analysis is clear. It is emphasised again that there does not appear to be consensus as to either of these issues in current geotechnical practice.

Adding the scaling of η_IL with b in the elastic model framework produces somewhat more nuanced outcomes. The effect of K₀ on initiation of failure is still present, but of a lesser magnitude, between K₀ = 0.5 and K₀ = 0.7. However, for K₀ from 0.7 to 0.8, the effect of K₀ (or, more fundamentally, Poisson's ratio) reverses slightly, making the overall variation in the prediction of failure less than when b scaling is not included. This outcome is a result of the combined effect of Poisson's ratio on K₀ and b (refer to Figure 11). For the models with a high Poisson's ratio, while this increases K₀, the significant effect this also has on b results in a ‘worsening’ of conditions in the context of comparing in situ η and the b-scaled value of η_IL. While the authors argue, based on previous discussions in this paper, that the evidence for incorporating b scaling of this kind is not definitive, the fact that such scaling would likely reduce some of the uncertainty in this process would be a fortuitous outcome were it clearly shown to be the case.

The inclusion of the SSR analyses with $su(yield)/σ′vo$ values of 0.20 and 0.25 to the comparison may suggest that inclusion of both b and α scaling of η_IL appears to result in excessively conservative predictions of failure for any value of BLG K₀. Similarly, as previously discussed, the tendency for most models with K₀ = 0.5 to produce excessively conservative estimates of failure conditions compared with the SSR models (regardless of b or α scaling assumptions) could be taken as being consistent with previous arguments (e.g. Jefferies and Been, 2015) against in situ values of K₀ as low as 0.5.

Despite the previous commentary, it must be acknowledged that the results of the models outlined here do not in themselves provide direct evidence as to whether η_IL should, in practice, be scaled based on b and α in such an analysis: given the interaction of the uncertainty related to the effects of α and b, along with that of K₀, no firm conclusion can be drawn. Indeed, the authors speculate that different conclusions could easily be drawn by different readers based on the results presented in Figure 14. It is therefore suggested that the following additional experimental and field data would be of significant benefit to engineering practice.

• Laboratory studies on the effect of α on η_IL under drained SLT stress paths such as CSD, where the consolidation stress path is carried out in a similar manner to that performed by Sivathayalan and Vaid (2002) such that drained rotation of principal stresses occurs prior to the subsequent CSD process.

• Further plane strain test programmes to indicate whether the effect of b on CSL elevation as observed by Wanatowski and Chu (2007) is seen in other soils, particularly tailings.

• Further attempts to carry out in situ stress measurements of tailings in order to provide better understanding of K₀ and b near to and below slopes, for example through increased use of self-boring pressuremeters.

In the meantime, it is also suggested that practitioners take cognisance of the implicit assumptions regarding K₀ that are, in essence, made when applying different constitutive models to problems such as that outlined in this paper.

To investigate the effects of K₀, b and α on the outcomes of SLT assessments, a series of numerical analyses was carried out, applying both the state-of-practice NorSand constitutive model and, through manipulation of the Mohr–Coulomb constitutive model, analyses to target specific K₀ values for comparison and parametric analysis. A CSD path on an idealised TSF geometry was used as the SLT condition in the numerical analyses.

These analyses demonstrated the significant variation in the prediction of SLT and resulting slope failure based on in situ K₀ that was assumed or that resulted from application of a particular constitutive model and input parameters. Examination of the results of different models against SSR screening analyses may agree with previous suggestions that in situ K₀ values in some cases could develop to levels higher than 0.5, consistent with some previous arguments on this topic. Finally, firm conclusions as to whether η_IL should be scaled on the basis of in situ b and α values in such screening analyses could not be drawn, and thus the need for additional experimental programmes to explore drained triggering and the effect of intermediate principal stress conditions on CSL elevation were outlined.

This work forms part of TAILLIQ (Tailings Liquefaction), which is an Australian Research Council (ARC) Linkage Project supported by financial and in-kind contributions from Anglo American, BHP, Freeport-McMoRan, Newmont, Rio Tinto, and Teck. The TAILLIQ project is being carried out at The University of New South Wales (UNSW), The University of South Australia, The University of Western Australia (lead university), and The University of Wollongong. The authors also thank Adrian Wightman for providing useful feedback to an early version of this paper.

To provide an illustration of the process of shear densification as predicted by NorSand, three element-test simulations were carried out using the same Cadia parameters as listed in Table 2. Each simulation commenced from an isotropic effective stress of 10 kPa and was taken to a final vertical effective stress of 1000 kPa, as summarised in Figure 15. The consolidation process of the three different simulations was: (a) one-dimensional (1D) compression, with no lateral strains, (b) isotropic compression and (c) anisotropic consolidation ramping from the initial isotropic state to a final K₀ value of 0.5. The anisotropic model with an intentionally decreased K₀ value (and thus increased stress ratio) tended toward the CSL, with a final state parameter of +0.035, like the below-slope behaviour shown previously in Figure 5. Alternatively, the 1D compression case exhibited behaviour more similar to that of the isotropic case, with the final K₀ of the 1D case being 0.84. This relatively high K₀ thus resulted in quite low stress ratios, again consistent with the observations in Figure 5. The propensity for NorSand to produce relatively high values of K₀ is discussed in more detail by Castonguay and Konrad (2016).

For the purpose of the NorSand models, a value of post-triggered strength (i.e. liquefied strength ratio) was required to assign to elements where SLT was detected. A relationship between the characteristic Ψ and $sr/σ′vo$ was developed based on the liquefaction case history data collected and analysed by Jefferies and Been (2015), as presented in Figure 16. Jefferies and Been (2015) developed Ψ values for each site, representing the estimated 80th percentile value.

It is important to note that the relationship used, consistent with the case history data, indicates a liquefied strength ratio $sr/σ′vo$ of approximately 0.11 for materials where Ψ = 0 (i.e. on the CSL). This contrasts with the typical experimental behaviour of TX-C samples from Ψ = 0, which would be highly unlikely to provide such a liquefied strength ratio (or, indeed, post-peak brittleness). The reasons for this disconnect between element-scale and field-scale behaviour remains unclear, although anisotropy and shear-localisation at the quasi-steady state are leading candidates. Regardless of the cause, the trend here is informed by field-scale behaviour as per Figure 16 rather than element tests or simulations.

For the parametric Mohr–Coulomb models, where Ψ is not relevant to the analysis, $sr/σ′vo$ = 0.10 was assumed, which is similar to the trend for Ψ ∼ 0.03 in Figure 16, consistent with the values seen in the triggered zone of the NorSand models.

The pseudo code for the NorSand models is shown in Figure 17 and the pseudo code for the Mohr–Coulomb models is shown in Figure 18.

Stability analyses using the SSR technique to find the phreatic conditions that produce a FoS of unity for the idealised TSF were carried out as part of this work. This enables the current work to be referenced more easily to state-of-practice yield stress ratio methods that follow such a general approach. The following sources were used to provide the 0.20–0.25 range of $su(yield)/σ′vo$ values used for the SSR analyses.

• A value of 0.20 adopted by the Cadia panel analyses (Jefferies et al., 2019).

• The typical range of back-analysed values for flow liquefaction failures (e.g. Olson and Stark, 2003).

• The comparison outlined in Figure 19 between the values of η_IL and $su(yield)/σ′vo$ developed through NorSand TX-C element simulations and the experimental data produced by the Cadia panel.

Finally, although the data presented here seem to support $su(yield)/σ′vo$ = 0.20–0.25, a value of 0.30 has been included in some comparisons, owing to the range observed by Olson and Stark (2003).

b

intermediate principal stress ratio (⁠$=(σ2−σ3)(σ1−σ3)$⁠)

c′

effective cohesion

CIU

isotropically consolidated undrained

CAU

anisotropically consolidated undrained

CK₀U

K₀ consolidated undrained

e_cs

void ratio on critical state

K

elastic bulk modulus

K_c

effective principal stress ratio (⁠$=σ′3σ′1$⁠)

K₀

geostatic stress ratio (⁠$=σ′hσ′v$⁠)

G

elastic shear modulus

H

plastic hardening

M_tc

critical state friction ratio, triaxial compression conditions

M_b

critical state friction ratio for the relevant value of b

N

volumetric coupling parameter

p′

mean effective stress ( = $σ′1+σ′2+σ′3/3$⁠)

q

deviator stress ( = $12[(σ1−σ2)2+(σ2−σ3)2+(σ3−σ1)2]12$⁠)

s_r

liquefied shear strength

s_u(yield)

yield stress shear strength

α

principal stress angle

η

stress ratio ( = q/p′)

η_IL

instability stress ratio ( = q/p′ at instability)

η_IL(TX-C)

instability stress ratio, triaxial compression conditions

ν

Poisson's ratio

σ_h

total horizontal stress

$σ′h$

effective horizontal stress

σ_v

total vertical stress

$σ′v$

effective vertical stress

$σ′vo$

effective vertical overburden stress

ϕ′

effective friction angle

χ

state dilatancy constant

Ψ

state parameter