A Two-Fluid Simulation of Tailings Dam Breaching

This paper presents the development and application of a dam breach model, EMBREA-MUD, which is suitable for tailings dams. One of the common failure modes for these structures is breaching due to overtopping, which together with the flow of liquefied tailings, is simulated by the proposed model. The model simultaneously computes the outflow of water and tailings from a tailings storage facility and the corresponding growth of the breach opening. Tailings outflows are represented by a separate non-Newtonian viscous layer, which together with a water layer, represent the two fluid components of the model. The third component represents dam material that can be eroded by the shear forces exerted by either water or mud. The water layer also exerts dynamic and erosional forces and can transport solids eroded from either the mud or dam layer. The model was verified against laboratory cases as well as two field cases reported in the literature, the failures of the Mount Polley tailings dam in Canada in 2015 and the Merriespruit dam in South Africa in 1994. The model results agreed well with the recorded narrative of the events, although in the latter case, careful calibration of one of the model parameters was necessary to obtain a good match.


Introduction
The consequences of a tailings dam failure can be catastrophic. A recent reminder of this is the failure of Dam I at the Córrego de Feijão Mine in Brazil in January 2019 ( Fig. 1), where approximately 270 people died (VALE 2020). In 2015, the Fundão dam failed, resulting in the worst ever environmental disaster in Brazil, with the release of over 30 million m 3 of tailings that polluted 670 km of watercourses on its way to the Atlantic Ocean, where it spread along hundreds of kilometres of the Brazilian coastline (IUCN 2018;Palu and Julien 2019). In 2014, the Mount Polley dam failure in Canada released 25 million m 3 of tailings and supernatant water into the environment (BCMEM 2015). There have been more than 30 tailings dam failures per decade in the period 1960-90 and around 20 per decade since 1990s. Although the number of tailings dam-related accidents has decreased since the 1990s, the number of severe failures (i.e. those that have released more than 100,000 m 3 of tailings and/or resulted in loss of life) has increased (Bowker and Chambers 2015).
An informed risk management strategy can minimize the probability and consequences of tailings dam accidents. An accident presents a considerable economic damage to the mine owner, danger to life and health of staff working at the facility, and potentially devastating consequences to communities and the environment downstream. The starting point for an assessment of the consequences is understanding the failure modes and the volumes of potentially released materials. Tailings dams can fail due to various reasons, such as slope failure, liquefaction of stored tailings, foundation failure, earthquake, and internal erosion. A common failure mode is also breaching due to overtopping, first by the supernatant water and later, as the breach opening grows, liquefied tailings. This paper presents the development of a two-fluid dam breach model, with the first fluid representing water and the second fluid representing tailings. The model outputs include time series of water and tailings outflows, which can be used to assess potential downstream impacts, thus contributing to risk management and emergency planning. The assessment of economic damages and health impacts, including loss of life downstream of the dam, requires other 1 3 models or assessments that are not discussed further in this paper, but are included in a paper by Lumbroso et al. (2020).

