Transient Seepage Analysis of a Tailings Storage Facility Constructed by On-Dam Cycloning Using FEFLOW

1. Introduction

On-dam cycloning is a common method of constructing tailings storage facilities (TSFs) whereby total tailings and/or cyclone overflow spigotted from the facility periphery settles out across one or more tailings beaches with decreasing grain size towards the tailings pond. The Kareerand TSF in Matlosana, South Africa is an example of a TSF constructed by this method [11]. To correctly assess slope stability of TSF dams, including those constructed by on-dam cycloning, it is of critical importance to analyze the seepage of water through TSFs [14]. For example, after breach of the Jagersfontein TSF in 2022, an investigation panel identified several factors contributing to failure including high phreatic surface level caused by the reprocessed tailings deposition method and absence of a decant facility [20]. This fact points to the importance of using seepage analysis during TSF design phase to conservatively predict the effect of phreatic surface level and pore water pressure on dam stability, and of using seepage analysis to explain deviations of as-built phreatic conditions from design phase predictions.

Important factors determining the location of phreatic surface level in a TSF dam are the length of the tailings beach between the dam crest and tailings pond, the gradient of hydraulic conductivity across the tailings beach, drainage, and the rate of TSF construction, which in the case of construction by on-dam cycloning controls the spigotting water recharge rate [21]. An appreciation for the qualitative influence of beach hydraulic conductivity gradient and spigotting recharge rate on dam phreatic surface level can be obtained using the Dupuit-Forchheimer approximation for seepage flow, in which seepage is assumed to be horizontal, and the phreatic surface height $\mathrm{h}(x)$ along the length of the dam satisfies the differential equation:

$${\displaystyle \frac{d}{dx}}\left(k_{s,\mathrm{h}}(x)\mathrm{h}{\displaystyle \frac{d\mathrm{h}}{dx}}\right)=-R,$$

where $x$ is the distance to the edge of the tailings pond, $k_{s,\mathrm{h}}(x)$ is the saturated horizontal hydraulic conductivity of dam material, and $R$ is the spigotting recharge rate [5]. Solutions to the Dupuit-Forchheimer equation subject to the boundary conditions $\mathrm{h}(0)=100$ [m] and $\mathrm{h}(300)=0$ [m] are shown in Figure 1 for constant and exponentially increasing hydraulic conductivities of ${10}^{-6}$ [m/s] and ${10}^{-8}·e^{0.01535x}$ [m/s], and recharge rates of 0 and 2.5 [m/year]. The exponentially increasing hydraulic conductivity is selected so that $k_{s,\mathrm{h}}(0)={10}^{-8} [\mathrm{m}/\mathrm{s}]$ and $k_{s,\mathrm{h}}(300)={10}^{-6} [\mathrm{m}/\mathrm{s}]$ [1]. Note that the phreatic surface level computed by Dupuit-Forchheimer approximation for the case of exponentially increasing hydraulic conductivity and non-zero recharge is not physically realistic, in that the phreatic surface rises

Figure 1 suggests an accurate approximation of how tailings beach recharge rate relates to beach hydraulic conductivity profile may have value during the TSF design process for accurately approximating phreatic surface level in factor of safety computations, particularly for dams constructed by on-dam cycloning with hydraulic segregation of tailings whose beach hydraulic conductivity increases from the tailings pond edge to the dam crest. Moreover, accurate simulation of seepage flow during dam staged construction may have value to assessing whether or not changes in as-built piezometer monitoring data for a tailings dam, such as rising phreatic surface level beneath the dam crest, are indicative of inadequate hydraulic segregation of beached tailings. Therefore, because dam staged construction occurs with partial saturation of beach material above the phreatic surface, the first purpose of this article is to solve Richards’ transient flow (RTF) equation in 1 and 2 spatial dimensions for idealized tailings beach staged construction processes to clarify what beach saturation profiles are anticipated based on numerical modeling, documenting any practical engineering advantage solution of the RTF equation offers over seepage analysis performed by solving the saturated transient flow (STF) equation [15]. The second purpose is to apply FEFLOW numerical solution of the RTF equation to analyze design-intended seepage flow through a TSF in 2 spatial dimensions for comparison with as-built seepage flow. The TSF design specifications are selected to be representative of a currently operational copper TSF constructed by on-dam cycloning, modified to allow for straight forward FEFLOW computation of phreatic surface, drain flows, and water balance.

