Open access

Application of Computational Fluid Dynamics (CFD) for Simulation of Acid Mine Drainage Generation and Subsequent Pollutants Transportation through Groundwater Flow Systems and Rivers

Written By

Faramarz Doulati Ardejani, Ernest Baafi, Kumars Seif Panahi, Raghu Nath Singh and Behshad Jodeiri Shokri

Submitted: October 16th, 2010 Published: July 5th, 2011

DOI: 10.5772/16927

Chapter metrics overview

5,250 Chapter Downloads

View Full Metrics

1. Introduction

Many environmental problems associated with the mining industry involve the understanding and analysis of fluid or gas flow. Typical examples include groundwater flow, transport of contaminants, heat transfer, explosions, fire development and dust movements. Both experimental work and numerical models can provide the necessary information for solution of any particular problem. The long-term pyrite oxidation, acid mine drainage generation and transportation of the oxidation products are noted to be the most important problems that can be modelled in order to predict the transport of the contaminants through groundwater and rivers flow systems, to interpret the geochemistry and achieve a better understanding of the processes involved.

The use of computational fluid dynamics (CFD) to simulate flow problems has risen dramatically in the past three decades and become a fairly well established discipline shared by a number of engineering and science branches. Associated with the widespread availability of high performance and advanced computers and computational methods, CFD is rapidly becoming accepted as a cost-effective design and predictive tool. Numerical solution methods may be used for CFD analysis for the simulation of fluid flow, heat and mass transport problems when it is expressed in terms of partial differential equations (PDEs). Recent improvement of CFD codes enables researchers to visualise easily the local velocity, temperature, concentrations of solutes and pressure fields in a domain by means of graphic facilities (Edwards et al., 1995).

CFD is a very powerful tool and applicable to a wide range of industrial and non-industrial areas, including aerodynamics of aircraft, automotive, pollution control, power plant, turbo machinery, electrical and electronic engineering, civil engineering, hydrology and oceanography, meteorology and medical science (Versteeg & Malalasekera, 1995).

CFD has been recently used in many applications relating to environmental studies as a way of estimating impacts and developing control strategies for the long-term impacts of mining on the environment and many other activities. The low response time and cost associated with the simulations compared to experiments are the main benefits provided that adequate accuracy may be obtained by the computer erally used in CFD can be classified as finite difference method (FDM), finite element method (FEM), finite volume method (FVM), or spectral methods. The choices of a numerical method and a gridding strategy are strongly interdependent to solving the CFD problems. For example, the use of finite difference method is typically restricted to structured grids (Harvard et al., 1999).

Some researchers used the finite element method to solve the PDEs for modelling of solute transport through groundwater flow systems (Pickens & Lennox, 1976; Pinder, 1973; Rabbani & Warner, 1994; Wunderly et al., 1996). Furthermore, some researchers used the finite volume method to solve the PDEs for modelling of pollutant transport through rivers flow systems (Ani et al., 2009; Kachiashvili et al., 2007; Zhi-Qiang & Hoon-Shin, 2009). Green and Clothier (1994) used the PHOENICS-code (Spalding, 1981) incorporating the finite volume method to simulate water and solutes transport into unsaturated soils. Edwards et al. (1995) have noted the applicability of the CFD analysis in mine safety and health problems such as methane control, gas or coal outbursts, dust suppression and explosions. Balusu (1993) developed a numerical model using a CFD code, FIDAP, to simulate airflow patterns and the respirable dust concentration at a longwall face in underground coal mines. The analysis of contaminant transport in groundwater systems using finite volume techniques have been carried out by Putti et al., (1990) and Binning & Celia (1996). Singh & Doulati Ardejani (2004) developed an one-dimensional numerical finite model using PHOENICS code to simulate long term pyrite oxidation, acid mine drainage generation and transportation of the oxidation products through the backfills of an open cut mine.

This chapter presents the application of computational fluid dynamics (CFD) using PHOENICS code to simulate pyrite oxidation, acid mine drainage generation and subsequent pollutants transportation through groundwater flow systems and Rivers. The simulations have been conducted by developing one- and two-dimensional models.


2. Computational fluid dynamics

Computational fluid dynamics (CFD) is defined as the analysis of systems involving fluid flow, heat transfer and mass transport and associated phenomena such as chemical reactions using computer-based simulation. To predict the way in which a fluid will flow for a given situation, a mathematical analysis of the fluid flow has to be made to formulate the governing equations of flow, and the CFD code enables users to calculate numerical solutions to these equations. To produce a solution, these equations have to be transformed into numerical analogues using discretisation techniques such as finite difference method (FDM), finite element method (FEM) and finite volume method (FVM).

Figure 1.

Main stages in a CFD simulation

CFD codes such as FLUENT, PHOENICS, FLOW3D and FIDAP are now widely available commercially, each with its own particular set of features to deal with fluid flow problems. Edwards et al., (1995) and Doulati Ardejani (2003) has given a comparison of some commercial CFD codes. In general, a flow analysis with CFD codes can be divided into three main phases:

  • Pre-processing phase includes the input of a flow problem to a CFD package using an operator-friendly interface and the subsequent transformation of this input data into an appropriate form for use by the solver. This phase mainly consists of the definition of the geometry of the problem of interest, mesh generation, specification of physical properties of the fluid and appropriate boundary conditions. An unknown flow variable such as velocity, pressure and temperature is solved at nodes inside each cell. The accuracy of a CFD solution is governed by the number of cells in the grid. The accuracy is improved by increasing the number of cells. Optimal meshes are often non-uniform. A finer mesh is constructed in areas where large variations occur from point to point and a coarser grid is used in regions with relatively small change (Versteeg & Malalasekera, 1995).

  • Simulation phase including solution of the governing equation for the unknown flow variable.

Post-processing phase including domain geometry and grid display, vector plots, surface plots, x-y graphs, line and shaded contour plots and animation for dynamic result display.

2.1. General stages of a CFD analysis

To produce a CFD simulation, a number of key steps should be followed in order to generate an exact picture of a particular problem. Figure 1 shows the main steps for a CFD analysis.

2.2. PHOENICS as a CFD software

The PHOENICS program has a few different modules to perform all these three phases of flow analysis, namely SATELLITE, EARTH, and post-processing facilities including VR VIEWER, PHOTON, AUTOPLOT and RESULT (Spalding, 1981). Figure 2 shows the solution performance in the PHOENICS package.


3. Finite volume method (FVM)

FVM (sometimes called the control volume method) is a numerical technique for solving governing equations of fluid flow and mass transport that calculates the values of the conserved variables averaged across the control volume. The calculation domain is broken down into a finite number of non-overlapping cells or control volumes such that there is one control volume surrounding each grid point. The partial differential equation is integrated over each control volume. The resulting discretisation equation containing the values of a variable for grid points involves the substitution of a variety of finite-difference-type approximations for the different terms in the integrated equation representing flow and transport processes such as convection, diffusion and source terms. The integral equations are therefore, converted into a system of algebraic equations. The algebraic equations obtained in this manner are then solved iteratively.

It is possible to start the discretisation process with a direct statement of conservation on the control volume. Alternatively we may start with the differential equation and integrate it over the control volume. The numerical solution aims to provide us with values of at a discrete number of points in the flow domain. These points are called grid points, although we may also see them referred to as nodes or cell centroids, depending on the method (Murthy, 2002).

Figure 2.

The procedure of solution in PHOENICS (CHAM, 2008)

FVM has become popular in CFD as a result, primarily, of two advantages. First, they ensure that the discretisation is conservative, i.e., mass, momentum, and energy are conserved in a discrete sense. While this property can usually be obtained using a finite difference formulation, it is obtained naturally from a finite volume formulation. Second, FVM does not require a coordinate transformation in order to be applied on irregular meshes. As a result, they can be applied on unstructured meshes consisting of arbitrary polyhedral in three dimensions or arbitrary polygons in two dimensions. This increased flexibility can be used to great advantage in generating grids about arbitrary geometries (Harvard et al., 1999).

