The failure of the Stava Valley tailings dams (Northern Italy): numerical analysis of the flow dynamics and rheological properties

Tailings dams are made up of mining residue deposits, and they represent a high risk, in terms of mechanical instability. In the event of collapse, the tailings in such dams may be released and flow over long distances, with the potential risk of extensive damage to property and life. The traditional geotechnical assessment of tailings facilities has mainly concentrated on the stability of tailings dams, while relatively few studies have investigated the flow of tailings released after dam failure. In this context, it is possible to state that, if the complex rheological behaviour of the tailings material is captured correctly during the flow, numerical modelling can be used to contribute to a better comprehension of the flow characteristics and for the assessment of the possible extension of the impact area. Considering the wide range of possible rheological behaviour that tailings flows can assume (from laminar to turbulent), this paper presents the new version of a computer model, which was designed to simulate the motion of rapid flow movements across 3D terrain. This new version integrates the existing rheological kernel (Frictional, Voellmy) with two new rheological laws (Bingham and Turbulent), and adds the possibility of changing the rheological properties of the flowing mass during the propagation process. The code has been applied to the disastrous flow that was caused by the failure of a pair of tailings impoundments in the Stava Creek Valley (Italy) in 1985. Since different interpretations on this flow behaviour already exist in literature, and since a large number of changes in the rheological values along the run-out path have been proposed to recreate its dynamics, new simulations, carried out with different rheological combinations, are presented and discussed here in order to obtain a better understanding of the flow dynamics and to identify the rheology that reproduces the phenomenon that occurred with the fewest possible changes in the rheological values along the runout path. The latter aspect is particularly important when numerical analyses are used for prediction purposes. The great rheological flexibility of the new code has allowed the Voellmy rheology and a combination of its parameters to be identified as the most suitable to describe the Stava flow, even where the run-out path presented critical characteristics.