The outline of the article is as follows. Chapter 2 provides MATLAB finite difference solutions of RTF equations for 1D tailings beach staged construction processes, and compares FEFLOW finite element simulation with MATLAB simulation results. Chapter 3 compares FEFLOW and solutions of the 2D RTF and STF equations for steady state seepage through a tailings beach at the final stage of a 2D staged construction process, and the final RTF steady state phreatic surface level with the transient staged construction phreatic surface level. Chapter 4 provides 2D FEFLOW simulation results for seepage flow through a TSF with two tailings dams constructed in stages by on-dam cycloning. These results include phreatic surface level, drain flows, water balance, and possible deviation of as-built seepage flow from design intended seepage flow. Chapter 5 concludes by summarizing the article results.

2. 1D Staged Construction

A staged construction process is considered in which tailings beach material is deposited annually in 5m lifts for 20 years on top of a drain at elevation 0m, seepage through a constructed column of material is analyzed both in 1D with a MATLAB finite difference solver and 2D with FEFLOW. Accuracy and stability of the MATLAB RTF finite difference solver has been verified against 1D analytic solutions for Gardner evaporative drying and Philip short-time infiltration with Van Genuchten unsaturated flow material parameters, and by comparing MATLAB and FEFLOW simulation results, a mesh size and time stepping procedure required for accuracy of the FEFLOW RTF equation solver applied to analyze 2D staged construction seepage flow is informed [6,12].

Transient Flow Equations

The STF equation is:

$$S_s{\displaystyle \frac{\partial h}{\partial t}}=\nabla ·\left({\mathit{k}}_s\nabla (h+z)\right)+Q,$$

where $h$ is pore water pressure head, ${\mathit{k}}_s$ is the saturated material hydraulic conductivity tensor, $S_s$ is material specific storage, and $Q$ is a fluid source/sink term. Typical parameter values assigned to different tailings dam construction materials are shown in Table 1 [5,17]. In this table, $k_{s,\mathrm{h}}$ and $k_{s,v}$ are the saturated horizontal and vertical material hydraulic conductivities.

The RTF equation describing unsaturated transient pore water flow is:

$$\left({\displaystyle \frac{d\theta }{dh}}+S_s\right){\displaystyle \frac{\partial h}{\partial t}}=\nabla ·\left(\mathit{k}(h)\nabla (h+z)\right)+Q,$$

where $u$ is pressure head, $z$ is elevation head, $\mathit{k}(h)$ is the unsaturated hydraulic conductivity tensor, $\theta (h)$ is the volumetric water content, and $Q$ is a fluid source/sink term. The Van Genuchten equation:

$$\theta (h)={\theta }_r+{\displaystyle \frac{{\theta }_s-{\theta }_r}{{\left(1+{\left(\alpha \left|h\right|\right)}^n\right)}^m}},$$

expresses the volumetric water content in terms of the residual volumetric water content ${\theta }_r$, saturated volumetric water content (i.e. unsaturated-flow porosity) ${\theta }_s$, air entry $\alpha$, and pore size distribution $n$, with $m=1-1/n$ taking a value between 0 and 1. The Mualem-Van Genuchten equation:

$$\mathit{k}(u)={\mathit{k}}_s{\mathrm{S}}_{\mathrm{e}}^{1/2}{\left(1-{\left(1-{\mathrm{S}}_{\mathrm{e}}^{1/\mathrm{m}}\right)}^m\right)}^2,$$

expresses the unsaturated hydraulic conductivity tensor in terms of the saturated hydraulic conductivity tensor and effective saturation:

$${\mathrm{S}}_{\mathrm{e}}={\displaystyle \frac{\theta (h)-{\theta }_r}{{\theta }_s-{\theta }_r}}={\displaystyle \frac{1}{{\left(1+{\left(\alpha \left|h\right|\right)}^n\right)}^m}}.$$

