Validation of computational liquefaction for tailings: Tar Island slump Open Access

0

initial condition

1, 2, 3

principal directions of stress or strain

c

critical state

tc

triaxial compression condition (θ = π/6)

e

elastic

p

plastic

p′

mean effective stress (⁠$=σ¯m$⁠) [FL⁻²]

q

triaxial deviator stress q = σ₁ − σ₃ (⁠$=σ¯q$⁠) [FL⁻²]

η

dimensionless distortional stress measure $η=σ¯q/σ¯m$ [−]

θ

Lode angle, $sin(3θ)=−13.5σ¯1σ¯2σ¯3/σ¯q3$ [Rad]

σ _{1, 2, 3}

principal stresses [FL⁻²]

$σ¯m$

mean effective stress $σ¯m=(σ¯1+σ¯2+σ¯3)/3$ [FL⁻²]

$σ¯q$

deviatoric stress invariant $σ¯q=[12(σ1−σ2)2+12(σ2−σ3)2+12(σ3−σ1)2]12$ [FL⁻²]

D

dilatancy, as ratio of strain rates dϵ_v/dϵ_q [−]

D ^p

plastic dilatancy, as ratio of strain rates dϵ_v^p/dϵ_q^p [−]

ϵ _{1, 2, 3}

principal strains (assumed coaxial with principal stresses) [−]

ϵ _q

distortional strain measure work conjugate with $σ¯q$ [−] $ϵq=13[(sinθ+3cosθ)ϵ1+2sinθϵ2+(sinθ−3cosθ)ϵ3]$

ϵ _v

volumetric strain ϵ_v = ϵ₁ + ϵ₂ + ϵ₃ [−]

e

void ratio [−]

ψ

state parameter, ψ = e − e_c [−]

M

critical friction ratio, equals η_c at the critical state; varies with Lode angle, value at triaxial compression (M_tc) taken as reference soil property [−]

N

Nova’s volumetric coupling coefficient in stress-dilatancy [−]

X

scaling factor (soil property) for state dilatancy [−]

Γ

reference void ratio on critical state line (CSL), defined at p′ = 1 kPa [−]

λ

slope of CSL in e−ln(σ_m) space for semi-log idealisation [−]

λ ₁₀

slope of CSL, but defined on base 10 logarithms (2.3 λ) [−]

Construction of a liquefaction-resistant dam at Franklin Falls, NH, USA, by the US Army Corps of Engineers in 1935 led to Casagrande’s canonical idealisation relating soil behaviour to void ratio, with the critical void ratio distinguishing between dilatant and contractive soil. In the subsequent decades, the critical void ratio developed from a simple idea guiding construction to a central tenet in the quantitative representation of soil using the mathematical theory of plasticity. In the last 5 years, critical state theory has gained prominence amongst practical engineers because of its role in understanding three large dam failures: Fundao in Brazil (failed 2015), Cadia in Australia (failed 2018) and Brumadinho in Brazil (failed 2019); critical state theory proved essential in understanding the sudden transition from drained and apparently stable conditions to an undrained liquefaction-driven slump (a transition that took about 5 s at Brumadinho).

Analyses of past failures (‘case histories’) have been constrained by both lack of data on the liquefying soil and the complexity of the failures. In four well-known instances of liquefaction at Fort Peck (MT, USA), Nerlerk (in the Canadian Beaufort Sea), Los Frailes (Spain), and Cadia, liquefaction was a consequence of a preceding foundation yielding/failure and the interplay between foundation movement, and subsequent liquefaction somewhat obscures clarity as validation cases. Similarly, Fundao was associated with complex yielding beneath the upstream raise. In the case of Aberfan in Wales, there are no data on the pre-failure conditions and failure was triggered by rising groundwater (a very different stress path).

One liquefaction case history stands out as a validation case: a static liquefaction slump at the Tar Island Dyke in Canada on 23 August 1974. The slump was a modest scale event during construction of an upstream raise to a hydraulic fill dam and was a construction inconvenience rather than a dam-stability issue; but the loading and geometry were simple: Plewes et al. (1989) describe the key aspects with further information from some original records of the time (Hardy, 1975). There is laboratory data on soil properties, in situ testing using seismic cone penetration test (CPT) and calibration of the CPT in a chamber using similar sand; overall, a possibly unique set of data. The development of this slump was analysed using commercial software incorporating the NorSand elasto-plastic critical-state model (Plaxis 2D; Bentley, 2021); a remarkable simulation of the slump is obtained using all information at face value.

The Athabasca oil (tar) sands, a vast resource comparable to the reserves of Saudi Arabia, are a mixture of crude bitumen within sand that forms part of the McMurray Formation in Alberta, Canada (Figure 1). Commercial development started in 1967 by strip mining, with subsequent extraction of the bitumen from the mined sand using hot water, a process that results in large quantities of sand, process water and residual bitumen: tailings (although these are not crushed rock often associated with that term).

The Tar Island Dyke is a tailing impoundment structure constructed for the Suncor oilsand mine. An aerial view of the mine and the tailing impoundment in the 1970s is shown in Figure 2 with the dyke being used to close-off an oxbow of the Athabasca River. Dyke construction started in the late 1960s with the associated tailing pond being decommissioned in 2006 (and now essentially reclaimed into a forest and small wetland), with Figure 3 illustrating the dyke in its final state at some 100 m high.

Oilsand tailings are predominantly fine sand and also provide the fill material for much of the dyke. Originally planned as a conventional (‘centreline’) compacted fill structure 23 m high, initial oil production found that the clayey material in the ore became modified during bitumen extraction such that the fine sand did not settle out as expected, thus an unanticipated need for tailing storage. The storage requirement made the dyke necessarily four times higher than initially envisaged and its construction changed to using compacted tailing sands placed in the modified upstream arrangement shown in Figure 4.

This modified upstream design, and construction procedures, resulted in three distinct zones within the dyke and adjacent tailings: compacted fill; hydraulically placed tailing sand deposited on beaches above the pond surface (‘beach above water’: BAW) and hydraulically placed sands deposited below the pond surface (‘beach below water’: BBW). Figure 5 shows an aerial view of the Tar Island Dyke from that era with a second photo illustrating the construction method.

The compacted fill used tailings that were sluiced into cells, 500–700 m long by 30–90 m wide, parallel to the dyke axis; the cells were formed using dozer-push sand bunds 2–3 m high. The tailings were discharged at one end of the cell with a spillway box at the far end in turn discharging water and unsettled solids to the pond. The sand within the cell was compacted by dozer tracking.