Although the finite element method is mainly used for the simulation of the contaminant transport problems through groundwater flow systems (see Barovic & Boochs, 1981; Bignoli and Sabbioni, 1981; Fried, 1981; Kerdijk, 1981; Nawalany, 1981; Putti et al., 1990), finite volume techniques have also been used by many researchers for the solution of groundwater flow and solute transport governing equations in saturated and unsaturated flow systems (see Binning and Celia, 1996; Green and Clothier, 1994; Putti et al., 1990; Svensson, 1997).

The basic principles and concepts of the finite volume formulation are very simple and easier to understand by engineers than other numerical techniques. One advantage of the FVM over finite difference methods is that a structured grid is not required although a structured mesh can also be utilised. Furthermore, the FVM is noted to be superior to other numerical methods due to the fact that boundary conditions can be applied non-invasively because the values of the conserved variables are located within the control volumes, and not at nodes or surfaces. FVM is particularly useful in dealing with coarse non-uniform grids. The conservation of any flow variable within a control volume is expressed as a balance equation between the various processes causing an increased or decreased value of (Versteeg & Malalasekera, 1995).


The quantities being balanced are the dependent variables such as the mass of chemical species, energy, momentum and turbulence quantities.

CFD packages such as PHOENICS, FLUENT, FLOW3D and STAR-CD contain discretisation techniques appropriate for dealing with the main transport processes such as convection (i.e. transport due to fluid flow), diffusion (i.e. random motion of molecules) as well as for the source terms (i.e. chemical reaction for energy or chemical species) and the rate of change of with time (i.e. accumulation within a cell).

3.1. The finite volume discretisation

Fundamental to the development of a numerical method (such as FDM, FEM and FVM) is the idea of discretisation. When the number of grid points is small, the departure of the discrete solution from the exact solution is expected to be large. A well-behaved numerical scheme will tend to the exact solution as the number of grid points is increased. The rate at which it tends to the exact solution depends on the type of profile assumptions made in obtaining the discretisation. No matter what discretisation method is employed, all well-behaved discretisation methods should tend to the exact solution when a large enough number of grid points is employed (Murthy, 2002).


4. Model verification

Two problems are described to show the applications of CFD analysis in the modelling of the environmental problems. The first problem considers oxygen transport, pyrite oxidation, acid mine drainage generation and long-term leaching of the oxidation products from the backfills of an open cut mine. The second problem deals with the modelling of metal pollutants transportation associated with acid mine drainage in Shour River at Sarcheshmeh copper mine, Iran. A multi-purpose computational fluid dynamics (CFD) package called PHOENICS has been modified to model these problems.

4.1. Problem 1- oxygen transport, pyrite oxidation, acid mine drainage generation and long-term leaching of the oxidation products from the backfills of an open cut mine

Acid mine drainage in the backfills of an open cut mine is initially generated by the direct oxidation of pyrite by oxygen and water producing Fe2+,SO42and H+(Singer & Stumm, 1970). The stoichiometric oxidation reaction is as follows:


The ferrous iron formed by Reaction 2 oxidises to produce ferric iron according to Reaction 3 (Rogowski et al., 1977):


According to Singer & Stumm (1970), Reaction 3 is slow under low pH conditions. It is catalysed by the bacteria Thiobacillus ferrooxidans under suitable environmental conditions (Singer & Stumm, 1970). The environmental conditions must be favourable, for example bacterial action is most important in solutions having a pH between 2.0 and 3.5 (Malouf & Prater, 1961). If conditions are not favourable, the bacterial influence on pyrite oxidation and acid generation will be minimal; for example, the bacteria become inactive at temperatures above 50 C (Cathles, 1975; Malouf & Prater, 1961). It has also been reported that the rate of ferrous iron oxidation may be accelerated in the presence of bacteria when oxygen is readily available (Jaynes et al., 1984a). The ferric iron produced by Reaction 3 can choose two path ways: It may react with pyrite and produce more ferrous iron, sulphate and hydrogen concentrations (Singer & Stumm, 1970). The stoichiometry of this process reduces to:


Alternatively Fe3+ may be precipitated as amorphous ferric hydroxide (Jaynes et al., 1984a; Walter et al., 1994a, 1994b) according to Reaction 5 as given below:


A number of numerical models of pyrite oxidation have been reported in the literature. They are mainly associated with mine tailings (Bridwell & Travis, 1995; Elberling et al., 1994; Walter et al., 1994a, 1994b; Wunderly et al., 1996), waste rock dumps (Cathles, 1979; Cathles and Apps, 1975; Davis 1983; Davis et al., 1986b; Lefebvre & Gelinas, 1995), open cut mines (Jaynes et al. 1984a, 1984b, Singh & Doulati Ardejani, 2004) and coal refuse piles (Doulati Ardejani et al., 2008b, 2010).

In this section a numerical finite volume model is presented for prediction of the long-term pyrite oxidation and subsequent pollutant transportation from backfilled open cut mines. This work is a further development of the model proposed by Jaynes et al., (1984a, 1984b) and is presented in both one- and two-dimensional forms with different numerical techniques. The results of a two-dimensional simulation are necessary for better understanding the mechanisms involved in pyrite oxidation and pollutant leaching related to open cut mines. The results of a two-dimensional simulation show how contaminants spread into groundwater systems and how the physical and chemical processes influence the transport of the pyrite oxidation products. Such results can be used for designing effective site remediation programs in order to minimise environmental effects arising from abandoned backfilled open cut coalmines.

The model incorporates oxygen transport, pyrite oxidation, enthalpy balance, groundwater flow and transport of the oxidation products. The model has been implemented in the general-purpose PHOENICS computational fluid dynamic package. By creating a Q1-file (PHOENICS input file) in PHOENICS input language, the necessary settings such as domain geometry, finite volume grid, boundary and initial conditions and fluid properties were specified and for all non-standard equations and specific source terms, the required coding in FORTRAN 99 language was supplied by developing a subroutine called GROUND. This subroutine was used by the PHOENICS solver module during the course of the solution. The conceptual model proposed in developing the model is given in Doulati Ardejani et al., (2002). Figure 3 shows a schematic diagram showing pollutant leaching in backfills of an open cut mine.

Figure 3.

A conceptual model of pyrite oxidation in the backfills of an idealized open cut mine

Figure 4 shows the physical and chemical processes affecting the fate and transport of pyrite oxidation products through the backfills of an idealized open cut mine.

4.1.1. Modelling governing equations

Groundwater flow

The governing equation for two-dimensional steady state groundwater flow is based on the well recognised Laplace equation. This equation is obtained by coupling the continuity equation and Darcy’s law. It is expressed as follows (Freeze & Cherry, 1979):



xand y = Cartesian coordinates;

Kxand Ky= hydraulic conductivities in x and y directions;

h = hydraulic head;

W = recharge or discharge rate per unit volume.

The components of the groundwater velocity are obtained by solving Equation 6 and are then used in the mass transport equation for each product of pyrite oxidation.

Figure 4.

Chemical and physical processes affecting the fate and transport of a pollutant

Pyrite oxidation mechanism

Pyrite oxidation can be expressed by a shrinking core model (Jaynes et al., 1984a, 1984b; Levenspiel, 1972). This model can be modified to take into consideration spherical pyritic grains and incorporating both surface reaction kinetics and oxidant diffusion into the particles containing pyrite (Cathles & Apps, 1975):



X = fraction of pyrite remaining within the particle (Kg/Kg);

t = time (s);

τC = total time required for complete oxidation of pyrite within a particle if the oxidation process is only controlled by the decreasing surface area of pyrite (s);

τD = total time required for full oxidation of pyrite within a particle assuming the oxidation rate is solely controlled by oxidant diffusion into the particles (s);

[O2] and [Fe3+]= oxygen and ferric iron concentrations.

In Equation 7, both oxygen and ferric iron participate in the pyrite oxidation process. τDand τC are expressed by Equations 8 and 9 (Levenspiel, 1972):