Typical Van Genuchten parameter values assigned to different tailings dam construction materials are shown in Table 2 [9].

Matlab Finite Difference Solution

The 1D RTF equation can be solved by finite differences using the implicit Euler method in MATLAB. For the purpose of simulating seepage flow during staged construction, the 1D RTF equation is solved in mixed head-saturation form:

$${\displaystyle \frac{{\theta }_j^{n+1}-{\theta }_j^n}{∆\mathrm{t}}}={\displaystyle \frac{1}{∆\mathrm{z}}}\left[{\mathrm{K}}_{\mathrm{j}+1/2}^{\mathrm{n}+1}\left({\displaystyle \frac{h_{j+1}^{n+1}-h_j^{n+1}}{∆\mathrm{z}}}+1\right)-{\mathrm{K}}_{\mathrm{j}-1/2}^{\mathrm{n}+1}\left({\displaystyle \frac{h_j^{n+1}-h_{j-1}^{n+1}}{∆\mathrm{z}}}+1\right)\right],$$

where $j$ indexes space and $n$ indexes time, so that fluid mass balance is preserved [4]. Between construction lifts, 0 hydraulic head and 0 flow boundary conditions are imposed at elevation 0m and the 1D column surface elevation. Total hydraulic head within newly added saturated material is initialized to be constant at the column’s new surface elevation except within a boundary layer of thickness $\mathsf{\delta }=50$[cm] between newly deposited saturated material and previously deposited unsaturated material at elevation $z_p$, where a pressure head gradient:

$$h=h_{dry}+\tanh \left({\displaystyle \frac{z-z_p}{\mathsf{\delta }}}\right)\left(h_{wet}-h_{dry}\right),$$

is defined to improve stability of the finite difference Picard nonlinear solver [8]. The finite difference grid spacing is $10$[cm], and the time step is adaptively increased from 0.1[s] immediately after new construction to a maximum value of 1.5 days by factors of 2 whenever the number of Picard iterations required for pressure head convergence within $1\times {10}^{-4}$[m] is less than 14.

Figure 2 shows pressure head and water content profiles for a 1D staged construction process after stages 14 and 20, computed by finite difference solution of the 1D RTF equation for a tailings beach material with $k_{s,v}=1\times {10}^{-7}$ [m/s], $S_s=1\times {10}^{-4}$ [1/m], $\alpha =0.5$, $n=2$, ${\theta }_s=0.5$ and ${\theta }_r=0.15$. This solution suggests that away from the top and bottom of the column of material, water content remains at a constant value ${\theta }_c\approx 0.42$ as construction proceeds beyond the first 10 stages of construction. Similar water content profiles are observed when saturated vertical hydraulic conductivity is varied between values $k_{s,v}=1\times {10}^{-8}$[m/s] and $k_{s,v}=1\times {10}^{-6}$[m/s], as indicated by Figure 3-left where ${\theta }_c$ at stage 20 is plotted against $k_{s,v}$. Figure 3-right shows a plot of annual average drain flow for stage 20 against $k_{s,v}$. For all conductivity values:

$$drain flow\approx v({\theta }_s-{\theta }_c),$$

where $\mathrm{v}$ is the rate of construction 5[m/year], and an upper bound for drain flow as $k_{s,v}$ is increased beyond $1\times {10}^{-6}$[m/s] is $v({\theta }_s-{\theta }_r)=5.55\times {10}^{-5}$[L/s$·{\mathrm{m}}^2$]. For the plotted range of saturated vertical hydraulic conductivities, the drain flow at $z=0$[m] is approximately constant throughout stage 20 at the value $k_{s,v}\left({\theta }_c\right)$, and approximately linear dependence of ${\theta }_c$ follows from solving the equation:

$$v({\theta }_s-{\theta }_c)=k_{s,v}\left({\theta }_c\right),$$

and using the fact that the Mualem-Van Genuchten hydraulic conductivity function satisfies:

$$\ln v=\ln \left({\displaystyle \frac{k_{s,v}\left(\theta \right)}{{\theta }_s-\theta }}\right)\approx \ln k_{s,v}+32.06\theta -13.06,$$

