Stability and characterisation of spoil heaps in European surface lignite mines: a state-of-the-art review in light of new data

The large amount of spoil material produced during the mining process imposes a significant economic and environmental liability on lignite producers. In this context, the present paper provides an overview of the geotechnical characteristics of European lignite mine spoil heaps and discusses their significance to the stability of the heaps. In order to achieve this, samples collected from spoil heaps of Polish, Czech and Greek mines are analysed and the results are compiled with data from the literature. A major conclusion drawn is that both physical and engineering properties of spoil heaps indicate a noteworthy variability, which is larger than typical in situ ground material. This is because of the additional factors affecting spoil heap deposition, such as the transportation and dumping method. Furthermore, failure mechanisms and case histories of large instabilities in lignite spoil heaps are critically discussed in order to better understand triggering failure mechanisms. It is concluded that classical assumptions made for natural soil slopes and relevant limit equilibrium models should be cautiously applied to spoil heaps. The challenges associated with numerical and probabilistic modelling of spoil heap stability, such as the inherent spatial variability of spoils and the time-dependent changes in their geotechnical properties, are also critically discussed. Finally, important research gaps in design and analysis of spoil heap stability, such as the absence of appropriate constitutive models developed specifically for spoil materials, are summarised.


Introduction
Lignite, also known as brown coal, is a low-rank coal (early stages of coalification) and one of the cheapest sources of energy worldwide. Lignite is basically a peat that has been transformed into a brown-black rock, which sometimes contains recognisable plant structures. Europe accounts for roughly 40% of global lignite reserves. The top European producers include Germany, Greece, Poland and the Czech Republic (BGR 2014;Euracoal 2017;EY 2014). Surface or opencast mining of lignite requires excavation of the overburden material that covers the coal seam. The ratio of the volume of the stripped overburden (known as spoil) to the mass of produced lignite is known as the stripping ratio, which is reported to range between 2.2:1 and 9.5:1 m 3 /t throughout Europe (EY 2014). As an example, the Bełchatów mine is one of the largest excavations in Europe and it annually produces 38.5 million tons of lignite and 100-120 million m 3 of excavated overburden materials (Bednarczyk 2016). The excavated overburden is usually transported by trucks or conveyer belts to an external dumping area (external heap) or may be used as backfill material within the mining pit (in-pit heap). The structure of a spoil material significantly differs from its original state (i.e. preexcavation) and the physical and engineering properties can be extremely variable, both spatially and temporally.
The large amount of spoil imposes a significant economic and environmental liability on lignite producers. The occurrence of movements or catastrophic failures in spoil heap slopes can lead to major interruption in production operations, in addition to potential disastrous environmental consequences. Moreover, post-mining reclamation is required so that excavated lands can be restored for different uses, such as forestry, agriculture, or leisure activities. Thus, the long-term stability of the slopes is of great importance to ensure a safe and sustainable reclamation and closure process. In order to understand both short-and long-term stability of spoil heaps, one needs to properly estimate the physical and mechanical characteristics of spoil materials. This paper provides an overview of the characteristics of materials typically found in spoil heaps, including new datasets obtained from European lignite mines and existing data available in the literature, with special attention given to the variability of geotechnical characteristics. In addition, typical failure mechanisms and case histories of large instabilities in lignite spoil heaps are discussed, and the advantages and shortcomings of slope stability analysis methods applied to spoil heaps are reviewed. Modelling techniques that are usually applied for stability of spoil heaps are also critically discussed based on a review of examples within the literature that use limit equilibrium analyses (Nguyen et al. 1984;Steiakakis et al. 2009;Ulusay et al. 1996), numerical models (Poulsen et al. 2014;Richards et al. 1981), and probabilistic analyses (Cheng and Usmen, 1987;Chowdhury et al. 1986;Nguyen et al. 1984). Finally, important research gaps and key shortcomings in the methods of design and analysis of spoil heap slopes in lignite mines are summarised.

Sites of interest: materials and methods
The present work includes new datasets obtained from the Bełchatów mine in Poland, the Vrsany mine in the Czech Republic, and the South Field mine/Kardia mine in Greece. The Bełchatów mine is one of the largest surface lignite mines in Europe. It is located in central Poland, about 150 km west of Warsaw, and it ranks first in the country in terms of the percentage of lignite extraction, accounting for about 65% of the total Polish lignite production (Widera et al. 2016). The Vrsany lignite mine is located in the North Bohemian coal basin near the town of Most, about 100 km north-west of Prague. The mine contributes a significant share of lignite in the Czech Republic's energy mix and its reserves (within existing mining limits) have the longest remaining life of any in the Czech Republic (Euracoal 2018). The South Field mine and the Kardia mine are two of the Ptolemais-Kozani lignite basin mines in Northern Greece, located more than 500 km north of Athens. In between these two mines, there is a huge disposal terrain, known as the Soulou spoil heap, that is about 5 km long, 150 m high, and several hundred metres wide (Zevgolis et al. 2018b). In addition to the above, data from the Turów mine in Poland (Borecka 2007), the As Pontes mine in Spain (Varela et al. 1993), and the Assarel mine in Bulgaria (Kaltchev and Germanov 2005) are compiled in the present study.
The Bełchatów mine samples were superficial materials collected from different levels of the external spoil heaps. Samples from the Vrsany mine and the Soulou spoil heap were cored from borehole samples taken from a variety of depths. Tests on samples from Bełchatów and Vrsany were conducted according to ASTM standards for: particle-size distribution (gradation), ASTM D6913; Atterberg limits, ASTM D4318; water (moisture) content, ASTM D2216; and density (unit weight), ASTM D7263. Tests on samples from the Soulou spoil heap were conducted according to Greek national standards, which are based on the previously mentioned ASTM standards, as well as: engineering classification, ASTM D2487; and consolidated undrained triaxial compression test for cohesive soils, ASTM D4767.