Beach sands were discharged from the upstream dyke crest, with the length of the above water beach depending on the pond level. Figure 5 illustrates the beach discharge, with variable accumulation of the beach apparent. The finer-sized tailings not dropping out of suspension on the beach travelled into the pond, forming an underwater delta-like front: BBW.

Dyke instability was experienced soon after upstream construction started (Mittal and Hardy, 1977), with four slumps into the impoundment as construction built out on recently placed tailings. The largest of the four slumps occurred on 23 August 1974 and is illustrated in Figure 6.

This largest slump was triggered by placing 12.8 m of compacted sand that stepped forward 36 m over a 9-m thick layer of BBW sand. The fill was placed in five lifts over the period 14 May to 22 August 1974, giving the average 0.1 m/day rate of rise shown in Figure 7; there was a clear acceleration of placement rate with lift 5 although a comparable placement rate had also been achieved for a short time earlier in mid-June. The front slope of the compacted raise sand was 1V:1H, a surprisingly steep slope in itself but real as illustrated in Figure 8. The retained pond level was at the toe of the compacted slope at the time of the slide.

The first four lifts apparently produced no response (e.g. no reported cracks along the dyke crest). Anecdotally, the slump initiated near completion of the last lift with a dozer working in the cell accelerating to its maximum speed to escape the developing subsidence.

The subsidence event comprised the whole 36 m (120 ft) width of the compacted cell dropping between 4.9 and 6.7 m with the compacted fill tilting slightly (as illustrated in Figure 6), but otherwise settling as a rigid block with minimal displacement upstream. Subsidence extended for some 450 m along the dike crest (almost the full length of the cell), so whether viewed from slope height or width of compacted fill, the aspect ratio of the slump was 13:1 or greater – essentially, a plane-strain event. Remarkably, soil movement extended about 250 m upstream from the toe of the raise, as illustrated in Figure 6, which is much more than might be expected with a ‘bearing capacity’ failure.

Liquefaction is assumed to have caused the slump because it arose suddenly during construction over loose tailing sands and with about 40% loss of undrained strength from peak (geometry at initiation of collapse) to residual (the post-slump geometry) and the slump developed rapidly. There were no piezometers in the slope to document pore pressures during the incident, as this was routine construction rather than being a test fill.

There was no investigation of the slide as such, say using CPT soundings. However, Plewes et al. (1989) report an envelope of steady-state strengths from triaxial tests. It is understood that the upper part of the envelope is characteristic of 5% fines while the lower part is representative of a much greater 30% fines.

An unusual aspect of the engineering at that time was the use of downhole density logging. This was done at a different location, but nevertheless provides unusual data on void ratio and with clear distinction between BAW and BBW tailings.

Many CPT soundings have been carried out in the impoundments since 1974 as part of routine engineering for ongoing dyke raising. Importantly, the construction methods have remained the same in the subsequent 25 years as the recommended practice described by Plewes et al. (1989). The tailings have remained from the same source. So, the recent CPT data reasonably represents the tailings and can be used to assess layering within BAW and BBW sands as well as to infer their respective in situ states.

Some of the CPT soundings have included geophones (‘seismic’ CPT) and were used to measure the elastic shear wave velocity, and thus the in situ elastic shear modulus was known.

Suncor is one of several mines in the area, and much relevant testing has been for these adjacent sites with their similar tailings derived from the same McMurray Formation. There are two relevant sets of data: Mildred Lake and Albian.

The Mildred Lake data comprises BAW sand that was recovered as a ten-tonne bulk sample for calibrating the CPT in a chamber and which, as well as providing a tailing-specific mapping of void ratio to CPT resistance, also provides triaxial tests, drained and undrained, over a range of void ratios and confining stress.

The Albian data comprises BAW sand that was used to assess tailing properties in the context of potential use of explosive compaction. There were both drained and undrained triaxial tests, over a range of void ratios and confining stress.

With the exception of the Tar Island triaxial tests and the Mildred Lake undrained triaxial tests, for which only paper records exist, the triaxial data was obtained using computer-controlled equipment and with the data available as digital records. The Mildred Lake triaxial testing used then (late 1980s) state-of-the-art GDS equipment; the Albian triaxial testing used comparable Trautwein equipment (early 2000s).

The Mildred Lake and Albion tests used ‘end-of-test-freezing’, with commensurate accuracy in void ratios to current best-practice standards (see Reid et al., 2020). The Tar Island tests may not have used this method of void ratio measurement.

The available CPT data are digital records and the tests were carried out using modern equipment by an experienced testing company.

The tailings comprise a uniform fine to medium sand with changing silt content depending on location. Plewes et al. (1989) give an overall range of 3–40% (Figure 9), but those fines distribute themselves during deposition with the cell and BAW having about 8% fines (Fear et al., 2014), with more fines being retained in the BBW.

Figure 9 shows the particle size distribution of the sands tested in the laboratory. The sand particles comprise almost entirely quartz and are angular, almost cubical in shape, as illustrated in the inset to Figure 9; the surface of the particles appears rough when viewed with a microscope. The specific gravity lies in the range 2.62–2.64, with G _s = 2.63 a reasonable characteristic value.

Geostatic stress was not measured at Tar Island but has been determined in other oilsand tailings using similar placement methods. Compacted cell and BAW show K ₀ ∼ 1 with BBW at K ₀ ∼ 0.6 (Shuttle et al., 2021).

The critical state locus was inferred from triaxial tests on loose samples. For Mildred Lake tailings, eight undrained triaxial tests were carried out to define the critical state line (CSL) with seven of these eight using a load-controlled protocol. The Albian tailings used the modern protocol of both drained and undrained tests (to more easily determine the CSL at high confining stress) and displacement control. Plewes et al. (1989) published trends in s _r, the shear strength in the critical state, and this has been converted to a CSL using p _c = 2 s _r/M. The inferred CSL are shown in Figure 10 using a semi-log idealisation that is a reasonable characterisation of the data:

The properties for each of the tailing sands are given in Table 1, and the void ratio scale used in Figure 10 approximately shows the accessible range for these sands.

There is limited dependence of the compressibility λ ₁₀ on fines content of the soil, with Figure 11 comparing the data on oil sand tailings with other natural soils. These oil sand tailings are within the range of usual experience with the exception of the 30%-fines data from Suncor which exhibits a λ ₁₀ about a third to a fifth of what might be expected from other soils of comparable fines content (the reason for this low λ ₁₀ appears not to have been investigated in studies at this site). In terms of engineering assessments, there is a case that oil sands tailings appear to display a near common λ ₁₀ ≈ 0.06 across their practically encountered range of in situ gradations.