for $0.34<\theta <0.48$.

As a verification of MATLAB simulation results, FEFLOW simulation of seepage flow during 1D staged construction was conducted using a 2D geometry of final height 100[m] and width 50[m] with no lateral variation of material parameters and an equilateral triangular mesh of side length approximately 1.44[m]. FEFLOW solution of the 2D RTF equation showed stage 20 pressure head and water content profiles similar to MATLAB finite difference profiles, with a maximum difference in ${\theta }_c$ values of less than 0.01 for the range of saturated vertical hydraulic conductivities $1\times {10}^{-8}<k_{s,v}<1\times {10}^{-6}$ [m/s].

3. 2D Staged Construction

RTF/STF Steady State Comparison

FEFLOW computed RTF and STF steady state solutions for the final construction stage of an idealized tailings dam in 2D have been verified against semi-analytic and analytic solutions, and used to compare RTF and STF phreatic surface level and transient decay to steady state. Figure 4-top shows RTF and STF quasi-steady state phreatic surface levels for a 70m high by 400m wide tailings sand embankment, computed using FEFLOW transient solvers to obtain stable convergence of the solution to steady state [19]. The upstream boundary condition is a constant hydraulic head of 70m, and the downstream boundary condition is a constant hydraulic head of 0m along a 100m long drain, with all other boundary conditions on the rectangular domain being Neumann 0 flow. Saturated horizontal hydraulic conductivity $k_{s,h}(x)$ of tailings beach material increases exponentially from 1e-7[m/s] to 1e-5[m/s] as distance from the upstream boundary increases from 0[m] to 300[m], and the vertical to horizontal saturated hydraulic conductivity ratio is 0.1. Material specific storage is 1e-4[1/m], and for RTF solution, Van Genuchten parameters $\alpha =0.5$, $n=2$, ${\theta }_s=0.5$ and ${\theta }_r=0.15$ are selected. The FEFLOW mesh used for solving the RTF and STF equations has 9483 nodes and 18576

The computed STF steady state solution $h_s(x,z)$ is roughly approximated by the analytic expression:

$$h_s(x,z)\approx 70e^{-{\displaystyle \frac{x\ln 100}{300}}}-z,$$

which solves the STF steady state equation on the rectangular domain with a non-zero Dirichlet boundary condition along the underdrain. The STF steady state drain flow computed with FEFLOW is approximately $0.65$[${\mathrm{m}}^2/\mathrm{day}$]. Semi-analytic solution for the RTF steady state pressure head is obtained for $\mathsf{\alpha }\ll 1$ by noting that in this limit:

$$(k_x,k_z)=(k_{s,h},k_{s,v}), h>0$$

$$(k_x,k_z)\approx (k_{s,h}(1-2{\alpha }^{n-1}{(-h)}^{n-1}),k_{s,v}(1-{2\alpha }^{n-1}{(-h)}^{n-1})), h\le 0$$

whereby the expansion:

$$h_r(x,z)\approx h_s(x,z)+{\alpha }^{n-1}{h}_1(x,z),$$

for the RTF steady state pressure head implies:

$${\displaystyle \frac{\partial }{\partial x}}\left(k_{s,h}{\displaystyle \frac{\partial h_1}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial \mathrm{z}}}\left(k_{s,v}{\displaystyle \frac{\partial h_1}{\partial \mathrm{z}}}\right)=2k_{s,h}{\displaystyle \frac{\partial {(-h_s)}^{n-1}}{\partial x}}{\displaystyle \frac{\partial h_s}{\partial x}}+2k_{s,v}{\displaystyle \frac{\partial {(-h_s)}^{n-1}}{\partial z}}\left({\displaystyle \frac{\partial h_s}{\partial z}}+1\right), h_s\le 0$$

Figure 4. Top: RTF and STF quasi-steady state phreatic surfaces computed using FEFLOW. Bottom: RTF and STF solution water content and average pressure head vs time triangular elements with side lengths between 1.2[m] and 2.5[m] and a typical side length of 1.75[m]. Observing solution phreatic surface with $\alpha$ values 1, 2, and 4 suggests it is necessary that the mesh size be at most 3.5$/\alpha$ to avoid visible spurious oscillations of the phreatic surface, although a mesh size of at most 1$/(5\alpha )$ is typically recommended for numerical accuracy.