Physical properties of spoils
The type, depth, and geological setting of the overburden material are mainly functions of the geological age and diagenetic origin of the lignite and may significantly differ from one place to another. The mineralogy of the spoils depends on the geological settings of the lignite-bearing area. For many basins, it is possible to consider a generalised stratigraphy of the area. For example, Piwocki and Ziembińska-Tworzydło (1997) provided a generalised stratigraphic column of the main lignite resources of Poland, as depicted in Fig. 1a. Figure 1b presents the general lithology of the largest lignite basin in Turkey (Ural and Yuksel 2004). Studies have reported the presence of magmatic and volcanic rocks (Kasiñski and Piwocki 2002), as well as carbonated metamorphic rock such as limestone in the overburden formation (Świtoniak et al. 2013).
The generalised stratigraphy can therefore help in identifying the general mineralogy of the spoils. The dominant mineral group in coal mine spoils are silicate aluminium clay minerals (e.g. kaolinite, illite, montmorillonite, smectite), while other common minerals are present, such as quartz, micas (biotite, muscovite), carbonates (e.g. calcite, dolomite), feldspars and sulphides (Li et al. 2014). In addition, small amounts of organic matter have been reported in spoils of different areas such as Poland and Turkey (Ediger et al. 2014;Kasiñski and Piwocki 2002). Table 1 presents a summary of the physical properties of spoil soils of lignite mines from the present study as well as data from other mines from around the world. It can be seen that there is significant variability in the reported physical characteristics of spoil soils. In fact, not only do the physical properties and mineralogy of samples from the various mines differ, but there is also a considerable variability within a single spoil heap. The spoils can contain a wide range of rock fragments and fines. This reflects the randomness in the dumping process and the fact that different methods have been used for excavation and transportation of the spoils. The inhomogeneous and heterogeneous nature of spoils can be observed in images of cores taken from the spoil heaps of the Turów (Fig. 2a, b), Vrsany (Fig. 2c), Bełchatów (Fig. 2d), and South Field/Kardia (Soulou spoil heap) (Fig. 2e) mines.
There is also great variability in grain size distribution of spoil materials reported in the literature. The grain size Fig. 1 a Generalised lithology of major Polish lignite resources: Adamów, Lubstów, Sieniawa, Turów and Bełchatów mines (after Piwocki et al. 1997;Widera 2013); b generalised lithology of Afsin-Elbistan lignite basin in Turkey (Ural and Yuksel 2004) Kaltchev and Germanov (2005) distribution of a number of samples from the Vrsany and Bełchatów mines is presented in Fig. 3a. It can be seen that the spoils from Vrsany contain a large proportion of fines with a clay content of 40-60 wt%, while samples from Bełchatów are mainly sands with occasionally low clay contents. The Atterberg limits of the spoils may be directly related to the presence of fines in the spoils, as evident in Table 1. The plasticity index of samples from Vrsany shows that both low plasticity and high plasticity clay are present ( Fig. 3b). In order to facilitate the interpretation of the material characteristics of spoils with regard to particle size, the United States Department of Agriculture (USDA) textural classification system was used. This classification does not provide a complete picture of the geotechnical characteristics of soils (it is mainly used in agricultural and forestry sciences), however it has been widely used and provides a good indication of the relative characteristics of the spoils. For example, Bełchatów spoils have been identified as mainly loam and clay loam (Świtoniak et al. 2013), while Turów mine spoil heap samples presented a wider range of textural classes, as depicted in Fig. 4a (Borecka 2007). On the other hand, it can be seen from Fig. 4b that spoils from different mines can have very different textural categories. Specifically, the samples from Vrsany mine were determined to fall within the clay and silty clay category. In order to visualise the extreme spatial variability within a single spoil heap, Fig. 5 presents the USCS (Unified Soil Classification System) classification of 100 samples taken from different depths of ten boreholes in the Soulou spoil heap. It can be seen that the spoil material varies significantly in both horizontal and vertical directions and no trend can be easily detected between the soil types at different locations. Nonetheless, it seems that silty material (MH, ML) is the most dominant. Furthermore, Table 2 provides statistical parameters of index properties of the Soulou spoil heap. As shown in the table, on average, fine-grained material (silts and clays) dominate (54.7%) within the spoil heap body. In terms of Atterberg limits, mean values of liquid limit and plastic limit were found to be 50 and 35%, respectively, with corresponding coefficients of variation (COV) equal to 26% and 34%. These values fall within the upper range, or even exceed, the usual COVs that are reported in the literature for in situ ground material. For instance, in general for Atterberg limits, Lacasse and Nadim (1996) reported COV values in the range of 3-20 (for clays), while Phoon and Kulhawy (1999) reported values in the range of 6-30 (for silts and clays). Table 2 also provides information for the bulk density and moisture (water) content. As far as unit weight is concerned (with a mean value of 16.8 kN/m 3 ), the obtained scatter is also larger than that usually reported in the literature for in situ soils. More specifically, while COV values typically range from 0 to 10% (Baecher and Christian 2003;Harr 1987), in the case of the Soulou spoil heap it was 14%. Larger than usual COV values were also found for water content; for example Phoon and Kulhawy (1999) reported COV values between 8 and 30 for clays and silts, whereas for the Soulou spoil heap the COV was computed as 46%.
Another source of extreme variability in spoils is the presence of materials that are deposited along with the excavated It is important to note that the presence of coal particles in spoil heaps have not been frequently reported, but Borecka (2007) indicated coal contents between 0 and 65 vol% in samples from spoil heaps of the Turów lignite mine. Borecka (2007) also reported a typical organic content of 1-10 vol% with some local anomalies of up to 65 vol%, while other studies have reported much smaller organic contents (e.g. Tripathi et al. 2012). These high coal and organic contents in the Turów mine spoil heaps are mainly because the waste materials from the Turów power station are also dumped on the heap along with the excavated overburden. The presence of power station wastes can have significant implications to slope instability due to their inherent moisture content, particle size, geomechanical properties, and geochemical/weathering reactivity.
The erosion and weathering processes are functions of several parameters, such as the heap's slope angle, the presence or lack of vegetation on the slopes, and the mineralogy of the soil combined with surface water and groundwater conditions. Slope angle has an obvious effect on its stability, but also affects vulnerability to erosion and weathering. It is also widely known that vegetation may have a significant influence on providing additional stability to soil slopes (Greenwood et al. 2004). Lastly, different minerals have varying solubility levels and hence basic or acidic environments can develop within the slope. For example, sulphide minerals can increase the acidity of the soil (Nyssen and Vermeersch 2010), while carbonate minerals act as a buffer, hence retarding the acidity (Rivas-Pérez et al. 2016;Schreck 1998). Studies have reported pH levels as low as 2.5 (Schreck 1998) which can speed up the mineral degradation of the spoil and hence endanger the long-term stability of the slopes. Rivas-Pérez et al. (2016) presented the results of a 20-year study on chemical characteristics of an afforested coal mine heap and showed that the acidity of the soil and groundwater may significantly fluctuate over time. They observed that the pH of groundwater generally decreased slowly over time, however in some cases, the level of  1 3 Page 7 of 18 505 acidity remained the same over the 20 year period. Acidification of the groundwater can lead to dissolution of carbonates and clay minerals, resulting in changes to hydro-mechanical properties of the spoil (e.g. porosity, permeability, and cohesion). Eventually, such an event could possibly impact the stability of slopes.

