1. Introduction

The construction of tailings deposits gives rise to displacements and pore pressure in soil mass. These factors may, depending on their magnitude, cause significant damage to the deposit and buildings nearby. To help ensure the safety of such work, it is of fundamental importance that engineers be able to predict hydro-mechanical behavior. With this ability, they can gather support for important decisions related to the placement of drainage systems and even decisions related to defining construction rates.

Among the numerical methods available in literature for this purpose, the one researchers have most widely used is the finite element method (FEM). To obtain reliable numerical results, however, it is necessary to adopt a constitutive model capable of adequately representing the stress and strain relationships of the materials. It is also necessary to possess an appropriate computer program to simulate the construction process while also taking into account the coupling between flow and deformation.

This study makes use of the computational program ANLOG - Nonlinear Analysis of Geotechnical Works (^{Nogueira, 1998}) to analyze the hydro-mechanical behavior of a tailings deposit from iron mining. The study analyzes this deposit in the plane strain state, in the field of small displacements and takes no account of the presence of the water level throughout the construction process. The materials are considered saturated and anisotropy of permeability. The geotechnical characteristics of the materials and the deposit geometry have been obtained from ^{Rezende (2013)} as described by ^{Braga (2016)}.

2. Numerical modeling of landfill construction

The numerical modeling by FEM of the landfill construction process involves defining the equivalent nodal force, due to the action of each new layer on those already built, and defining the stiffness of the new set of elements.

This study has adopted the technique known as "turn on gravity." By this technique, the effect of a new layer is simulated by applying an equivalent nodal force due to body forces (**F**_{b}) on the elements that belong to this layer. This technique still adopts the procedure of dynamic mesh. In this case, a single finite element mesh-containing all layers to be built-has previously been defined. These elements have a "turn off" status (that is, the respective degrees of freedom are initially restricted). At each construction stage, the group of elements to be built has the status changed to "turn on" (resulting in the release of their respective degrees of freedom). Hence, the “turn on,” or active, elements are allowed to contribute to the overall stiffness of the system (^{Teixeira, 2012}).

Applying FEM to solve the mechanical problem of static equilibrium, related to the landfill construction process, coupling the flow and the deformation lead to an algebraic equation system in terms of displacements (**u**) and pore pressure (**p**). This equation should be solved over time for each stage of construction.

Thus, from an equilibrium configuration *n*, where the states of effective stress (**σ**’), pore pressure (**p**) and strain (**ε**) are known, an incremental solution, associated with a given time increment (Δt) is obtained (^{Nogueira et al., 2009}). In case of problems involving physical nonlinearity, the predicted solution should be corrected by successive iterations until reaching a new equilibrium configuration *n+1*. In this strategy, the solution of the problem is obtained, updating the global nodal vector of unknown variables (**Û**^{T} = [**û p̂**]) at every new balance configuration, as follows:

In which **J**_{t} is the Jacobian global matrix; and, Δ**R**_{ext} and Δ**R**_{int}
^{k} are, respectively, the increment of the global nodal equivalent external and internal force, given by:

where iter is the number of iterations needed to achieve convergence in the current step. The matrix **J**, called Jacobin global matrix, is composed of the global matrixes of stiffness (**K**_{t}), coupling (**C**) and flux (**H**) of the active elements in a given construction stage.

The coupling matrix of each active element does not depend on the material properties. It depends only on the finite element approximation. The flux matrix of each active element is obtained as a function of the permeability matrix (**k**), defined, in the context of the tow dimensional flow in saturated soil media, as follow:

where k_{1} is the major principal permeability; β is the angle formed between the x direction and major principal direction, and F_{k} is the anisotropy factor.

The stiffness matrix of each active element is obtained as a function of the tangent constitutive matrix (**D**_{t}) which depends on the constitutive model.

In the field of nonlinear elastic models, the tangent constitutive matrix is equal to the tangent elastic constitutive matrix (**D**_{e}), which depends on the tangent elastic deformability modulus (E_{t}). In this context, the hyperbolic constitutive model proposed by ^{Duncan (1980)} adopts the following tangent elastic deformability modulus and the Bulk modulus (B) written in terms of the effective minor principal stress level (σ’_{3}) and the deviatory stress level (σ_{1} - σ_{3}):