Figure 4. Top: RTF and STF quasi-steady state phreatic surfaces computed using FEFLOW. Bottom: RTF and STF solution water content and average pressure head vs time triangular elements with side lengths between 1.2[m] and 2.5[m] and a typical side length of 1.75[m]. Observing solution phreatic surface with $\alpha$ values 1, 2, and 4 suggests it is necessary that the mesh size be at most 3.5$/\alpha$ to avoid visible spurious oscillations of the phreatic surface, although a mesh size of at most 1$/(5\alpha )$ is typically recommended for numerical accuracy.

Figure 5. 2D plot of the pressure head perturbation ${\mathrm{h}}_1(x,z)$ upon substitution in the RTF equation. Solving this differential equation with a MATLAB 2D finite difference solver gives the plot of ${\mathrm{h}}_1(x,z)$ shown in Figure 5. The RTF quasi-steady state drain flow computed with FEFLOW is approximately $0.32$[${\mathrm{m}}^2/$day], in close agreement with the prediction of Dupuit-Forchheimer theory.

Figure 5. 2D plot of the pressure head perturbation ${\mathrm{h}}_1(x,z)$ upon substitution in the RTF equation. Solving this differential equation with a MATLAB 2D finite difference solver gives the plot of ${\mathrm{h}}_1(x,z)$ shown in Figure 5. The RTF quasi-steady state drain flow computed with FEFLOW is approximately $0.32$[${\mathrm{m}}^2/$day], in close agreement with the prediction of Dupuit-Forchheimer theory.

Figure 4-bottom shows plots of total water content vs time and average pressure head vs time for FEFLOW quasi-steady and steady state solutions of the RTF and STF equations initiated from STF steady state and RTF quasi-steady state solutions. The STF solution fit shows exponential decay of average pressure head to a constant value with time constant 14.28[days], in approximate agreement with the longest decay time constant of 13.70[days] computed in MATLAB by solving a Helmholtz linear eigenvalue problem with mixed homogeneous Dirichlet and Neumann boundary conditions on the solution domain. The RTF solution fit shows approximately exponential decay of total water content to a constant value over the first 20 years of simulation, with fitting time constant 8.69[years], before continuing to increase at a rate in disagreement with exponential decay. When a different initial condition with a phreatic surface above the RTF steady state phreatic surface is selected, exponential decay of total water content with a fitting time constant 2.19[years] is observed over the first 5 years of simulation. Deviation of RTF transient decay time constants from STF transient decay time constants is not explained by analyzing transient corrections to the STF steady state solution of order ${\alpha }^{n-1}$, since these transients decay with STF decay time constants. A possible explanation for disagreement of long term RTF transient decay with exponential decay is that for Van Genuchten/Mualem material parameters, the RTF long term description of the dry capillary fringe is approximately:

$$S_s{\displaystyle \frac{\partial h}{\partial t}}\propto \nabla ·\left({|h|}^{-(n-1)(0.5+2/m)}{\mathit{k}}_s\nabla h\right),$$

which has similarity solutions of the form:

$$h(x,z)=t^a\overline{h}(x{t}^b,z{t}^b),$$

$$a={\displaystyle \frac{1+2b}{1+(n-1)(0.5+2/m)}},$$

with $a>0$ and $b<0$ whose effective saturation ${\mathrm{S}}_{\mathrm{e}}$ decreases to 0 in time by a power law in $t$ rather than an exponential.

RTF Staged Construction Solution

Figure 6-top shows the end of construction RTF transient phreatic surface level and saturation computed throughout 14 stages of 5[m] height increase each year. Figure 6-bottom shows the final stage quasi-steady state phreatic surface level and saturation. From Figure 6 it is evident that the RTF transient phreatic surface is slightly higher than the RTF steady state phreatic surface, unlike the STF transient phreatic surface which quickly decays to steady state within 365 days of new construction. It is also clear from the figure that the water content of unsaturated material within the tailings beach is different between RTF staged construction and quasi-steady state solutions.