Geotechnical properties of spoils
Hydraulic and geomechanical properties are of great importance when considering the stability of spoil heaps. Table 3 presents a summary of the data from this study and those available in the literature on permeability and shear strength of the spoil materials of lignite mines. It is evident that strength parameters (cohesion and friction angle) vary considerably even among samples from the same site. As far as post-yield strength parameters are concerned, these can be significantly reduced from peak values (residual cohesion about 50% of the peak value). This strain-softening behaviour often characterises brittle soils and indicates a potential for progressive failure types. The extensive dataset from the Soulou spoil heap provides an opportunity to consider some interesting statistical parameters from samples obtained from a single spoil heap, as presented in Table 4 and based on 27 consolidated undrained triaxial compression tests with pore pressure measurement (CU-PP) (PPC 2010). Results refer to the mean and median value, the standard deviation and coefficient of variation, and the so-called characteristic value (Eurocode 7). The latter is taken as half a standard deviation below the mean, as suggested by Schneider (1997) and widely adopted within the Eurocode 7 framework. In terms of friction angle, mean and standard deviation values were computed as 25.9° and 7.2°, leading to a characteristic value equal to slightly larger than 22°. The COV of the friction angle was 28%. This value is very high compared to coarse-grained soils, and it is towards the upper limit of COVs reported for finegrained soils (Baecher and Christian 2003; Lacasse and Nadim 1996). Cohesion has an extremely large COV (143%) and a large difference between the mean (14.2 kPa) and the median (6.6 kPa) values. These results are due to the fact that nearly half of the measurements of cohesion (11 out of 27) were found equal to or very close to 0. From a design point of view, the characteristic value of cohesion is 4.1 kPa.
Studies have shown that the strength parameters and permeability of soil can be related to the proportions of sand, silt, and clay. Knoll and Knight (1994) studied the permeability and porosity of a sand and clay (kaolinite) mixture and showed that as clay content increases, the porosity and permeability decrease up to a limiting value, then increase as further clay is added. This limiting clay content determines the boundary between clayey sand and sandy clay (see Fig. 6a, b). Monkul and Ozden (2007) conducted a series of direct shear and oedometer tests to investigate the effect of kaolinite content on compressional behaviour of soils. The results of the direct shear tests showed that adding to the clay content does not change shear strength up to a particular point, after which shear strength decreases (Fig. 6c). The point at which the shear strength starts to decrease represents the transition from clayey sand to sandy clay. Their oedometer tests showed that the compression indices of soil increase as clay content increases; mainly because as clay content increases, the number of contacts between coarse grains decreases (Fig. 6d). In addition, Rybicki and Wozniak (2010) . 6 The effect of clay content on geotechnical properties of sand: a porosity (Revil et al. 2002); b permeability (Revil et al. 2002); c shear strength (Monkul and Ozden, 2007); d compression index (C c ) and granular compression index (C c−s ) (Monkul and Ozden 2007) discussed the significance of the lumpy structure of spoils in their hydro-mechanical behaviour and concluded that spoils are in fact dual porosity materials (inter-lump and intra-lump porosities). Studies have shown that the hydraulic and geomechanical properties of spoil are also related to their age and the reclamation activity undertaken. Varela et al. (1993) reported that permeability of superficial soils of spoil heaps increased from 8 to 20 mm/h over a period of 5 years (porosity changed from 30 to 59%). In addition, Tripathi et al. (2012) reported an increase in hydraulic conductivity of spoil soil from 0.31 to 1.21 mm/h over a period of 12 years after afforestation, while porosity increased from 12 to 18.2%. Tripathi et al. (2012) also showed that cohesion of the spoils increased by 100-150% over 12 years of vegetation, while the angle of friction increased by 2-3 degrees. These observations may suggest that, although consolidation of the spoil heaps increases the strength of soil by reducing the porosity, other processes (e.g. weathering, mineral dissolution) can deteriorate the texture or structure of the soil, which in turn results in significant increases in porosity and permeability.