ρPy = density of pyrite in the particle (mol/m3);

R = particle radius (m);

b = stoichiometric ratio of pyrite to oxidant consumption (mol/mol);

De[Ox] = effective diffusion coefficient of oxidant in oxidised rim of the particle (m2/s);

K[Ox] = first-order surface reaction rate constant (m/s) ;

αPyrock = surface area of pyrite per unit volume of particle (m1);

λ = thickness of the particle in which pyrite oxidation occurs (reaction skin depth) (m);

C[Ox] = concentration of oxidant in the water surrounding the particle (mol/m3).

Enthalpy model

A simple enthalpy balance is employed to model the heat transfer using Equation 10 (adopted from Cathles & Apps, 1975):



T = temperature (C);

ρb = bulk density of the spoil (Kg/m3);

ρS = molar density of pyrite in spoil (mol/m3);

ϕ = porosity of spoil;

KT = thermal conductivity (J/mCS);

CP = specific heat capacity of the spoil (J/KgC);

ρW = density of water in the spoil (Kg/m3);

uxjW= water velocity(m/s);

CW = heat capacity of water (J/KgC);

B = moles of pyrite consumed in pyrite-oxygen reaction per heat generated(mol/J);

B = moles of pyrite consumed in pyrite-Fe3+ reaction per heat generated(mol/J).

Oxygen balance

Oxygen has an important role in pyrite oxidation within the spoils of an open cut mine. Oxygen is transported within the spoils by the process of diffusion. It is consumed by the pyrite oxidation reaction, chemical oxidation of ferrous iron and Thiobacillus ferrooxidans bacteria (Jaynes et al., 1984a). The governing equation of oxygen transport incorporating the volumetric consumption terms reduces to:



ϕa = air-filled porosity of the spoil(m3/m3);

u = gaseous concentration of oxygen in the spoil pore spaces(mol/m3);

De = effective diffusion coefficient (m2/s);

SKPyO2,SKFe2+O2,SKBO2= volumetric oxygen consumption terms (mol/m3s) by the pyrite oxidation reaction, chemical oxidation of ferrous iron and oxygen consumption by bacteria. The oxygen consumption terms can be found in Jaynes et al. (1984a).

Transport equation

This model assumes that the oxidation products are transported through the groundwater flow system after leaching from the spoils. The governing transport equation can be expressed as follows (Freeze & Cherry, 1979; Rubin & James, 1973):



Ci = solute concentration in aqueous phase(mol/m3);

C¯i = solute concentration in adsorbed form(mol/Kg);

ρb = bulk density of the flow medium(Kg/m3);

qre = surface recharge(m/s);

S = sink and source terms representing the changes in aqueous component concentrations due to the chemical reactions(mol/m3s);

qj = vector components of the pore fluid specific discharge(m/s);

D = hydrodynamic dispersion coefficient(m2/s);

nc= number of dissolved species.

The H+- spoil interaction

The interaction between H+ produced by oxidation reactions and spoil was incorporated by a simple empirical relationship (Jaynes et al., 1984a):



H+ = Hydrogen ion in solution (mol/m3);

HGen+ = H+ generated through chemical reactions (mol/m3);

GA = an empirical constant that defines the buffer system.

Equation 13 was slightly modified to calculate the H+ consumption per cubic metre of spoil per unit time as given below:



Δt= simulation time step (s).

Equation 14 was incorporated as a sink term in the governing transport equation for the hydrogen ion. In this particular case the solution pH was maintained above 2.5.

4.1.2. Model parameters

The model input data were taken from Cathles & Apps (1975), Lide, (1998) and Jaynes et al., (1984b). A one-dimensional simulation was conducted to verify modelling accuracy. A two-dimensional simulation was then carried out when the reasonable output results were obtained in comparison with those results obtained from the POLS finite difference model developed by Jaynes et al., (1984a, 1984b). The model input parameters are given in Table 1.

Model parameter

Pyrite fraction
Kg pyrite/ Kg spoil
Jaynes et al., 1984b

Fraction of coarse particles containing pyrite

The bulk density of spoil

The porosity of the spoil

Water-filled porosity

Surface area of pyrite per unit volume of spoil
Cathles & Apps, 1975

Radius of particles
Jaynes et al., 1984b

Diffusion coefficient of oxidant in water

Molar density of pyrite within particle

First-order rate constant for Fe3+

First-order rate constant for oxygen

Recharge value

constants for chemical oxidation of Fe2+
K1=1.3x10-10 K2=1.7x10-9

Effective diffusion coefficient for oxygen transport

Enthalpy of Reaction (1)
Lide, (1998)

Enthalpy of Reaction (3)
Lide, (1998)

Moles of pyrite consumed per heat generated in FeS2-O2 reaction

Moles of pyrite consumed per heat generated in FeS2-Fe3+ reaction

Thermal conductivity

Specific heat capacity

Pyrite fraction0.0025Kg pyrite/ Kg spoilJaynes et al., 1984b
Fraction of coarse particles containing pyrite75%
The bulk density of spoil1650Kg/m3
The porosity of the spoil0.321-
Water-filled porosity0.225-
Surface area of pyrite per unit volume of spoil80cm-1Cathles & Apps, 1975
Radius of particles2cmJaynes et al., 1984b
Diffusion coefficient of oxidant in water1.0x10-11m2/s
Molar density of pyrite within particle23mol/m3
First-order rate constant for Fe3+4.4x10-5m/s
First-order rate constant for oxygen8.3x10-10m/s
Recharge value0.5m/yr
constants for chemical oxidation of Fe2+K1=1.3x10-10 K2=1.7x10-9mol2/(m3)2(s)
Effective diffusion coefficient for oxygen transport1.0x10-7m2/s
Enthalpy of Reaction (1)-1.45x106J/molLide, (1998)
Enthalpy of Reaction (3)-1.78x104J/molLide, (1998)
Moles of pyrite consumed per heat generated in FeS2-O2 reaction6.9x10-7mol/J
Moles of pyrite consumed per heat generated in FeS2-Fe3+ reaction5.6x10-5mol/J
Thermal conductivity2.1J/m°CS
Specific heat capacity2500J/Kg°C

Table 1.

Parameters used for simulation (Cathles & Apps, 1975; Jaynes et al., 1984b; Lide, 1998)

4.1.3. Model calibration

A one-dimensional simulation was carried out to calibrate the model and verify its capability for simulation of pyrite oxidation and leaching process associated with open cut coal mines. To perform the simulation, a 10-metre spoil column was simulated as 20 equal size control volumes. The influx and background concentrations of dissolved species and pH are listed in Table 2.

Aqueous components Concentration (mol/m3)
Influx dataBackground data

Table 2.

Influx and background chemistry (Jaynes et al., 1984b)

4.1.4. Boundary conditions

A first type boundary condition was assigned at the top of the spoil profile for the oxygen transport model equal to its atmospheric concentration (8.9mol/m3). The spoil column initially contained no oxygen concentration. A zero-gradient boundary condition was considered at the profile base.

The spoil temperature was initially set at 15C for the enthalpy model. The top surface of the spoil profile was specified as a first type boundary condition equal to15C.

A constant recharge of 2×108m/swas maintained at the profile top surface. A dispersivity of 0.005 m was applied for the transport model in order to achieve consistency with the POLS model. The main reason of selecting a small value of dispersivity is that the mechanism of the leaching process used in the POLS model is different from the present finite volume model. The POLS model did not consider the processes of diffusion and hydrodynamic dispersion effects.

The following boundary and initial conditions were specified for the transport of the pyrite oxidation products:

  1. The top boundary of the model was assigned as first type with respect to concentration of solutes

  2. The background solute concentrations were simulated by specifying an initial boundary condition

  3. A zero concentration gradient boundary condition was specified at outflow boundary.

4.1.5. One-dimensional simulation results