Based on 1D MATLAB finite difference simulations the initial FEFLOW simulation time step is selected as 0.1[s] with maximum adaptive time step increase of $\times 2$ up to a maximum time step of 1.5[days]. With these settings, and $\alpha =0.5$, the FEFLOW solver shows no indication of Picard iteration instability, meaning the residual rate balance:

$$(\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w})-(\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w})-(\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e})$$

remains less than the FEFLOW default specified error tolerance of 0.0001[${\mathrm{m}}^2/$day]. However, when the Van Genuchten parameter $n$ is decreased from 2 to 1.5, FEFLOW indicates Picard iteration instability during stage 11, due to increased water retention in unsaturated material increasing the residual rate balance. If the parameter $n$ is not too close to 1, this problem can be solved by increasing the FEFLOW default number of Picard iterations, as is the case for $n=1.5$ if the number of iterations is increased from 12 to 30.

4. 2D TSF Staged Construction

Feflow Model Definition

Tailings Storage Facility Geometry

TSF layout, 2D cross sectional geometry, and 3D geometry are shown in Figure 7. The TSF is a cross valley copper tailings impoundment with two tailings dams, West and East, providing containment. Each dam is constructed on top of a 30[m] thick sand and gravel foundation with clay lenses, and contains a foundation underdrain, starter dam, cycloned sand downstream shell, tailings beach, and tailings slimes. The top of the foundation is assumed to coincide with the local groundwater table at the start of TSF construction, at elevation 0[m], except where 2.5[m] thick underdrains have been constructed with excavation of foundation material so that the top of the underdrains is located at elevation 0[m]. Note that this geometry implies underdrain material is always saturated, thereby avoiding unphysical oscillatory switching of underdrain hydraulic conductivity during staged construction that may otherwise occur with simulation of unsaturated flow. At the end of 2025 the dam crest elevations were 183[m], and the elevation of the tailings pond was 173[m]. Sloping sidewalls of the valley create a longitudinal dam cross section that is 750[m] wide at elevation 173[m] and 250[m] wide at the bottom of the TSF foundation at elevation -30[m]. The effect of sidewall sloping on cross sectional seepage flow is not considered for the purpose of 2D TSF seepage analysis.

The TSF is reactivated, in that the TSF was inactive for 15 years 1997-2011 before being used for another 14 years 2012-2025 starting from a dam crest elevation of 113[m]. The original TSF, used for 25 years before 1997, was built up to a dam crest elevation of 50[m] by centreline construction, followed by upstream construction until 1997. The reactivated TSF was built up to elevation 183[m] by centreline construction. West and East Dam tailings beaches of the reactivated TSF have been maintained at widths of 525[m] and 560[m]. West and East Dam tailings beaches existing prior to 1997 are not shown in Figure 7 but are assumed to have been maintained at widths 525[m] and 560[m] throughout all stages of construction.

Initial Conditions

To specify phreatic conditions within the TSF existing at the time of reactivation, the effect of 15 years of TSF inactivity on phreatic conditions is simulated using the FEFLOW RTF solver starting from a condition of complete water saturation of the TSF for which the pressure head of all materials in the TSF above elevation 0[m] is initialized to 0[m] and the hydraulic head of the underdrains and foundation is initialized to 0[m]. Seepage collection ponds with surface elevation 0m are present downstream of both West and East Dams, whereby 0[m] hydraulic head boundary conditions are assigned to the left edges of the West Dam underdrain and foundation and right edges of the East Dam underdrain and foundation, and no flow boundary conditions are assigned elsewhere under the assumption no tailings pond is present during deactivation. Depth dependent void ratio and hydraulic conductivity of tailings slimes are assigned based on 1D self weight consolidation simulation, and depth dependent tailings beach hydraulic conductivities are assumed to decay exponentially with distance from the dam crest to tailings slimes values at the material boundaries. TSF material parameters used to simulate phreatic conditions within the TSF during 15 years of inactivity are shown in Table 1 and Table 2.

Staged Construction Boundary Conditions