Background
Some several thousand million tonnes of mining waste is produced each year. The wastes resulting from the chemical and mechanical processes of mining extraction, called tailings, are mixtures of crushed rock and processing fluids from mills, washeries or concentrators that remain after the extraction of metals, minerals, mineral fuels or coal from the mine resource (Kossoff et al. 2014).
Tailings are often accumulated in upstream/downstream valleys, or in ring impoundments and, in all these cases, are retained by a dam. If, for any reason, the dam wall should break, the accumulated tailings can escape the impoundment, thus giving rise to a flow of material that is characterized by a high energy content, with ensuing serious socioeconomic and environmental consequences.
Guaranteeing the stability of such impoundments is one of the most challenging tasks in mine waste management. These storage basins are in fact particularly vulnerable to failure for the following reasons (Rico et al. 2008a): (1) dams made of locally derived filling material (soil, coarse waste, overburden material from mining operations and tailings); (2) dams subsequently raised as solid material coupled with a severe increase in effluent (plus runoff from precipitation); (3) lack of regulations on specific design criteria; (4) their extremely large area makes punctual criticality identification difficult, and would require an extensive and continuous monitoring, which, apart from the very high costs, is also time consuming, and (5) the high cost of remediation works, following the closure of mining activities.
Even when a tailings dam has been designed correctly, there is always a risk of its operational failure. Despite this, even though 1) knowing how far the tailings can travel and 2) being aware of the extent of damage could help in assessing the risk associated with tailings collapses, the traditional geotechnical assessment of tailings facilities has mainly been concentrated on the stability analysis (Vick 1990), and relatively few studies have investigated the flow of tailings released following dam failure.
Many numerical dam-break analysis models have been developed to analyze the breakage of water-storage dams. However, these models cannot be used directly to estimate the run-out distance of a flowing mass in the case of tailings dam collapse, especially because of the differences in the rheological nature of the involved mass. Breaching floods are in fact usually composed of highly water-saturated, oozy sediments, which can exhibit various kinds of fluid behaviour, ranging from debris flows to muddy flood water (Rico et al. 2008a). If the outflow volume versus run-out distance is considered, these floods can be grouped into two main categories: (1) floods with high-viscosity spilled mine waste and (2) floods with large volumes of water inside tailings dams.
Because of this complexity, Hunter and Fell (2001) recommend not using empirical correlations for flow slides in tailings dams. Predictions based on numerical modelling are considered more appropriate, and, since debris and sediment flows exhibit similar behaviour to tailings flows, numerical models implemented for the analysis of rapid landslide run-out are suitable for the modelling of tailings dynamics. The selection of the flow rheology is important for both situations, and can determine the accuracy of the method.
The above mentioned mathematical models may be divided roughly into: (i) one-phase models, which describe the flow resistance behaviour of either the water and fine material slurry, or of the entire fluid-solid mixture(e.g. O'Brien et al. 1993;Pirulli 2005;McDougall 2006;Pastor et al. 2009); and (ii) two-phase models, which consider both a fluid phase and a solid phase separately (e.g. Bozhinskiy and Nazarov 2000;Iverson and Denlinger 2001;Pudasaini 2012). The discussion is here limited to one-phase approaches since the here presented model was designed to facilitate both a systematic comparison of single-phase flow resistance relations, such as made by Naef et al. (2006) and Rickenmann et al. (2006) for other case studies, and to identify the rheology that is able to back analyse a tailings dam failure in the most accurate way with the fewest possible changes in rheological values along the runout path.
In fact, because of the wide range of rheological behaviour that tailings flows can assume, a new version of the single-phase continuum-mechanics based RASH 3D code (Pirulli 2005), which was originally designed to study flow-like landslide run-outs, has been implemented and is here described briefly. The introduced upgrade code integrates the existing rheological kernel (Frictional, Voellmy) with two new rheological laws (Bingham and Turbulent), and adds the possibility of changing the rheological parameters of a flowing mass, in a GIS (geographic information system) environment, during the propagation process.
The new version of RASH 3D has then been used to simulate the well documented case of a pair of tailings dams whose failure, on 19 July 1985, generated a catastrophic flood in the area of Stava, Northern Italy (Berti et al. 1988;Chandler and Tosatti 1995).
With the aim of contributing to the definition of rheological laws and their range of parameters that can be used for prediction purposes, the back-analysis of the Stava event has been carried out in an attempt to try to understand whether a single rheology is able to simulate the whole process, from triggering to deposition, without any changes in its rheological values along the run-out path, or with a minimum number of changes. Numerical analyses with Turbulent, Voellmy, Frictional and Bingham rheologies have therefore been carried out by combining different sets of rheological values. The obtained results are discussed and compared with Takahashi's (1991) 1) velocity estimates of the flowing mass attained at various points along its path, and 2) physical measurements of flow super-elevations measured in channel bends.

The post-failure rheological behaviour of tailings
The liquefied material that originates from the failure of tailings has fluid-like characteristics, and the flow behaviour has been reported to be mainly laminar on relatively flat slopes and turbulent on steeper slopes (Jeyapalan et al. 1983a;Blight 1997). Jeyapalan et al. (1983a) suggested that the Bingham plastic rheological model is a good approximation model for the flow of liquefied tailings, when the flow is laminar. Instead, the use of flood routing computer programs is recommended to model the turbulent flow of tailings. Jeyapalan et al. (1983b) applied a Bingham plastic model (TFLOW computer program) to both the Aberfan case, in South Wales (1966), and to the Gypsum Tailing Dam incident, in East Texas (1966), and laminar flow behaviour of the released tailings was pointed out. The computation results fitted well with the actual observations, even though the model was limited to one dimension, and was therefore unable to take into consideration the spreading of the flow. A turbulent flow analysis (GVFP 1978 andFread 1978 computer programs) was instead made, through the use of the Manning resistance relation, for the Buffalo Creek Dam failure, which occurred in the USA in 1972, due to the presence of a very fluid material made up of a mixture of waste and large quantities of water. The computation results were able to represent the actual flow data rather well.
In most flow studies on tailings deposits, the flow will be laminar in nature, and the liquefied tailings will behave approximately like a Bingham plastic fluid, exhibiting both a yield stress and viscous flow characteristics, once this stress has been exceeded.
However, the entrainment of free water in the sliding mass can significantly alter the flow properties of the fluid. In laminar-type flows, the entrainment of the water in the lower portion of the flow can significantly reduce the yield strength, at the base of the sliding mass, thus resulting in much longer travel distances. Large volumes of water were reported to have been stored within the impoundments of Stava (Italy), Bafokeng (South Africa), Merriespruit (South Africa) and Buffalo Creek (USA) at the time of their failure.
Furthermore, in several cases, the confinement or partial confinement of the flow in running streams and rivers resulted in the development of hyper-concentrated stream flows that travelled for great distances on relatively flat slopes. This phenomenon becomes particularly significant where the tailings enter an actively flowing watercourse. It is considered that the flows of the Stava (Italy), Buffalo Creek (USA) and Bafokeng (South Africa) events transformed into turbulent stream-flows as a result of their confinement in flowing watercourses (Hunter and Fell 2001).
It can be determined whether a flow will be laminar or turbulent by using the Hanks and Pratt (1967) or the Takahashi (2007) criteria. Hanks and Pratt (1967) made a detailed analysis of a large number of published experimental data, and proposed a chart to determine the transition conditions for Bingham plastic fluids. In this chart, the critical Reynold number for transition from turbulent to laminar flow is expressed in terms of the Hedstrom number (Jeyapalan et al. 1983a). Takahashi (2007) defined two kinds of debris flows in a wider sense: one is the quasi-static-debris flow, in which Coulomb friction stress dominates, and the other is the dynamic debris flow. In order to transmit quasistatic Coulomb friction stress, the solids concentration (C) has to be large enough to guarantee that the particles are always in contact, even though their relative position changes continuously. Bagnold (1966) stated that this condition is fulfilled when C is larger than 0.51 for natural beach sand, but this value would depend on the size of particles. Under such a densely concentrated condition, the other stresses become small, and the motion would be quasi-static. Dynamic debris flows can instead be subdivided into three further subclasses: when grain collision stress dominates, the debris flow becomes a stony-type flow, when turbulent mixing stress dominates, it becomes a turbulent-muddy-type flow, and when viscous stress dominates, it becomes a viscoustype flow. Transition from one type of dynamic flow to another depends on the solid concentrations and on the relative depth of the flow, h/d (where h is the depth of the flow and d is the particle size), on the Bagnold number (ratio of the characteristic shear stresses due to grain collisions and liquid viscosity) and on the Reynolds number (ratio of the characteristic shear stresses due to grain collisions and solid friction without the solid friction coefficient) (Iverson and LaHusen 1993). When the Bagnold number is large, and the relative depth is small, a stony-debris flow occurs. When the Bagnold number and the Reynolds number are small, a viscous-debris flow occurs. When the relative depth and the Reynolds number are large, a turbulent-muddy debris flow occurs.

Theoretical aspects of RASH 3D
The RASH 3D numerical code is based on a one-phase continuum mechanics approach, and on depth-averaged St. Venant equations. The real heterogeneous mass is replaced with an incompressible equivalent fluid, whose behaviour is described by the depth-averaged balance equations of mass and momentum: where v x ; v y denote the depth-averaged flow velocities in the x and y directions (z is normal to the topography), h is the fluid depth, τ zx , τ zy are the shear resistance stresses, σ xx ; σ yy are the depth-averaged normal stresses, ρ is the mass density and g x , g y are the projections of the gravity vector in the x-and y-directions, respectively.
The governing equations in RASH 3D (1) are solved in an Eulerian framework, on a triangular finite element mesh, through a kinetic scheme that is based on a finite volume approach (see Audusse et al. 2000;Bristeau and Coussin 2001;Mangeney-Castelnau et al. 2003).
The rheology of the material is modelled by a single term, which describes the basal shear stress that develops at the interface between the moving mass and the sliding surface. Two new rheologies have been implemented into the code, according to Jeyapalan et al. (1983a) and Takahashi (2007), in order to run numerical analyses concerning the collapse of tailing dams. These rheologies are the Bingham one, which is used in the case of a laminar flow, and the Turbulent one, which is used in the case of a turbulent flow. Furthermore, a new geographic information system (GIS) integrated function, which makes it possible to change the type of rheology and/or the rheological parameter values along the run-out path, has been introduced into the improved RASH 3D version, to allow changes to be made to the flow characteristics during flow propagation.

Rheological kernel
A flow may switch to a turbulent regime, which is characterized by intense mixing, at relatively high inertial to viscous stress ratios. The turbulent basal shear resistance is proportional to the square of the depth-averaged flow velocity, and it can be calculated using the Manning equation: where n is the Manning roughness coefficient, and the subscript i = x,y, respectively.
A commonly used alternative to eq. (2) is the Chézy equation: where C is the Chézy coefficient, which is related to the Manning coefficient (n) by C = h 1/6 /n.
One disadvantage of this approach is that it cannot reproduce the cessation of motion on gently sloping surfaces. Nevertheless, Costa (1997) and Jin and Fread (1999) showed that the flow depth and the velocity of a channelized flowing mass can be simulated reasonably well after calibration with Manning or Chézy coefficients. However, this limitation can be overcome by including a term in the rheological formulation that describes the stopping of the flow on a sloping surface (e.g. Hungr and McDougall 2009;Naef et al. 2006;Rickenmann et al. 2006). In this regard, the Voellmy flow relation, which consists of a turbulent term, ξ (= C 2 ) that accounts for velocity-dependent friction losses, and a Coulomb or basal friction term, μ (= tanδ), which is used to describe the stopping mechanism, where the basal friction angle δ is generally only a fraction of the Coulomb angle φ (McDougall and Hungr 2005), had already been included in the RASH 3D rheological kernel (Pirulli and Marco 2010).
On the other hand, the Bingham resistance model combines plastic and viscous behaviour. A so-called Bingham fluid behaves like a rigid material below a given threshold yield strength, but like a viscous material above this threshold. The basal shear resistance can be determined by solving the following cubic equation: where τ y is the Bingham yield stress and μ B is the Bingham viscosity.
The third-order polynomial has been solved and implemented in RASH 3D using the polynomial economization technique proposed by Pastor et al. (2004).

The collapse of tailings dams at Stava
The Stava Creek, with a drainage area of approximately 21 km 2 , is located in the Dolomite Mountains, North-Eastern Italy, in the upper valley of the Avisio River (Fig. 1).
On July 19, 1985, the fluorite tailings dams of the Prealpi Mineraria mine,located on the Porcellini Creek, a small tributary of the Stava Creek, at an elevation ranging from 1330 to 1380 m asl., failed (Genevois and Tecca 1993). The tailings dam was built on a hill-slope, and consisted of two basins (Fig. 2). The failure occurred when the up-slope basin collapsed into the lower one. The inflow of this material caused overtopping of the lower basin and its subsequent collapse.
The tailings that spilled out of the failed dams flowed downstream at a high speed. According to a survey that was conducted after the collapse, the total volume of sediment that had flowed out was estimated to be 88.300 m 3 (12.000 m 3 of fine sand and 76.300 m 3 of silt). Moreover, after an estimation of the volumes of the supernatant water and the interstitial water, a total volume of 185.000 m 3 of mud debris, with a volumetric solids concentration of 0.476, was calculated to have been released as a muddy debris flow (Takahashi 2007).
The villages of Stava and Tesero, located along the stream channel and at its end, respectively, were completely wiped out or buried. The flowing mass killed 268 people and completely destroyed three hotels, eight bridges, 53 homes and six industrial buildings. Nine more buildings were seriously damaged (Van Niekerk and Viljoen 2005).
The released mass surged down the 600 m long 10°m ountain slope, which was covered by grass in the upper part and by forest in the lower part, and the flow direction then changed perpendicularly after clashing into the cliff on the left side of the Stava Creek (Fig. 3b). At that moment, it destroyed many hotels and buildings. It then flowed channelized along the Stava Creek for about 3800 m, weaving right and left, and finally stopped at the junction with the Avisio River (Takahashi 2007) (Fig. 3).
After the event, a layer of dense mud, between 200 and 400 mm thick, covered an overall area of 43.5 ha over a length of 4200 m (Genevois and Tecca 1993;Van Niekerk and Viljoen 2005). At the inquest that followed the disaster, it became apparent that the tailings dams had not been subjected to any detailed stability checks for over a period of more than 20 years.
The possible causes of instability included the following (Van Niekerk and Viljoen 2005; Takahashi 2007): overloading, due to the additional elevation of the embankment that was under construction; liquefaction of the slime, due to water seeping from the slope; toe failure of the upper embankment, caused by seepage; failure of the upper embankment, due to an incorrect installation and blockage of the drainage pipes; insufficient separation between the pooled water and the embankment. Furthermore, the increase in seepage water might have been affected by the heavy rainfall that had occurred two days before the collapse (Takahashi 2007), and by the fact that the embankments had been constructed using the up-stream method. The up-stream method is the oldest and most economic construction method, which begins with the construction of a starter dam at the downstream toe. The resulting tailings generally have a low relative density and a high water saturation degree (Vick 1990).

Definition of the flow characteristics
Although the flowing mass had an intensive destructive power, as well as a high fluidity, the Stava Creek channel itself did not suffer from much erosion or deposition (Berti et al. 1997). Takahashi (1991Takahashi ( , 2007 stated that no erosion had occurred, because the solid fraction inside the water-sediment mixture was so high (estimated to be about 0.48) that the flow could not have become denser because of erosion. In fact, different debris flow models are based on the assumption that a debris flow can erode and incorporate material from the bed until the solid concentration of the flow is as large as the equilibrium concentration for the channel gradient (e.g. Takahashi et al. 1992). This concept is also supported by Hungr et al. (2005), who stated that flows with lower volumetric sediment concentrations can be expected to be more erosive than flows with larger sediment concentrations. Finally, the laboratory experiments of Rickenmann et al. (2003) pointed out that the erosion rate initially increased as the volumetric sediment concentration increased, and then decreased as the sediment concentration increased.
As far as the grain size distribution is concerned, ISMES (1986) analysed the material collected on site, and found a clean separation between the materials that seemed to be part of the embankment and those finer ones inside the basin. The two types of material were almost completely differentiated by the grading curves, which defined them as silty sand and clayey silt. Takahashi (1991) underlined that the involved material was so fine that it had a relative depth of about 10 5 . In this condition, the resistance to flow is similar to that of a plain water flow, and the Manning equation can be applied. As a consequence, the Manning roughness coefficients were obtained by Takahashi (1991) in each section by means of reverse calculation from the data on the velocity computed with the Lenau formula (1979), which was in turn used to estimate the flow velocity at the super-elevations at the bends (Table 1).
Both theoretical and experimental data Takahashi 1983, 1986;Takahashi 1991;2007) have indicated that the velocity near the bed increases and the velocity in the upper layer approaches a uniform state for either decreasing concentrations or increasing relative depth, thus supporting Takahashi's assumption; the velocity distribution forms approach those of a plain water flow. The dual dependence on the concentration and relative depth indicates that there is no a single value of the sediment concentration that can be considered representative of a given type of phenomenon, but this value depends on the combination of the concentration and the relative depth value, for example (Takahashi 2007).
Nevertheless, Berti et al. (1997) subdivided the Stava process into two phases: flow slide (from the tailings dams to the Stava Creek) and debris flow (along the Stava Creek valley to the confluence of the Stava Creek with the Avisio River), and evidenced, through simple scaling analysis (based on the work of Iverson and LaHusen 1993), that the friction effect had prevailed over the behaviour of the Stava debris flow.
As far as velocity estimation of the flowing mass at Stava is concerned, Takahashi (1991) made simple onsite measurements of the flow cross-section, at the positions indicated by the numbers in Fig. 4, and, by measuring the flow super-elevation on the outerbank versus the internal bank where the flow entered the curves, estimated the flow velocity (Table 1) with the Lenau formula (1979). Lenau (1979) proposed a method for estimating the height of crests along an outer edge, when a bending flow passes through a trapezoidal cross-section. The formula, rearranged by extrapolating v, is where E max is the maximum elevation of the outer edge from the surface of the undisturbed flow, m is the reciprocal of the side slope of the channel cross-section, h 0 is the centre line depth of the undisturbed upstream flow, b is the base of the trapezoidal cross-section, v is the velocity of the undisturbed flow, and r 0 is the centre line radius of the bend. This formula can also be applied to a rectangular channel, as indicated by Takahashi (2007).
Since a seismogram at the nearby Cavalese station (3.7 km from the collapsed tailings dams) had recorded the event, Takahashi used the available data to interpret the large vibrations as collisions of the flowing mass against particular obstacles along its path (Takahashi 1991). On the basis of records of the event, it was estimated that the whole event, from the flowing into the Stava Creek (Fig. 4, section1) to the flowing out to the Avisio River, lasted about 425 s, and that velocities of up to 37 m/s were reached. In particular, the forefront arrived at the Romano Bridge (Fig. 4, section 13) about 225 s after it had passed section 1 (Fig. 4).
The propagation rate of the forefront along the Stava Creek, according to the shock records at Cavalese, is in good agreement with the velocity values computed using the Lenau formula (1979).  Takahashi (1991) applied the kinematic wave theory to simulate the flowing process of the Stava tailings dams, but this approach requires a quasi-one-dimensional flow, and it completely neglects the effects of extreme bends, irregular geomorphological slopes and cross-sections of the channel. These important assumptions are satisfied in the central part of the Stava Creek, due to the narrow shape of the gorge, but a more complex and fully threedimensional approach is necessary to simulate the initial change in direction of the flowing mass and the final deposition at the bottom of the valley. Three-dimensional analyses (i.e. without simplifying the topography) were carried out with the RASH 3D code (Pirulli 2005). The event interpretation was extended outside of the previously studied Stava Creek area, that is, from the tailings dams triggering area to the Stava village and from the Stava Creek to the junction with the Avisio River. This has allowed the following to be investigated:

Methods
the propagation of the front, after the dam rupture, that overflowed onto the open terrain, which was partly covered by grass and partly by trees, before the impact with the village of Stava; the channelized flowing of the mass along the Stava Creek, using the real irregular geomorphology of the channel; the interaction of the flowing mass with the Romano bridge, and the final deposition at the intersection with the Avisio River.
The complex geomorphology of the terrain was reproduced accurately by downloading the 10 m digital elevation model (DEM) from the Public WEB GIS of the Trentino Geocartographic Portal of the Autonomous Province of Trento (http://www.territorio.provincia.tn.it/portal/server.pt/community/lidar/847/ lidar/23954). The tailings dams were reconstructed using the available cartography and preserving the original event volume. Even though a post-event DEM was used in the analyses, it was considered representative of the 1985 morphology, since, according to different Authors (e.g. Takahashi 1991;Berti et al. 1997), the Stava Creek channel did not suffer from much erosion or deposition. On this basis, the erosional process was neglected in the numerical simulations.
Four rheologies were tested: Frictional, Voellmy, Bingham and Turbulent. Taking advantage of the GIS function that is now implemented in RASH 3D , analyses were also run changing the rheological values along the runout path.
The obtained results were analysed, in terms of mass distribution (cross sections, plan distribution and longitudinal profile), and compared with Takahashi's flow velocity estimation (Table 1). Furthermore, the Δv value, that is, the difference between the RASH 3Dcomputed values and Takahashi's estimated velocity values, were also represented for each analysis to facilitate a rapid identification of the best simulations. A negative Δv value indicates an underestimation of the computed velocity with respect to the estimated one.

Results and discussion
According to Takahashi (1991), Manning's equation was adopted in the first numerical analysis (T-t1) to describe the resistance to flow with the combination of parameters reported in Table 1. Since the suggested rheological values only concern the run-out path between section 2 and section 13 (Fig. 4), the first value of Manning's coefficient (n = 0.04) was extended upstream, while the last value (n = 0.12) was extended downstream to allow RASH 3D to perform a complete simulation of the phenomenon (from triggering to deposition) ( Table 3, T-t1). A rather good agreement can be observed from the obtained results between the RASH 3D computed values and Takahashi's estimated velocity values, along the whole propagation path (Fig. 5a), with the exception of sectional reach 6′-7, where Δv = -12.5 m/s (Fig. 6a, T-t1).
In order to reduce the number of different rheological values (i.e. Manning's coefficient) to use in a single analysis, and to investigate the possibility of obtaining a better approximation of the mean velocity at sectional reach 6′-7, two new analyses were run: T-t2, with n = 0.04, and T-t3, with n = 0.02 (Table 2). It was observed that, compared with the velocity values calculated by Takahashi, the T-t3 results overestimated the mean velocities, with the exception of sectional reach 6′-7 (Δv = -0.7 m/s), while theT-t2 results approximated Takahashi's mean velocities in a rather good way, except for sectional reach 6′-7 (Δv = -12.5 m/s), and after section 10 (Figs. 5 and 6).
Although one of the aims of these analyses was to find the smallest number of different rheological values necessary to use in a single analysis, the T-t4 analysis was carried out by combining two values of Manning's coefficient: n = 0.04 from the triggering point to section 10, and n =0.12 in the remaining portion of the run-out path (Table 3). However, as can be observed in Figs. 5 and 6, the results did not improve.
On the basis of these results, a different rheological law was tested, still looking for the possibility of reproducing the event with the combination of parameters that requires the smallest number of changes of rheological values along the runout path. A Voellmy rheology, where friction coefficient, μ, and turbulent coefficient, ξ, have to be calibrated (see Eq. 4), was then applied. The rheological values suggested by Golder Associates Ltd (1995) were initially used with a single In the B-t3, B-t4 and B-t5 tests, the flowing mass stopped before reaching sectional reach 12′-13 combination of rheological parameters for the whole propagation path (Table 2: V-t1, V-t2, V-t3 and V-t4). Figures 5 and 6 point out that the V-t1 results simulated the flow velocity after section 10 in a rather good way, but they greatly underestimated the velocities before section 10 of the run-out path. The opposite conclusion can be drawn for the results of the V-t2, V-t3 and V-t4 analyses.
As a consequence, a combination of Voellmy rheological parameters was adopted in the V-t5, V-t6 and V-t7 simulations, using a first set of values up to section 10, and a second set from section 10 onwards (Table 3). A good approximation of the mean velocity was obtained for the different sectional reaches, except for the V-t5 analysis, even though the value for sectional reach 6′-7 remained critical (Fig. 5).
The possibility of simulating the event with a reduced number of changes in the rheological values along the runout path was highlighted by moving from the Turbulent (T) to the Voellmy (V) rheology. This is an important Fig. 6 Flow velocity difference, Δv, defined as the difference between the computed (RASH 3D ) and measured (Takahashi 1991) values. a T: Turbulent rheology; b V: Voellmy rheology; c B: Bingham rheology; d F: Frictional rheology. -t: indication of the test number. In the B-t3, B-t4 and B-t5 tests, the flowing mass stopped before reaching sectional reach 12′-13 aspect, since it is difficult to foresee and define a large number of local rheological changes along the runout path when forward-analyses are run.
Since Berti et al. (1997) stressed that friction played the most significant role in determining the behaviour of the Stava debris flow, a further set of analyses was run adopting a simple Frictional rheology, where the only parameter that needed to be calibrated was the basal friction angle (δ). As for the previous rheologies, the first analyses used a single rheological value for the whole run-out path (Table 2: F-t1 with δ = 10°, F-t2 with δ = 8°, F-t3 with δ = 7°). Figure 5 shows how it was impossible to simulate the flow velocity in the whole propagation process correctly with these three analyses. Although F-t1 gave rather good results for sectional reaches after section 10, it underestimated the velocities before this section; F-t2 and F-t3 gave a better approximation of some of the velocities before section 10, but overestimated the velocities in the other sectional reaches. The F-t4 analysis was carried out in an attempt to combine different friction angle values. However, the obtained results were unsatisfactory, and underlined the difficulty in using the Frictional rheology to simulate the dynamics of this type of phenomenon, even when back-analysis is conducted. The use of a combination of frictional values in a forward analysis would be even more complex.
The Bingham rheology had never previously been applied to simulate the Stava event, on the basis of the interpretations of several Authors (see the "Tailings postfailure rheological behaviour" section). However, this rheology was adopted in the last proposed set of analyses in order to test its eventual capacity to simulate the flow velocity along the run-out path.
The first analysis conducted with the Bingham rheology (see Eq. 5) was carried out with the combination of rheological values that had been defined by Pastor et al. (2004) for a gypsum tailing impoundment in East Texas, that is,τ y =1 kPa and μ B =0.05kPa⋅s (Table 2: B-t1). The results of this run only gave a rather good approximation of the velocity between sections 7′ and 10 of the flowing channel (Fig. 5). It can be observed, in Table 2, that the combination of rheological parameters used for the B-t2 analysis approximated the velocities in  a rather good way, but only in the lower part of the propagation path (downstream from section 7), while the B-t3 results simulated the upper part of the flow path rather well, but largely underestimated the velocities in the lower part, and showed that the mass in fact stopped before reaching section 9 (Figs. 5 and 6). However, combining the rheological values of B-t2 and B-t3 in analysis B-t4 (Table 3) led to results that were even less satisfactory than those obtained for a single combination of parameters for the whole run-out path (Figs. 5 and 6). After considering the previous results, the last Bingham analysis (B-t5) was carried out with the combination of parameters indicated in Table 3. Although a large number of rheological values were considered in this single analysis, the results were not able to reproduce the complete phenomenon, which would seem to confirm that this flowing mass had different rheological behaviour.
Considering all of these analyses, two critical points pertaining to the Stava run-out path have clearly emerged: a change in the flow direction took place when the mass reached the Stava village; sudden increases in the channel cross-section occurred downstream from sections 10-10′, and this was combined with a great resistance to flow, which was induced by the houses on the flood plain and the backwater effects of the Romano bridge.
Although the Voellmy rheology was considered with the V-t6 and V-t7 simulations, was capable of reproducing the whole Stava event dynamics with just a few changes in the rheological values along the run-out path, a single couple of rheological parameters (e.g. V-t2 and V-t4) was only able to reproduce the change in the flow direction and the channelized process, but a different combination of rheological values was always necessary from section 10 onwards. This pointed out a change in the flow behaviour from sections 10-10′ onwards.
The possibility of simulating the event in a satisfactory way with the Voellmy rheology, by reducing the number of changes in the rheological values, compared to the number of changes adopted by Takahashi (1991), has evidenced the combined turbulent and frictional nature of the Stava flow mass. This work has then combined the turbulent interpretation of the event given by Takahashi (1991) with the prevalent frictional process pointed out by Berti et al. (1997). Figure 7 compares the computed heights (H) of the flood marks (super-elevation) on the right and left banks of the flow, with respect to the lowest level at each cross-section, obtained from the V-t6 and T-t1 numerical analyses, with the values indicated by Takahashi Fig. 7 Super-elevation on a the left bank and b the right bank at the cross-section locations shown in Fig. 4. Comparisons between Takahashi's (1991) measured values and the computed RASH 3D results (1991,2007). The V-t6 numerical simulation was the best simulation, in terms of flow velocity and a limited number of changes in the rheological values along the run-out path, and the T-t1 numerical simulation reproduced the event with a 3D numerical model using Manning's values as suggested by Takahashi (1991) The super-elevations computed with the RASH 3D code for both banks mirror the height trend observed by Takahashi (1991) in a satisfactory way. However, as far as the left bank is concerned (Fig. 7a), a certain discrepancy can be observed between the computed and observed values, especially in sections 3 and 10′. As already mentioned, these sections were located where the mass initially changed flow direction, and where the channel-section suddenly increased. The numerical results pertaining to the right bank (Fig. 7b) follow the general height trend, but they indicate a certain overestimation with respect to the observed values from sections 2 to 9.
However, it should be pointed out that the positions of the cross-sections analysed by Takahashi (1991), and to which reference has been made in the present paper, were obtained from a printed image. As a consequence, errors were certainly introduced by stretching and rotating the image to place it in the "right" position on the GIS map. Although this could have partially contributed to the discrepancy between the computed and measured flow velocities and the super-elevations, the obtained numerical results reflect the trend of the real event.
Finally, considering the analysis of the entire event process, in terms of plan distribution of the debris along the whole run-out path, Fig. 8 visualizes a time sequence of the numerical event, up to the final deposition, for the previously identified best simulation (i.e. V-t6). A fairly good agreement with the in situ observed run-out area can be observed.

Conclusions and further developments
This paper presents an updated version of the one-phase continuum-mechanics-based RASH 3D code, which integrates the existing rheological kernel (Frictional, Voellmy) with two new rheological laws (Bingham and Turbulent), and adds the possibility of changing the rheological properties of the flowing mass during the propagation process, using an integrated GIS function.
The upgrade of the code originated from the intent to investigate the complex rheological behaviour of tailings flows released following a dam failure, which till now has been the subject of relatively few studies, despite the serious consequences that the collapse of a tailings dam can induce.
The collapse of the Stava tailings dams was chosen as a case study, due to the amount of information that was available, and to the existence of different interpretations of the flow dynamics of this event in literature (e.g. Takahashi 1991 andBerti et al. 1997).
A large number of numerical analyses with Turbulent, Voellmy, Frictional, and Bingham rheologies have been carried out by combining different sets of rheological values, with the aim of reproducing the dynamics of the past event, in terms of flow velocity, depth and plan distribution, using the fewest possible changes in rheological values along the runout path. The idea of identifying a single rheology that would be able to simulate the whole process, from triggering to deposition, without changing its values along the run-out path or with the fewest number of changes, could in fact contribute to obtaining rheological indications that could be used for forward-analyses. Takahashi (1991) interpreted the Stava flow as turbulent, Berti et al. (1997) considered the the flow as being mainly frictional, but all the literature on this flow agree on the inapplicability of the Bingham rheology to this case.
The results obtained from the present analyses have indicated that the Voellmy rheology, which combines a frictional and a turbulent term, is the most appropriate to simulate the Stava event, with the simplest combination of rheological values. The Turbulent rheology, suggested by Takahashi (1991), required Manning's value to be changed many times along the run-out path. The same can be stated for the friction angle of the Frictional rheology. In the case of the Bingham rheology, the combination of a large number of values did not lead to satisfactory results, thus confirming the inapplicability of this rheology to the Stava event.
The best results of the Voellmy rheology have been obtained from test T-t6, in terms of both the mean velocity in the sectional reaches, and the super-elevation on the left and right banks of the whole run-out path.
In reality, the complex process of the collapse of tailings dams involves different mechanisms, any of which may dominate at different locations. However, considering the current limited understanding of these mechanisms, it has been shown that simplification of the phenomenon, by means of an appropriate rheological law and rheological parameters, should not preclude the successful simulation of the bulk characteristics of a real event.
Finally, it can be stated that the new model significantly expands the possibility of simulating the dynamic behaviour of flow-like landslide events of different natures. With continued back-analysis of real events and investigations on the correlation between the characteristics of the event and the most suitable rheology, important indications can be obtained on the rheological choices necessary to make forward-analyses.
Further developments of this research will certainly include the use of more sophisticated methods (e.g. twophases models). However, the generalized use of more complex numerical models for the back-and forward-analyses of real events still remains limited, since information on individual sites and the characteristics of events are usually limited.