where K_{i}, n, c (cohesion), f (frictional angle), R_{f}, K_{b} and m are the model’s parameters. The atmospheric pressure (p_{a}) is used in order to normalize the stress level. Using this model version it is possible to obtain a null volumetric change post peak.

In the field of elastoplasticity model, the tangent constitutive matrix is given by:

where **b** and **a** are, respectively, the gradient of the plastic potential function (G) and the gradient of the yielding function (F); and, H is the hardening module.

The Lade-Kim model used in this work is written in terms of the first and third invariant of effective stress tensor (I_{1} and I_{3}) and in terms of the second invariant of deviator stress tensor (I_{2D}) by adopting a nonlinear elasticity and a non-associate plasticity with hardening-softening law based on the plastic work. The tangent elastic deformability modulus used on the elastic constitutive matrix (**D**_{e}) can be obtained, according to^{Lade and Kim (1995)}, by doing:

or, according the ^{Lade and Kim (1988a)}, by doing:

where *ν* is the Poisson coefficient, K_{ur} and n; and, M and *λ* are elastic parameters.

The yielding function is defined according to ^{Lade and Kim (1988b)} as follows:

In which f ’_{p} is a stress level function and f ’’_{p} is the hardening-softening plastic work function. The stress level function is defined as:

In which

where α and h are yield parameters; and η_{1} and m_{1} are failure parameters. The function q is null during the isotropic compression and the unity at failure.

The hardening-softening work plastic function is defined according to the stress level. During the hardening (before failure) the work plastic function is defined as:

where c and p are hardening parameters obtained during the isotropic compression. During the plastic flow the yield function increases isotopically until it reaches the failure surface (^{Lade and Kim, 1988b}). The relationship between the yield function and the plastic work is described by a monotonic growth function in which its gradient decreases while the plastic work increases. During the softening (after failure) the work plastic function is defined as:

where A and B are evaluated at the failure (or peak) stress state condition by doing:

b’ is the softening parameter, which is considered positive and constant. The softening parameter equal to zero means the model present prefect plasticity. The magnitude of softening increasing with the parameter b’ as presented by ^{Yang (2009)}.

The potential function can be obtained, according to ^{Lade and Kim (1988b)}, by doing:

where ψ_{2} and m are the plastic parameters.

The hardening modulus depends on the hardening/softening function and is given by:

The vector Δ**R**_{ext} (Equation 3) depends on the increment of the external force equivalent to the body force of the layer that is being built (Δ**F**_{b}) and the parcel Δt**Hp**_{n} that represents, for the case where there is no water level, the volume variation applied in the current step. The vector Δ**F**_{int} (Equation 4) represents the global arrangement of the increment of internal force vector on each element (Δ**F**^{e}
_{int}) evaluated in each iteration k and defined as:

In which

where Δ**σ**’ and Δ**ε**^{k} are, respectively, the increment of effective stress and strain vectors, evaluated in each iteration at the active element level. The matrix **B**_{u} is the cinematic matrix.

At the end of each iterative cycle, a convergence state of the solution is verified by using a criterion that relates the Euclidian norm of the unbalanced force vector (**ψ**^{k} - Equation 1) with the Euclidian norm of the external force vector. Thus, for a given tolerance and at each increment, the iterative scheme ensures the overall balance by satisfying the compatibility conditions, boundary conditions, and constitutive relationships.

3. Numerical simulation

Figure 1 presents a schematic illustration of a tailings deposit, which was built by the upstream method in the form of drained stacking.

An initial compacted fill, called a starter dike (156 m of base length, 16 m top length, 40 m high, upstream slope of 1V:1.5H and downstream slope of 1V: 2H) was installed in a valley where the geological profile indicates a saprolite soil foundation of 40 m thick. The drainage system of this dike consists of a draining mat on its base.

After reaching the final starter dike’s height, the deposit was increased by using a 14 dikes (5m height each), called a raised dike, following the same downstream slope inclination of the starter dike. The deposit drainage system consists of a 120 m-long draining carpet built at an elevation of 36 m.

Figure 1 also shows some control points in order to present some numerical results. The points EM and EJ, in the starter dike, are adopted to accompany the excess of pore pressure generated during the process. A vertical section S1 through the D1 point on the crest of the 5th raised dike is used to evaluate the horizontal displacements throughout the deposit.