Overall, there is considerable similarity in the CSL with the differences possibly reflecting changes in gradation between the sands (although the difference for the sands with fines in the range 3–8% is also within cross-laboratory precision (see Reid et al., 2020).

Frictional strength in the critical state, the soil property M, can be estimated from tests on loose samples, drained or undrained, but this is imprecise because the samples usually bulge despite the use of lubricated end plattens – the ‘area’ correction adopted has a surprising effect on the computed stresses and the resulting imprecision becomes significant after 10% axial strain.

The alternative procedure is to derive the strength properties from drained tests on dense samples by plotting peak strength (η_max) against minimum dilatancy (D_min, minimum because of the compression positive convention). This procedure uses data at axial strains of commonly less than 5% axial strain and before any localisation; the calculation of stress from measured load remains quite accurate. There is only data to do this for the Mildred Lake and Albian sands, the data being shown in Figure 12. Both tailings are fitted using the strength–dilatancy equation:

where M is the critical state friction ratio and N is the volumetric coupling coefficient introduced by Nova (1982). The Mildred Lake data is particularly clear with a well-defined trend line corresponding to M_tc = 1.29 and N = 0.45. There are only two points for Albian sand, with the trend line through that corresponding to M_tc = 1.39 and N = 0.33.

Soil strength is controlled by the limiting dilatancy, D_min in Equation 2, and this limiting dilatancy in turn depends on the state parameter through the soil property X:

Again there is no data for Suncor sand but, as shown in Figure 13, there is a very clear trend with the Mildred Lake data and moreover the single test on dense Albian sand lies within the bandwidth of the Mildred Lake data around the inferred trend. The slope of the trend corresponds to X = 4.1; this is a ‘not unsual’ value for predominantly quartz sands with few fines and, as shown below, nicely validates when computing the complete stress–strain behaviour.

There were no bender element tests during Tar Island’s triaxial testing nor were any local strain transducers or the like used (this was testing for practical engineering). Subsequently, a large amount of vertical seismic profiling has been carried out in similar tailings placed using the same procedures, with examples of this data shown in Figure 14.

Elastic modulus in soils is nicely fitted by the idea that compressibility tends to zero as the void ratio tends to its minimum value (a little denser than measured using the ASTM procedure). Further, Hertzian contact theory leads to the expectation that elastic modulus will depend on confining stress. These ideas are captured, for shear modulus G_max, by:

where p_ref is a reference pressure (by convention taken as 100 kPa) for simplicity so that the elastic properties A, e_min and b are all dimensionless and where p and G are in the same units as p_ref. Two trend lines are shown in Figure 14, using Equation 4, with the properties of Erksak sand, a hard quartz sand with a few per cent fines (A = 260, e_min = 0.35, b = 0.5; estimated by minimising the root mean square error when fitting Equation 4 to geophysically measured G_max data and thus comparable to the Tar Island measurements); the void ratio of ‘loose’ and ‘dense’ Erksak trends in Figure 14 are only to illustrate the effect of void ratio on elastic modulus and should not be used to infer void ratios for the in situ tailings.

Comparing the oilsand tailings data with that of conventional natural sand, the compacted cell and BAW show elastic moduli generally consistent with quartz sand albeit with significant scatter (attributed to variable compaction, freeze–thaw cycles and desiccation). The loose BBW, with its higher fines contents, suggests less than the usual increase in modulus with increasing stress.

True elastic moduli are independent of strain level, and thus geophysically measured values should be applicable in any elasto-plastic model. However, the authors weighted the assessed elasticity to the low bound of the geophysical measurements based on their experience that geophysical values form the upper bound to elastic shear moduli measured by unload–reload cycles in high-resolution pressuremeter tests. The concern was that there might be an effect of frequency with the geophysical tests at ∼1 kHz versus the rather slow loading from construction, a concern amplified by the need to use less than geophysially measured G_max to calibrate elasto-plastic models to undrained triaxial tests. The G_max trend used in modelling the BBW is illustrated in Figure 14 and corresponds to G_max = 50 MPa at σ′_v = 100 kPa and with an exponent of nG = 0.5 for the effect of confining stress on the elastic modulus.

Table 1 summarises the soil properties. No constitutive model has been involved, with properties derived from triaxial test data by appropriate plotting. These properties are focussed on ‘strength’ and ‘compressibility’. The missing aspect is the soil’s plastic shear stiffness (typically termed ‘hardening’) for which there is no standard property; this will be a model-specific property when calibrating the numerical approach.

NorSand (NS) is an elasto-plastic model derived from the axioms of critical state theory. It has similarities to the familiar Original Cam Clay (Schofield and Wroth, 1968) but with two crucial differences: (i) the yield surface generally does not intersect the CSL, only moving to the CSL with distortional strain; and (ii) the yield surface hardening is constrained by Equation 3 such that accurate dilatancy is captured regardless of the soil’s void ratio. Any admissible form of CSL may be used, the other properties being M, N, X discussed above together with one NS-specific property: the dimensionless plastic hardening modulus, H. Key ideas in NS are illustrated in Figure 15.

NS has evolved through several developments. The derivation (Jefferies, 1993) from the axioms of critical state soil mechanics (CSSM) was in the context of triaxial compression with N controlling the shape of the yield surface and constant X for all soils (which was based on the tests on sands by Been and Jefferies, 1985). The triaxial derivation was extended to general three-dimensional (3D) conditions (Jefferies and Shuttle, 2002) by introducing the Resende and Martin (1985) strain invariants for work conjugacy, imposing strict honouring of the work-dissipation postulate in 3D, and allowing X to be a soil-specific material property as opposed to a ‘universal’ value. The further evolution (Jefferies and Shuttle, 2011) was based on the idea introduced by Dafalias and co-workers (Li et al., 1999; Manzari and Dafalias, 1997) that M should depend on ψ; this idea has the effect of moving M into an ‘operating’ friction ratio M _i and with yield surfaces now all having constant-shape ‘bullet’ form; how the proportion of intermediate principal stress affects M was also modified at the same time for better computational performance. These modifications simplify NS as well as giving greater accuracy and were made standard from 2011 onward.

The version of NS used for the present analysis is the standard implementation in the Plaxis code, the kernel being given in Table 2. This implementation has two modifications from the NS version documented in Jefferies and Been (2015: appendix C).

First, the derivation of NS applied the axioms of CSSM to the yield surface, which then related the limiting dilation to the ‘image’ condition involving a transformation of Equation 3. A recent theoretical derivation (Jefferies, 2021) indicates that Equation 3 should be used directly as opposed to being mapped onto the yield surface; this change leads to slightly better numerical stability with some property combinations. Compared to prior use of NS, there is a small change in the numerical value of the calibrated H because the internal cap – which forms the ‘target’ of the hardening rule – now evolves slightly more than in the 2015 version of NS.

Second, while Equation 2 is an accurate characterisation of the strength of dense soils, it becomes ambiguous for soils looser than the CSL because D _min is then the end-point rather than a transient ‘peak’ condition. Jefferies and Shuttle (2011) made M _i symmetric around ψ = 0, in essence mapping Figure 11 of Been and Jefferies (1985). However, Equation 2 is a modification of the ‘interlocking’ model of soil behaviour suggested by Taylor (1948) and formalised by Bishop (1950); and Bishop (1950) presented data suggesting (in the current paper’s notation) M _i = M for ψ > 0. As Bishop’s idea fits some loose soils rather well, the NS implementation now provides both idealisations of M _i for the user to choose based on their calibration to the soil being modelled.

NS uses the CSL as best-fit to laboratory data and the properties M, N, X as determined above. Elasticity is likewise determined from data, either in situ or laboratory; the assumption is that modern geophysical methods will be used; thus, NS adopts the elastic shear modulus G _max and Poisson’s ratio. Poisson’s ratio is rarely measured in geotechnical engineering, and a ‘not unreasonable’ ν = 0.2 is commonly used based on Bellotti et al. (1996). This leaves a single NS-specific property, the plastic hardening modulus H

Conceptually, H could be mapped to (say) initial stiffness in drained triaxial compression. A better approach, however, is to use NS to simulate the entire stress–strain behaviour adjusting H for the best-fit: a constrained optimisation. This optimisation is carried out for each drained triaxial test, allowing minor change in the input void ratio during fitting because laboratory accuracy is Δe = ±0.02 even when using the ‘end of test freezing’ method. Figure 16 illustrates examples of fits achieved using constrained optimisation and the properties given in Table 1 and the Dafalias–Bishop option for the zero-dilatancy stress ratio M _i. Constrained optimisation often shows that H depends on ψ ₀; this is illustrated in Figure 17. Thus, the hardening modulus is represented by two properties:

With the NS properties established from drained tests, undrained tests are then simulated. Experience is that the geophysically measured elastic modulus must be reduced to get the best fit to undrained triaxial tests (both stress paths and stress–strain behaviour). Figure 18 illustrates an example fit, in this case to one of the triaxial tests reported by Plewes et al. (1989) from the original investigation of the slump, using a reduced elastic moduli from the geophysically measured trend to obtain the calibration.

CSSM uses void ratio as the kernel variable, the soil properties in themselves not giving stress–strain behaviour (or even strength) without knowing ψ. In the laboratory, the additional information about void ratio comes from its direct measurement. In situ, something else is needed.

It is possible to recover ‘undisturbed’ samples of tailings using thin-wall tubes without swaged ends, but any transport of the samples away from the drillrig causes densification. And, extruding an ‘undisturbed’ sample is impossible with current technology. For these reasons, assessing conditions of tailings within an impoundment must focus on in situ testing and the CPT has become the reference test.

Many CPT soundings have been carried out at Suncor since 2000, with testing campaigns exploring how soil stratigraphy evolves beneath cells and with distance into the impoundment from the dyke axis. It is this data that is used here; this involves assuming that more recent measurements reflect what would have been the case 30 years earlier, but that is not overly onerous as Suncor continued to construct and operate the Tar Island tailing facility in the manner developed during the 1970s and as described by Plewes et al. (1989).

The tailing densities from downhole nuclear logging techniques (Plewes et al. 1989) are summarised in Table 3. The density logging data was only reported at depth. Density logging data are only used for assigning the bulk weights in the analysis, the state parameter being taken from the CPT.

An example of a CPT sounding working from tailing surface in BAW is shown in Figure 19. The water table is just below ground surface, and the sounding is largely drained; where the sounding departs from drained, there is also a large increase in the friction ratio, a situation that would usually be associated with finer grained soils but in this instance, as annotated on the figure, attributed to the clay-bitumen residuals found in these tailings. These residuals appear, surprisingly, dilatant with B _q ∼−0.1.

Recent CPT investigations have been carried out along transects perpendicular to the dyke crest and extending from the cell across the BAW. However, no soundings have been carried out from barges so the data only extends upstream as far as the BAW was trafficable to the CPT equipment. Figure 20 shows data from one such transect, chosen as a plausible characterisation of the tailings involved in the slump shown in Figure 6. There are four soundings, the spacing being shown on the figure, and the distance across these soundings is comparable to the distance slumping propagated (see Figure 6); the sounding denoted ‘101’ on this figure is that of Figure 19.

The boundary is shown between compacted cell and tailings over which the cell was constructed. Thicker instances of clay-bitumen-contaminated layers are annotated and do not appear continuous over distances of the order of 100 m. The water table falls within the body of the dyke because of its engineered drainage.

Stratification within the tailings is more complex than the schematic of Figure 6, with the CPT soundings indicating stochastic variation overlying depositional trends. The spacing of the CPT is such that ‘loose pockets’ cannot be discerned at a horizontal extent of about 80 m but may exist at shorter distances. As indicated by the dashed-line AA on the figure, denser layers do appear to extend over the horizontal distance of the scale of the slump movements; the slope of AA is consistent with the beach slope and is plausibly relic BAW from an earlier stage of the impoundment. Taking this line AA as a reference, the looser ground above is then taken as BBW for the analysis. The characteristic states of these two strata will be assessed after discussing the calibration of the CPT.

The relation between void ratio and CPT resistance is based on calibration chamber testing, where the CPT is pushed into what amounts to a very large triaxial sample of known void ratio and under controlled confining stress (see figure 4.8 in the book by Jefferies and Been (2015) for details of the chamber used to test Mildred Lake tailings). Repeated tests provide a mapping between q _t, e, and p′ (or σ′_v), such mappings being treated as the reference standard within the industry. However, there is no unique mapping from one soil to another because soil properties affect the CPT response as well as void ratio. Been et al. (1987) provided the first systematic approach to evaluating the CPT while allowing for the soil properties giving the basic relation (which is expected from cavity expansion theory):

where k, m are coefficients that depend on soil properties; Been et al. (1987) suggested k, m primarily depend on λ ₁₀.

Shuttle and Jefferies (1998) investigated the effect of soil properties on k, m using a large-strain finite-element analysis of cavity expansion with NS as the constitutive model. This was then extended by Ghafghazi and Shuttle (2008) to show the methodology recovered the worldwide calibration chamber data. Equation 6 was found accurate, but apart from the soil properties λ ₁₀, M, N, X it was found that k, m also depend on the elastic shear rigidity (I _r = G _max/p′). This finite-element model is colloquially known as the CPTwidget (it is a public domain open-source code). The results from this program using the Mildred Lake soil properties are shown in Figure 21.

The Mildred Lake tailings were also used to physically calibrate the CPT in a large chamber (Horsfield and Been, 1987). That physical calibration data is shown in Figure 22 alongside the numerical calibration developed from Figure 21 for the test conditions (e ₀, p′₀) of each physical calibration. The ‘scatter’ seen in the numerical calibration is the effect of I _r. The greater scatter with the physical calibration arises in part because assessing the characteristic tip resistance involves judgment and in part because the measured void radio is not overly accurate as some 2 tonnes of soil was involved in each single calibration test; Ghafghazi and Shuttle (2008) discuss these testing issues in further detail. The CPTwidget is optimised across the set of published chamber test data and thus preferred in developing CPT inversion coefficients.

The CPT data for tailings shown in Figure 20 was processed to compute in situ ψ. This processing used the computed k, m (Figure 21) and was for: (i) for tailings with F <1.3%, established by inspection of the data with greater F implying fines/bitumen contamination making the k, m values invalid (invalid points are not plotted); and (ii) using K ₀ = 0.7 (the theory is for normalisation using mean, not vertical, stress). The G _max used was the low-bound trend to ‘loose’ BBW results discussed earlier. The results of the processing are shown in Figure 23, with three soundings on the same plot to illustrate the characterisation of tailing state within the depositional idealisation. The two upstream CPT soundings were carried out from tailing surface and the third from newly placed cell fill (Figure 20), with the processed data plotted as a common depth below the tailing surface at the toe of the cell. Each point plotted corresponds to a digital ‘scan’ at 25 mm depth intervals for these soundings.

All soil deposits, including ‘uniform’ man-made fills, show variability of the in situ state. This variability, and its pattern, can be characterised and then stability and/or seepage assessments made using stochastic models; notable examples of such an approach in the context of liquefaction are Popescu et al. (1997) and Hicks and Onaspherou (2005). However, data to substantiate a stochastic approach is rarely available, and stochastic modelling is rarely used for design. Rather, it is usual to choose a parameter value (e.g. undrained strength or ψ) from the distribution such that the behaviour of the soil mass with this chosen parameter value is consistent with what would be found from a stochastic analysis; this chosen parameter value is called ‘characteristic’ and denoted by the subscript ‘k’. In the case of liquefaction, the characteristic value appears to be within the lower 80–90 percentile band of the parent distribution (see the book by Jefferies and Been (2006) for extended discussion). Practically, it is usual to select ψ _k from judgment and allowing for possible continuity or extent of loose zones; that judgment approach is followed here as now described. The estimates for ψ _k are shown as annotated lines in Figure 23.

The results from surface to a depth of about 5 m are BAW deposition and suggest layered strata, but for the purpose of analysing the Tar Island slump, that level of detail is not needed as the BAW amounts to a non-liquefying layer that largely just increases the vertical stress on the underlying tailings. A lightly dilatant value of ψ _k = −0.10 is a ‘not unreasonable pick’ for analysing the event.

The BAW is underlain by a transition zone that appears to be an intermingling of BAW and BBW as opposed to a uniform trend of becoming progressively looser. This is modelled as an uncertain thickness of BAW rather than a linear trend in ψ _k.

The transition is underlain by extensive BBW, on this transect to a depth of about 19 m. The median state is about ψ ₅₀ ≈ +0.05. The apparent loose zones from clay-bitumen inclusions have been filtered in the data processing, but some of the very loose points appearing in Figure 20 likely reflect the filtering not eliminating all such instances; thus, the true loose tailings appear denser than about ψ ₁₀₀ < +0.12. The characteristic value must lie between these two limits; by eye, which essentially weighs the data by number of equally valued instances, ψ _k ≈ +0.06 for BBW.

Underlying the BBW is a denser layer about 3 to 5 m thick, below line A–A in Figure 17, and where all three CPT show similar results. This layer is inferred to be earlier BAW that is now submerged, and where the characteristic state has become looser lying in the range −0.02 < ψ _k < +0.00.

The tailings appear looser and more variable below the A–A denser layer. These deeper tailings are from an earlier stage of the impoundment when the dyke crest was lower and thus these locations were further from the spigot points. Instances of reasonably dense soil are seen, but equally there is a greater amount of yet looser soil than the overlying BBW. Such behaviour has been seen at other tailing facilities, where hydraulic segregation develops even though the slurry discharged is narrowly graded; a migration to looser states accompanies the segregation to finer-grained soils. A characteristic ψ _k ≈ +0.065 is reasonable as a starting point, but a more contractive situation also appears possible.

The domain used for the analysis is illustrated in Figure 24 and extends to the underlying natural ground at the base and some 450 m upstream to ensure the boundary did not constrain movement.

The compacted dyke embankment is dense sand and idealised as a drained Mohr–Coulomb soil (orange in Figure 24); in reality, its properties do not matter as it does not participate in the slump (as observed at the time).

The compacted cell is also idealised as a drained Mohr–Coulomb soil (violet in Figure 24), and which was placed on site in five discrete stages (Figure 6). The finite-element model used incremental ‘gravity turn on’ for each layer in turn to replicate the weight of the accumulation of sand during construction. A constraint on the soil properties was the reported 1:1 outer slope of the cell (Figure 8). If the fill is weak, collapse of the cell itself will develop long before the underlying BBW. Thus, the reported geometry constrains the compacted cell to c′ = 5 kPa, φ′ = 45°; a tension cut-off of 10 kPa was also used. The common G = G_max/3 characterisation for shear stiffness in Mohr–Coulomb was adopted, giving a constant G = 33 MPa based in Figure 14.

The crest settlement documented at the time (see the sketch in Figure 6) suggests that a tension crack developed at the back scarp of the slump, but it was not possible to introduce tension cracking within the present finite-element model without constraining the deformations in the tailings below the crack. The influence of a tension crack was approximated by including a thin weak layer between the cell being raised and the dyke itself so that the cell would not ‘hang up’ on the interface; this weak layer is represented by Mohr–Coulomb elements with c′ = 2 kPa, φ′ = 5°.

The cell used to raise the dyke was constructed on a mat. Whilst Plewes et al. (1989) do not give details of the mat, practically this should have been BAW tailings and with the water table at least 0.5 m below surface before allowing construction equipment onto the mat to build the confining ‘bunds’ to retain the tailings to form the cell. A mat thickness of 1 m was assumed. In the absence of other information, the idealisation has taken the denser upper layer seen in the two upstream CPTs of Figure 20 as representative of the mat: thus, taking ψ_k = −0.10 as the starting point for cell construction.

The Plewes et al. (1989) description has recent ‘beach sand overboarded into water’ overlying ‘previously overboarded sand’, Figure 6. That description does not match ground conditions encountered by the CPTs, Figure 20. Because the CPT soundings are consistent with beach development (Figure 5), the site has been idealised primarily looking to the CPT: a discrete, continuous submerged BAW stratum is idealised as 3 m thick and oriented parallel to the surface of the tailings. The continuity of this submerged BAW is open to question, particularly when the tailing deposition illustrated in Figure 3 suggests a delta-like situation of channels, and whose location will change as the spigot locations are changed; this topic of continuity is revisited below when developing scenarios for assessing the event.

Both BAW and BBW were represented using NS; the Mildred Lake properties were chosen. The only difference between the BAW and BBW was ψ_k. The data shows a range of in situ states (Figure 23); scenarios were chosen to investigate the effect of alternative judgement regarding ψ_k as detailed below. The ‘base case’ for scenario development was taken as shown in Figure 24.

There are no self-bored pressuremeter (SBP) data at this site. However, such SBP data does exist at other oilsand tailing facilities and these indicate K₀ ∼0.6 for normally consolidated tailings (Shuttle et al., 2021). The idealisation adopted was to equalise the finite-element model to K₀ = 0.6 prior to cell construction.

Plewes et al. (1989) report a pond level coincident with the cell toe, but that cannot be correct as construction equipment could not operate on the mat with such a water table. It has been assumed that the mat was raised 0.5 m above the pond. Upstream of the dyke, the hydrostatic pore pressure has been idealised below pond level, which is consistent with the CPT soundings. The water table falls within the dyke, because of engineered drainage, and that is found in the CPT soundings; hence, the water table idealisation shown in Figure 24. This is a slight simplification because the water table within the dyke implies horizontal movement of pore water whereas the hydrostatic-from-pond implies no such movement. However, this is a small inconsistency within wider issues as to the continuity and thickness of the submerged BAW stratum.

Plewes et al. (1989) assessed the BBW behaviour in the context of ‘Point A’, which was mid-depth in the upper BBW and central under the footprint of the cell. As illustrated in Figure 25, Point A has been retained and then a further two points, B and C, added at the same elevation. Point B is located in tailings that are expected to experience greater deviatoric stress, or at least at an earlier stage, than A. Point C is added to provide results to help understand how liquefaction propagates. Results from these three locations are shown later.

The discretisation used the standard 15-node triangular element of Plaxis, with the mesh density weighted to where the soils were going to experience the most intense distortion. Figure 25 illustrates the discretisation in the part of the model where liquefaction will develop.

These simulations were carried out using the ‘Dafalias–Bishop’ option for M_i. ‘Soft Flag’ was set for cap softening in the hardening law, this being an accelerator for when excess pore pressure is changing rapidly.

An objective in analysing this case history is to understand static liquefaction within procedures based on plasticity theory. Such liquefaction events result in strength loss, usually considerable; in this case history, the constructed cell settled about 5 m of its constructed 12 m height because of this strength loss. Large strain simulations were needed to track the interplay between crest settlement and evolving strength during the event. Plaxis provided such a mode and it was used for the simulations reported here.

While there were no piezometers to document lack of excess pore pressure in the BBW tailings, it is evident from Plewes et al. (1989) that the dyke’s engineers were aware that the rate of rise needed to be controlled. Further, the recent instance of liquefaction at Brumadinho had involved piezometer arrays which showed no excess pore pressure in the period before the event (Robertson et al., 2019).

Taking a representative drainage path for consolidation of about 15 m, and a hydraulic conductivity using Hazen’s correlation to grain size, a time factor T = 1 (which corresponds to full consolidation) implies a loading period of about a month. This overall cell-raise took place over 3 months; individual steps were placed over about a week and then followed 2 weeks or more of no fill placement.

Steps 1–4 of the cell construction were therefore simulated as developing drained conditions in both the BBW and submerged-BAW tailings. Step 5 was broken into two sub-steps, with the lower half of the raise applied drained. The second half of the raise was applied undrained, testing whether sufficient strength had developed by the tailing consolidation to resist a ‘rapid perturbation’ (and which appears to have been the operation of a bulldozer with associated addition of fill). The nature of the tailing response to this undrained perturbation is discussed below after presenting the results of the analyses.

The base case is the situation where everything is taken at face value, which is the central thread of the present discussion of site conditions. Looking beyond the base case, there are four significant uncertainties within this case history, and which are modelled as alternative scenarios (Table 4).

The recent CPT showed quite complicated stratification which was consistent with the aerial photograph of tailing spigotting (see Figure 5(b)). The base-case scenario blended the CPT results with the description of Plewes et al. (1989), but there was a clear possibility that the submerged BAW had sufficient meanders that it was not a continuous stratum. While a stochastic scheme could be developed, there was minimal data on which to assess the length scale of continuity; even the frequency of denser pockets would be uncertain. Small and unconnected dense pockets will largely not affect the overall response of the BBW strata. Therefore, as an ‘end-member’, the submerged-BAW layer was given the same state as the BBW with all else the same: Scenario 1.

As usual, the in situ state showed stochastic variability (see Figure 23). The stability of the raise was marginal, with the possibility that the assessed ψ _k = +0.06 was too contractive and thus Scenario 2 adopted ψ _k = +0.04, with all else unchanged from the base case. The thickness of the possible dense layer was also uncertain, so Scenario 3 was introduced where this layer was treated as 5 m thick and continuous – in all likelihood, a too-strong situation, but explored to assess if such a feature would cause toe movement at greater distance.

Finally, the authors explored the effect of elastic stiffness. In an elastic-plastic representation of soil behaviour, the elastic bulk modulus partially controlled the post-peak brittleness during the development of liquefaction. The above scenarios used an elastic bulk modulus based on the measured geophysical G _max and a ‘not unreasonable’ Poisson’s ratio (ν = 0.2); this parameter combination normally gives too stiff elastic bulk modulus when calibrating NS to laboratory tests on reconstituted samples. The final scenario was therefore to reduce G _max to a value giving reasonable representations of laboratory behaviour (e.g. Figure 18): Scenario 4.

Results are first presented in terms of strain increments and displacements to illustrate how the slump mechanism develops in each scenario and highlight the similarities and differences between scenarios; these plots also allow judgment of how close a scenario comes to honouring the description of the slump given by Plewes et al. (1989). Following discussion of slumping mechanisms, the results are then presented and discussed in terms of the stress paths and stress–strain behaviour of the tailings at the ‘marker points’ A−C.

The development of liquefaction is illustrated in Figure 26 using contours of distortional strain increment, these strain increments giving the clearest picture of how liquefaction initiates and propagates.

Figure 26(a) shows the incremental deviatoric strain increment contours during the first half of load step 5, with the cell nearly at its ultimate height. The tailings are drained; there is little strain localisation and no obvious collapse mechanism developing.

Figure 26(b) shows the same contour diagram, but now for the situation 10 days later as placement of the last lift concludes and collapse triggers. The load increment between the prior drained situation and insufficient undrained strength was 0.64 m of additional fill.

Figures 26(c) and 26(d) show the contours as liquefaction propagates, initially propagating both above and below the submerged BAW (i.e., the layer A–A in Figure 20) as as the cell slumps (Figure 26(c)), before propagation stabilises at the post-liquefaction configuration (Figure 26(d)); no additional fill is placed, and the change from Figure 26(b) to 26(d) that occurs rapidly is purely caused by loss of tailing strength with the finite-element algorithm following the yielding as the soil loses strength.

Overall, a clear slump within the upper BBW is computed which ‘bottoms out’ on the underlying, somewhat stronger, submerged BAW. However, once strength loss develops, it propagates quickly and the final slump mechanism punches through the BAW to form a log-spiral mode extending upstream some 150 m.

It is difficult to appreciate just how dramatically static liquefaction can develop from Figure 26. The software’s animation function has been used to develop a video of the computed liquefaction, which is provided as online supplementary material, showing how the liquefaction develops and propagates – there are considerable similarities to the propagation of liquefaction during the Brumadinho failure (compare this computed animation with the widely available video record of B1-CAM1-Barragem).

The movements generated during the liquefaction are shown in Figure 27. The cell settles some 3.5 m during the liquefaction, in a manner that most observers in such an event might describe as ‘uniform’, which is what was reported at the time. The tailing movements are some 500 mm at 120 m upstream, also consistent with the reported observations of the time. Overall, this base case computes the reported mechanism with a small bias to lesser settlement of the raise and corresponding reduced extent upstream. The point of instability is nicely captured.

The full post-liquefaction mechanism for the scenario of the submerged BAW being discontinuous, and thus not affecting the collapse mode, is shown in Figure 28. There is considerable similarity between Scenario 0 (Figure 26(d)) and Scenario 1 (Figure 28(a)), which is unsurprising since the collapse mode punches through the continuous BAW in Scenario 0. The absence of the stronger layer makes the collapse mode slightly deeper and extending further upstream in Scenario 1. The crest settlement during liquefaction is now close to that reported, at 4.5 m.

Since there were no reported slips or other warning signs during steps 1–4, the possibility of less contractive BBW was simulated to investigate if the CPT evaluation was biased towards too-contractive. The results of Scenario 2 are illustrated in Figure 29. Compared to Scenario 0, the base case, the less contractive BBW of Scenario 2 has little effect on the collapse mechanism in itself although post-liquefaction displacements are now less with a computed crest settlement of 3.3 m.

Scenario 3 investigates a thicker relic BAW than the base case, a scenario developed from Figure 20 where there is uncertainty in stratification around the A–A line, and these CPTs were measured on similar tailings placed the same way, not the actual ground conditions of the slump. Scenario 3 retained the features and parameters of the base case except that the base of the relic BAW was lowered to give a 5-m thickness for the relic BAW. The results are illustrated in Figure 30.

Compared to base case Scenario 0, the thicker relic BAW of Scenario 3 is now sufficient to limit the mechanism at depth (see Figure 30(b)); but rather than the consequence being a greater distance of toe movement, a much-reduced crest settlement of only 1.1 m develops during the liquefaction.

Scenario 4 retains the features and parameters of Scenario 0 (the base case), in particular, the characteristic ψ_k in the BBW and relic BAW. The central feature of Scenario 4 is to reduce the elastic moduli (both bulk and shear) in the BBW and BAW; such reduction is usually needed to get the best fits to undrained triaxial tests on loose reconstituted samples as shown in Figure 18. The reduction factor used was G = G_max/3, and the results of Scenario 4 are illustrated in Figure 31.

The reduced elastic modulus of Scenario 4 has a dramatic effect on the computed collapse mode, with collapse mechanism becoming a local ‘bearing capacity’ rotational slip rather than a flowslide. The cell does not settle as a whole, with the upstream crest of the cell settling only 1.5 m, which is much less than reality.

The stress paths and stress–strain behaviour of the three locations A–C are shown in Figure 32 for the base case scenario. Considering point B first, this was located at the toe of the raise (Figure 25) in the area where soil liquefaction developed fastest with correspondingly large η/M. The drained response as lifts 1–5a are placed takes the stress state at B above its instability locus such that the small load of lift 5b, applied undrained, causes an immediate loss of strength – a negative modulus to load increment.

The loss of load capacity at B then has to be shared by that load being carried by adjacent soil, which can be seen at point A where there is a small capacity for additional load, with the stress path moving to increasing σ_q, but that in turn then takes A beyond its peak strength and also now into rapid loss of strength. Collapse propagates. Point C is in an area liquefying towards the end of the propagation and shows useful undrained strength before being eventually overloaded by the developing slump.

The nature of the instability locus follows directly from NS hardening and involves no additional features or physics. The NS yield surface is characterised by its mean effective stress at the image condition, p_i. Distortional strain increments initially cause p_i to increase, reflecting the underlying formation of contacts between the soil particles and with the soil developing greater strength as the yield surface expands, but this hardening only goes so far as further yield surface evolution depends on the state parameter once the hardening limit p_mx is reached, and the yield surface may even have to contract. Thus, there is a universally stable soil response, drained or undrained, while p_i < p_mx, but once that condition is reached the response then depends on how p_mx evolves; the evolution of p_mx depends on the soil properties, not just the evolution of state parameter. If p_i > p_mx, the soil is vulnerable to a ‘rapid perturbation’ (small transient load applied quickly, forcing a short-term undrained response) that will produce a negative modulus to the stress–strain behaviour. The issue then is how that load capacity reduction from the negative modulus can be shared by adjacent soil – and with the speed at which that load must be redistributed to adjacent soil corresponding to the speed of propagation of a plastic shear wave for dead load situations, such as here, and hence the almost explosive propagation of liquefaction in the animation (and in reality at Brumadinho).

The mobilised strengths at onset of liquefaction and after liquefaction has fully developed are shown in Figures 33(a) and 33(b), respectively, for Scenario 0, the base case. The mobilisation is plotted as the strength ratio η/M as that ratio allows direct comparison on an evolving situation in plane strain, where the Lode angle changes, to triaxial compression data (which underlies perceptions about soil behaviour); η/M = 1 corresponds to soil in its critical state.

Comparing Figure 33 with the corresponding deviator strain increments, Figure 26, it is apparent that the localisation of the strains – and thus the collapse mechanism – is from the tailings shearing into their critical state. Once that occurs, subsequent movement is driven by tailings flowing until cell settlement reduces the driving forces (in limit equilibriums terms) to match the now reduced tailing strengths. This mechanism is illustrated by the animation provided as online supplementary material, with the relation between crest settlement and propagation of liquefaction being clearly evident.

All three of the points A–C developed their critical state strengths during liquefaction, evident from the stress–strain curves at each location (Figure 32), but there was not a single post-collapse strength ratio because each location experienced differing amounts of ‘shear-induced densification’ during its drained loading as the cell was constructed. So, although starting from the same assessed characteristic ψ_k in the greenfield condition, significant divergence developed in ψ between the three locations.

Taking the mean effective stress p_5a at the end of constructing cell lift 5a as a normalising stress, the post-liquefaction strength ratios realised were s_r/p_5a = 0.101 at A (mid-point beneath cell), a slightly larger s_r/p_5a = 0.143 at B (toe area, with most densification pre-collapse) and s_r/p_5a = 0.042 at C (a reasonable characterisation of tailings still in the greenfield state).

Although there is a tendency in the literature to regard a ‘case history match’ as being proof of how a situation/event developed, that can never be the case: there will always be uncertainties about subsurface conditions and rarely sufficient instrumentation measuring what happened (e.g. data-logged piezometers in the failing soil). These uncertainties manifest themselves in a range of scenarios being plausible matches to observations. Four central uncertainties (Table 4) were explored in the numerical simulations around the central estimate (base case).

A striking result of this exploration of uncertainties is the importance of honouring true elastic moduli (easily measured geophysically in present practice). Simulations involving reduced moduli, which are needed to fit the response seen in laboratory tests on reconstituted samples (e.g. Figure 18), simply do not develop into the observed collapse mode at this site. In contrast, simulations using G _max from the seismic CPT all come close. This raises questions about whether ‘aging’ (e.g. time-dependent creep) leaves reconstituted samples in the laboratory, the standard procedure for testing silts and sands, as unrepresentative of the same soil’s stiffness in situ – wider use of self-bored pressure meter testing in tailings, with its inherent testing of stress–strain behaviour in situ, would seem helpful.

In terms of the submerged-BAW stratum, the range explored was based on credible idealisations of ground conditions revealed by the CPT. These conditions were for the past 30 years, and perhaps the depositional environment was different, but it seems that denser layers in the tailings would have little effect on the collapse mode – there is too much of the BBW and that BBW controls.

The closest match of the computed results to the reported observations from the time was Scenario 1, which treats the denser zones in the tailings as discontinuous and gives a very close match to settlement and upstream extent of movement.

Liquefaction can be triggered through loading sufficiently quickly so that the soil response becomes undrained. If the soil has been stressed above its instability locus that undrained increment will be small – in the simulations presented, a load increment of just 11 kPa triggered the flowslides.

The engineers of the day involved with Tar Island developed an empirical rule, based on four further liquefaction slumps, relating raise height to its rate of construction (Plewes et al., 1989: figure 7). Implicit in such a rule is the concept of drainage. Permeameter test on the Mildred Lake tailings revealed a hydraulic conductivity k ∼ 1 × 10⁻⁶ m/s, which is much less than would be expected from the tailing gradation (Horsfield and Been, 1987). The Mildred Lake samples, despite being ‘clean’ smelled strongly of bitumen, with the reduced conductivity being attributed to residual hydrocarbon traces affecting the available pore space for water movement; this hypothesis is further supported by the end-of-test water contents using the freezing method being clearly offset from what would be expected based on volumetric strain during drained tests. The implication was that even though the upstream raise was extending over sands, there was potential for a much slower dissipation of construction pore pressures than the 1-month expectation discussed earlier. Treating the second half of the step 5 raise as ‘rapid’ is a reasonable perturbation.

The engineers of the day did not find precursor indicators for the flowslide (e.g. excessive settlement or tension cracks). Most recently, the Brumadinho failure had comprehensive piezometric monitoring and surveillance, with all involved reporting no precursor indicators, and that lack of indicators was exactly what was computed at Tar Island.

Recall, that all loading was drained up to and including lift 5a at which time there was 11.4 m of fill. Even if piezometers had been deployed at the time, no excess pore pressures would have been measured: the same as reported at Brumadinho.

In terms of using settlement as an indicator for an observational method approach, Figure 34 shows the computed settlement as would be observed using a settlement indicator placed on the mat at the start of cell raising. Practically, it would be difficult to discern between the differing scenarios for foundation conditions using settlement. Most importantly, no scenario shows accelerating settlement prior to liquefaction; indeed, there is a perceptible stiffening of the settlement caused by additional load once about 7 m of fill is in place – again, exactly aligned with reports of no apparent distress prior to collapse.

Tar Island is arguably the simplest full-scale liquefaction event to be found in the literature and is a good case history to test liquefaction assessment procedures: validation. Using triaxial and CPT data to calibrate soil properties and in situ state at ‘face value’, finite-element modelling using Plaxis/NS (which implements critical state soil mechanics using the state parameter) gives an excellent match to the reported observations about this liquefaction failure. What stands out from the numerical work is the importance of honouring geophysically measured elastic moduli: an aspect that is not discussed or mentioned in any of the current liquefaction literature. The elastic shear modulus influences, arguably controls, the redistribution of load after soil is stressed past its instability locus and it is this redistribution that produced a large difference in the computed flowslides outcomes for this case history. This importance of using ‘G_max’ is possibly the most significant conclusion from this study, albeit very reassuring to also find critical-state-based procedures so accurate.

The CPT data shown were obtained by Conetec working under the direction of Klohn Crippen Berger, and the authors appreciate access to this data. Particular thanks are owed to Ed McRoberts, who took over as Engineer-of-Record for the Tar Island Dyke and provided photographs from the time as well as information on aspects of how the structure was built.