Mechanisms of instability associated with spoil heaps
The mechanism of slope failure in any soil depends on its physical, hydraulic, and mechanical characteristics, as well as climate conditions and slope geometry. As shown in Tables 1 and 3, physical and mechanical properties of spoil vary, which can influence the type of instability and its consequences. Spoils are initially composed of loose rock fragments, sand, silt, and clay particles, but their physical, hydraulic and mechanical properties may change over time, mainly due to the changes in loading and geometry of the slope, rainfall infiltration, weathering, vibrations and earthquakes, and post-mining reclamation activities. It is also important to note that the stability of both external and internal spoil heaps depends on the geotechnical performance of the foundation and the mine pit. Thus, it is important to identify the modes of expected failures during both the short-and long-term, so that suitable analyses and design methodologies can be carried out. Varnes (1978) identified five general instability types in slopes: falls, topples, slides, spreads, and flows. In light of this classification, the instabilities of spoil heaps in lignite mines can be better reviewed. Table 5 presents an overview of those instabilities reported in the literature. Unfortunately, the slope instabilities of spoil heaps are not generally welldocumented and there are not many comprehensive datasets, but some general conclusions can be made using the data presented in the table, as well as others reported in the literature. Almost all of the reported instabilities indicated slide and spread mechanisms; mostly rotational and occasionally translational (or their combination). Although there are some reports of flow mechanisms in the literature (e.g. Wichter 2007), no data could be found to confirm them. The size of the reported sliding wedges varies significantly, even within the same heap, which is believed to be due to the inhomogeneity and anisotropy of the hydro-mechanical properties of spoils. The moving mass in the Turów mine (Poland) was 1300 m long and 750 m wide (Bednarczyk 2016), while in the Southfield mine (Greece) it was 1100 m long, 550 m wide, and 30-90 m deep (Steiakakis et al. 2009).
Important considerations of spoil heap instabilities are their total displacement and the rate of displacement, as these indicate the possible adverse consequences. Total displacements of up to 600 m are presented in Table 5, some of which resulted in severe damage to mining facilities and disruption of mining activities. Only a few studies presented the movement rate for the whole period of the failure process, mainly due to lack of monitoring instrumentation. The data reported by Steiakakis et al. (2009), which seem to be the most comprehensive movement data available in the literature, was obtained by installation of monitoring stations after  Nyssen and Vermeersch (2010) the first signs of instability in the spoil heaps (see Fig. 7a).
Large displacement rates were observed in the initial stages of the failure, which slowly reduced until the sliding wedge completely stopped. However, it should be noted that all these displacements were recorded during or after mitigating measures were put in place, without which the displacement rates could have been larger (especially at the later stages). The time required for the wedge to stop sliding is notable (approximately 2 months of continuous sliding, Fig. 7a). As seen in Table 5, failures in spoil heaps usually occur after periods of heavy rainfall. An example of a rainfallinduced failure is from Central Anatolia in Turkey, where cumulative rainfall was almost at a peak when a spoil instability occurred (Fig. 7b). Kasmer et al. (2006) showed that the mean precipitation in May 2001, before the failure occurred, was considerably higher than the mean precipitation of previous years (nearly 30%). The drainage mechanisms in spoil heaps affect water seepage and consequently the effective stress distribution, which is vital when considering stability. However, the heterogeneity of materials and the possibility of impermeable layers within the heaps make evaluation of drainage characteristics very difficult. Some reports have mentioned the existence of drainage holes in spoil heaps, however there are no specific standards for drainage design in spoil heaps. Moreover, the hydraulic properties of the formation beneath the spoil heaps are of great significance. Steiakakis et al. (2009) identified that the seasonal water from adjacent karstic carbonate rocks infiltrated into the aquifer beneath the spoil heap before failure. They estimated a water pressure of 1.0 MPa was built up during the heavy rainfall period which transmitted into the heap and led to a large instability.
The geometry of spoil heaps can also play an important role in susceptibility of the slope to failure. The role of slope angle is rather intuitive, however the effect of height with respect to slope angle on the susceptibility to failure and instability mechanisms is not well-defined. Ulusay et al. (1995) presented a series of field data from numerous spoil heaps at the Eskihisar mine where they plotted the height of the heaps versus the slope face angle in relation to their instability mechanisms, as depicted in Fig. 8. For this specific mine, the data indicate a boundary between stable and unstable heaps. Although generalisations should be applied with caution, if the same type of data were to be obtained for other mine sites, some practical guidance of critical height Fig. 7 a The recorded movements of sliding mass in spoil pile of South Field mine in Greece (Steiakakis et al. 2009); b cumulative precipitation during the 6 months immediately before a spoil heap failure in Central Anatolia, Turkey (Kasmer et al. 2006) Fig. 8 The stability of the spoils considering the slope angle and height of the piles (adapted from Ulusay et al. 1995) could be defined based on the angle of the slopes (as well as other potential factors), which could assist with decisionmaking and land management processes. Nyssen and Vermeersch (2010) reported numerous spoil heap instabilities in Belgium that occurred following events of spontaneous combustion. There have been relatively few investigations on the effect of small-and large-scale spontaneous combustion events on the stability of spoil heaps. In countries such as Australia and South Africa, where dry and hot summers are common, spontaneous combustion is an expected occurrence (Carras and Leventhal 2000), whereas the relatively lower coal content and cooler climate in European lignite mines do not generally facilitate selfburning processes. However, a number of cases in European coal mines are reported in the literature (e.g. Hajra et al. 2009;Jelínek et al. 2015;Nyssen and Vermeersch 2010). The occurrence of spontaneous combustion is related to the content of coal and iron sulphides. Thus, the risk of combustion-induced instabilities in spoil heaps is directly related to the source of the material being dumped. In some cases, the tailings of coal washing plants and furnace slags of power stations are dumped alongside the spoils, which leads to a higher risk of spontaneous combustion and slope failure. For example, as discussed earlier, Turów mine spoils have a relatively high coal content, which may make them susceptible to spontaneous combustion if subjected to hot and dry seasons.