After reactivation, the TSF impoundment is increased in height by 5[m] each year. This height increase is modeled in FEFLOW by activating a new 5[m] layer of the TSF geometry each year. Within each newly activated layer, the hydraulic head at each location at time of deposition is initialized to a constant equal to the height of TSF at that location, under the assumption consolidation of newly deposited material occurs immediately, except at the interface between newly and previously deposited material where pressure heads are matched with a gradient smoothing function. No flow boundary conditions are specified on all TSF boundaries except at the location of the tailings pond where a constant hydraulic head equal to the height of the TSF is imposed, at the left edges of the West Dam underdrain and foundation where 0m hydraulic head boundary conditions are imposed, and at the right edges of the East Dam underdrain and foundation where 0[m] hydraulic head boundary conditions are imposed. Using the FEFLOW Python Interface Manager (IFM), depth dependent void ratio and hydraulic conductivity of tailings slimes are updated with activation of each new construction stage based on 1D self weight consolidation simulation, and depth dependent tailings beach hydraulic conductivities are updated with each new construction stage to decay exponentially with dam crest distance from dam crest values to tailings slimes values. The transient effect of undrained construction loading on TSF pore water pressure is not considered in the analysis, whose primary objectives are to compute phreatic surface level, underdrain flows, and TSF water balance.

Consolidation Water

While undrained construction loading occurs suddenly without change in underlying material porosity, consolidation is the gradual process by which the porosity of material decreases in delayed drainage response to construction loading. Authors of the Fundão investigative report used 1D FSConsol simulation of tailings sands and slimes consolidation to update the porosity and consolidation water production of tailings dam structural elements at each time step of the FEFLOW STF solver with an external program tracking changes in effective stress within the dam. Here, to avoid numerical instability of the FEFLOW RTF solver, a different 1D approximation of tailings slimes consolidation is used, whereby 1D self weight consolidation of tailings slimes is simulated by solving the 1D finite strain consolidation equation by finite differences, and the resulting time dependent porosity profile is used to adjust the porosity of tailings slimes elements each time a new tailings lift occurs [16,18]. So that fluid mass balance is preserved, a consolidation water source is introduced for each tailings slimes element with a water production rate specified by the annual change in element porosity. This approximation accounts for consolidaton water production in TSF transient seepage analysis without directly coupling RTF and consolidation equation solutions. However, it should be noted that because the geometry of newly activated layers is fixed in FEFLOW, this approximation does not account for consolidation

1D consolidation equations expressing void ratio and saturated vertical hydraulic conductivity in terms of vertical effective stress:

$$\mathrm{e}=\mathrm{A}{{\mathsf{\sigma }}_{\mathrm{v}}^{\prime }}^{\mathrm{B}}+\mathrm{M}$$

$${\mathrm{k}}_{\mathrm{s},\mathrm{v}}=\mathrm{C}{\mathrm{e}}^{\mathrm{D}},$$

are used with copper tailings slimes constants $(\mathrm{A},\mathrm{M},\mathrm{C},\mathrm{D})=(2.15,0.5,1.8\times {10}^{-8} [\mathrm{m}/\mathrm{s}],3.8)$ to simulate self weight consolidation of tailings slimes [2]. A maximum void ratio of 1.46 is assigned to slimes at the surface of the TSF. Figure 8 shows 1D simulation results for tailings slimes impoundment height across the history of TSF staged construction, and end of 2025 vertical profiles for tailings slimes saturated vertical hydraulic conductivity, void ratio, and excess pore pressure. 1D consolidation water flux out of the top and bottom of the impoundment is computed to be 1.6[m/year] and less than 0.1[m/year] at the end of 2025.

Feflow Simulation Results

Quasi-Steady State Phreatic Surface Level and Water Balance

As a point of reference for staged construction simulation results, the end of 2025 TSF quasi-steady state phreatic surface level and drain flows have been computed using a 15 year approach to steady state. Figure 9 shows the computed phreatic surface level and material saturation, using a mesh with 26462 nodes and 52240 triangular elements.