For one-dimensional simulation, the effective diffusion coefficient for oxygen transport was selected to be 1.0×107m2/sand air-filled porosity of 0.1 was set for the simulation. The temperature rise due to the pyrite oxidation reactions was ignored and a constant temperature equal to 15C was considered for modelling purposes. The iron-oxidising bacteria were allowed to be active in the spoil profile where the conditions were favourable. The mathematical expressions describing the role of bacteria were taken from Jaynes et al., (1984a) and relevant FORTRAN codes were supplied in the GROUND subroutine accessible by the user of the PHOENICS package.

Figure 5 shows the total iron-discharging rate as a function of time predicted in the water leaving the spoil profile for two different cases. In Case 2 the iron- oxidising bacteria are active and the interaction between the spoil and H+ is incorporated. In Case 1 no bacteria are allowed to be present but in this case the effective diffusion coefficient for oxygen transport is four times greater than that in Case 2. A comparison was made with those results predicted by the POLS model (dots). The iron-leaching rate showed a similar pattern with time for both cases and the maximum rate occurred between 1750 to 2100 days. As Figure 5 shows, the leaching rate in Case 1 is greater than in case 2 due to bigger in the effective diffusion coefficient.

The H+- spoil interaction increased the pH of the solution. It caused an increase in the bacterial activity therefore more ferric iron was produced. Consequently, the pyrite oxidation rate increased considerably over time (Figure 6).

Figure 5.

Total iron leaching rate versus time in the water leaving the spoil profile for Cases 1 and 2

As Figure 6 shows, about 36 % of the pyrite was oxidised after 10000 days of the simulation. For this time, the POLS model predicted that 39 % of pyrite was consumed. The difference between the finite volume modelling predictions and those predicted by the POLS model can be explained in that unlike the POLS model, the ferric complexation reaction was not incorporated in the finite volume model. The complexation reaction increases the ferric concentration. Consequently the rate of pyrite oxidation increases.

The mole fraction of oxygen within the spoil profile versus depth was illustrated in Figure 7. Oxygen decreased linearly with depth. As Figure 7 shows, oxygen diffusion was limited in the surface elements of the spoil profile up to about 3 metres from the spoil surface. Bacterial activity (SKBO2in Equation 11) was the main reason for this sharp reduction of oxygen concentration within the surface layers of the spoil profile, consuming most of the oxygen over this part of the profile (Case 2). In Case 1, oxygen diffused over entire profile because no bacteria were active to consume the oxygen.

Figure 6.

Comparison of the POLS results (dots) with the PHOENICS results (solid lines) for the pyrite oxidised vs. time over entire spoil column with and without H+-spoil interaction

Figure 7.

Comparison of the POLS modelling results (dots) with the finite volume simulated results (solid lines) for oxygen mole fraction as a function of depth for 5-year period

4.1.6. Two-dimensional simulation

Two-dimensional simulations were also performed to demonstrate the capability of the finite volume model for prediction of the long-term pyrite oxidation and transportation of the oxidation products from backfilled open cut coal mines. Except for few, all other model parameters are almost similar to those input data used for one-dimensional simulation (Table 1). Chemical species considered in the two-dimensional simulations are the same as those in the one-dimensional model. Slight modifications were made to their influx and background concentrations (Table 3).

Aqueous components Concentration (mol/m3)
Source Background

Table 3.

Source and background concentrations of aqueous components used for two-dimensional simulations

The two-dimensional cross-sectional dimensions are 50 m horizontally by 20 m vertically and this domain is discretised into 40×20 control volumes of size 1.25 m horizontally × 1 m vertically. The groundwater flow system was assumed to be steady. The difference head between the left and the right boundaries was maintained at 0.1 m. An average recharge value of 0.5 m/yr was considered for the upper boundary (spoil surface). For the simulation it was assumed that reactive pyrite was contained only in a 12.5-m-wide segment of the unsaturated zone of the spoil (Figure 8).

Figure 8.

Two-dimensional simulation cross-section showing velocity vectors and the segment of the spoil where oxidation reactions take place

The upper 4 m of the grid was assumed to be unsaturated and the remainder fully saturated with a constant porosity of 0.321. For simplicity the horizontal component of the velocity was ignored in the unsaturated zone and flow was only assumed to be vertical in this zone. Horizontal and vertical hydraulic conductivities of 1.40×105 m/s and 1.40×106 m/s were used for the simulation. Horizontal and vertical dispersion coefficients of 7.0×109 m2/s and 5.0×109 m2/s were specified for the model. The steady state flow system in terms of velocity vectors is shown in Figure 8.

The spoil temperature was assumed constant at 20 C and the temperature rise in pore water due to the oxidation reactions was predicted using the enthalpy balance. An effective diffusion coefficient of 9×108m2/s was assigned for the oxygen transport model. The bacterial action and the interaction between H+and the spoil were also considered for two-dimensional simulation. The spoil surface was maintained as a first-type boundary condition for oxygen equal to its atmospheric concentration (8.9 mol/m3 0.21mol/mol). First – type boundary conditions were specified above the water table for the oxygen transport model. The spoil solution initially contained no oxygen. To avoid non-linearity problems, no ferric complexation reactions were allowed to take place. An oxidation period of 10000 days (27 years) was considered for the simulation. 100 time steps were considered for the simulation.

The pH of the system was not allowed to drop below 3.0. The interaction between H+and the spoil and also maintaining the pH closer to the optimum pH required for the maximum activity of the iron-oxidising bacteria (3.0), caused more ferric iron to be produced by the bacterially mediated oxidation of ferrous iron. Consequently more pyrite was oxidised.

Figure 9 shows the oxygen concentrations for simulation times of 5 and 10 years. In the segment where chemical reactions take place, more oxygen was consumed due to bacterial activity. As time progresses, more oxygen is being diffused to the reaction site.

Figure 9.

Oxygen concentration after: (a) 5 years; (b) 10 years of simulation

The ferrous iron concentrations for simulation times of 5, 10, 12 and 15 years are illustrated in Figure 10. After 5 years of the simulation (Figure 10a), pH reached about 3.3, a favourable pH required for bacterial activity. The bacteria catalysed the ferrous iron oxidation reaction and more Fe2+converted toFe3+. Therefore, the ferrous iron concentration decreased in the unsaturated zone whereas the ferric iron concentrations increased considerably in this zone.

At 10 years (Figure 10b), Fe2+is almost removed from the unsaturated zone due to ferrous- ferric oxidation reaction and a downward recharge water flow. Therefore, An Fe2+ peak (greater than 1.05mol/m3) reached to the saturated zone. In the saturated zone, Fe2+spreads by groundwater flow. By 12 years (Figure 10c), the Fe2+ peak moved further down into the saturated zone.

After 15 years of the simulation (Figure 10d), pH is still in the range that is favourable for bacterial activity; therefore, ferric concentration is dominant in the unsaturated zone and Fe2+ peak moved further down into the saturated zone.

Figure 11 shows the solution pH for simulation times of 5, 10, 12 and 10 years. As pyrite oxidation takes place the pH dropped significantly (average pH is 3.3 in the unsaturated zone) but it never drops below about 3.2 due to taking into account the reaction between the spoil andH+. By 15 years the average pH is about 3.5 in the unsaturated zone due to the transport of hydrogen ions downward by the water flow.

Figure 10.

Formula: Eqn137.wmf>concentration after: (a) 5 years; (b) 10 years; (c) 12 years; (d) 15 years of simulation

Figure 11.

pH after: (a) 5 years; (b) 10 years; (c) 12 years; (d) 15 years of simulation

Figure 12 shows the ferric iron concentrations after 5, 10, 12 and 15 years of simulation. The concentrations of Fe2+are high in the unsaturated zone where oxygen is available and the conditions for maximum bacterial activity are favourable. The ferric concentrations are limited to the unsaturated zone, but as time progresses, some is being transported below the water table. In the case where pyrite is present below the water table, Fe3+is converted back into Fe3+ through the Fe2+- pyrite reaction.

Figure 12.