The raised dike and its respective layers of tailings were considered as a single stage of the landfill construction. The construction of the deposit was performed in 20 steps, as illustrated in Figure 1. The first stage (ET1) took 365 days and consisted of building the starter dike. The stages ET2 to ET5 (31, 61, and 90 days long) consisted of building the tailing deposit base and installing the 183 m-long draining carpet. Stage 6 (ET6), 22 days long, consisted of the construction of a downstream rockfill at the starter dike. Stages ET7 to ET20 consisted of constructing the raised dikes and their tailings layers; these were built respectively in 100, 122, 121, 92, 92, 123, 120, 122, 123, 120, 92, 153, 212, and 65 days.

The following parameters were used by materials. Starter dike: *γ* =17.52kN/m^{3}, k =5.3x10^{-10}(m/s), F_{k} =0.2; Rockfill: *γ* =25.5kN/m^{3}, k =10^{-2}(m/s), F_{k} =1.0; Draining carpet: *γ* =20kN/m^{3}, k =10^{-3}(m/s), F_{k} =1.0; Raised dike: *γ* =20.29kN/m^{3}, k =2.5x10^{-6}(m/s), F_{k} =0.30; Tailings: *γ* =17.54kN/m^{3}, k =3.3x10^{-6}(m/s), F_{k} =0.3.

In the context of linear elasticity, the following parameters were used by materials. Starter dike: E =15MPa, n =0.33; Rockfill: E =50MPa, *ν* =0.33; Draining carpet: E =60MPa, *ν* =0.33; Raised dike: E =30MPa, n =0.33; Tailings: E =22MPa, *ν* =0.33.

A set of CID (consolidated isotropically drained) conventional triaxial compression (CTC) tests was held with undisturbed samples obtained from the raised dike and the tailings deposit. The tests were carried out for confining pressures equal to 75 kPa, 150 kPa, 300 kPa, and 550 kPa.

Based on the results of these tests, ^{Braga (2016)} obtained the following constitutive parameters for the hyperbolic model (^{Duncan, 1980}): Raised dike: K_{i} =158.23, n =1.11, c =36.9kPa, f = 39.5º, R_{f} =0.47, K_{b} =223.92, m =0.9; Tailings: K_{i} =98.22, n =0.6, c =23.9kPa, f = 31.4º, R_{f} =0.75, K_{b} =173.94, m =0.28.

For the Lade-Kim model, ^{Braga (2016)} obtained the following constitutive parameters: Raised dike: M =46.67, *λ* =0.56, *ν* =0.33, m_{1} =0.27, *η* =99.7, c =0.00358, p =0.43, y_{2}=-3.22, *ψ* = 2.39, h=0.18, a =0.39, b’ =0.5; Tailings: K_{ur} =240.20, n =0.67, n =0.30, m_{1}=0.23, h =39.15, c =0.107, p=2.79, y_{2} = -3.17, *ψ* = 2.47, h =0.96, *α* = 1.13.

The laboratory test results of the raised dike material had shown softening behavior clearly for high level of confinement stress (Figure 2). The laboratory test results of the tailings material had shown a slightly softening behavior for a high level of confinement stress and large strain. So, only for the raised dike material was there adopted a softening behavior whereby the b’ parameters adopted were the ones which provided the best laboratory curve fitting.

Figure 2 presents the comparison between the laboratory tests’ curves and the results of a numerical simulation of an axisymmetric and uncoupled triaxial test, adopting a displacement control strategy for the solution's algorithm and considering the material parameters listed above. Using this strategy, it is possible to numerically obtain the post peak behavior of the CTC test.

Regarding the raised dike material (Figure 2a), the results of numerical simulation with the hyperbolic model provide good agreement in terms of stress, at least until the axial strain level corresponding to the maximum point of the stress’s curves. In terms of volumetric strain, the results are poor since the hyperbolic model takes no account of either the dilatancy or the softening effect. The results obtained with the Lade-Kim model, however, provide good agreement in terms of stress and volumetric strain, at least until approximately 8% of axial strain. Figure 2 shows the softening and dilatancy effect captured by the Lade-Kim model.

With relation to the tailings material (Figure 2b), the results obtained with the hyperbolic model provide a better approximation, in terms of stress, of the laboratory results than the results obtained using the Lade-Kim model. In terms of volumetric strain, the numerical results were not well represented by any of the models. It is important to note that the laboratory results for the tailings material do not represent the standard behavior of compressive materials. This fact led to a great difficulty in getting the constitutive parameters and certainly influenced the low quality on numerical simulation for these triaxial tests. In the absence of a larger set of laboratory testing, these were adopted and are still used in obtaining the constitutive parameters.

