Consolidation behavior of various types of slurry tailings co-disposed with waste rock inclusions: a numerical study

The co-disposition of mine tailings and waste rock in tailings storage facilities (TSFs) could contribute to increase the consolidation rate and decrease long-term settlement of tailings. Non-linear change of tailings properties during the filling process and interaction between tailings and waste rock inclusion (WRI) are critical to mechanical analysis but can, however, be complicated to simulate. The question of net volume gains or losses of tailings was also raised. In this study, fully coupled analysis which considered continuous variation of hydraulic conductivity and stiffness of tailings were performed to assess the evolution of consolidation of various tailings types in the presence of WRI. Effects of volume ratio of tailings over WRI on the net volume was investigated. Finally, effect of several practical aspects such as filling rates, and instantaneous filling assumption were considered. Results indicated that WRI could increase by 3.3 times the rate of consolidation of tailings. The zone of influence of WRI on tailings consolidation varied for each tailings. Using updated properties showed significant differences compared to models with constant values. The application of WRI can reduce volume available for the storage of tailings and net volumetric change due to settlement of the tailings with or without WRI could be estimated explicitly. Equations predicting evolution of net volume with the changes in the volume ratio of tailings and WRI were accordingly proposed. WRI effects were more pronounced with the increase of the filling rate. Finally, instantaneous filling assumption had little effect on the simulated rate of consolidation.


