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}{\mathit{dx}}}\left(k_{s,\mathrm{h}}(x)\mathrm{h}{\displaystyle \frac{d\mathrm{h}}{\mathit{dx}}}\right)=-R,$$

where $x$ is the distance to the edge of the tailings pond, $k_{s,\mathrm{h}}\left(x\right)$ 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}\left(0\right)=100$ [m] and $\mathrm{h}\left(300\right)=0$ [m] are shown in Figure 1 for constant and exponentially increasing hydraulic conductivities of ${10}^{-6}$ [m/s] and ${10}^{-8}·{\mathrm{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}}\left(0\right)={10}^{-8} [\mathrm{m}/s]$ and $k_{s,\mathrm{h}}\left(300\right)={10}^{-6} [\mathrm{m}/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 above the level of the tailings pond at 100m as the distance to the tailings pond increases from 0m. Therefore, more

accurate description of seepage through saturated tailings beach material deposited near the pond edge is required to obtain physically meaningful results. 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 benchmarked FEFLOW finite element numerical modeling to analyze design-intended and as-built seepage through a hydraulically segregated TSF in 2 spatial dimensions. The hydraulically segregated TSF design specifications are selected to be representative of a currently operational copper TSF, 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 steady state phreatic surface level with the transient phreatic surface level obtained by solving the 2D RTF equation throughout the entire staged construction process. Chapter 4 provides results of simulating staged construction of a hydraulically segregrated TSF with two tailings dams in 2D with FEFLOW, including phreatic surface level, drain flows, and water balance. 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\right)+Q,$$

where $\mathrm{h}$ is pore water total hydraulic head, $k_s$ is the saturated material hydraulic conductivity tensor, ${\mathrm{S}}_s$ is material specific storage, and $\mathrm{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 u}{\partial t}}=\nabla ·\left(\mathit{k}\left(u\right)\nabla (u+z)\right)+Q,$$

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

$$\theta \left(h\right)={\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}\left(u\right)={\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 \left(h\right)-{\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}{∆z}}\left[{\mathrm{K}}_{\mathrm{j}+1/2}^{\mathrm{n}+1}\left({\displaystyle \frac{h_{j+1}^{n+1}-h_j^{n+1}}{∆z}}+1\right)-{\mathrm{K}}_{\mathrm{j}-1/2}^{\mathrm{n}+1}\left({\displaystyle \frac{h_j^{n+1}-h_{j-1}^{n+1}}{∆z}}+1\right)\right],$$