Methods of slope stability analysis
Continuum mechanics methods are generally used to evaluate the safety level of a slope. In other words, a spoil heap is evaluated, at least at a preliminary design stage, as a typical embankment. The limit equilibrium method is perhaps the most widely used method for slope stability analysis. Many limit equilibrium procedures are used, including Bishop (1955); Janbu (1954); Morgenstern and Price (1965); Spencer (1967), and Generalised Limit Equilibrium (Lam and Fredlund 1993) methods. Most of these are based on the well-known method of slices.
The main reason for the use of limit equilibrium methods in slope stability analysis is their relative simplicity and the confidence gained in their use over the years. Safety level is determined by means of a factor of safety (FoS), which provides a quantitative evaluation that is directly associated with an assumed failure surface (Duncan et al. 2014). Typically, FoS is defined simply as the ratio of shear strength available along a slip plane to the shear stress mobilised along the slip plane. Nguyen et al. (1984) proposed a limit equilibrium model for spoil heaps based on a two-wedge planar sliding failure mechanism. Based on their results, they provided useful slope stability charts associated with the most commonly observed failure mechanisms in Australian spoil heaps. Similar bi-planar wedge limit equilibrium models were developed for the case of Turkish spoil heaps which were subsequently utilised to investigate possible remedial measures for the stabilisation of identified incipient failure mechanisms (Ulusay et al. 1996). In Greece, Steiakakis et al. (2009) used limit equilibrium methods to conduct a back-analysis of a large-scale failure in an external waste dump of the South Field lignite mine in northern Greece.
The limit equilibrium framework has several disadvantages which limit its legitimate applicability for more complex cases. In particular, limit equilibrium disregards the developing stress-strain distributions inside the slope mass. As a result, different types of slopes formed by different construction procedures give the same FoS provided they are comprised of the same material and have identical geometric features (Abramson 2002). Pasternack and Gao (1988) found that a limit equilibrium analysis often yields a satisfactory FoS in terms of catastrophic failure, but unsatisfactory FoS in terms of lateral movements. In addition, time-dependent changes in hydro-mechanical processes (e.g. rainfall infiltration in unsaturated slopes and soil degradation) can significantly contribute to the development of unstable stress-strain states. Moreover, incorporation of extreme material variability, as observed in spoil materials, in limit equilibrium analysis (e.g. Low 1997;Pascoe et al. 1998) has been limited to Monte Carlo simulations of uniform slopes, which does not fully represent the variable and anisotropic nature of spoil materials.
The inherent limitations of limit equilibrium methods can be overcome by advanced numerical methods, such as the Finite Element Method (FEM) and/or the Finite Difference Method (FDM). The former is possibly the most versatile and broadly used numerical method in geotechnical engineering and is widely used for slope stability analysis. FEM can readily simulate the formation process of a slope, either via the sequential removal (cut slope) or the successive placement (embankment slope) of material into a finite element mesh domain. The actual stress history of the soil mass that composes the slope can be credibly represented if an appropriate selection of soil constitutive models and relevant parameters is used. Moreover, no assumption needs to be made about the shape, direction, and location of slip surfaces: failure occurs "naturally" through zones on which the soil's shear strength is unable to sustain the applied shear stresses (Diederichs et al. 2007;Griffiths and Lane 1999). Ideally, finite element analyses should be accompanied by a verification of results against experimental or case-history data. However, even in the absence of such data, FEM models that are carefully developed by experienced users may provide useful information for the evaluation of the stability of spoil heaps.
Traditionally, results of numerical analyses have been expressed in terms of displacements of the soil masses. However, over the last 15-20 years, the so-called shear strength reduction (SSR) technique has become an integral part of many numerical analysis programs. In an SSR analysis, a so-called strength reduction factor (SRF) is calculated by progressively reducing the shear strength of the material in order to bring the slope to the verge of failure (i.e. until non-convergence occurs within a specified number of iterations and tolerance in a numerical solution scheme). The ratio of the soil's actual shear strength at equilibrium to this reduced shear strength is the SRF. If φ′−c′ and φ′ r −c′ r are the original and the reduced parameters, respectively, then the SRF can be written according to (Dawson et al. 1999;Matsui and San 1992): As indicated by several researchers (Griffiths and Lane 1999;Zheng et al. 2009), the SRF is equivalent to the FoS in limit equilibrium methods. The greatest advantage of the SSR technique is that it provides the familiar quantitative index (i.e. the SRF) of safety and stability while working with FEM or FDM models. On the other hand, some shortcomings of the general application of the SSR technique include the use of a Mohr-Coulomb constitutive model (instead of advanced models that are able to capture more complex behaviour of spoil heaps) and the fact that resulting deformations are fictitious and have no physical meaning. The latter disadvantage is often overcome by also conducting regular finite elements analyses to evaluate displacements.
Bearing in mind all of the above, one can conclude that the unique properties and behaviour of spoils cannot be easily captured using limit equilibrium methods and that numerical methods need to be considered. There have been a number of efforts to numerically model the stability of spoil heaps. Richards et al. (1981) used a finite element model to simulate the progressive failure of a spoil heap, taking into account residual shear strength of soil, ensuing from soil strain-softening and wet-softening phenomena. They concluded that slope failure in the Goonyella mine spoil heap in Australia was initiated by weakening of a moisture-sensitive material at the base of the slope, and the subsequent cracking of bulk spoil material was due to settlement and compaction movements. Poulsen et al. (2014) utilised finite element models and a strength reduction technique to back-analyse a failure case and concluded that residual values of shear strength dominated the stability of the spoil heap.
Spoils are very complex geo-materials exhibiting significant anisotropy, both in terms of strength and stiffness. They are exposed to numerous environmental factors over a long period of time and therefore may sustain severe degradation through softening of inter-particle bonds, a phenomenon known as de-structuration (Ulusay et al. 1995). Moreover, time-dependent and strain-rate processes are usually present in spoil heaps due to continuous variations in the loading states and the geometric boundary conditions during the heap construction period (Holubec 1976). Load hysteresis and strain-hardening phenomena accompanied by high nonlinearity in small strain domains have also been associated with the spoil material (Rai et al. 2012). The uneven and sometimes unknown deposition processes, as well as lack of compaction effort during the placement of the spoils, lead to further complexities and uncertainties in relation to the behaviour of spoil dumps (Zevgolis 2018).
Consideration of these complex characteristics of spoil materials necessitates the need for advanced constitutive models that may not be necessary for classical engineered embankments. Various advanced models have been developed over the years that can adequately capture the complex features of specific material behaviour (Adachi and Oka 1982;Dafalias and Taiebat 2014;Martindale et al. 2013;Rezania et al. 2018;Sivasithamparam et al. 2015;Yin et al. 2002). Constitutive models strictly adapted to spoils are not found in the literature. Development of specialist models for spoils could enhance the prediction of spoil heap performance, and ultimately provide more refined and costeffective slope designs.
Another important direction of research could be the development of coupled hydro-mechanical models where the physico-chemical processes associated with spoils can be simulated. While many recent studies have considered some hydro-mechanical processes in slopes (e.g. Borja and White 2010;Choo et al. 2016;Davies et al. 2014;Hu et al. 2011Hu et al. , 2018Liu et al. 2017;Qi and Vanapalli 2015;Wu et al. 2016;Yang et al. 2017), they are more suited to natural slopes subjected to rainfall. For spoil dumps, the long-term stability is vital to the feasibility of land rehabilitation and reclamation plans. Therefore, slope stability models for spoils need to account for the long-term physico-chemical processes, as well as their impact on hydro-mechanical characteristics of saturated and partially saturated spoils, as discussed in Sect. 2.2 of this paper. Also, models should account for environmental conditions, such as cycles of dry/hot and wet/ cold seasons, and how these fluctuations may lead to local or global failures of spoil heaps.
In addition to the development of advanced constitutive models and coupled hydro-mechanical modelling approaches, emphasis should also be put on the value of numerical methods themselves, and their capability to simulate the complex engineering processes and failure mechanisms occurring in spoil heaps. For instance, Chowdhury et al. (2009) stated that in the lower parts of spoil heaps, near the foundation base, spoils are often weakened due to high levels of shear strains experienced during the construction period, a phenomenon which can lead to progressive slope failure (initiating from the toe and propagating to the crest). This progressive failure mode is often caused by the continual saturation of the spoil heap, due to seasonal rainfall and water infiltration within the spoil dump, which is made worse by inadequate or non-existing drainage infrastructure. Richards et al. (1981) observed similar progressive failure mechanisms in spoil heaps, pointing out that the strength of basal materials of heaps has a greater influence on the stability of the structure than that of the bulk of the spoil heap. In addition, due to the high levels of heterogeneity associated with spoil materials, more than one failure mode can occur within the same heap (Okagbue 1984). Hence, these different failure modes that may be either strongly or weakly related to each other, contribute individually to the development of an overall resultant failure mechanism (Chowdhury and Xu 1994). These cases of composite engineering response-selected among several others-can only be approximated by sophisticated computational methods, such as FEM. These issues again highlight the importance of developing advanced constitutive models for spoils which can be incorporated into FEM and FDM models.