The pond inflow rate is 0.67[${\mathrm{m}}^2/$day], and the total drain/foundation outflow rate is 1.90[${\mathrm{m}}^2/$day], indicating that after 15 years of approach to steady state and stabilization of the phreatic surface, water storage in unsaturated parts of the dam remains far from steady state.

Staged Construction Phreatic Surface Level

TSF staged construction simulations use a mesh with 8235 nodes and 16016 triangular elements staged. The simulation maximum time step is reset to 0.1[s] after each new lift with maximum adaptive time step increase of $\times 1.04$ up to a maximum time step of 1.5[days]. Figure 10 shows end of 2025 TSF phreatic surfaces for 4 different material parameter settings:

• I: Horizontal hydraulic conductivities of tailings beach at dam crest and entire uncompacted cycloned sand shell are ${10}^{-6}$[m/s] and ${10}^{-5}$[m/s]
• II: Horizontal hydraulic conductivity of tailings beach at dam crest decreased from ${10}^{-6}$[m/s] to ${10}^{-7}$[m/s].
• III: Horizontal hydraulic conductivity of uncompacted cycloned sand shell existing before TSF reactivation decreased from ${10}^{-5}$[m/s] to ${10}^{-6}$[m/s].
• IV: Both settings II and III.

Setting II is imposed to simulate the effect of constructing the tailings beach with some mixture of total tailings and cyclone overflow, rather than total tailings alone, whereby the ratio of saturated horizontal hydraulic conductivity of

Staged Construction Water Balance

Table 3 presents annually averaged rates [${\mathrm{m}}^2/$day] of underdrain and foundation outflow, pond inflow, pond outflow, consolidation water inflow, and TSF total water content change including change in specific storage. The table demonstrates that for working underdrains, the West+East underdrain outflow is more than an order of magnitude greater than West+East foundation outflow bypassing the drains, and that the majority of seepage outflow is attributable to change in water storage of the dam shells and tailings beaches and consolidation water inflow rather than direct seepage from the tailings pond.

The 2D seepage analysis results, extruded to widths 250m and 750m, can be used to set rough lower and upper bounds on design-intended 3D underdrain, foundation, and pond seepage flows. For 2025, the lower and upper bounds on total underdrain and foundation flow are (700[${\mathrm{m}}^3/$day], 2099[${\mathrm{m}}^3/$day])-West and (640[${\mathrm{m}}^3/$day], 1920[${\mathrm{m}}^3/$day])-East, or (8.1[L$/$s], 24.3[L$/$s])-West and (7.4[L$/$s], 22.2[L$/$s])-East, in accordance with Table 3. This result implies measured foundation flows an order of magnitude greater than these upper bounds are indicative of as-built TSF structural problems such as internal erosion. 2D simulations performed with reduced underdrain hydraulic conductivity can further be used to analyze drain clogging problems that may be detected by piezometers in practice [14].

5. Conclusion

1D staged construction seepage flow computations have been presented as approximations of how tailings beach hydraulic conductivity influences water content of unsaturated material above the phreatic surface, and as a demonstration of how finite difference solver stability requirements can be used to inform finite element simulation mesh size and time stepping. In 2D, steady state solutions of the RTF and STF equations for a single stage of tailings dam construction have been verified against semi-analytic and analytic solutions, demonstrating significant ~10[m] differences in RTF and STF phreatic surface level are possible, and two orders of magnitude difference in transient decay time to steady state are possible. FEFLOW simulation results demonstrate that for a TSF constructed by on-dam cycloning with long tailings beaches and a permeable foundation, 2D RTF phreatic surface, drain flows, and water balance can be computed with an accuracy useful for initial TSF dam design screening and diagnosing as-built structural problems, such as high dam crest piezometer readings and unintended foundation seepage. In particular, 2D simulation results demonstrate TSF dam crest piezometer readings are more sensitive to hydraulic conductivity of the cycloned sand downstream shell than the hydraulic conductivity gradient of the tailings beach, and place upper and lower bounds on design intended 3D foundation flows. While simulation results have been obtained for an idealized 2D TSF geometry, it may be of value to obtain 3D results for a real TSF topography to clarify computational requirements for accurate and stable solution of the RTF equation.