where $\mathrm{j}$ indexes space and $\mathrm{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], ${\mathrm{S}}_s=1\times {10}^{-4}$[1/m], $\alpha =0.5$, $\mathrm{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 }_{\mathrm{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 ${\mathrm{k}}_{\mathrm{s},\mathrm{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 }_{\mathrm{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 $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 }_{\mathrm{c}}\right)$, and approximately linear dependence of ${\theta }_{\mathrm{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 }_{\mathrm{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,\mathrm{h}}$ 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 STF steady state solution ${\mathrm{h}}_s(x,z)$ for pressure head can be computed precisely, and is roughly approximated by the analytic expression:

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

which solves the STF equation on the rectangular domain with a non-zero Dirichlet boundary condition along the underdrain. Semi-analytic solution for the RTF steady state pressure head is obtained for $\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 \left(k_{s,h}\right(1-2{\alpha }^{n-1}{(-h)}^{n-1}),k_{s,v}(1-{2\alpha }^{n-1}{(-h)}^{n-1}\left)\right), 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 z}}\left(k_{s,v}{\displaystyle \frac{\partial h_1}{\partial z}}\right)=2k_{s,x}{\displaystyle \frac{\partial {(-h_s)}^{n-1}}{\partial x}}{\displaystyle \frac{\partial h_s}{\partial x}}+2k_{s,z}{\displaystyle \frac{\partial {(-h_s)}^{n-1}}{\partial z}}\left({\displaystyle \frac{\partial h_s}{\partial z}}+1\right), h_s\le 0$$

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.

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. 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 quasi-steady state solution capillary fringe is equivalent to a porous medium equation description whose transient decay to steady state is power law rather than 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. The FEFLOW mesh used for solving the RTF equation in both cases has 9483 nodes and 18576 triangular elements with side lengths between 1.2m and 2.5m. 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 evident that the water content of unsaturated material within the tailings beach is much different between RTF transient and quasi-steady state solutions.

4. 2D TSF Staged Construction

FEFLOW Model Definition

Tailings Storage Facility Geometry

The 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 equilibrium groundwater table at the start of TSF construction, at elevation 0[m], and except where 2.5[m] thick underdrains have been constructed with excavation of foundation material and 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 assignment of Van Genuchten parameters to underdrain material. 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 this 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 are simulated using FEFLOW STF and RTF solvers 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. 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 model boundaries between beach and slimes. 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 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 and underdrain flows.

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 of 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 change in element porosity. This approximation accounts for consolidaton water production in TSF transient seepage analysis without directly coupling RTF and consolidation equation solutions.

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}}^{\mathrm{\text{'}}}}^{\mathrm{B}}+\mathrm{M}$$

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

are used with copper tailings slimes constants $(A,M,C,D)=(2.15,0.5,1.8\times {10}^{-8} [\mathrm{m}/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

Phreatic Surface Level

Figure 9 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 cyclone underflow to adjacent tailings beach is increased from 10 to 100 [13]. Setting III is imposed to simulate the reduction in void ratio of the 2011 downstream cycloned sand shell due to reactivated TSF construction [7]. Figure 9 shows the end of 2025 phreatic surface through the TSF for each of the 4 material parameter settings, demonstrating that condition III has a greater influence on dam crest phreatic surface level than condition II, as may be relevant to diagnosing changes in dam crest piezometer readings across the history of TSF staged construction.

Water Balance

Table 3 presents annually averaged rates $[\mathrm{mL}/\mathrm{m}·s$] of underdrain and foundation outflow,

pond inflow, and TSF total water content change. For the purpose of writing this table, 2D transient seepage analyses were conducted with consolidation water production turned off, under the assumption consolidation water outflow from the top and bottom of the tailings slimes impoundment is more accurately quantified by 1D finite difference simulation than 2D finite element simulation results. To explain this assumption, it is noted that for each year, due to limitations of the RTF solver accuracy (i.e. residual rate error), fluid volume imbalance is approximately a quarter of a percent of the total TSF fluid volume at 350$\left[{\mathrm{m}}^3/\mathrm{m}·\mathrm{year}\right]$, or $11.4\left[\mathrm{mL}/\mathrm{m}·s\right]$, while consolidation water seeping out of the tailings impoundment into the foundation is estimated to have a volume of at most $1000\left[\mathrm{m}·0.1\mathrm{m}/\mathrm{year}\right]$, or $3.2\left[\mathrm{mL}/\mathrm{m}·s\right]$. 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 cyclone underflow and spigotting water infiltrating the dam shells and tailings beaches rather than direct seepage from the tailings pond.

Without direct RTF solution of 3D TSF transient seepage, the 2D solution, 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, these bounds are (0.6[L/s],1.8[L/s])-West and (0.7[L/s],2.1[L/s])-East in accordance with Table 3, implying 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 as initial indicators of 3D drain clogging problems detected by piezometers in practice [14].

5. Conclusions

1D staged construction unsaturated seepage flows 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 requirements. In 2D, steady state solutions of the RTF and STF equations for a single stage of idealized tailings dam staged 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, 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 and effects of geometric irregularity on RTF solution numerical stability. It may also be of value to extend simulation results to model the TSF with a partially saturated foundation and compute the effect of this partial saturation on RTF predicted drain flows.