Modeling groundwater flow and salinity evolution near TSF Żelazny Most . Part I – groundwater flow

Tailings which are by-product of the extraction of various metals (copper, gold, silver, molybdenum, etc.) are often stored in so called Tailings Storage Facilities (TSF), where they are deposited as a soil-water mixture by spigotting. In many cases the water discharged together with tailings to the TSF is rich in salts and other chemical compounds imposing negative pressure to the groundwater environment. Even in the case of total or partial lining of such facilities and well-developed drainage systems to control leaching, some portion of contaminated water often seeps either through the surrounding dams or the bed into adjacent groundwater bodies. Numerical models can be very helpful tools to assess the extent of the contamination and particularly to predict its potential development in the future. This paper and the companion one describe such a numerical model developed for Żelazny Most Tailings Storage Facility (south-west Poland), one of the world's largest tailings sites. In the first part general information about the facility is provided and a 3D hydrogeological numerical model of the structure is described. Groundwater flow pattern near the facility obtained from numerical simulations is confronted with the measurements from a comprehensively developed monitoring system. Part II will be focused on the modelling of chloride transport in groundwater.


INTRODUCTION
TSF Żelazny Most is a depository for post-flotation tailings, which are by-product of copper mined in three mines: Rudna, Lubin and Polkowice located in south-west Poland, operated by KGHM company.It belongs to one of the largest facilities of its kind in the world.Since this is only place to store tailings produced by the Polish copper mines it is a key element in this production.The tailings are transported to the TSF in the form of low concentration slurry which is discharged into the facility by spigotting along the whole perimeter of the object.Semi-liquid tailings are retained within the structure by earth dams surrounding the depository.During the spigotting the tailings undergo natural process of segregation which causes that the coarser material deposits near the dams whereas the finer one is gravitationally transported with flowing water towards the central part of the TSF, where the water is clarified forming a water pond (Figure 1).In order to accept new supplies of tailings (the annual production of which amounts 28 mln tonnes) the facility has to be continuously developed.The Żelazny Most tailings facility is being raised using the upstream construction method.Only at early stage the so-called starter dams were constructed from borrow earthen material.Afterwards the dams have been raised by 5 m high embankments utilizing coarse tailings deposited near the dams.Total area occupied by the facility amounts 14 ha with the length of perimeter 14.3 km.Due to the natural terrain relief the present height of the surrounding dams varies: 41 m for the south dam, 47 m for the north one, 57 m for west dam and over 70 m for the east dam, at the crest elevation 185 m a.s.l.Total volume of tailings stored in the TSF exceeds 600 mln m 3 and according to the future plans of copper production in KGHM it should be ready to accept next 300 mln m 3 .The dams as well as the bed are not sealed which means that the water discharged with tailings can almost freely infiltrate.In order to contain infiltration an extensive system of drainage has been developed.It consists of various elements, including several floors of circumferential drainage, the aim of which is to intercept the seeping waters at the highest possible elevation in order to control the position of phreatic surface within the dam body.At present there exist four floors of circumferential drainage, constructed every 10 m in height starting from the elevation 152 m a.s.l. and installed in the course of the facility development.Moreover, there is pipeline drainage of starter dam and finally the dich drains surrounding the facility at the dams' toe (Figure 2). Figure 2. Schematic cross-section through the TSF Żelazny Most dam.Despite the comprehensive system of horizontal drainage some waters still seeps into the subsoil and next downstream, infiltrating into the groundwaters.The waters discharged with tailings to TSF come from dewatering of mining fields therefore they are highly contaminated, mostly by salts.In order to reduce the negative impact of the waters seeping downstream a system of drainage wells (approximately forty wells) located around the facility outside the dams' toe has been installed.The second purpose of these wells is to reproduce the groundwater flow conditions prior to the construction of the TSF.The drainage wells continuously pump out the groundwater from the small local aquifers preventing the contaminated water from flow downstream, however small amounts of leachate still pass the hydraulic barrier.In order to recognize the existing flow paths of saline water, verify the efficiency of existing drainage system, but first of all to predict the development of salinity front with the raise of the dams and tailings, a numerical groundwater model of the TSF Żelazny Most and surrounding area has been developed.