4. Results and discussion

Several numerical simulations were conducted by ANLOG using the finite element mesh illustrated in Figure 3 (2945 nodal points, 914 quadrilateral elements Q8Q4, isoparametric displacements quadratic and superparametric pore water pressure linear). The study evaluated the influence of the starter dike’s drain, the anisotropy of permeability, the construction rate and the constitutive model in generating an excess of pore pressure and displacements.

Initially, a linear elastic analysis was adopted to investigate the influence of the starter dike’s drain but considering isotropy of permeability (k_{x} = k_{y}) for all involved materials. The presence of the drain was simulated by assuming a null excess of pore pressure in the nodal points shown in Figure 3.

Figure 4 illustrates the variation of excess pore pressure on points EM and EJ along the stages of construction depending on drain presence. According to Figure 4a, at point EM, there is no significant influence. The important fact at this point is the high level of averaged excess pore pressure around 75% of the vertical stress due to self-weight. At point EJ, however, the presence of the drain affects the magnitude of the excess pore pressure, both during the starter dike construction stage (ET1) as well as during the tailing deposit construction stages (ET2-ET20).

As the tailings deposit’s construction without the starter dike’s drain has provided larger values of excess pore pressure, the numerical simulation of the construction process without using this drain was adopted in order to study of the effect of the permeability anisotropy.

To evaluate the influence of the permeability anisotropy, a linear elastic analysis was also performed. Figure 5 presents the variation of pore pressure excess on points EM and EJ throughout the building process as a function of the permeability anisotropy. Anisotropy has increased the magnitude of the excess of pore pressure at both points analyzed.

To investigate the influence of the construction rate, a linear elastic analysis was run (without the dike’s starter drain and considering the permeability anisotropy) by using three different rates: v_{1}, v_{2}, and v_{3}. Rate v_{2} corresponds to the construction rate previously mentioned. The rates v_{1} and v_{3} are defined as 0.5v_{2} and 1.5v_{2}, respectively. Figure 6 shows the excess of pore pressure for these rates at points EM and EJ. As expected, the excess of pore pressure increases with the increasing of construction rates.

The excess of pore pressure at points EM and EJ (Figure 6), at the end of the starter dike construction, do not depend on the construction rates, since for the permeability order of magnitude of the dike material, an undrained response is expected (^{Nogueira et al., 2009}). At the end of the construction process, these excesses are even more evident and vary by approximately 100 to 120 kPa.

Figure 7 shows the iso-curves of excess pore pressure at the end of the construction process for the analyses carried out with the models-linear elastic, hyperbolic and Lade-Kim. In these analyses, the study adopted as the permeability anisotropy, the construction rate of v_{2}, and the starter dike without the drain. In all three analyses the excess of pore pressure inside the deposit never exceeded 60 kPa. This value is considered low compared to the excess of pore pressure accumulated inside the starter dyke. Due to the high saturated permeability of the tailings (0.284 m/day), it can be said that the deposit demonstrated draining behavior throughout all of the construction process.

Figure 8 shows the variation of the horizontal displacements throughout the vertical section S1 at the points on the surface and 30 m depth along the constructive process. The analyses with the linear elastic and Lade-Kim models presented a very similar behavior. The biggest difference between the two analyses is around 15cm at the end of process and at 30 m depth (Figure 8b). The results obtained by using the hyperbolic model provided the highest horizontal displacements both at the surface and at a depth of 30 m.

5. Conclusion

From the results presented in herein, we can state that the presence of the dike’s starter draining carpet is essential to reducing the level of excess pore pressure and that the permeability anisotropy leads to a larger amount of water in the body of the starter dike. Due to the order of magnitude of the permeability of the tailings and the starter dike material as compared with the construction rate analyzed, what is observed is practically a draining behavior for the deposit and no draining for the starter dike. The numerical results highlight two matters of importance-having a more detailed set of CTC tests in order to facilitate the obtaining of constitutive parameters and having in place an efficient drainage system.

Moreover, it is important to highlight the necessity of taking into account the water level variation in the hydro-mechanical behavior of numerical simulations of landfills that are hydraulically constructed.