Rehabilitation of tailings dams for final closure is one of the biggest challenges facing engineers. Not only is there the potential for producing leachate, but stability and settlement must also be considered.
The water balance of a tailings dam is unique because it usually contains a pond that causes a phreatic level in the dam. The position of the phreatic surface defines the saturated and unsaturated zones in the impoundment, which of course varies spatially and temporally.
Predictive modelling is difficult, as numerical models capable of analysing the combined saturated/unsaturated zones are not adequately refined to accurately solve the flux boundary problem for infiltration at the surface of the tailings.
The use of a flux function to accurately describe and predict the surface fluxes through the top unsaturated tailings cross-section is a great benefit for the tailings engineer.
SoilVision Systems has developed 3D modelling software, SVFlux, to verify a net infiltration flux function to bridge the gap between accurate surface flux boundary calculations and regional flow modelling software, and to simplify the 3D model, without sacrificing accuracy.
The software was used to model a gold tailings dam in Queensland, Australia. The climate here is tropical, with annual rainfall of 702mm, falling mostly as high intensity showers between November and March. The potential evaporation of 1650mm a year means the annual site water balance is negative.
The tailings dam was built in a stream valley using hydraulically placed tailings behind a 5.8km long engineered embankment encircling 70% of the dam. The tailings area covers 310ha including a 100ha pool.
An additional 104ha catchment impacts on the dam because it is built against a hillside. The embankment height (and thus tailings depth) varies from 32m to less than 1m.
Distribution of surface fluxes on the tailings dam surface Both a saturated and an unsaturated zone exist in the tailings dam because of the pool. The shape of the phreatic line is governed by the tailings' properties and the exit is determined by drains in the embankment walls.
In a typical cross-section through the dam (Figure 1), there is an unsaturated zone of tailings varying in thickness from the embankment end to the pool end.
The top tailings dam surface (beach profile) along a cross-section would be subject to the water balance components of precipitation (P), evapotranspiration (ET), infiltration (I), runoff (R), recharge (Re), and seepage (S).
However, spatial variation in the magnitude of these components would be expected between the embankment and the pool. This is because of the availability of moisture in the profile, which is governed by the phreatic level.
Therefore, at a point close to the embankment, evaporation would be expected to be a minimum, and towards the pool the evaporation should increase until it reaches a maximum (potential evaporation) at the pool edge. Similarly, infiltration would be expected to be a maximum close to the embankment, decreasing towards the pool (see graph in Figure 1).
Flux boundary modelling
Estimating tailings dam water balances has always been an important issue, be it for operational or closure water management. Most of the saturated zone water balance components are relatively well understood and can be estimated or measured with relative ease and with a high degree of confidence. However, the same cannot be said for the surface flux components above the unsaturated zone.
The measurement of these fluxes is difficult, expensive and time consuming, and as a result engineers look towards numerical modelling to provide the answers. Important advances have been made, with the development of codes like SoilCover, HELP, UNSAT-H, SWACROP, HYDRUS, and SWIM, to name but a few.
These models all attempt to calculate the surface flux components using a multitude of methods and assumptions. The most important single component is accurate calculation of evaporation, and the most accurate way of calculating this is using the modified Penman formulation as proposed by Wilson et al (1994). The only known model that uses the modified Penman formulation is SoilCover, making it the most appropriate tool.
Due to the detailed field data required for accurate use of a model like SoilCover, and because it is only a 1D model, the surface flux is often over-simplified using coarse recharge numbers when complex problems such as the gold tailings dam are modelled.
Commonly, water balance problems are solved using saturated/unsaturated groundwater modelling software. However, these seepage modelling packages do not allow for the actual calculation of the surface fluxes, but require some form of estimated recharge input. To obtain this recharge value the modeller will calibrate towards a known parameter, typically a phreatic level, and as such the most suitable recharge value might not represent the surface flux situation accurately.
To overcome this problem, and to bridge the gap between the two modelling systems, a system has been developed that allows for the calculation of flux functions that accurately describe the surface fluxes for any tailings dam cross-section. These flux functions can then be used as an actual boundary condition in regional 2D or 3D numerical modelling.
Calculation of the flux functions rests on a principle of selecting a typical tailings dam cross-section and solving the surface fluxes along this cross-section using the SoilCover model.
The top physical boundary of the typical cross-section consists of the tailings dam beach profile. Due to the hydraulic deposition of tailings, particle segregation takes place and an exponential expression can be developed to define this shape, verified by accurate survey data and physical tailings property testing.
The phreatic line, the location of which is determined by the tailings properties as well as the presence of any drains and the pond elevation, determines the base of the unsaturated zone to be modelled. Using observational techniques with piezometers it is possible to develop a mathematical function that will describe the shape of this boundary.
To carry out SoilCover modelling on the typical tailings dam crosssection the section is divided into a number of zones (the width of each zone is determined based on the section shape). For the study described here 13 zones, each 50m wide, were selected. An individual SoilCover model is run for each zone, and by integrating the computed surface fluxes over the entire tailings cross-section a good estimate of the cumulative result would be obtained. This approach therefore allows for a 2D solution using the 1D SoilCover model.
The material properties required for the SoilCover modelling were determined from an extensive field and laboratory testing programme, consisting of 66 grain size distribution and 25 soil water characteristic curve tests. The tailings material varies from well graded sands, with an average D 50 of0.26mm, to fine sands, and an average D 50 of0.03mm.
The soil water characteristic curves indicated a saturated volumetric water content ( u s) ranging between 34% and 56%, an Air-Entry Value (AEV) ranging between 1.5 and 11kPa and a residual suction ( c r)ranging between 2.5kPa and 600kPa. The database of tested tailings properties was used to define three main tailings types for modelling purposes: coarse, intermediate and fine.
Using a single averaged material type for the entire tailing dam cross-section would not be accurate, as measured data had shown that the tailings becomes progressively finer away from the embankment towards the beach (due to natural particle segregation). The choice of three material types was also based on work done by Kealy & Busch (1971), who modelled seepage from mill-tailings, and found that three tailings types best describe the model.
The three soil water characteristic curves for these tailings types were respectively selected based on the 75%, 50% and 25% tile values of the u s, AEV and c rmeasured values (Table 1). The steepness of the curves caused modelling instability and these curves had to be flattened in the high matric suction range.
The positions of the transitions between the tailings types down the beach profile was determined by solving of the quasi-2D SoilCover model. Kealy & Busch (1971) used a 6:6:1 ratio for their three tailings types, moving from the embankment towards the pool. For the typical tailings dam cross-section in this study a final tailings zone ratio of 5:5:3 was adopted as a good transition for the tailings types.
To define the vertical saturated hydraulic conductivity, k s(m/s), on the tailings beach profile for each of the zones, a theoretical function for saturated hydraulic conductivity was developed. This was verified using measured field data consisting of 29 laboratory saturated permeability tests, 62 Guelph permeameter tests, eight double-ring infiltrometer tests and 14 rainfall simulator tests. The function is described by the following expression:
ks= 1.94 x 10 -5 . e (-0.00977.H ) [Equation 1] Where H (m) is the distance along the beach profile as measured from the embankment, and e is the base of the natural log. Any vertical tailings profile is not homogeneous, and the physical and hydraulic properties of each of the horizontal tailings layers can vary significantly. However, physical tailings characterisation and model calibration supported the assumption that for the purpose of the modelling described here homogeneous vertical profiles would suffice.
The flux functions
Solving of the individual SoilCover runs, using each of the chosen tailing types, as well as the final optimal combined solution, each presents different spatial distribution functions for the surface flux components. Figure 2 presents the distribution of the actual evaporation/potential evaporation ratio (AE/PE), combined with the net infiltration ratio (NIR). The net infiltration ratio (NIR) is used to present the net infiltration (NI) data on the same basis as the AE/PE ratio and is defined as:
NIR = 1 - NI z [Equation 2] NI tIn equation 2, NI zis the zonal net infiltration for each of the 13 modelled zones, and NI tis the total net infiltration integrated over the entire typical tailing dam cross-section. The NI is defined as:
NI = P - R - ET [Equation 3] It is evident that the AE/PE ratio is the least close to the embankment and gradually increases to a value of 1 near the edge of the pool.
This is consistent with the proposed hypothesis in Figure 1. Similarly the NIR trends in Figure 2 follow the spatial infiltration hypothesis, with the maximum net infiltration occurring at the embankment end, and the least happening at the pool end of the tailing dam.
The application of the above flux functions allows for the calculation of a water balance for the typical tailings dam cross-section.
Expressing all data in terms of annual precipitation suggests that 40% runoff occurs, 114% evapotranspiration and -55% net infiltration. The negative net infiltration indicates a net negative water balance from the system, which in this case would imply the lowering of the phreatic level in the long term.
3D seepage modelling
To prove that the spatial flux functions in Figure 2 are a reasonable approximation of the actual surface fluxes, they had to be used as an input in regional seepage modelling software, and the seepage rates from the drains of the tailings dam predicted.
If the seepage rates were a good match, the flux functions could be considered to have fulfilled their function. Modelling of seepage in the tailings dam presented a unique challenge. First, the complexity of the model dictated that a 3D seepage model be used. A 2D cross-section of the tailings dam would not capture the essence of the problem nor would it yield accurate flow rates.
The 3D modelling of the flow regime through the tailings pond presented a significant challenge in itself. The model requirements included complex, irregular geometry, unsaturated flow, irregular flux sections, and highly complex boundary conditions.
For practical purposes, development of the 3D model also had to be done within a reasonable time. The SVFlux model developed by SoilVision Systems was selected based on its ability to model complex, highly irregular 3D problems with a relatively short learning curve.
3D model soil properties and boundary conditions The selected material properties were identical to those used in developing the surface flux functions described earlier, resulting in the use of three soil-water characteristic curves and a saturated surface hydraulic conductivity function.
There was detailed survey data (2000 points) for the base and surface topography of the problem and this was used to describe the model.
Furthermore, 13 zones had been isolated containing flux boundary conditions. However, the detailed data was simplified to a geometry grid of 90 points and three flux zones to allow solution of the model within a reasonable time frame.
The simplification of the geometry grid allowed reasonable representation of the model topology. The flux zones were averaged in such a way as to allow the total volume of flux to remain the same. SVFlux allowed the flux boundary conditions to be entered as free-form equations for each soil region. A water table formed the bottom boundary condition and was entered according to actual piezometer data.
3D model solution
A steady state model was first run to provide initial conditions for the transient state model. The automatic mesh generation and refinement of SVFlux allowed a steady state solution in 18s on a PIII 866 with 2889 nodes, 1368 cells, a single element error of 8.1x10 -5 m and a maximum problem error of 8.8x10 -4 m. The steady-state model solution is shown in Figure 3.
The transient state model was run for 365 days (the same as for developing the flux functions). The flux functions were added on the top boundary to simulate the combined effects of precipitation and evaporation, ie net infiltration. The transient model solved in 65mins on a PIII 866. An example of the flux function for one of the zones is presented in Figure 4.
Combined flux functions caused a net loss of 139,305m 3of water from the system over one year. Upward flow through the unsaturated zone is significant, as the net flow is negative. The negative flux boundary was countered by the water source, the pond at the centre of the tailings (the level of which was kept constant for the model run duration).
Seven seepage flux surfaces were placed in the 3D model to monitor seepage outflow from the tailings dam. A total volume of 23,452m 3ofwater exited the system over one year with an average flow rate of 4.3 litres/s. Actual seepage flow rates from the tailings dam measured over one year and four months over the life of the mine suggest an average annual seepage rate of 6 litres/s.
Considering the inaccuracies involved in the actual seepage measurements, which includes 15% over-estimation due to surface runoff intercepted in the seepage drains, combined with the complexity of the problem as a whole, the flux function seems to provide an excellent solution.
The combined advantage of using this function as a direct input into 3D seepage modelling software to assist in long-term water balance calculations makes it a worthwhile effort. Another great benefit of this study is the way the 3D model could be simplified using the same boundary conditions and material properties used in developing the flux function, effectively eliminating the guesswork normally associated with setting up complex 3D problems.
References Blight GE (1997). Interactions between the atmosphere and the earth. Geotechnique, Vol 47(4), pp 715-767.
Fayer ML and Jones TL (1990). UNSAT-H version 2: Unsaturated soil water and heat flow model, PNL-6779, Pacific Northwest Laboratory, Richland, Washington, USA.
Feddes RA, Wesseling JG and Wiebing R (1984). Simulation of transpiration and yield of potatoes with the SWACROP-model. 9th Tri-annual Conference of the European Association of Potato Research (EAPR), Interlaken, Switzerland, July 2-6.
Holtz RD and Kovacs WD (1981). An introduction to geotechnical engineering. PrenticeHall, Inc, Englewood Cliffs, New Jersey, USA.
Kealy CD and Busch RA (1971). Determining seepage characteristics of mill-tailings dams by the finite element method. Report of Investigations 7477, US Department of the Interior, Bureau of Mines, January, 51pp.
Ross PJ (1990). SWIM - a simulation model for soil water infiltration and movement. CSIRO Division of Soils, Davies Laboratory, Townsville, Queensland, Australia.
Rykaart EM (2001). Spatial infiltration to the Kidston tailing dam. PhD Thesis (in preparation), University of Saskatchewan, Department of Civil Engineering, Saskatoon, Saskatchewan, Canada.
Schroeder PR, Lloyd CM and Zappi PA (1994). The hydrological evaluation of landfill performance (HELP) model user's guide for version 3, EPA/600/R-94/168a, USA.
Simunek J, Huang K and van Genuchten MTH (1998). The HYDRUS code for simulating the one dimensional movement of water, heat, and multiple solutes in variably-saturated media, Version 6.0. Research Report No 144, USA Salinity Laboratory, USDA, ARS, Riverside, California, 164pp.
SoilCover (1997). SoilCover User's Manual. Unsaturated Soils Group, Department of Civil Engineering, University of Saskatchewan, Saskatoon, Saskatchewan, Canada.
SoilVision Systems (2001). SVFlux 3D: saturated/unsaturated automated 3D seepage modelling. SoilVision Systems.
Staley, R (1957). Effect of depth of water table on evaporation from fine sand. MSc thesis, Colorado State University, Fort Collins, Colorado, USA.
Wilson GW, Fredlund DG and Barbour SL (1994). Coupled soil-atmosphere modelling for soil evaporation. Canadian Geotechnical Journal, Vol 31, pp 151-161.