STRUCTURE OF 3D HYDROGEOLOGICAL NUMERICAL MODEL
The 3D model, which was originally developed at the beginning of the current decade and successively updated, reproduces flow conditions in the area of TSF Żelazny Most including the following elements (Świdziński et al., 2011): -infiltration of saline waters through the bed of water pond, the beaches and next through the mass of tailings into the groundwater, -downstream flow of groundwater and pollution migration, -surface watercourses and the change of their quality, -interaction between ground and surface waters, -drainage system of the facility, -drinking water intake wells.
The 3D model has been developed with GMS (GroundWater Modelling System) commercial software package, initially using version 6.5 and more recently using version 8.3.The groundwater flow equations are solved by MODFLOW software integrated with GMS.The model covers the area of the current impact of the facility determined by the field measurements and the predicted area of long-term future impact.The borders of the model are defined by watercourses or local watersheds.The model domain also includes Retków-Stara Rzeka water intake located within the area of Main Groundwater Reservoir (GZWP -314) (Figure 3).Steady-state groundwater flow is considered.The bottom boundary of the 3D model was determined based on the assumption that the most important role in the impact of the facility on groundwater is played by the shallowest multi-aquifer formations i.e. quaternary and upper tertiary aquifers.The latter is well isolated from the lower inter coal (middle and early Miocene) and under coal (Oligocene) aquifers.The elevation of the model bottom was assumed at -70,0 m a.s.l., i.e. the lowest point of the tertiary aquifer formation.The model domain has 154.8 km 2 (12.9 x 12.0 km).It was subdivided into a grid of squares and rectangles with alternate sizes.Due to essential gradients in elevation and resulting hydraulic gradients, the area of very facility was covered by square mesh with the size 25 x 25 m, zone 1 km distant from the dams toe by squares 50 x 50 m, whereas the rest of the model area with squares 100 x 100 m (Figure 3).Such discretization allowed an appropriate modelling of the work of drainage wells as well as good reproduction of the inclination of downstream slopes of the facility.The orientation of the mesh has been chosen in such a way so that the mesh line follows main direction of flow and migration of saline waters (Świdziński et al. 2011).The grid for a single calculation layer is built of 289 rows and 282 columns.Geological and hydrogeological conditions of the modelled area are extremely complex since they were significantly impacted by three glaciations which passed over Żelazny Most in Pleistocene.The various ice sheets, which are believed to have been at least 1000 m thick, have induced widespread glacio-tectonic phenomena causing many thrusts of Pleistocene deposits into the initially horizontally bedded freshwater Pliocene sediments (Jamiolkowski 2014).As a consequence it was very difficult to identify and separate single aquifers and aquitards in the subsoil (see Figure 2) as it is done in a standard approach in the schematisation process and generalisation of hydrogeological conditions.Moreover, despite a huge number of geological profiles (several thousands) obtained from various field investigations (standard boreholes, CPTUs, geophysical surveying) only small portion of them could have been used for identification of hydrogeological conditions since the majority of the profiles were located near the dams of the facility and a lot of other ones were too short to be included in the interpretation.Furthermore, very differentiated geological conditions which were changing very often and rapidly in plane and depth, prevented to use standard GMS tools to construct individual hydrogeological layers of similar hydraulic properties.Therefore, it was decided to construct 3D model of the subsoil in a non-standard way.First of all, for the generalization purposes four basic hydrogeological layers have been assumed: -quaternary layers of permeable formations -k = 10 -4 m/s, -quaternary layers of impermeable formations -k = 10 -9 m/s, -tertiary layers of permeable formations -k = 5x10 -4 m/s, -tertiary layers of impermeable formations -k = 10 -9 m/s.geometric average.This procedure was applied for each hydrogeological borehole taken into account.Distributed averaged values of hydraulic conductivities were next interpolated between neighbouring boreholes for each model layer using the natural neighbour method.
In such a way a fuzzy distribution of hydraulic conductivities was obtained covering the full and continuous range of permeabilities i.e.: from impermeable soils, via semi-pervious ones, and finally to fully permeable soils.The model of the foundation has been supplemented by the system of rivers and surface watercourses running within the model as well as along its borders, drainage wells around the facility, groundwater intake and farm wells in this area.
In turn, hydrogeological model for the facility also consists of 10 calculation layers, the lowest three of which have differentiated thickness so as to have horizontal top of the third lowest layer at the elevation of 160 m a.s.l.The subsequent 7 higher layers of the depository have equal thickness of 5 m reaching final elevation of 195 m a.s.l.Each of the calculation layers consists of external embankment and the tailings filling the depository with decreasing permeability going towards the pond.In order to reflect the change of the fines content in deposited tailings, for each calculation layer 9 zones with the same permeability have been distinguished following the decrease of permeability coefficient as a function of the distance from the embankment.Finally, based on accessible digital maps and half-tone screen of the terrain of concern the system of surface watercourses as well as all drainage elements installed in the depository together with water intakes were implemented in the model.For the majority of river sections flowing along the borders of the model the first-type (Dirichlet) boundary condition was assumed (specified head).At the north-west corner, the only artificial border of the model, the so-called specified-flux boundary was assumed which in fact is a second-type boundary condition.For all watercourses located inside the model area the third-type boundary condition (Cauchy) was specified.The work of any wells (drainage, relief and intake wells) was modeled by the second-type (Neumann) boundary condition (constant discharge).Recharge was assumed based on the long-term mean annual precipitation at the level of 597 mm reduced respectively to the effective infiltration depending on the type of the surface soils.With regard to the facility, the infiltration of saline waters throughout the mass of tailings under the pond was simulated by first-type boundary condition assuming elevation of water head whereas the infiltration of waters discharged on the facility beaches during spigotting process was simulated by third order boundary condition using MODFLOW RIVER package.

Figure 3 .
Figure 3. Finite difference mesh of 3D model for TSF Żelazny Most.

Figure 4 .
Figure 4. Observed and calculated elevations of piezometric head in a) drainage wells b) in observation points downstream the facility.