Failure Modes and Breach Outflow Classification
There are four primary failure modes of tailings dams: slope failure, foundation failure, internal erosion and surface erosion (James et al. 2017;Liu 2018): • Slope failure (or instability): This most commonly occurs after saturation of the dam embankment due to heavy rain (or snow) or poor surface water (pond) control, which causes a rise in the phreatic surface. Slope failures can also be caused by increased pressure on the dam from liquefaction of stored tailings or excessively high rates of dam rise (Liu 2018). A recent example of this type is the 2018 failure of the Cadia dam in Australia (Jefferies et al. 2018); however, owing to prompt action and the site location, there were no deaths, unlike the disastrous Stava dam slope failure in Italy in 1985, where more than 200 people died in two downstream villages that were flooded (van Niekerk and Viljoen 2005). • Foundation failure: Sudden or excessive loading may cause the foundation to deform if it is not sufficiently strong, which may lead to a local or an overall failure of the dam. An example of this mode is the failure of the Mount Polley tailings dam in Canada in 2014 (BCMEM 2015). • Internal erosion (Seepage and piping): A phreatic surface that is too high may lead to seepage and piping eroding a hole through the embankment that grows in size as the flow progressively erodes the surrounding material.
This can lead to the collapse of the dam above the piping hole and local and general dam failures.  (Lyu et al. 2019) and the Padcal dam in the Philippines in 2012 (AGHAM et al. 2013). Overtopping can also occur after the dam initially partially failed due to another failure mode that locally lowered the crest level, enabling water outflow from the pond. An example of this is the Mount Polley tailings dam failure (BCMEM 2015), where the dam initially experienced a foundation failure. The resulting outflow of more than 10 million m 3 of process water further eroded and enlarged the breach opening, which enabled tailings to flow out as a mudflow. In many cases, this type of failure results from inappropriate water management and/or inadequate beach length (i.e. the distance between the tailings dam and the pond).
In addition to these failure modes, a failure can also occur due to an earthquake (ICOLD 2001;James et al. 2017). Recent examples of this type are the failures of the Las Palmas tailings storage facility after an 8.8 magnitude  (Lyu et al. 2019) and of the Kayakuri tailings dam after the Tohoku Earthquake hit Japan in 2010 with a magnitude of 9 (Ishihara et al. 2015). This type of failure is due to the development of excess pore water pressure in the tailings during an earthquake, ultimately leading to liquefaction and collapse of the tailings dam.
Liquefaction occurs when a soil-like material undergoes continued deformation at a low constant residual stress or with no residual resistance, due to build up and maintenance of high pore-water pressures that reduce the effective confining pressure to very low values (Puri and Kostecki 2013). Tailings can liquefy if earthquake loading creates a breach or by the loss of confinement resulting from breach of the dam by another mechanism. The potential for liquefaction of tailings is a function of several factors, but most importantly, their density and stress state. Thus, the potential for liquefaction varies during the life cycle of the impoundment and generally decreases with time due to consolidation and ageing (James et al. 2017).
Two factors that have been shown to have an important influence on the flows from tailings dam breaches (James et al. 2017): the presence of surface water near the breach and the potential for liquefaction of tailings near the breach. Depending on the status of each factor, there are four possible breach outflow classes: • Class 1A: surface water present near crest and liquefaction of tailings: Dam break. Flow from the breach will consist of water, possibly with eroded tailings, and liquefied tailings • Class 1B: surface water present near crest and no liquefaction of tailings: Dam break. Flow from the breach will consist of water with eroded tailings. • Class 2A: surface water far from crest or no pond, and liquefaction of tailings: Slope failure with tailings released as debris flow or mud flow • Class 2B: surface water far from crest or no pond, and no liquefaction of tailings: Slope failure without outflow, only displaced tailings.
The simplest type of assessment of stored tailings outflow are empirical equations that assess the volume of released tailings and their runout distance (Larrauri and Lall 2018;Rico et al. 2008). They were derived from a database of previous tailings dam failures and use readily available parameters such as the volume of stored tailings and the height of the tailings dam. They are appropriate for use in scoping studies and to identify dams for which more in-depth risk assessments are required.
A higher level of assessment uses event trees. These methods incorporate more information about dam and tailings properties, as well as user judgement, and pass these input data through a failure event tree. Once the failure event and mode is confirmed through the tree, the input data can also be used to produce outflow hydrographs for the water and tailings (O'Brien et al. 2015).
Finally, numerical models have the advantage of using fundamental principles of physics to simulate the release of tailings and water. There are a number of models that simulate both the breaching of the dam and the spreading of the water following a dam break (flow class 1B). The flow of tailings can be simulated in a similar way. Models that simulate the spreading of mudflows (flow class 2A) must consider how tailings, unlike water, behaves as a non-Newtonian viscoplastic material (mudflow). These models are typically used in combination with assumptions of breach parameters and the released volume of tailings (Moon et al. 2019).
A very common case of tailings dam failures includes the flow of both water and mud. Models for the simultaneous simulation of the how both water and mud flow have been developed (Chen et al. 2007), but are less common than water-only or mud-only models. Breach growth, i.e. erosion of the dam, must also be simulated to accurately predict the outflow of water and mud from a tailings storage facility (TSF). This paper presents an attempt to develop such a model, a type of which, to the best of our knowledge, has not been available previously.

The EMBREA Dam Breach Model
The model presented in this paper is an extension of the EMBREA model for simulation of a dam breach. The model was first developed by Mohamed et al. (2002) under the name HR BREACH. A series of European funded collaborative research projects allowed further testing and refinement. These projects included the CADAM, IMPACT, FLOODsite, and FloodProBE projects. In particular, the IMPACT project incorporated large-scale field and laboratory data for model performance validation (Morris et al. 2005), with five field and 22 laboratory tests performed. EMBREA was also tested against three historical dam failure cases: Oros in Brazil, Banqiao in China (Morris 2011), and Teton in USA (Mohamed et al. 2002).
EMBREA predicts the growth of a breach in embankment dams due to erosion caused by outflowing water and the quantity of water released in the form of a hydrograph (Mohamed et al. 2002). Both processes, breaching and outflowing, are simulated simultaneously based on the characteristics of the dam material, predicting the evolution of the breach opening without the need to make assumptions regarding the dimensions of the breach. EMBREA is a onedimensional and a one-fluid (water) model. In the next section, we present the development of its extension to two-fluid modelling (water and mud). The developed model is called EMBREA-MUD.

Model Components and Interactions Between Them
The model describes dam breach dynamics with the following layers: • Fluid 2: Water, a fluid with Newtonian behaviour. This layer represents water but also includes suspended solids (eroded tailings and dam material). • Fluid 1: Mud, a fluid with a visco-plastic non-Newtonian behaviour. This layer represents liquefied tailings including eroded dam material. • Solid (layer 0): Dam material (can only be eroded).
Fluid layers are characterised by its depth (h) and flow velocity (u). The interaction between layers that occurs as a result of shearing due to different velocities is described through shear stress (τ), as shown in Fig. 2. To achieve this interaction, the values of variables are passed between each layer at the end of each computational time step.
Two shear forces between the layers are simulated by the model. The first is the shear stress τ 1 between the water and mud layers, or between the water and dam material if no mud is present. This is a dynamic force acting on water, as well as mud if mud is present; and an eroding force for the mud layer, or dam material if no mud is present. The second is the shear stress τ 0 between the mud and dam layers. This is a dynamic force acting on mud and an eroding force for dam material.
Shear stress τ 1 is related to water (Newtonian fluid) and is computed from the Manning equation for flow resistance: where γ 2 is the specific weight of Fluid 2 (water), n is the Manning coefficient, R 2 is the hydraulic radius of water layer, and u 1 and u 2 are the flow velocities of Fluids 1 (mud) and 2 (water), respectively. Shear stress τ 0 is related to mud (non-Newtonian fluid) and is computed using to the Herschel-Bulkley fluid model (Chen et al. 2007): where K is the viscosity coefficient and m is the flow index, a measure of the degree to which the fluid is shear-thinning (m < 1) or shear-thickening (m > 1). The well-known Bingham model is a special case of the Herschel-Bulkley model where m = 1. The symbol R' = R 1 (1 − τ y / τ 0 ) and R 1 is the hydraulic radius of the mud layer. As in EMBREA (Mohamed et al. 2002), the quasi-steady approach is used to solve the dynamic equations for flow depth h 2 and velocity u 2 for water layer. On the other hand, fully unsteady flow equations are solved to compute mud depth (h 1 ) and velocity (u 1 ), according to the LHLL numerical scheme described in detail in Chen et al. (2007).

Spatial Discretisation
EMBREA is a single point (0-D) model for the water storage behind the dam and a one dimensional (1-D) model for the dam area, where computational points are represented by dam cross sections. The assumption of one water level representative for the whole storage is justifiable for water, but not for mud, which is a visco-plastic material. Therefore, in EMBREA-MUD, the TSF is also discretised as a set of 1-D computational points. Many TSFs can be well approximated as having a rectangular shape. Therefore, in the present version of the model, the tailings storage facility is represented as a rectangle with a width W and length L.

Erosion and Breach Growth
The erosion rate E i for the layer i is determined by its erodibility coefficient K D,i and critical shear stress τ c,i . Index i corresponds to either of the two erodible layers, i.e. i = 1 for mud layer and i = 0 for dam. If the shear stress is less than the critical value, there is no erosion (E i = 0). When the acting shear stress is higher than the critical value, the erosion rate E i is computed from the equation: In computational cross sections where mud is not present but there is water, the shear stress used to compute dam erosion (τ 0 ) equals the shear stress exerted by the water, τ 1 . After the erosion rate is calculated with Eq. 3, it is also taken into account that this value is a nominal erosion rate for a cross section, while its local value along the wetted perimeter varies. This variation, as well as modifications of the cross section resulting from a varying erosion rate, are done using one of the options in EMBREA (Mohamed et al. 2002), where variation of shear stress along banks is proportional to depth and sidewall collapse is assumed to occur at the same instant as any scour at the toe (Fig. 3).
In EMBREA, all the eroded material enters and is transported by the water column. In EMBREA-MUD, the fate of the eroded material is simulated differently, depending on where it is eroded form: • Bank material eroded from below the top of the mud layer is added to the mud layer and assumed to behave as mud once it is eroded, • Bank material eroded from above the top of the mud layer is added to the water layer and removed from the system in suspension, and • Material eroded from the bed is added to the mud layer if there is mud in the present cross section; otherwise, it is added to the water layer

Model Input Data and Parameters
The model accounts for a number of user-specified parameters and initial conditions for the dam, TSF, tailings, and water. In addition to those, a boundary condition of water inflow into the facility can be optionally specified. The dam is assumed to have a trapezoidal section initially: the user must specify the dam height, crest width, and length and embankment slope on the dam's downstream and upstream side. In the case of multiple embankment sections, the user should select one for which the failure is simulated. The geometric properties of the TSF are described by stage-area or stage-volume curve. The initial conditions are the level of water and tailings. The water level is the same throughout the facility (where this level is higher than the level of tailings; otherwise, the water depth is zero). The tailings surface can vary along the facility at a user-prescribed slope, different above and under the water surface. There are also material and flow parameters used throughout the simulation. There are a number of dam-related parameters, as in the original EMBREA model: most importantly, the dam erodibility coefficient K D,0 , and the critical shear stress τ c,0 . Further details are given in Mohamed et al. (2002) and Morris (2011). For tailings, in addition to the erodibility coefficient K D,1 and the critical shear stress τ c,0 for erosion, there are three parameters that describe mud flow: K, the viscosity coefficient; τ y , the yield stress; and m, the flow index. For water flow, the Manning n roughness coefficient must be specified, per Eq. 1.
The initial failure that may occur due to causes other than erosion may be described by the depth and width of an initial gap in the dam. Otherwise, a small gap from where the erosion will be initiated once the water level reaches its lowest point can be prescribed (for example, this could be 30 cm deep and 3 m wide). The user can also prescribe the maximum erosion depth and width, if physical restrictions to dam erosion exist.

Results and Discussion
The model was validated against laboratory experiments and compared against observations made during failure of two actual tailings dams.

Validation of EMBREA-MUD with Experimental Laboratory Results
Dam-break type experiments found in the literature usually consider only one (viscous) layer. For our purpose of testing the numerical behaviour of the newly added mud layer, this was sufficient. The water layer and dam erosion had already been tested during the original developments of the EMBREA model (Mohamed et al. 2002;Morris et al. 2005), which did not include a mud layer.
A series of dam break experiments with a viscous fluid were conducted by Jeyapalan et al. (1983) in a 0.305 m wide flume. The dam failure was simulated by a quick removal of a barrier representing a dam. Behind the barrier, there was a tank filled with a viscous fluid. In the tests, oil was used as a viscous Bingham fluid with the following properties:  density ρ 1 = 900 kg/m 3 , viscosity K = 3.9 Ns/m 2 , and yield stress τ y ~ 0 (0.01 N/m 2 ). Here we present the numerical results of the model against the observed longitudinal profile obtained for Test 2. In this test, the tank was 0.61 m long and the depth of fluid in the tank was 0.152 m initially. A comparison between the numerical and observed profiles is shown in Fig. 4. The match between the two is satisfactory, indicating that the newly added mud component is capable of simulating a dam break scenarios for viscous fluids.

Simulation of the Mount Polley Tailings Dam Failure
The  (2015), is shown in Fig. 6.
Following the foundation failure, water and tailings started flowing out of the TSF. The chronology of events relevant for assessment of the breaching event is given in Table 1.
The volume of outflowing water was significant, with two channels scoured by water (Fig. 4): a smaller and shorter channel scoured in the right hand side (looking in the direction of outflow), and a bigger and longer one stretched across most of the length of the pond and turning towards to the further corner (about 2 km long). It appears that the movement of mud was mainly limited to these two channels. For modelling purposes, it was important to predict the width of the channel(s). The median value of the crest length after the collapse (360 m) and two estimates based on an empirical relation for flushing channels was used. Drawdown flushing is a controlled measure for removing sediment from a water storage reservoir. White (2000) proposed the following equation for the width of these channels: where W fc is the channel width in [m] and Q f is the flushing discharge given in [m 3 /s]. The first width estimate was based on the overtopping discharge (630 m 3 /s), estimated from the head above the crest and the width after the foundation failure, resulting in a width of 320 m. The second estimate was based on the average flow through the breach (about 1200 m 3 /s, see results in Fig. 7) once it was scoured, which produced a width of 440 m. The median value of the crest length and the two predictions was used for modelling (360 m).
The main parameters considered in the simulation of Mont Polley failure are given in Table 2. The model estimated an outflow tailings volume as mudflow of 15.06 million m 3 and an outflow in suspension of 510,000 m 3 , totalling 15.57 million m 3 . This is slightly above the observed outflow of tailings and interstitial water (estimated between 11 and 15 million m 3 ); however, it includes the eroded dam material as well. Hydrographs of outflows from the dam are shown in Fig. 7.
For the comparison of the model results to the narrative of the observations in Table 1, it was assumed that the failure occurred instantaneously at around midnight, though to 100 m) is probably a result of the model's assumption that the sides were stable until vertical, which is likely too steep. With respect to the outflow of water and mud, the report mentioned that the pond was empty after five h, but that a high outflow of muddy water continued. If the pond had been completely free from surface water, there should have been no water outflow, only mudflow. This description could therefore be interpreted as water having been limited to the channels scoured in the tailings, while the rest of the pond area was dry. The EMBREA-MUD model predicted that water became limited to channels at 4:48 from the start of simulation, while water release ended at 5:04 (the pond was empty of water at this time). Mud outflow was at its highest (between 1500 and 2580 After a pump at the Perimeter Seepage Collection Pond (PSCP) was turned on earlier (a normal procedure) to decrease water level in the PSCP, the water level at this time was no longer decreasing, indicating external uncontrolled inflow (equal to the pump capacity) to the pond, a possible indicator that overflow had already occurred by this time 4 Aug. 2014 0:30 to 1:00 Mine personnel observed the water level in TSF is dropping. They drove to the corner near to the breach location and heard "roaring like a 50 feet waterfall" 4 Aug.  The simulated outflow of mud fell to about 130 m 3 /s after 6.5 h. This fits the observation of a much smaller outflow, compared to the peak outflow rate, which was 20 times higher. The predicted mud outflow then decreased further; between 10 and 14 h, it ranged between a few hundred L/s and a few m 3 /s. The outflow dropped to 5 L/s after 14 h and to 2 L/s by 16 h. The observation record mentions that at 12:00 (roughly 12 h after the breach), there was still some outflow, which had ceased by 16:00.
The observed final breach width at the narrowest section was 40 m (BCMEM 2015); at an average cross section, it was 92 m. The simulated width ranged between 56 and 76 m with an average of 65 m. A possible explanation for these differences is that the model did not account for variability in the tailings dam's properties, using an average erodibility coefficient instead. A general view of the final breach geometry is shown in Fig. 9.
While some differences between the observed and predicted dimensions can be noted, the overall volumes and chronology of the event were reproduced reasonably well by EMBREA-MUD. No particular calibration of model parameters was undertaken in this case.   (Hanson and Hunt 2007) Embankment critical shear stress 10 N/m 2 Estimated as 10 times mud critical stress. In a dam breach event, the shear stress is much higher than the critical shear stress, and the value of the latter is not of significant importance Clear water started to flow through the village as early as 19:00, indicating that overtopping started 2 h before the dam failed. The primary cause of the failure of this tailings dam was overtopping. The probable sequence of events leading to failure was summarised as follows (van Niekerk and Viljoen 2005;Wagener 1997): • Rainfall falling on the upper compartment (see Fig. 10), with an area of 0.60 km 2 , flew through a gap in the dividing wall into the lower compartment (with an area of 0.72 km 2 ). The water accumulated in a pool (with an area of 0.15 km 2 ) near the northern wall and started overflowing • Water flowing through a small breach at the top of the tailings dam embankment eroded material from the outer slope face. • As the material was removed, the slope became steeper; the shear strength was inadequate to maintain these steeper slopes, leading to local slumping failures. • This eroded material was carried away by water released from the tailings pond, preventing stabilisation of slope to be reached. • At one point, this combination of retrogressive failure and water erosion exposed tailings that were in a metastable state in situ and instead of merely slumping, flowed as a liquid.
Unlike the channel-type scars that remained in the storage facility after Mount Polley failure, the failure of Merriespruit dam resulted in a half-funnel shape scar near the breach location (Fig. 10). This difference seems to be due to the amount of water being insufficient to scour channels in the Merriespruit case; it was the geotechnical failures that defined the extent of the mudflow. The dimensions of the scar area were about 300 m wide and extended about 200 m into the TSF (Figs. 10 and 11).
After the flow stopped, the reported slope of the tailings surface were between 2° (Blight and Fourie 2003) and 4.8° (Wagener 1997) in the breach section, while the ground slope was about 1.5° (Blight and Fourie 2003). From the graphical material presented in Blight and Fourie (2003) and Wagener (1997), the depth of tailings at rest in the breach section was about 5 m. Slopes were much steeper around the perimeter of the scar, from 10° to 20° (Blight and Fourie 2003). Liquefaction might not have occurred in this area. There were two main challenges in modelling of the breaching of the Merriespruit dam with EMBREA-MUD. First, breaching occurred not only due to water erosion but also due to slope failures. To simulate this additional effect, not explicitly included in the model, the value of dam erodibility coefficient was increased. Second, the liquefaction pattern was complex and varying. The parameter controlling whether tailings movement occurs or not is the yield stress τ y . The value of τ y can be calculated from the balance of forces at rest as: where g is gravitational acceleration (taken as 9.8 m/s 2 ), ρ is density, h 1 is depth, and S 1 is slope; the subscript 1 refers to mud (liquefied tailings). Density does not have a great impact on the simulation results as the initiation of mud movement depends on the yield stress per unit mass. Based on the aforementioned observed depth values (h = 5 m) and slope in the breach section, the values of τ y could range from 1800 (for S = 1.5° = 0.025 m/m) to 2500 (for S = 2° = 0.035 m/m) and 6000 Pa (for S = 1:12 = 0.083 m/m). The value of τ y would be significantly higher near the scar perimeter. The values of the parameters relevant for modelling are given in Table 3.
The size and bed slope (10°) of the modelling domain for mudflow were based on the surveyed data shown in Fig. 11. Tailings below the bed were considered not liquefied and therefore not included in the model. The results for various simulations are given in Table 4.
A comparison of different runs show that the model results show some sensitivity to the selected input parameters. The time to peak is shorter if the dam erodibility K D,0 is higher, and the two values are almost inversely proportional. Peak outflow rates are also markedly higher when erodibility (5) y = 1 g h 1 S 1 is higher, as would be expected. The same effect can also be observed on the final breach width, although the sensitivity of this output is less. The output parameter most sensitive to yield stress τ y is peak outflow rate, while time to peak and final breach width are less sensitive. Tailings outflow volume was not very sensitive to τ y or K D,0 because in all cases, most of the liquefied tailings were discharged, a consequence of a relatively steep underlying surface at 10°.
The predicted outflow volumes by the model were between 490,758 and 594,315 m 3 , which is lower but fairly close to the estimated outflow of tailings (600,000 m 3 ; Wagener 1997). This good match is a result of the model geometry, which was based on the observed final dimensions of the scar that developed on the surface of the tailings storage facility as a result of the failure. In practical applications, where these dimensions are not known, they have to be estimated. They could be obtained from the dam height and an assumed slope of post-failure profile. In the Merriespruit case, this slope in the flow direction was about 10°.
Other model outputs that could be compared to available information were also examined, in particular the final breach width and the chronology of outflow (i.e. the duration of the water outflow and the time of collapse). The final observed breach width was 150 m (van Niekerk and Viljoen 2005), which from a photograph taken by the Virginia Publicity Association (Duvenhage 1998) appears to refer to the width at the top of breach. The same photograph shows the side slopes of the breach to have a steepness of about 1:1.5 to 1:2. At the bottom, the width would therefore be around 50 m and around 100 m on average over the depth. The modelled widths ranged from 52 to 80 m, which is in the same range as the observed widths.
According to the narrative of the events recorded in the literature, the failure occurred rapidly at 21:00 on 22 Feb. 1994, two hours after water discharging from the facility was noticed. During this time, clear water, presumably from the dam, ran through the streets of the village (Blight and Fourie 2003). Hydrographs of outflow of water and mud for both runs with τ y = 2500 Pa are shown in Fig. 12.
The model results for both runs predicted an initial slow rise in water flows until they reach a rate between 20 to 30 m 3 /s, about 10 min before the peak in tailings outflow is reached. At that point, tailings start flowing out of the facility as mud, while before this time they were only present in smaller quantities as solids suspended in water. Mudflow increases dam erosion and both water and tailings outflow start rising rapidly. This point in time could be considered the moment of failure. The peak water discharge is reached five minutes later and the peak tailings outflow five minutes later. A similar pattern was observed in other model runs. These short time intervals corresponds to the observed sudden rupture of the dam in that there was no time to warn the inhabitants of the village once the dam breached (Wagener 1997). Between all runs, the simulated outflow peak times ranged between 1:35 and 3:04, i.e. the failure predicted by the model occurred in the range of between 1:25 and 2:54 after the start of water release, which is close to when the actual failure was observed (two hours after the start of water release).