Formula: Eqn142.wmf>concentration after: (a) 5 years; (b) 10 years; (c) 12 years; (d) 15 years of simulation

Figure 13 shows Fe3+concentrations for time periods of 5, 10, 12 and 15 years. Both oxygen and Fe3+ react with pyrite and produceSO42.

Figure 13.

Formula: Eqn146.wmf>concentration after: (a) 5 years; (b) 10 years; (c) 12 years; (d) 15 years of simulation

Although not illustrated, early in the simulation the Fe3+peak occurs in the unsaturated zone. Pyrite-oxygen reaction is recognised to be the only important source for SO42 production at the beginning of the simulation.

At 5 years (Figure 13a), an SO42 peak (greater than 58SO42) occurred in the saturated zone due to a downward recharge water flow. In the saturated zone, SO42spreads by groundwater flow. As time progresses (Figure 13b-d), the SO42 peak moved further down into the saturated zone whereas mol/m3 concentrations in the unsaturated zone began to decrease. No significant concentrations of SO42 were produced in the saturated zone by the reaction between ferric iron and pyrite, because, for two-dimensional simulations, it was considered that only pyrite in a 12.5 m wide segment of the unsaturated zone participates in the oxidation reaction.

4.2. Problem 2- modelling of pollutants transportation associated with AMD in Shour River at Sarcheshmeh copper mine, Iran

During the past years there was a major public interest concerning the management of water resources. A hot spot of this field is the availability of water of acceptable quality under circumstances of pollutant discharge in rivers. The water quality management in such situations requires fast decisions based on knowledge related to the distribution of pollutant concentration along the river downstream of the releasing source.

The second problem deals with the modelling of transportation of metallic pollutants associated with AMD in the Shour River at Sarcheshmeh copper mine, Iran. This case attempts to numerically model the water pollution problems in the Shour River passing near the Sarcheshmeh copper mine. This river has been seriously affected by AMD emanated from the Sarcheshmeh mine, sulphide wastes and acidic leachates from heap leaching facilities in the study area.

Sarcheshmeh is a large open pit mine in the Kerman Province of Iran and considered to be one of the largest copper deposits in the world. Sarcheshmeh complex contains substantial amounts of molybdenum, gold and other rare metals. Approximately 40,000 tons of ore at 0.9 % Cu and 0.03 % molybdenum is extracted per day in Sarcheshmeh mine (Banisi & Finch 2001). The Sarcheshmeh copper mine is located about 160 km southwest of Kerman and 50 km south of Rafsanjan. Sarcheshmeh ore body, situated in the central part of Zagros ranges, consists of folded and faulted early tertiary volcano- sedimentary rocks. The mine is recognized with latitude of 29˚ 58΄ and a longitude of 53˚ 55΄. The average annual precipitation at the site varies from 300 to 550 mm. The temperature varies from +35 ˚C in summer to -20 ˚C in winter (Bani Asadi et al., 2008; Doulati Ardejani et al., 2008a). Production units of Sarcheshmeh Copper Complex involve the mine itself, concentrator, smelter, refinery, foundries and heap leaching facilities. Different aspects of the Sarcheshmeh copper deposit have been previously described by Aftabi & Atapour (2000), Atapour & Aftabi (2007), Bani Asadi et al., (2008), Banisi & Finch (2001), Bazin & Hubner (1969), Doulati Ardejani et al., (2005), Doulati Ardejani et al., (2008a), Ghorashizadeh (1978), Majdi et al., (2007), Shayestehfar et al., (2007).

Mining operation has placed many low grade waste dumps and has posed many environmental problems. AMD is a serious problem in Sarcheshmeh copper mine which has a detrimental impact on groundwater and surface water quality. Sulphide minerals oxidation in particular pyrite and subsequence discharge of acidic drainages containing toxic metals into Shour River and groundwater aquifer is the major source of water pollution.

AMD is generally formed due to the oxidation of sulphide minerals in particular pyrite. When pyritic minerals are exposed to air, their rapid oxidation occur which produces acidic drainage (Atkins & Pooley, 1982; Ricca & Schultz, 1979; Rubio & Del Olmo, 1995).

The impacts of AMD on the quality of the surface and groundwater were investigated by sampling in the Shour River and analysing them for hydrogeochemical parameters in particular toxic metals. The concentration of SO42- in both surface and groundwater samples were high. The pH of the water samples in Shour River varies from 2 to 3.9. However, the pH of groundwater samples (6.3-7.2) was found within the permissible limit (Bani Asadi et al., 2008; Doulati Ardejani et al., 2008a).

Numerical modelling of metallic pollutants transport is necessary to predict their migration through the Shour River. Modelling the metallic pollutants fate and transport can help in the design of a mining operation to minimize the various effects on the environment during the activities. The effects of AMD may remain for long times after mine closure, hence careful environmental management is a necessary task during mining operations and a basic mine closure strategy has to be developed to prevent the subsequent generation of AMD. Simulation of long-term pyrite oxidation and the transportation of oxidation products in particular metallic pollutants provide useful information which can be used for the design of efficient remediation programs.

4.2.1. Sampling

Water samples were collected from the Shour River. The Shour River is the main stream flows in the study area. The Sarcheshmeh mine drainages enter the Shour River. Figure 14 shows this river and the sampling locations. The water samples were collected in February 2006. Sampling locations consist of upstream of the Shour River, coming from the Sarcheshmeh mine in the upstream of heap facility, in the entrance of acidic leachate of heap structure to the river, run-off of leaching solution into the river and river downstream (Table 4). The water samples were immediately acidified to less than pH=2 using pure HNO3. Analysis results of samples are shown in Table 5. The pH of the samples was measured using a portable pH meter. Concentrations of dissolved metals in the water samples were determined using an atomic adsorption spectrometer (AA220) in water Lab of the National Iranian Copper Industries Co (Bani Asadi et al., 2008).

Concentrations of heavy metals and SO42- in water decrease with distance from upstream to downstream of the Shour River due to a dilution effect by uncontaminated waters, adsorption process with suspended fine particles in river water and riverbed sediments and also due to precipitation effect resulting from increase in pH. Such processes can be well described by mathematical models for transport and fate of pollutants in rivers (see Ani et al., 2010).

pH of the Shour River water plays an important role in the transportation of heavy metals. In the Shour River, Fe is the most rapidly depleted from aqueous phase. pH varies from 2.0 (leaching solution run-off) to 4.3 (Shour River diluted water). This range of pH is lower than the acceptable values of pH from any surface mine site (Stiefel & Busch, 1983). Due to low pH, the heavy metals remain in dissolved state in river water for long distances. Subsequently, this leads to potential environmental damages in the study area. Furthermore, high concentrations of heavy metals and low pH cause the death of aquatic life.

Figure 14.

Sampling locations in Shour River, from Sarcheshmeh porphyry copper mine

StationY (m)X (m)Sampling place
S13315792339809leaching factory overflow (leaching solution run-off)
S23315734389737river water (downstream of the leaching factory)
S33315095390234acidic leachate of heap No.2
S53315422390013river water mixed with heap leachate
S63315320390013river water (crusher downstream)
S73314578339765river water (crusher upstream)

Table 4.

Sampling details from Shour River. Sample S4 was removed from Table 1 due to analysing error

SampleDistance (m)Cu Mn Fe CdPbCrMoZnSO4pH

Table 5.

Chemical analysis of water samples from Shour River

Variations in pH strongly control the mobility of dissolved heavy metals with a typical inverse relationship (Dinelli et al., 2001). The heavy metals are significantly elevated in water sampled at the nearest points of river to the heap structure. Concentrations of heavy metals in water decrease with distance from the heap facility due to a dilution effect by uncontaminated waters and due to an increase in pH. This shows that high concentrations of metals are greatly influenced by leachates from heap structure. The concentrations of most metals are significantly higher than their standard limits found in potable water (WHO, 2008).

4.2.2. Governing equation of pollutant transport through the Shour River