Introduction
Tailings generated from mining activities are typically deposited hydraulically (i.e., as a slurry) in tailings storage facilities. Consolidation of mine tailings can take a very long time because of their initially high-water content and their relatively low hydraulic conductivity (Vick 1990;Blight 2010). Progressive deposition of slurry tailings in TSFs might lead to the build-up of excess pore water pressure (PWP) that has been linked to various geotechnical instabilities issues, including failure of tailings impoundments and static/dynamic liquefaction of tailings (Azam and Li 2010;Morgenstern et al. 2016;Robertson et al. 2019). A novel approach of tailings and waste rocks co-disposal was developed for surface tailings storage facilities that could contribute to accelerate tailings consolidation. In this method, waste rocks are placed as linear inclusions within the tailings and raised at the same time (Fig. 1). This technique can contribute to accelerate the dissipation rate of excess PWP generated during the continuous accretion of tailings (and thus speed up the consolidation of tailings) by promoting horizontal movement of pore water (Bolduc and Aubertin 2014). WRI can also provide other advantages such as reducing the size of waste rock piles (i.e., WRI could contribute to completely remove the necessity to build a waste rock piles for underground mines) (Aubertin 2013(Aubertin , 2018Saleh-Mbemba and Aubertin 2021a). Strength and stiffness of tailings can consequently be improved more quickly. The increase in the consolidation rate and magnitude of tailings in the presence of WRI were experimentally evaluated using physical laboratory tests (Saleh-Mbemba and Aubertin 2021a, b), while consolidation behavior of tailings in the vicinity of WRIs was also numerically investigated at the large scale (Jaouhar et al. 2013;Bolduc and Aubertin 2014;Boudrias 2018). WRIs can also help to provide mechanical reinforcement to the tailings impoundment enhancing both static and dynamic stability of the retaining structures (Pépin Nicolas et al. 2012;Ferdosi et al. 2015). Effect of WRI configuration was evaluated by changing various values of width of the inclusions and center-to-center spacing between these inclusions (Ferdosi et al. 2015). In general, results indicated that downstream slope deformation tended to decrease with the increase of the number of inclusions and to increase with the increase of spacing between inclusions (Ferdosi et al. 2015).
Critical concerns when studying consolidation of tailings with WRI include capturing properly the evolution of tailings consolidation characteristics and the acceleration of dissipation of excess PWP in tailings. A balance should also be found between maximizing the quantities of tailings deposited and ensuring a good performance of WRI by reducing the distance between them. WRI design optimization, thus, requires to consider net volume gains or losses due to the addition of waste rocks within the tailings and thus volume available for the storage of tailings afterwards (Bolduc and Aubertin 2014; Saleh-Mbemba and Aubertin 2021b). Several factors can pose uncertainties to the efficiency of this technique including the filling rate and hydro-mechanical properties of tailings and waste rocks (Bolduc and Aubertin 2014). Finally, comparative research on the behaviors of different types of tailings applying co-disposal technique can accordingly contribute further to the application of this technique to management operations of various types of tailings and somewhat improve optimization of co-disposal technique as a function of tailings properties.
The objective of this paper was, therefore, to evaluate numerically the effect of WRI on various types of tailings sequentially filled in a simplified TSF and considering the evolution of consolidation properties. Potential differences of models with constant and updated properties when applying WRI were subsequently investigated. Effects of volume ratio between tailings and WRI on the volume gains or losses of tailings in the impoundment were studied. Influences of the filling rate on the consolidation rate of tailings were examined. Finally, influences of the assumption of instantaneous filling in simulations (i.e., no consolidation of tailings occurs during the filling process) compared to the progressive filling scheme in practice (i.e., consolidation occurs during the filling process) was studied. Limitations of the models were also discussed in this study.

Methodology
Many multi-dimensional codes have been developed to investigate the settlement of slurry tailings under sequential filling (Liu and Znidarčić 1991;Fredlund et al. 2015) and the explicit finite volume code Flac3D (Itasca Int. 2021) is one of the most frequently used to simulate multi-dimensional and non-linear consolidation of tailings materials. For example, Shahsavari and Grabinsky (2015) used Flac3D to investigate the effect of non-zero boundary condition on the Fig. 1 Conceptual configuration of a simple WRIs system in a tailings impoundment with a only tailings and water can only move upward and b tailings codisposed with WRI and water can move either vertically or horizontally due to the presence of WRI, reducing the length of the drainage pathway (modified after James (2009)) consolidation of the cemented paste backfill layers, while Zhou et al. (2019) simulated mine tailings consolidation in an impoundment with the presence of wick drain systems. In this study, a series of simulations were performed using Flac3D to model tailings placed sequentially in a simplified TSF with WRI. Common practice for TSFs design usually considers 1-D geometries (Fredlund et al. 2015). 1-D models are fast and particularly suited for cases where impoundments have small depths compared to the width and the length of the retaining structures. In these cases, the fluid flow direction and the consolidation are essentially vertically oriented (Priestley 2011;Fredlund et al. 2015). However, such assumption may be not applicable when highly permeable regions and horizontal drainage systems (such as WRIs or wick drains) are installed in TSFs. In these cases, 2-D or 3-D analysis are required to capture the principal effect of lateral flow and strain (Coffin 2010;Fredlund et al. 2015).

Conceptual models
A 100 × 5 m section was modelled to simulate the consolidation of tailings in the TSF and the effect of the WRI (Fig. 2). The dimensions of the model were chosen to reduce the mesh size and computation time, yet remaining representative of field conditions. This width of the TSF was assumed to be far enough from the dams so they would not affect simulated settlements. The influences of the volume ratio between tailings and WRI on the calculated results were also investigated. The simulated WRI represented a simple central inclusion whose crest width was around 12 m. WRI geometry was simplified as parallelepiped instead of the typical successive trapezes usually constructed on real mine sites. This assumption has shown no influence on the results but reduces numerical instabilities (Bolduc and Aubertin 2014). Saturated tailings were filled into the TSF in successive 10 layers over a period of 10 years. Each layer was 3 m thick for a total height of 30 m. Tailings deposition was assumed instantaneous (i.e., a new 3 m-thick layer was added every year at once and left to consolidate before a new layer was added one year later). Four types of tailings were modelled in this study (see below). WRI was raised sequentially at the same rate as tailings (similar to field conditions).
Displacements at the bottom of the TSF were fixed in all directions, and an impervious bottom was considered. Side boundaries of the models were allowed to move vertically only and also considered impervious (i.e., axis of symmetry). Zero PWP was assigned at the surface of the model to simulate groundwater table and to allow upward movement of water in the tailings. The deposition of a new layer of tailings induced the generation of excess PWP, followed by the dissipation of PWP (towards the surface and the WRI) and the consolidation of the materials.
A parametric analysis was also conducted to evaluate the impact of filling rate and effect of the filling scheme.

Material properties
Four different tailings were considered in this study: 2 types of gold mine tailings (tailings A and B), uranium mine tailings (tailings C) and bauxite mine tailings (tailings D). Tailings properties (i.e., unit weight, relative density, friction angle, and relations of void ratio-effective stress and hydraulic conductivity-void ratio) of tailings A were obtained from experiments by authors, while those for other mine tailings were obtained from literature (Somogyi 1976;Bhuiyan et al. 2015;Lévesque 2019). Unit weight was highest for tailings A (e.g., around 19.7 kN/m 3 ), while that for tailings B, C, D was around 18.7, 13.0 and 17.5 kN/m 3 , respectively. Specific gravity of those materials was 2.61, 2.82, 2.70 and 3.27, respectively (the relative high value of tailings D was attributed to the presence of iron). These corresponded to dry density of 1.55, 1.39, 0.46 and 0.98, respectively (Das 2010), which indicated high initial void ratios of uranium and bauxite tailings relative to those of gold tailings. Friction angles were assigned as 37°, 36°, 30° and 39° for those tailings, respectively. All four tailings were cohesionless as mine tailings usually are (Blight 2010). The relations between void ratio and effective stress, hydraulic conductivity and void ratio, and Young's modulus and effective stress were derived from column test for tailings A, B and C. There were some differences in terms of dimensions of samples and load applied in the tests for tailings C: the dimension of the columns used was around 10 × 10 cm, and tailings were underwent incremental loads (from 1 to 8 kPa) while those for tests on tailings A and B might reach around 300 kPa (Bhuiyan et al. 2015;Lévesque 2019;Nguyen and Pabst 2020). Details on how to derive these relations from compression column tests can be seen in Essayad and Aubertin (2021) and Lévesque (2019). Those relations were obtained from conventional consolidation tests in conjunction with hydraulic conductivity measurements for tailings D (Somogyi 1976).
The same waste rocks were considered for WRI in all models. WRI were assigned a simple linear elastic constitutive model with a high Young's modulus to represent the very high stiffness of the material (i.e., E = 500 MPa) and a Poisson's ratio of 0.277 (Bolduc and Aubertin 2014). WRIs were assigned a hydraulic conductivity of 3 × 10 −6 m/s (around 20 times higher than that of tailings) to facilitate the horizontal drainage of water in tailings. In practice, hydraulic conductivity of WRI might be greater but a smaller value was chosen to reduce computational time. This is because the presence of the two materials with contrastive permeabilities in the model can significantly increase the computational time of the model (Itasca 2021).

Constitutive models and update of tailings properties
Tailings properties vary with time as their void ratio decreases and the assumptions of constant properties might not be representative of actual and non-linear behaviours of slurry materials (Schiffman 1982;Townsend and McVay 1990;Morris 2002). Several studies have shown discrepancies between tailings consolidation simulated with assumption of constant properties and field observations. For instance, simulated consolidation rates were around 3 times faster than observation in the case of gold tailings disposed in pit (McDonald and Lane 2010). Various mathematical functions representing the relations effective stress-void ratio and void ratio-hydraulic conductivity for slurry tailings have been proposed. These include power function (Somogyi 1980), extended power function (Liu and Znidarčić 1991), logarithmic function (Bartholomeeusen et al. 2002), Weibull function and modified Kozeny-Carman equation (Mbonimpa et al. 2002;Chapuis and Aubertin 2003). Among these, the power function is often used for its simplicity and good representability (Priestley 2011;Agapito and Bareither 2018), and was chosen in the present study to represent variation of tailings hydraulic conductivity with void ratio (Table 1). From hydraulic conductivity relations (Table 1), tailings A and B exhibited a reduction of around 2.5 times after the first 3 loading before continuing to decrease slightly for the next loading steps (Lévesque 2019; Nguyen and Pabst 2020), while those for tailings C and D decreased by around 3 and 4 times (Somogyi 1976;Bhuiyan et al. 2015). Young's modulus of those tailings remained relatively small at low level of loads and non-linearly increased with the increased of applied loads, and values of tailings A and B were usually somewhat double that of tailings D (Table 1).
Linear elastic-perfectly plastic Mohr-Coulomb model was used to simulate tailings. The non-linear relations between e-σ', k-e and Young's modulus and effective stress of the four types of tailings followed an approach applied in Flac3D to modify these properties and were based on relations from Table 1. During the analysis, the value of void ratio was updated at every iteration based on the value of effective stress, and this new value of void ratio was then used to conduct an automatic update of hydraulic conductivity and stiffness in FLAC3D. The validity and precision of continuous update of tailings properties in Flac was presented by Zhou et al. (2019). This approach was compared with simulations using a constant saturated hydraulic conductivity.

Effect of WRI on the dissipation of excess PWP
Simulation results indicated that WRI significantly contributed to the horizontal dissipation of PWP and therefore to the acceleration of consolidation. For example, the PWP iso-contours immediately after the addition of the 7 th layer for tailings A (i.e., t = 7 years) were curved close the WRI, indicating a higher dissipation rate of PWP (Fig. 3). In other words, the dissipation of PWP close to the WRI could occur both vertically (to the surface of the TSF) and horizontally (to the WRI). Also, PWP always reached hydrostatic equilibrium at the end of each consolidation period (even for the cases without inclusions), therefore suggesting that no excess PWP remained before adding new tailings layers. Similar trends were observed for tailings B and C. The presence of WRI also contributed to accelerate tailings consolidations of tailings D, but excess PWP did not completely dissipate during the considered period and their effect was more limited compared to those for other tailings. For example, maximum degree of consolidation, U, after the addition of the 4 th layer was around 93% at 2 m from the WRI, while that reached only around 80% at 6 m from the WRI and excess PWP thus cannot fully dissipate. Such limited effect can be attributed to the relatively low hydraulic conductivity of these tailings (Newson et al. 2006;Matt 2015;Zhang et al. 2021). Results, therefore, indicated that WRI effectiveness is strongly affected by the tailings' hydraulic conductivity and that WRI for clay-like tailings might be significant in terms of improvement of the overall strength of the impoundment (Ferdosi et al. 2015) rather than improvement of consolidation rate of tailings. Simulation results showed that the effect of WRI on consolidation decreased rapidly with the distance from WRI (Fig. 4). For example, after the deposition of the 5 th layer, the degree of consolidation (at the base of the TSF, i.e., z = 0 m) for tailings A reached 95% in 13 days at 2 m from the inclusion, 34 days 14 m from the inclusion and 43 days 29 m from the inclusion (value for the case without inclusion was 45 days) (Fig. 4a). Consequently, the rate of dissipation of excess PWP was increased by over three times, but only in the direct vicinity of the inclusion. The decrease in the effect of WRI on the consolidation rate with the increase of the distance for tailings B, C, and D were also demonstrated for points at 6 m and 29 m after the placement of the 5 th layer (Fig. 4b, c, d). Accordingly, the radius of influence of the WRI was estimated as the location where time t 90 for the cases with and without WRIs were less than 5% different. For example, the radius of influence for tailings A was around 54 m when the tailings thickness was 30 m, 44 m for the thickness of 24 m and 34 m for the thickness of 18 m. The ratio between WRI's radius of influence and the thickness of tailings could be then estimated for each type of tailings (e.g., equal around 1,8 for the tailings A) (Fig. 5). This ratio for tailings B and C were around 1.8 and 2.6 at the end of the filling process (i.e., the 10th layer), respectively (Fig. 5). The WRI's zone of influence for tailings D could not be estimated from t 90 values as excess PWP did not fully dissipate, and this value was thus estimated as the distance where the differences of U between models with and without WRI were smaller than 5%. The radius of influence in that case was around 3.4 times the thickness of tailings (Fig. 5).
The radius of influence was different for each type of tailings, but they all exhibited a similar trend, that is a rapid and strong decrease after the filling of the first four layers before becoming relatively stable (Fig. 5). The ratio between WRI's radius of influence and the thickness of tailings after the placement of the 2 nd layer was around 4 for tailings A and B, 8 and 10 for tailings C and D, respectively, before decreasing to somewhat stable values when the tailings thickness was higher than 12 m (i.e., after the placement of the 4 th layer) (Fig. 5). The decrease of the radius of influence during the deposition of the first few layers can be attributed to the decrease of the hydraulic conductivity of the tailings under the applied loading and the decrease of void ratio under tailings settlement. These decreases in hydraulic conductivity of tailings can lead to the change in the hydraulic gradient and therefore possibly in the radius of influence (Hansbo et al. 1981;Indraratna et al. 2012;Han 2015;Saleh-Mbemba and Aubertin 2021b). Several endeavors with models of different Simulations also indicated using updated hydraulic conductivity values (changing with void ratio during consolidation) or constant values had a strong influence on the results. For example, the time required to reach 90% consolidation (or t 90 ) for tailings A 65 m from the inclusion Ratio of radius of influence of WRI over the tailings thickness for various types of tailings. The influence radius was different for each type of tailings and varied during the filling process: it exhibited a significant decrease after the placement of the first few layers before tending to a relatively constant value after the thickness reached around 12 m Fig. 6 Comparison of t 90 at the base of the TSF (tailings A) after the placement of the 10th layer as a function of the horizontal distance from the WRI and for various estimations of the hydraulic conductivity. A value of 1.3 × 10 −7 m/s generated somewhat similar results to those obtained with the model using updated hydraulic conductivity, while a value of 3.5 × 10 −7 m/s could yield results 3 times lower than those for model with updated hydraulic conductivity after the placement of the 10th layer was 130 days for model with updated properties, 41 days with k sat = 3.5 × 10 −7 m/s, and 162 days with k sat = 1.0 × 10 −7 m/s (Fig. 6). Using updated properties, however, requires to measure hydraulic conductivity in the laboratory for various void ratios to fit a descriptive models (Bhuiyan et al. 2015), and to input the hydraulic conductivity functions in the code, which is not always convenient. In practice, using a constant value k sat = 1.3 × 10 −7 m/s yielded somewhat similar results in terms of consolidation rate which corresponded to a void ratio of 0.62 (Fig. 6). The equivalent hydraulic conductivity for tailings B and C was 1.2 × 10 −7 and 1.1 × 10 −7 m/s, corresponding to the void ratio of 0.87-3.09, respectively.

Tailings settlement
Settlements occurred as pore pressures dissipated and varied depending on tailings properties, time and distance to the WRI. For example, the total settlement of the tailings body at the end of the placement process (t = 10 years) for tailings A was around 1.67 m, around 1.65, 0.5 and 1.94 m for tailings B, C and D, respectively. The increase in displacements of layer 1 reduced with the increase of tailings thickness, which indicated that the stiffness of the materials increaed under staged filling. For example, displacements in layer 1 for tailings A one year after the filling of the 3rd layer increased by around 47% but only 5% after the placement of the 10th layer (Fig. 7). Similar trends were observed for the other tailings with the increase of displacement after the filling of the last layer was generally less than 5%. Models thus represented very well the change in the stiffness of the materials. Settlement of tailings tended to be smaller closer to the WRI, which was attributed to the different in the strength and stiffness of WRI and the tailings as loads from tailings were partially transformed to the WRI. This depends on the modulus ratio of the WRI and surrounding tailings (van Eekelen et al. 2013;Han 2015) and the interface elements used in the models (Li and Aubertin 2009).
The volume change of tailings due to the presence of WRI can then be estimated based on the settlement evolution of tailings.

Volume change of tailings with the presence of WRI
The amount of volumetric compression of the tailings body resulting from settlement of tailings (i.e., equal to the average settlement of the tailings body times the area of the TSF, which was automatically calculated using a Fish function). Simulated results indicated that volume change always increased with time as more tailings were disposed (Fig. 8a). For example, volume change after the placement of the 2nd layer for tailings A without WRI was around 62 m 3 and increased to around 132 m 3 after the placement of the 10th layer. The changes in the slope of the curves also indicated a slower rate of volume change as tailings height increased (Fig. 8a). Tailings C had the smallest volume change (around only 40 m 3 ) among simulated tailings, which Fig. 7 Displacements of the first layer 1 year after the addition of the 2nd, 3rd, 9th and 10th layer (noted as L2, L3, L9 and L10) for 4 types of tailings Fig. 8 a Evolutions of average volume changes as a function of time resulting from tailings settlement for models with and without WRI for tailings A and b differences of volume change (expressed in percentage) between models with and without WRI for various types of tailings was explained by the smallest settlement of this type of tailings (session 3.2).
The presence of WRI can increase the consolidation rate of tailings, but they also occupy some space in the TSF. In practice, results have shown that most of the simulated tailings were fully consolidated after around 6 months. Using WRI therefore did not have a significant influence on the volume gain, but rather contributed to decrease the space available for tailings deposition. In other words, the volume of tailings that could be deposited in the TSF with the presence of WRI was always smaller than the case without WRI for all types of tailings (Fig. 8a). For example, the volume change of tailings A without WRI after the placement of the 10th layer was around 132 m 3 /per 5 linear meter, i.e., around 14% more than that of model with WRI (113 m 3 / per 5 linear meter) (Fig. 8a). The differences of volume change (DV) between models with and without WRI was estimated as where V t the volume change of model with tailings only (m 3 ); V w the volume change of model with WRI (m 3 ).
DV tended to become relatively constant after the deposition of the first fifth layers (Fig. 8b). For instance, DV of tailings A decreased from 25% after the placement of the 2nd layer to 16% after the filling of the 5th layer, while that of tailings C exhibited a reduction of nearly 17% (Fig. 8b). DV for tailings C and D tended to be greater than those of tailings A and B. For example, DV for gold tailings after the filling of the last layers were around 15%, while that for tailings C were somewhat higher, at around 20% (Fig. 8b). These differences could be explained by the difference in the stiffness evolution during compression of other types of tailings.
The additional height of the TSF required to dispose the same amount of tailings with WRI presence as without WRI can be accordingly estimated. This height can be estimated as the sum of volume occupied by WRI and total differences of volume change of tailings body between models with and without WRI divided by the area of the TSF (Table 2). These heights were around 2.22, 2.17, 2.06 and 2.45 m for tailings A, B, C and D, respectively. In other words, the increase of thickness of the TSF was less than 8% total thickness of tailings (Table 2). This was very meaningful for practitioners to briefly calculate the additional height of TSFs required when applying co-disposal technique (at least at the primary design stage). Ratio between volume of tailings and WRI disposed in this study was around 15.7. In practice, this ratio could widely change depending on the practical considerations of each mine sites. For example, this could depend on the primary goal of mine waste management of each mine site. Some might prioritize the volume of WRI that can be disposed to eliminate the need for the construction of waste rock piles. Coarse waste rocks might, indeed, contain sulfide minerals that can oxidize and generate acid mine drainage, and the placement of waste rock as inclusions inside the TSF could therefore efficiently contribute to prevent potential AMD generation (Jahanbakhshzadeh et al. 2019; Bussière and Guittonny 2020). The volumes of tailings and waste rock produced by mining activities can also vary widely in practice depending on the ore body characteristic and production technique (Blight 2010). The ratio of tailings and waste rocks volume was thus changed by changing the width of WRI in the model (i.e., from 6 m to 10, 12, 18 and 24 m corresponding to the volume ratio of 15.7, 9.0, 7.3, 4.6 and 3.2, respectively) to estimate the volume loss as a function of the quantity of waste rocks used in the TSF as inclusions. In general, the more volume of WRI placed in the pit, the less space dedicated for tailings that can be stored in the TSF (Fig. 9). For example, the DV values of tailings A for model with WRI width of 24 m was around 30% which was nearly 2 times higher than that of model with WRI width of 6 m. The same trend was recorded for tailings B and C (Fig. 9).  . 9 Evolution of DV with the ratio of volume of tailings and WRI in the TSF. Calculations assumed that the displacement at the boundary interface is zero and the width of the tailings part was large enough to reduce the effect of TSF sides The trend of tailings A and B was somewhat similar to each other, while DV for tailings C seemed to be higher than those of tailings A and B. Accordingly, equations presenting the evolution of DV with the changes in the volume ratio of tailings and WRI were derived as y = 49.6x −0.48 (R 2 = 0.97) for tailings A, y = 55.6x −0.58 for tailings B (R 2 = 0.99) and y = 56.8x −0.39 (R 2 = 0.980) for uranium tailings. These results can, thus, provide practitioners with a brief estimation (at least for the primary design stage) the potential volume of tailings that might be loss due to the presence of WRI for various volume ratios.
Effect of the filling rate TSF depth and deposition rates vary with mine size, production rate and type of operations (surface or underground) (Blight 2010;MEND 2015). Filling rate for TSF at Malartic mine is, for example, around 3 m per year (Boudrias 2018), and between 2 and 14 m per year at Rabbit Lake mine (MEND 2015). Additional simulations corresponding to 6 and 9 m/year filling rate were, therefore, carried out for tailings A to evaluate effects of tailings thickness on the evolution of tailings consolidation with the presence of WRI. Filling rate seemed to have a negligible effect on the volume change between the cases with and without WRI, which remained around 14% at the end of the filling process for all models regardless of layer thickness (Fig. 10a). A ratio R t90, i.e., the ratio of t 90 between models with and without WRI, was introduced to represent the effect of filling rate on the consolidation rate of tailings (Fig. 10b). The effect of WRI on the consolidation rate was increased for greater filling rates, which indicated that the use of WRI would be more beneficial with higher filling rate in practice. R t90 at 10 m from the WRI after the deposition of the 6th layer was 1.45, 1.70 and 1.83 for the 3-m, 6-m and 9-m thick models, respectively (Fig. 10b). The ratio also increased with the number of tailings layers. This ratio was, for example, around 1.46 and 1.54 for the 6-m thick model after the placement of the 2nd and 4th layer (Fig. 10b).
Models with 3 m/year but with different filling scheme instead of an instantaneously adding (i.e., layer of 1.5 m per 6 months and layer of 1 m per 4 months, respectively) were performed for tailings A to evaluate the potential effect of the assumption of instantaneous filling scheme on the rate of consolidation of tailings. Simulated results indicated that the effect of instantaneous filling on the rate of consolidation of tailings was insignificant. It took essentially the same amount of time with such filling intervals to dissipate 90% of excess PWP. The same trend was also observed by Boudrias (2018).

Discussion
WRI can contribute to accelerate the rate of consolidation of tailings but also reduces the volume available for tailings storage at the same time. The geometry of the TSF in this study and the configuration of the WRIs were simplified, and the results can be applied for the relatively simple geometry TSFs. This paper only dealt with one row of inclusion, while there might be several orthogonal inclusions (not only parallel inclusions) installed for the typical design of WRIs in TSFs (Aubertin et al. 2016). Other simulations on the more complex geometries of the TSFs and various permeable conditions (i.e., only central WRI was considered in this study) can be further performed to evaluate their influences on the evolution of tailings consolidation.
The models carried out with various type of tailings highlighted the importance of the constitutive relations describing the changes in the consolidation properties of tailings, which are unique for each type of tailings. For example, the relations for the tailings C were derived from large strain consolidation tests with pressure up to 8 kPa (Bhuiyan et al. 2015). In practice, the loads resulting from the filling of Fig. 10 a DV values for models with various filling rates and b R t90 values for the point located at 10 m from WRI at the TSF base after the placement of the 2nd, 4th and 6th layer for models of 3 m, 6 m and 9 m-thick layers for tailings A tailings might be significantly higher, and the change in the hydraulic conductivity and stiffness of tailings might become less non-linear when the applied loadings increase. Tests covering the actual range of pressures that are expected in the field are therefore recommended to improve the accuracy of constitutive relations.
Models with several values of constant Young's modulus of tailings were performed for tailings A to test the effect of stiffness differences on the influence of WRI on the settlement of tailings materials. In general, the higher the stiffness of tailings, the smaller zone of influence of WRI on tailings settlement. Also, tailings settlement simulated close to the WRI may be affected by the interface element used in this numerical study. The interface elements introduced in the simulations might not necessarily represent the real contact between WRI and tailings in practice. It is recommended that stiffness and frictions angles of the contact interface between these two materials should be estimated to improve the precision of the calculation. Despite these boundary effects, general behaviour are expected to remain similar (Li and Aubertin 2009;Yang 2016).
Finally, field-scale investigations would be required to validate or calibrate the numerical results and the upscaling of the approach in this study.

Conclusion
An investigation on the effect of a drainage system composed of WRI on the consolidation behaviors of several types of tailings in a simple TSF were numerically performed using Flac3D. Simulations were conducted considering continuously updating properties of tailings under sequential filling. Several practical considerations related to the operational aspects that could have effect on the settlement of tailings were also evaluated. Results including excess PWP dissipation rate and magnitude of consolidation were compared, bringing the following conclusions: WRI had a significant effect on the acceleration of PWP dissipation and thus on tailings consolidation rate. The effect of WRI depended on the distance and the thickness of the tailings. PWP always reached the hydrostatic value at the end of each filling stage for tailings A, B and C, while excess PWP in tailings D were not fully dissipated after one year because of their relatively low hydraulic conductivity.
Significant differences (up to 3 times) were observed when using a model with updated properties or with constant properties. Equivalent constant values of hydraulic conductivity were derived for these types of tailings, which can bring about an acceptable tolerance of t90 compared to the model with updated properties, reducing the computational efforts and the necessity of laboratory tests on non-linear evolution of tailings' properties.
The zone of influence of WRI was different for each type of tailings, decreasing significantly during the filling of the first few layers before remaining relatively constant after the thickness of tailings reached 12 m. The zone of influence was usually a function of the thickness of the tailings, around 1.8 times tailings thickness for tailings A and B, 2.6 for tailings C and 3.4 for tailings D, respectively. Settlement of tailings C was smallest, while largest value was observed for the tailings D. Changes in the tailings settlement was less pronounced after the placement of few last layers which demonstrated very well the increases in the stiffness of the materials during consolidation process.
The volume change due to settlement of tailings of model with only tailings was always greater than that of model with WRI for all types of tailing. In other words, the use of WRI decreased the volume available for the storage of tailings in the TSF, despite the acceleration of their consolidation. Additional height required when applying WRI was less than 8% of total thickness of tailings. DV for tailings A, B and D after the filling of the last layers were around 15%, while that for tailings C was around 20%.
Space available for tailings to be stored in the TSF was lower with the increase of volume of WRI placed in the TSF, and formulations estimating evolution of DV with the changes in the volume ratio of tailings and WRI were proposed.
Increase of the filling rate led to a more pronounced influence of WRI on the consolidation rate of tailings but had a negligible effect on the DV values.
3D effects because of complex drainage paths and permeable conditions are being investigated in an on-going project in the group. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.