Summary and Conclusions
EMBREA-MUD, a physically-based numerical model for simulation of tailings dam breaching is described in this paper. The model predicts the outflow rates of water and tailings and growth of the breach opening by simulating interactions between three layers. The water layer corresponds to the supernatant water stored above the tailings and is modelled as a Newtonian fluid. The mud layer corresponds to liquefied tailings and is modelled as a non-Newtonian fluid. The third layer is the dam itself, which is subject to erosion by the other two components. For modelling the water flow and dam erosion, the model shares the functionality with the extensively tested EMBREA model (Mohamed et al. 2002). For the non-Newtonian component, testing against laboratory experiments showed good agreement between the EMBREA-MUD predictions and observations for dambreak type flow. The model was then applied to back-analyse two different tailings dam failures reported in the literature, where overtopping and flow erosion played an important role and could therefore be modelled. The first case was the Mount Polley failure where a large amount of supernatant water was present in the tailings storage facility. The second case was the Merriespruit failure, with a smaller pond but increased water levels after a rainfall event.
These two cases provided valuable conclusions regarding approaches to modelling tailings dam breaches with EMBREA-MUD. Two different approaches to the definition of the model domain and selection of model parameters were used. In the case of Mount Polley, overtopping occurred after an initial foundation failure of the embankment reduced a section of the crest by more than 3 m. Following this initial collapse, a breach formed within this section and grew as a result of water and mud erosion. The discharging water was also able to scour a channelized path through the tailings dam through which the release of tailings occurred. To estimate the width of this channel, it was assumed that a formula for scour channel width that develops in water reservoirs during drawdown flushing could be used. This assumption, in combination with values of mud and erosion parameters collected from the relevant literature, produced a reasonably good agreement with the observations in terms of general chronology of the event and the amount of released tailings.
In the Merriespruit case, overflow started after rainwater increased the water level in the TSF. The failure process, however, was a result of both flow erosion (due to water and mud) and slope failures. To accommodate both factors within a model that does not simulate slope failures, the dam erodibility coefficient was increased. A sensitivity analysis with two values between the erodibility of the dam material and the erodibility of the tailings was performed. A sensitivity on the yield stress (within the values deduced from various reports on slope and depth of tailings "at rest") was also performed. With respect to the two most reliable observed parameters, the simulated time to peak and the final breach width, the results showed some differences; however, they were not greatly affected by varying these parameter values. The model was able to predict the sequence and duration of the events reasonably well. The area from which the tailings were released, however, had to be estimated based on the post-event survey. The main determining factor appeared to be slope failures. Hence, for modelling dam breach flows, the assumptions of model area should be determined based on geotechnical considerations.
For application of EMBREA-MUD to other tailings dams, the modelling approach can be selected depending on their hydrological and geotechnical similarity with one or the other case presented in this paper. If the dam is an active one, storing a significant amount of process water, the failure and outflow from the tailings storage facility is largely driven by water flow and the modelling approach described for Mount Polley can be used. If the tailings dam is an inactive one, with little ponded water, failure might be more related to slope failures in the dam and tailings, and the approach described for Merriespruit is more suitable. The Ouƞlow (m3/s)

Time (hours:minutes)
Water Mud selection of the approach is easier in the case of back analysis than if the model is used for prediction. Furthermore, if the model is used for predictions, we recommended that you perform a sensitivity analysis for the model's parameters (in particular, the dam the erodibility coefficient and the mud yield stress), within their expected ranges, in order to put confidence limits on the model outputs.
The ability of the model to predict breach outflow based on a range of possible physical properties of tailings is an advantage compared to simpler regression models. Other advantages include modelling water -tailings interactions (including tailings erosion), inclusion of bed slope configuration in the model (the proportion of discharged tailings will be higher from a TSF in a steep valley than from one on a flat ground), and outflow timing.
The model has proven to be a very useful tool to estimate outflow hydrographs, which could be used as input for other studies to simulate the spreading of water and mud flows downstream and therefore, to support an understanding of their economic, social and environmental impacts, contributing to better manage tailings infrastructure. The assessment of impacts could be used, for example, at the feasibility stage when selecting the location for the dam, for an existing dam to assess its impacts, or to assess the impact of dam raising in the (unlikely) case of dam failure.