In this section, a CFD model incorporating the finite volume discretisation scheme is presented to simulate pollutants transportation processes through the Shour River at Sarcheshmeh copper mine. It is assumed that the adsorption is the main process occurred during the course of the physical transportation of the pollutants in river. This process may take place between dissolved pollutants and suspended fine particles in river water and riverbed sediments. The model incorporates the Langmuir isotherm expression to describe the non-linear adsorption process in the Shour River. The governing equations of the model were numerically solved using PHOENICS software. The field data presented in this study were compared with those results predicted by the numerical model.

The Langmuir isotherm can be expressed as follows (Doulati Ardejani et al., 2007; Gokmen & Serpen, 2002; Hachem et al., 2001; Shawabkeh et al., 2004; Sheng & Smith, 1997):


where,SO42 is the adsorbed concentration, C is the equilibrium concentration, Q0 denotes the maximum absorption capacity and SO42is the Langmuir constant.

The second-order partial differential equation that expresses the pollutant transport at the surface rivers and groundwater aquifer and consist of both physical (advection and dispersion) and chemical (adsorption) processes can be rewritten as:



q= Cartesian coordinatesKl;

ϕCt+ρbqt=D2C2xjvjCxj±S= dispersion coefficientxj;

(m)= porosity;

D= vector components of the specific discharge (m/s);

(m2/s)= time (s);

ϕ= bulk density (kg/m3);

vj= sink and source term representing the change in aqueous component concentrations due to the chemical reactions.

The concentration of any chemical species in the aqueous system can be related to its sorbed concentration in the solid phase using the following equation:


Substituting Equation 17 into Equation 16 yields the following equation:



Differentiating the Langmuir isotherm expression with respect to the concentration S gives:


Substituting Equation 19 into Equation 18 yields the following equation:


Equation 20 can now be rewritten as:


The rearrangement of Equation 21 yields:


where, terms in ϕCt+ρbQ0KL(1+KLC)2Ct=D2Cxj2vjCxj±S have been added to both sides of Equation 21.

For a continuous system and assuming that the adsorption is the only mechanism taking place during pollutant transportation, Equation 22 reduces to:


In Equation 23, Ct=D2Cxj2vjCxj+(1ϕ)CtρbQ0KL(1+KLC)2Ct±Sand C/t

Equation 23 contains one additional source term, which can be evaluated by FORTRAN99 codes located in GROUP 13 of GROUND routine in PHOENICS software. Given appropriate boundary conditions and an initial condition, the non-linear partial differential Equation 23 is solved numerically using PHOENICS software for each dissolved pollutant.

In the case of a single-phase problem, the partial differential equation solved by PHOENICS has the following general form (Spalding, 1981):



ρb=1 = any of the dependent variable;

= time;

t(ρψ)+xj(ρujψΓψψxj)=Sψ= PHOENICS-term for density;

ψ= velocity component in the t direction;

ρ= diffusive exchange coefficient foruj;

xj= source rate ofΓψ.

4.2.3. Modelling setting and input data

In order to model the transport process, a one-dimensional modelling was performed using PHOENICS software. The Shour River was simulated as a finite volume model with a length, a width and a depth of 1200, 1 and 1 m respectively. The model was divided into 1500 equal size control volumes and considered as a continuous system in which the pollutants transport take place. An average water velocity of 1 (m/s) was considered. 30 time steps with a power law distribution (power=1.5) were assigned. The modelling time was 1200 seconds. Total iteration of 3500 was assigned to the simulation. A grid distribution with geometric pressure was used with a power/ratio of 1. An upwind differencing scheme was considered (Figure 15).

A hydrodynamic dispersion equal to 1.5 ψ was assigned. Table 6 gives the Langmuir isotherm constants and bulk density of the flow medium used in pollutant transport model.

In this paper, the results predicted by the finite volume model were compared with the field data. Figures 16 to 25 compare the numerical modelling results to the field data for Cu2+, Mn2+, Fe2+, Cd2+, Pb2+, Cr2+, Mo2+, Zn2+ and SO42- concentrations distributions and pH in the Shour River. A close agreement between the results of numerical model and field data explains that the adsorption process is well described by the Langmuir isotherm in the Shour River.

Figure 15.

Finite volume model for Shour River

Metallic pollutantKl(L/mg)Q0(mg/g)ρb(kg/m3)References
Cu2.1038.5981Doulati Ardejani et al., (2011)
Mn0.3472.3811Doulati Ardejani et al., (2011)
Fe1.907.151Not published data
Cd2.053.651Not published data
Pb1.455.631Not published data
Cr0.653.301Not published data
Mo2.06.201Not published data
Zn2.256.151Not published data
SO41.658.751Not published data
H+1.935.751Not published data

Table 6.

Langmuir isotherm constants and bulk density of the flow medium used in the simulation

Figure 16.

Comparison of model predictions (▬) and field data (■) for Cu2+ concentrations versus distance in the Shour River

Figure 17.

Comparison of model predictions (▬) and field data (■) for Mn2+ concentrations versus distance in the Shour River

Figure 18.

Comparison of model predictions (▬) and field data (■) for Fe2+`concentrations versus distance in the Shour River

Figure 19.

Comparison of model predictions (▬) and field data (■) for Cd2+`concentrations versus distance in the Shour River

Figure 20.

Comparison of model predictions (▬) and field data (■) for Pb2+ concentrations versus distance in the Shour River

Figure 21.

Comparison of model predictions (▬) and field data (■) for Cr2+ concentrations versus distance in the Shour River

Figure 22.

Comparison of model predictions (▬) and field data (■) for Mo2+`concentrations versus distance in the Shour River

Figure 23.

Comparison of model predictions (▬) and field data (■) for Zn2+`concentrations versus distance in the Shour River

Figure 24.

Comparison of model predictions (▬) and field data (■) for SO42- concentrations versus distance in the Shour River

Figure 25.

Comparison of model predictions (▬) and field data (■) for pH versus distance in the Shour River

4.2.4. Reactive transport modelling of Fe (II) and pH including ferrous iron oxidation and ferric iron hydrolysis reactions

The design of water management policies and water use technologies as well as management of mine drainages can only be conducted properly when it is based on appropriate mathematical models (Kaden & Luckner, 1984). For more accurate prediction of Fe (II) and pH fate and transport in Shour River, it is necessary to incorporate hydrolysis of Fe (III) reaction. Sψproduced by the oxidation of ferrous iron may be hydrolysed and precipitated as amorphous ferrihydrite according to Equation 25:


Equation 25 is important in the Shour River and it entirely depends on pH. Reduced conceptual water quality models have been developed for the concentrations of [m2/s] and [Fe3+] in mine water treatment plants, rivers and remaining pits (Hummel et al., 1985). The typical chemical reactions for all models in the one-phase system ‘water’ are the oxidation reactions of Fe (II) and the hydrolysis of Fe (III). According to Hummel et al., (1985), the following assumptions are necessary to be considered when such models are used for surface water bodies.

  1. The chemical reactions are considered as non-equilibrium reactions with complete stoichiometric turnover of the initial substances.

  2. Oxygen is enough for oxidation processes in the surface water bodies and the partial pressure is constant (Fe2++14O2+52H2OFe(OH)3+2H+).

  3. The transport processes are one dimensional.

  4. All the ferrous hydroxide formed is restricted within the reaction time; no mathematical modelling is therefore necessary to reflect the sedimentation processes.

The kinetics of the chemical oxidation of Fe2+ have been previously studied and the necessary rate expression was found. Singer & Stumm (1970) proposed the following rate law:



PO2=0.21barrate constant (Fe2+);

rFe=[Fe2+]t=k[OH]2PO2[Fe2+]=kPO2Kw[Fe2+][H+]2Concentration of hydroxyl ions;

k=Concentration of ferrous ions;

mol2L2min2bar1Concentration of hydrogen ions.

Equation 26 can be rewritten as follow:


[H+]=has the unit of inverse timerFe=[Fe2+]t=k1[Fe2+] and it varies betweenk1=kPO2Kw[H+]2=k*[H+]2 and k1 (quoted in Hummel et al., 1985). Substituting Equation 28 into Equation 27 yields the following equation:


Based on the Equation 25, the stoichiometric ratio between protons and ferrous mass formation rate, 1.6×1013, is:


Considering Equation 30, the rate expression for rFe=[Fe2+]t=k*[Fe2+][H+]2 may be expressed as:


Incorporating the kinetics expressions 29 and 31, the following transport equations can then be written to describe reactive transport of KFe=[H+][Fe2+]=3.58×102(molH+grFe2+) of H+ in the Shour River.


The values of H+ and [Fe2+]t=Dx2[Fe2+]x2v[Fe2+]xk[H+]2[Fe2+] are [H+]t=Dx2[H+]x2v[H+]xk*kFe[H+]2[Fe2+]andkFe, respectively.

Equations 32 and 33 can now be solved by defining the appropriate initial and boundary conditions. Figures 26 and 27 compare the results of numerical model with the field data. The hydrolysis reaction of Fe (III) was incorporated.

Figure 26.

Comparison of model predictions (▬) and field data (■) for Fe2+`concentrations versus distance in the Shour River. The hydrolysis reaction of Fe (III) was incorporated

Figure 27.

Comparison of model predictions (▬) and field data (■) for pH versus distance in the Shour River. The hydrolysis reaction of Fe (III) was incorporated

As Figures 26 and 27 show, the results predicted by the finite volume model for k and 0.0385 agreed closely with the field data. The reason is that, the oxidation of Fe (II) and the hydrolysis of Fe (III) reactions removed ferrous iron from the solution phase in the direction of the water flow.


5. Conclusion

In this chapter the general features of the CFD analysis and the fundamentals of the finite volume numerical technique were described. PHOENICS as a CFD package, its features, capabilities and applications were briefly reviewed. Finite volume discretisation of the surface and groundwater flow and transport equations is presented. The basic principles of the finite volume method can be easily understood and it is particularly strong on coarse non-uniform meshes. One important advantage of the finite volume method over other numerical methods is its strong capability in dealing with boundary conditions and specific source terms. This feature of the finite volume method is very important for acid mine drainage cases associated with production (source) and consumption (sink) terms for many chemical components and the specification of different boundary and initial conditions. The applications of CFD analysis in the modelling of the environmental problems have been considered by the simulation of two different problems. A numerical finite volume model using PHOENICS as a general-purpose CFD package has been developed.

The first problem deals with the simulation of long-term pyrite oxidation and pollutant generation within the spoil of an open cut mine. In this first modelling, the following conclusions were made:

  1. Gaseous diffusion process is an important source for supplying oxygen to the reaction site.

  2. Both oxygen and ferric iron produced by bacterially mediated oxidation of ferrous iron participate in the oxidation of pyrite.

  3. It was found that a pyrite oxidation model in the backfilled site of an open cut mine should take into consideration the role of iron oxidising bacteria.

  4. Bacterial activity caused a sharp depletion of oxygen in the vadose zone.

  5. In the absence of bacteria, oxygen is the only important oxidant of pyrite and the oxidation rate is highly dependent on the effective diffusion coefficient.

  6. The oxidation of only a small fraction of pyrite is enough to generate an acid mine drainage load. Considerable pyrite presence creates a long-term source of acid mine drainage.

  7. The lowering of pH in the range between 2.5 to 3.5, results in the bacterial oxidation of pyrite being enhanced. Subsequently the bacterial action produces more2×1013, Fe2+, H+, andFe3+.

  8. The results of two-dimensional simulations clearly indicate that the most of dissolved ferric iron remains above water table while a ferrous iron peak appears below it, in the saturated zone.

  9. Sulphate initially generated by the pyrite-oxygen reaction and subsequent production from the SO42- pyrite reaction due to bacterial action mainly appears below the water table.

The second problem deals with the modelling of transportation of metallic pollutants associated with acid mine drainage in the Shour River at Sarcheshmeh copper mine, Iran. Acid mine drainage production and heavy metals release and transportation have created many environmental problems. The results show that the mechanical dispersion and advection are the main pollutant transport mechanisms through river water. The results further show that the adsorption is the main mechanism for metal removal from river water. Moreover, this process fits well the Langmuir non-linear isotherm. A better agreement was achieved between the field data and model predictions when ferrous iron oxidation and hydrolysis of ferric iron processes were considered.

The numerical models presented here can help to design, optimise and predict the performance of field remediation and polluted waters treatment programs. The results obtained from this study can be used for developing an appropriate environmental management program in order to monitor surface and groundwater pollution problems and to have a better understanding of the pollution transport mechanisms.



The authors would like to thank the Faculty of Mining, Petroleum and Geophysics, Shahrood University and Technology, School of Civil, Mining and Environmental Engineering, University of Wollongong, Australia and National Iranian Copper Industries (N.I.C.I.Co.) for supporting this research.


  1. 1. AniE. C.HutchinsM. G.KraslawskiA.AgachibP. S. 2010 Assessment of pollutant transport and river water quality using mathematical models. Revue Roumaine de Chimie, 55 4 285291
  2. 2. AniE. C.WallisS.KraslawskiA.SerbanAgachi. P. 2009 Development, calibration and evaluation of two mathematical models for pollutant transport in a small river. Environmental Modelling & Software, 24 11391152
  3. 3. AtkinsA. S.PooleyF. D. 1982 The effects of bio-mechanisms on acidic mine drainage in coal mining. International Journal of Mine Water, 1 3144
  4. 4. BalusuS. R. 1993 Design and development of a multi-scrubber dust control system for longwall faces: Experimental and modelling studies. Ph.D. Thesis, University of Wollongong, 188246 .
  5. 5. BaniAssadi. A.DoulatiArdejani. F.KaramiG. H.DahrAzma. B.AtashDehghan. R.AlipourM. 2008 Heavy metal pollution problems in the vicinity of heap leaching 3 of Sarcheshmeh porphyry copper mine, 10th International Mine Water Association Congress, 355358 , Karlovy Vary, Czech Republic, 2-5 June 2008
  6. 6. BanisiS.FinchJ. A. 2001 Testing a floatation column at the Sarcheshmeh copper mine. Mineral Engineering, 14 7 785789
  7. 7. BarovicG.BoochsP. W. 1981 Two- and Three-dimensional mathematical models of contaminant movement in groundwater. Studies in Environmental Science, 17 849855
  8. 8. BignoliG.SabbioniE. 1981 Long-term prediction of the potential impact of heavy metals on groundwater quality as a result of fertilizer use. Studies in Environmental Science, 17 857862
  9. 9. BinningP.CeliaM. A. 1996 A finite volume Eulerian-Lagrangian localized adjoint method for solution of the contaminant transport equations in two-dimensional multiphase flow systems. Water Resources Research, 32 1 103114
  10. 10. CathlesL. M.AppsJ. A. 1975 A model of the dump leaching process that incorporates oxygen balance, heat balance and air convection. Metallurgical Transactions B, 6B 617624
  11. 11. CHAM 2008 The PHOENICS On-Line Information System,
  12. 12. LideD. R. .Ed. 1998 Handbook of chemistry and physics, 78th edition, CRC press LLC, Boca Roton, New York
  13. 13. DinelliE.LucchiniF.FabbriM.CortecciG. 2001 Metal distribution and environmental problems related to sulfide oxidation in the Libiola copper mine area (Ligurian Apennines, Italy). Journal of Geochemical Exploration, 74 141152
  14. 14. DoulatiArdejani.F).2003 Hydrogeological investigation of backfilled Surface of Coal Mine Site, PhD Thesis, University of Wollongong, Australia, 435 p.
  15. 15. DoulatiArdejani. F.BadiiKh.YousefiLimaee. N.MahmoodiN. M.AramiM.ShafaeiS. Z.MirhabibiA. R. 2007 Numerical modelling and laboratory studies on the removal of Direct Red 23 and Direct Red 80 dyes from textile effluents using orange peel, a low-cost adsorbent. Dyes and Pigment, 73 2 178185
  16. 16. DoulatiArdejani. F.KaramiG. H.BaniAssadi. A.AtashDehghan. R. 2008a Hydrogeochemical investigations of the Shour River and groundwater affected by acid mine drainage in Sarcheshmeh porphyry copper mine, 10th International Mine Water Association Congress, 235238 , Karlovy Vary, Czech Republic, 2-5 June 2008
  17. 17. DoulatiArdejani. F.JodieriShokri. B.MoradzadehA.SoleimaniE.AnsariJafari. M. 2008b A combined mathematical-geophysical model for prediction of pyrite oxidation and pollutant leaching associated with a coal washing waste dump. International Journal of Environmental Science and Technology, 5 4 517526
  18. 18. DoulatiArdejani. F.JodieriShokri. B.BagheriM.SoleimaniE. 2010 Investigation of pyrite oxidation and acid mine drainage characterization associated with Razi active coal mine and coal washing dumps in the Azad Shahr-Ramian region, northeast Iran. Environmental Earth Sciences, 61 15471560
  19. 19. DoulatiArdejani. F.MarandiR.SoleimanifarH. 2011 Biological extraction of copper and manganese ions from acid mine drainage using the fungal species Aspergillus niger. Journal of Environmental Science and Technology (Persian version), accepted for publication.
  20. 20. EdwardsJ. S.RenT. X.JozefowiczR. 1995 Using computational fluid dynamics (CFD) to solve mine safety and health problems, Technical Proceedings, Application of Computers and Operations Research in the Minerals Industries, APCOM xxv 1995, 4147 , Brisbane, Australia, 9-14 July 1995
  21. 21. FriedJ. J. 1981 Groundwater pollution mathematical modelling: important or stagnation. Studies in Environmental Science, 17 807822
  22. 22. GokmenV.SerpenA. 2002 Equilibrium and kinetic studies on the adsorption of dark colored compounds from apple juice using adsorbent resin. Journal of Food Engineering, 53 221227
  23. 23. GreenS. R.ClothierB. E. 1994 Simulation water and chemical movement into unsaturated soils. The PHOENICS Journal of Computational Fluid Dynamics & Its Applications, 7 1 7692
  24. 24. HachemC.BocquillonF.ZahraaO.BouchyM. 2001 Decolourization of textile industry wastewater by the photocatalytic degradation process. Dyes and Pigments, 49 117125
  25. 25. HarvardL.ThomasH. P.DavidW. Z. 1999 Fundamentals of Computational Fluid Dynamics, NASA Ames Research Center, 267p.
  26. 26. HummelJ.FischerR.LuscnerL.KadenS. 1985 Submodels of water quality for the analysis of regional water policies in open-pit lignite mining area, Mine Water Proceeding of the second International Congress, 2 673684 , Granada, Spain, 1985
  27. 27. JaynesD. B.RogowskiA. S.PionkeH. B. 1984a Acid mine drainage from reclaimed coal strip mines, 1, Model description. Water Resources Research, 20 2 233242
  28. 28. JaynesD. B.RogowskiA. S.PionkeH. B. 1984b Acid mine drainage from reclaimed coal strip mines, 2, Simulation results of model. Water Resources Research, 20 2 243250
  29. 29. KachiashviliK.GordezianiD.LazarovR.MelikdzhaianD. 2007 Modeling and simulation of pollutants transport in rivers. Applied Mathematical Modelling, 31 13711396
  30. 30. KadenS.LucknerL. 1984 Groundwater management in open-pit lignite mining area, International Symposium, Montreal, Canada, Proc. Vol. I, 6978 , 21-23 May 1984
  31. 31. KerdijkH. N. 1981 Groundwater pollution by heavy metals and pesticides from a Dredge spoil dump. Studies in Environmental Science, 17 278285
  32. 32. LevenspielO. 1972 Chemical reaction engineering, John Wiley, New York
  33. 33. MurthyJ. Y. 2002 Numerical Methods in Heat, Mass, and Momentum Transfer, School of Mechanical Engineering Purdue University, 196p.
  34. 34. NawalanyN. 1981 Observation of aquifer pollution state. Studies in Environmental Science, 17 985990
  35. 35. PatankarS. 1980 Numerical heat transfer and fluid flow, Taylor & Francis, USA
  36. 36. PickensJ. F.LennoxW. C. 1976 Numerical simulation of waste movement in steady groundwater flow systems. Water Resources Research, 12 2 171180
  37. 37. PinderG. F. 1973 A Galerkin-finite element simulation of groundwater contamination on Long Island, New York. Water Resources Research, 9 6 16571669
  38. 38. PuttiM.YenW. W. G.MulderW. A. 1990 A triangular finite volume approach with high-resolution upwind terms for the solution of groundwater transport equations. Water Resources Research, 26 12 28652880
  39. 39. RabbaniM. G.WarnerJ. W. 1994 Shortcomings of existing finite element formulations for subsurface water pollution modelling and its rectification: One-dimensional case. SIAM Journal on Applied Mathematics, 54 3 660673
  40. 40. RiccaV. T.SchultzR. R. 1979 Acid mine drainage modelling of surface mining. Mine Drainage, Proceedings of The First International Mine Drainage Symposium, G.O. Argall, Jr. C.O. Brawner (Eds.), 651670 , Miller Freeman Publications, Inc., U.S.A, 1979
  41. 41. ShawabkehR.Al-HarahshehA.Al-OtoomA. 2004 Copper and zinc sorption by treated oil shale ash. Separation and Purification Technology, 40 251257
  42. 42. ShengD.SmithD. W. 1997 Analytic solutions to the advective contaminant transport equation with non-linear sorption, Research report 158 12.1997, 0-72591-019-4 of civil, surveying and environmental engineering, The University of Newcastle, Australia
  43. 43. SingerP. C.StummW. 1970 Acidic mine drainage: The rate determining step. Science, 167 11211123
  44. 44. SinghR. N.DoulatiArdejani. F. 2004 Finite volume discretisation for solving acid mine drainage problems. Archives of Mining Science, 49 4 531556
  45. 45. SpaldingD. B. 1981 A general purpose computer program for multi-dimensional one- and two-phase flow. Mathematics and Computers in Simulation, 23 3 267276
  46. 46. StiefelR. C.BuschL. L. 1983 Surface water quality monitoring, In: Surface Mining Environmental Monitoring and Reclamation Handbook, L.V.A. Sendlein, H. Yazicigil and C.L. Carlson (Eds.), 189212 , Elsevier Science Publishing Co., Inc., New York
  47. 47. SvenssonU. 1997 Modeling ground- water flow on a regional scale. The PHONICS Journal: Computational fluid Dynamics and Its Applications, 10 4 442450
  48. 48. U.S. Environmental Protection Agency. 1994 Acid mine drainage prediction. Technical document, office of solid waste, special waste branch, EPA530 -R-94-036, NTIS, PB94-201829, 48p.
  49. 49. VersteegH. K.MalalasekeraW. 1995 An introduction to computational fluid dynamics, the finite volume method, Prentice Hall, UK
  50. 50. WHO, 2008 Guidelines for drinking-water quality, 3rd edition, 1 Recommendations, 3rd edition incorporating 1st and 2nd addenda.
  51. 51. WunderlyM. D.BlowesD. W.FrindE. O.PtacekC. J. 1996 Sulfide mineral oxidation and subsequent reactive transport of oxidation products in mine tailings impoundments: A numerical model. Water Resources Research, 32 10 31733187
  52. 52. Zhi-QiangD.Hoon-ShinJ. 2009 Scaling dispersion model for pollutant transport in rivers. Environmental Modelling & Software, 24 627631

Written By

Faramarz Doulati Ardejani, Ernest Baafi, Kumars Seif Panahi, Raghu Nath Singh and Behshad Jodeiri Shokri

Submitted: October 16th, 2010 Published: July 5th, 2011