Probabilistic analysis of slopes
Regardless of the adopted method, there will be uncertainties associated with any slope stability analysis. Due to their extreme variability and heterogeneity, the level of uncertainty involved in the design of spoil heaps is much higher than for conventional embankments. In practice, FoS evaluation is performed deterministically, resulting in a single value that is meant to reflect the complex evaluation of risk associated with slope failure. These evaluations are often subjectively influenced by regulations, company practices, or an individual's judgment. The FoS value alone cannot adequately express the inherent risk of slope failure, because even a large value of FoS does not mean that the probability of failure is zero (Baecher and Christian 2003). Deterministic slope stability analyses do not explicitly account for quantified uncertainty and instead rely on conservative parameters and designs to deal with uncertain conditions. However, seemingly conservative designs have not always been immune to failure (El-Ramly et al. 2002). In addition, the economic and safety impact of such subjective and conservative design approaches cannot be directly evaluated. In contrast, probabilistic slope stability analyses allow for the uncertainty of input variables to be quantified and incorporated rationally in the design process. As a result, probabilistic analyses provide, in addition to the conventional safety factor, outcomes that can contribute to design optimisation.
Using probability theory, the reliability of a spoil heap (i.e. the probability of a specified performance threshold) may be obtained, hence providing more assurance and confidence to the geotechnical engineer (Ang and Tang 1984). However, consideration of geotechnical variability, and interpretation of respective results, should be made with caution. Probabilistic analysis should not be viewed as a substitute to the well-established deterministic models in geotechnical engineering, but as a complementary element. Knowing both an approximate value of slope safety factor, and an approximate value of its failure probability, is better than knowing more precisely either one alone (Duncan 2000). Comprehensive risk analyses may also be conducted by associating the probability of unsatisfactory performance with potential accompanying consequences. In such cases, risk analysis may be utilised for optimisation of design, securing the sustainable operation of the engineering project (Zevgolis et al. 2018a).
Various probabilistic methods have been proposed in the literature, each with particular advantages and disadvantages. The most commonly employed probabilistic methods are the first-order second moment (FOSM), the second-order second moment (SOSM), the first-order reliability method (FORM), the point estimate method (PEM), and the Monte Carlo simulation (MCS). Geotechnical uncertainty and variability is explicitly incorporated in the above methods by considering critical geotechnical parameters as random variables with assigned probabilistic and statistical parameters. Depending on the adopted method and the available data, parameters such as the mean value and standard deviation (or coefficient of variation), the probability density function, the spatial correlation length, and the coefficient of cross-correlation (between different random variables) may be taken into consideration. A performance function is then defined, usually in terms of limit states (i.e. any condition that would lead to unsatisfactory performance of the slope, including slope failure and/or excessive deformations). Ultimately, one of the above methods is employed in order to evaluate the performance function in terms of probability of failure (or probability of unsatisfactory performance). A detailed description of the individual characteristics of each method is not within the scope of this paper; only a concise outline is given here.
The FOSM employs the first terms of a Taylor series expansion of a performance function to assess both the mean value and the variance of the performance function. The variance is the highest order statistical moment possible to be used and derived by the method. The main shortcoming of the method is the need for the evaluation of the performance function's partial derivatives, which in cases of complex problems could be cumbersome or even impossible. The SOSM utilises the second-order terms of a Taylor series expansion of a performance function, so it is more precise than FOSM. However, due to the added complexity of calculations, the method has not been widely employed. The PEM is an approximate, yet often sufficiently accurate, method for the evaluation of the lower order statistical moments of a performance function. Among the most prominent advantages of PEM are its simplicity, robustness, explicitness, convenience, and accuracy for a moderate number of random variables. In addition, the method can be easily used in combination with numerical models. The FORM is an enhanced version of the FOSM, evaluating a performance function's reliability index by means of the geometric representation of the safe distance between the expected performance and critical states. It is generally considered more accurate than the FOSM, and the greater sophistication accompanied does not necessarily comprise a large constraint for its utilisation in practical geotechnical problems. Last, the MCS is a stochastic technique that generates a large number (usually on the order of tens of thousands) of repeated simulation processes (realisations). Each realisation is based on the generation of a series of values of one or more random variables. The procedure requires complete definition of the random variable probability distributions, but, through simple computations, it provides empirical outcomes of numerically simulated random realisations of the performance function. Statistical estimates can then be obtained for the probability of failure and the safety factor.
Research thus far conducted on the employment of probabilistic analysis methods in the evaluation of spoil heap slope reliability is limited. Nguyen et al. (1984) used a twowedge limit equilibrium model and evaluated the reliability of several spoil heaps in Queensland, Australia, employing the Monte Carlo simulation (MCS) method and the point estimate method (PEM). Excellent agreement was obtained between the two methods, and the particular value of PEM for practical design purposes was provided, mostly due to its great efficiency and satisfactory accuracy. Chowdhury et al. (1986) modified their two-wedge limit equilibrium model to simulate progressive failure phenomena in spoils. Cheng and Usmen (1987) employed PEM and MCS to account for the large variability of geotechnical parameters in spoils to assess the safety of an embankment consisting of coal refuse material. The major drawback common to these early works is the use of very simplified forms of the performance function based on limit equilibrium methods. This could be improved by the use of more sophisticated numerical methods and advanced constitutive models, implemented within a probabilistic framework. Efficient probabilistic methods, such as PEM, have been conveniently coupled with numerical methods for slope stability analysis (Ahmadabadi and Poisel 2015;Ng et al. 2014).
Statistics of geotechnical parameters can be evaluated based on laboratory and/or field data. However, this is largely a manual and often laborious process, especially for tasks related to the verification of data consistency or identification of discrepancies or trends within the data. The implementation of data within a probabilistic analysis framework is even more complex and requires significant time and computational resources. In addition, practitioners have been reluctant to adopt probabilistic techniques, mainly due to a lack of familiarity with probabilistic theories and the absence of a link between conventional deterministic and probabilistic assessments (El-Ramly et al. 2002). While researchers have demonstrated the benefits of using probabilistic analyses for slope stability problems (e.g. Cho 2007; Griffiths and Fenton 2004;Griffiths et al. 2009;Griffiths and Lane 1999;Huang et al. 2017;Ji et al. 2012;Jiang et al. 2015;Srivastava et al. 2010;Suchomel and Mašı´n 2010), more work could be done to develop tools to help professionals with statistical processing of data and with the implementation of a probabilistic framework.
A major area of interest relating to the improvement of probabilistic analysis of slope stability is the implementation of spatial variability. Usually, the random field theory (Vanmarcke 1983) is employed to consider the spatial variability of soil properties. An auto-correlation function is used to describe the spatial distribution of the properties, where the value of a property at every spatial point is assigned based on the value of adjacent spatial points (i.e. conditional random field). While the use of auto-correlation functions for a single random geotechnical property is relatively straightforward (e.g. Andrade et al. 2008), difficulties arise when two or more geotechnical properties are varied. It is well known that geotechnical properties are related to each other, and thus it is not appropriate to assign independent random fields for each property. On the other hand, the relationship between properties through linear or non-linear regression is rather simplistic if not incorrect. As a result, a cross-correlation structure between random properties needs to be considered. Slope stability analyses have adopted this approach (e.g. Griffiths et al. 2009;Javankhoshdel and Bathurst 2015), where simplified cross-correlation techniques are used to generate random fields of cohesion and friction angle considering only lognormal distributions for both parameters. For extremely variable spoils, however, the dependency of the parameters can be much more complex, as three or more random dependent variables need to be generated. Therefore, application of more complex cross-correlation techniques in multivariate geotechnical random fields (e.g. Zhu et al. 2017) seems necessary in probabilistic analysis of spoil heaps. Such an approach would not only improve the accuracy of the probabilistic analyses, but also increase their efficiency as a smaller number of simulations are required when simultaneous random variables are used, as opposed to using only one random variable at a time.

Concluding remarks and identified research gaps
The problem of slope stability in spoil heaps is a significant challenge for lignite mining operations. Instabilities of lignite mine spoil heaps results in significant environmental and economic liabilities and thus prediction and prevention of such problems must be a priority. However, prediction of these instabilities can be very difficult owing to the significant spatial and temporal variability of the spoil properties. This paper provided an overview of the spoil properties and the mechanisms by which these properties can change. The main conclusions drawn from this study are as follows.
• The physical and geotechnical properties of spoils are variable, mainly due to anthropogenic impacts, i.e. the excavation and dumping methods. The results of index testing showed a wide range of material properties that can significantly influence the stability of slopes. The extensive dataset from the Soulou spoil heap demonstrated the randomness in spatial distribution of material types and the extent of variability within a single spoil heap, with coefficients of variation of key material parameters all being larger (and often significantly larger) than those reported for in situ soils. • Design guidelines should ideally address the effects of the variability of spoil materials. A probabilistic approach may be best suited to ensure that the spatial variability of spoil characteristics is properly considered. However, probabilistic analyses should be able to account for the inter-dependency of geotechnical properties of spoils within a spatial distribution; that is multivariate random geotechnical fields. • The long-term stability of spoil heaps is vital to safe restoration activities and therefore design guidelines should consider the time-dependent changes in hydro-mechanical properties of spoil materials. As such, the mineralogy of spoil is important because geochemical reactions with groundwater can lead to significant changes in the mechanical properties of spoils. While the consolidation of spoils can enhance the stability, chemical changes can deteriorate strength. On the other hand, the consolidation process leads to a lower permeability which in turn can increase the potential for pore pressure build-up and endanger the stability of slopes. A better understanding of the time-dependent chemo-hydro-mechanical processes within spoil materials would enable better longterm predictions of spoil heap stability, as indicated in the case of the lignite mine located in As Pontes, Spain (Rivas-Pérez et al. 2016). • Spoils are a special kind of anthropogenic soils and the classical assumptions made for natural soil slopes should not be directly applied to spoil dumps. This becomes more significant given the time-dependent changes in geotechnical properties of spoils and their potential impact on long-term processes. Therefore, the use of simplified limit equilibrium analyses to assess the reliability of spoil dumps may be inappropriate and lead to erroneous results. Advanced numerical models, in which the characteristics of spoils are adequately implemented, need to be employed in order to ensure higher confidence levels and better predictions. Although various advanced models of soil behaviour have been developed over the years, (which adequately capture the complex features of specific material behaviour), no constitutive models specifically adapted to spoil materials are found in the literature. • Review of reported instabilities in lignite mine spoil heaps showed that heavy rainfall and resulting increases of excessive pore pressures within/beneath the heaps was responsible for the majority of instabilities. Specification of drainage requirements within spoils heaps that could be implemented during and/or after construction could have a beneficial effect on heap stability. In addition, slope stability analyses should be able to account for the impact of weather conditions by simulating rainfall infiltration and/or groundwater seepage in and around the heterogeneous and anisotropic spoil materials within an unsaturated soil mechanics framework.