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

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).

$$\left(\begin{array}{l}\text{Net}\text{flux}\text{of}\text{}\varphi \text{}\text{due}\text{to}\\ \text{convection}\text{and}\text{diffusion}\\ \text{into}\text{the}\text{control}\text{volume}\end{array}\right)\pm \left(\begin{array}{l}\text{Rate}\text{of}\text{production}\text{or}\\ \text{consumption}\text{of}\varphi \text{}\\ \text{inside}\text{the}\text{control}\text{volume}\end{array}\right)=\left(\begin{array}{l}\text{Rate}\text{of}\text{change}\text{of}\text{}\varphi \text{}\\ \text{in}\text{the}\text{control}\text{volume}\\ \text{with}\text{respect}\text{to}\text{time}\end{array}\right)$$E1

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 $$F{e}^{2+},$$$$S{O}_4{}^{-2}$$and $$H^{+}$$(Singer & Stumm, 1970). The stoichiometric oxidation reaction is as follows:

$$Fe{S}_2+\frac{7}{2}{O}_2+H_{2}O\to F{e}^{2+}+2S{O}_4{}^{2-}+2H^{+}$$E2

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

$$F{e}^{2+}+\frac{1}{4}{O}_2+H^{+}\to F{e}^{3+}+\frac{1}{2}{H}_{2}O$$E3

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 $$\text{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 $${}^{\circ }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:

$$Fe{S}_2+14F{e}^{3+}+8H_{2}O\to 15F{e}^{2+}+2S{O}_4{}^{2-}+16H^{+}$$E4

Alternatively $$F{e}^{3+}$$ may be precipitated as amorphous ferric hydroxide (Jaynes et al., 1984a; Walter et al., 1994a, 1994b) according to Reaction 5 as given below:

$$F{e}^{3+}+3H_{2}O\to Fe{\left(OH\right)}_3+3H^{+}$$E5

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):

$$\frac{\partial }{\partial x}\left(K_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)+W=0$$E6

where,

$$x$$and $$y$$ = Cartesian coordinates;

$$K_x$$and $$K_y$$= 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):

$$\frac{\partial X}{\partial t}=\frac{-3X^{{}^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}}{6{\tau }_{D{}_{\left[O_2\right]}}{X}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\left(1-X^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\right)+{\tau }_{C{}_{\left[O_2\right]}}}+\frac{-3X^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}{6{\tau }_{D{}_{\left[F{e}^{3+}\right]}}{X}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\left(1-X^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\right)+{\tau }_{C{}_{\left[F{e}^{3+}\right]}}}$$E7

where,

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

t = time (s);

$${\tau }_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);

$${\tau }_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);

$$\left[O_2\right]$$ and $$\left[F{e}^{3+}\right]$$= oxygen and ferric iron concentrations.

In Equation 7, both oxygen and ferric iron participate in the pyrite oxidation process. $${\tau }_D$$and $${\tau }_C$$ are expressed by Equations 8 and 9 (Levenspiel, 1972):

$${\tau }_D=\frac{{\rho }_{Py}{R}^2}{6b{D}_{e\left[Ox\right]}{C}_{\left[Ox\right]}}$$E8

$${\tau }_C=\frac{{\rho }_{Py}R}{b{K}_{\left[Ox\right]}{\alpha }_{Py}^{rock}\lambda C_{\left[Ox\right]}}$$E9

where,

$${\rho }_{Py}$$ = density of pyrite in the particle ($$mol/m^3$$);

R = particle radius (m);

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

$$D_{e\left[Ox\right]}$$ = effective diffusion coefficient of oxidant in oxidised rim of the particle ($$m^2/s$$);

$$K_{\left[Ox\right]}$$ = first-order surface reaction rate constant (m/s) ;

$${\alpha }_{Py}^{rock}$$ = surface area of pyrite per unit volume of particle ($$m^{-1}$$);

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

$$C_{\left[Ox\right]}$$ = concentration of oxidant in the water surrounding the particle ($$mol/m^3$$).

Enthalpy model

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

$${\rho }_b{C}_P\frac{\partial T}{\partial t}=K_T\left(\frac{{\partial }^{2}T}{\partial x_j{}^2}\right)-{\rho }_W{C}_W{u}_{x_{j}W}\frac{\partial T}{\partial x_j}+\frac{3\left(1-\varphi \right)B^{-1}{\rho }_S{X}^{\frac{2}{3}}}{6{\tau }_{D_{O_2}}{X}^{\frac{1}{3}}\left(1-X^{\frac{1}{3}}\right)+{\tau }_{C_{O_2}}}+\frac{3\left(1-\varphi \right){B^{\prime }}^{-1}{\rho }_S{X}^{\frac{2}{3}}}{6{\tau }_{D_{F{e}^{3+}}}{X}^{\frac{1}{3}}\left(1-X^{\frac{1}{3}}\right)+{\tau }_{C_{F{e}^{3+}}}}$$E10

where,

$$T$$ = temperature ($${}^{\circ }C$$);

$${\rho }_b$$ = bulk density of the spoil ($$Kg/m^3$$);

$${\rho }_S$$ = molar density of pyrite in spoil ($$mol/m^3$$);

$$\varphi$$ = porosity of spoil;

$$K_T$$ = thermal conductivity ($$J/m{}^{\circ }CS$$);

$$C_P$$ = specific heat capacity of the spoil ($$J/Kg{}^{\circ }C$$);

$${\rho }_W$$ = density of water in the spoil ($$Kg/m^3$$);

$$u_{x_{j}W}$$= water velocity$$(m/s)$$;

$$C_W$$ = heat capacity of water ($$J/Kg{}^{\circ }C$$);

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

$$B^{\prime }$$ = moles of pyrite consumed in pyrite-$$F{e}^{3+}$$ 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:

$${\varphi }_a\frac{\partial u}{\partial t}=D_e\left(\frac{{\partial }^{2}u}{\partial x_j^2}\right)-S_{K_{Py-O_2}}-S_{K_{F{e}^{2+}-O_2}}-S_{K_{B-O_2}}$$E11

where,

$${\varphi }_a$$ = air-filled porosity of the spoil$$\left(m^3/m^3\right)$$;

$$u$$ = gaseous concentration of oxygen in the spoil pore spaces$$\left(mol/m^3\right)$$;

$$D_e$$ = effective diffusion coefficient ($$m^2/s$$);

$$S_{K_{Py-O}{}_2},S_{K_{F{e}^{2+}-O}{}_2},S_{K_{B-O}{}_2}$$= volumetric oxygen consumption terms ($$mol/m^{3}s$$) 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):

$$\varphi \frac{\partial C_i}{\partial t}+{\rho }_b\frac{\partial {\overline{C}}_i}{\partial t}=D\frac{{\partial }^2{C}_i}{\partial x_j{}^2}-q_j\frac{\partial C_i}{\partial x_j}+q_{re}{C}_i\pm S$$

$$i=\mathrm{1,2,...,}{n}_c$$E12

where,

$$C_i$$ = solute concentration in aqueous phase$$(mol/m^3)$$;

$${\overline{C}}_i$$ = solute concentration in adsorbed form$$(mol/Kg)$$;

$${\rho }_b$$ = bulk density of the flow medium$$(Kg/m^3)$$;

$$q_{re}$$ = surface recharge$$(m/s)$$;

$$S$$ = sink and source terms representing the changes in aqueous component concentrations due to the chemical reactions$$(mol/m^{3}s)$$;

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

$$D$$ = hydrodynamic dispersion coefficient$$(m^2/s)$$;

$$n_c$$= 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):

$$\frac{H_{}^{+}}{H_{Gen}^{+}}=\left\{1.0-EXP\left(G_A-pH\right)\right\}$$E13

where,

$$H^{+}$$ = Hydrogen ion in solution ($$mol/m^3$$);

$$H_{Gen}^{+}$$ = $$H^{+}$$ generated through chemical reactions ($$mol/m^3$$);

$$G_A$$ = 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:

$$\frac{\partial H_{}^{+}}{\partial t}=-\frac{H_{Gen}^{+}}{\Delta t}\exp \left(G_A-pH\right)$$E14

where,

$$\Delta 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; Value; Unit; Reference; Pyrite fraction; 0.0025; Kg pyrite/ Kg spoil; Jaynes et al., 1984b; Fraction of coarse particles containing pyrite; 75; %; The bulk density of spoil; 1650; Kg/m3; The porosity of the spoil; 0.321; -; Water-filled porosity; 0.225; -; Surface area of pyrite per unit volume of spoil; 80; cm-1; Cathles & Apps, 1975; Radius of particles; 2; cm; Jaynes et al., 1984b; Diffusion coefficient of oxidant in water; 1.0x10-11; m2/s; Molar density of pyrite within particle; 23; mol/m3; First-order rate constant for Fe3+; 4.4x10-5; m/s; First-order rate constant for oxygen; 8.3x10-10; m/s; Recharge value; 0.5; m/yr; constants for chemical oxidation of Fe2+; K1=1.3x10-10 K2=1.7x10-9; mol2/(m3)2(s); s-1; Effective diffusion coefficient for oxygen transport; 1.0x10-7; m2/s; Enthalpy of Reaction (1); -1.45x106; J/mol; Lide, (1998); Enthalpy of Reaction (3); -1.78x104; J/mol; Lide, (1998); Moles of pyrite consumed per heat generated in FeS2-O2 reaction; 6.9x10-7; mol/J; Moles of pyrite consumed per heat generated in FeS2-Fe3+ reaction; 5.6x10-5; mol/J; Thermal conductivity; 2.1; J/m°CS; Specific heat capacity; 2500; J/Kg°C	Value	Unit	Reference
Pyrite fraction	0.0025	Kg pyrite/ Kg spoil	Jaynes et al., 1984b
Fraction of coarse particles containing pyrite	75	%
The bulk density of spoil	1650	Kg/m^3
The porosity of the spoil	0.321	-
Water-filled porosity	0.225	-
Surface area of pyrite per unit volume of spoil	80	cm^-1	Cathles & Apps, 1975
Radius of particles	2	cm	Jaynes et al., 1984b
Diffusion coefficient of oxidant in water	1.0x10^-11	m^2/s
Molar density of pyrite within particle	23	mol/m^3
First-order rate constant for Fe^3+	4.4x10^-5	m/s
First-order rate constant for oxygen	8.3x10^-10	m/s
Recharge value	0.5	m/yr
constants for chemical oxidation of Fe^2+	K1=1.3x10^-10 K2=1.7x10^-9	mol^2/(m^3)^2(s); s^-1
Effective diffusion coefficient for oxygen transport	1.0x10^-7	m^2/s
Enthalpy of Reaction (1)	-1.45x10^6	J/mol	Lide, (1998)
Enthalpy of Reaction (3)	-1.78x10^4	J/mol	Lide, (1998)
Moles of pyrite consumed per heat generated in FeS2-O2 reaction	6.9x10^-7	mol/J
Moles of pyrite consumed per heat generated in FeS2-Fe^3+ reaction	5.6x10^-5	mol/J
Thermal conductivity	2.1	J/m°CS
Specific heat capacity	2500	J/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 data	Background data
Fe2+	0	0.5
Fe3+	0	0
SO42-	50	50
pH	5	5

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.9$$mol/m^3$$). 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 $$15{}^{\circ }C$$ for the enthalpy model. The top surface of the spoil profile was specified as a first type boundary condition equal to$$15{}^{\circ }C$$.

A constant recharge of $$2\times {10}^{-8}m/s$$was 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\times {10}^{-7}{m}^2/s$$and 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 $$15{}^{\circ }C$$ 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 ($$S_{K_{B-O_2}}$$in 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
Fe2+	0.6	0.6
Fe3+	$$2.1\times {10}^{-5}$$	0
SO42-	55	55
pH	4	5

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$$\times$$20 control volumes of size 1.25 m horizontally $$\times$$ 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\times {10}^{-5}$$ m/s and $$1.40\times {10}^{-6}$$ m/s were used for the simulation. Horizontal and vertical dispersion coefficients of $$7.0\times {10}^{-9}$$ $$m^2/s$$ and $$5.0\times {10}^{-9}$$ $$m^2/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 $${}^{\circ }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\times {10}^{-8}$$$$m^2/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/m^3$$ $$\approx$$ 0.21$$mol/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 ($$\approx$$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 ($$\approx$$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 $$F{e}^{2+}$$converted to$$F{e}^{3+}$$. 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), $$F{e}^{2+}$$is almost removed from the unsaturated zone due to ferrous- ferric oxidation reaction and a downward recharge water flow. Therefore, An $$F{e}^{2+}$$ peak (greater than 1.05$$mol/m^3$$) reached to the saturated zone. In the saturated zone, $$F{e}^{2+}$$spreads by groundwater flow. By 12 years (Figure 10c), the $$F{e}^{2+}$$ 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 $$F{e}^{2+}$$ 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 and$$H^{+}$$. 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 $$F{e}^{2+}$$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, $$F{e}^{3+}$$is converted back into $$F{e}^{3+}$$ through the $$F{e}^{2+}$$- 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 $$F{e}^{3+}$$concentrations for time periods of 5, 10, 12 and 15 years. Both oxygen and $$F{e}^{3+}$$ react with pyrite and produce$$S{O}_4{}^{2-}$$.

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 $$F{e}^{3+}$$peak occurs in the unsaturated zone. Pyrite-oxygen reaction is recognised to be the only important source for $$S{O}_4{}^{2-}$$ production at the beginning of the simulation.

At 5 years (Figure 13a), an $$S{O}_4{}^{2-}$$ peak (greater than 58$$S{O}_4{}^{2-}$$) occurred in the saturated zone due to a downward recharge water flow. In the saturated zone, $$S{O}_4{}^{2-}$$spreads by groundwater flow. As time progresses (Figure 13b-d), the $$S{O}_4{}^{2-}$$ peak moved further down into the saturated zone whereas $$mol/m^3$$ concentrations in the unsaturated zone began to decrease. No significant concentrations of $$S{O}_4{}^{2-}$$ 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 SO₄^2- 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 HNO₃. 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 SO₄^2- 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

Station	Y (m)	X (m)	Sampling place
S1	3315792	339809	leaching factory overflow (leaching solution run-off)
S2	3315734	389737	river water (downstream of the leaching factory)
S3	3315095	390234	acidic leachate of heap No.2
S5	3315422	390013	river water mixed with heap leachate
S6	3315320	390013	river water (crusher downstream)
S7	3314578	339765	river water (crusher upstream)

Table 4.

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

Sample	Distance (m)	Cu	Mn	Fe	Cd	Pb	Cr	Mo	Zn	SO4	pH
(ppm)
S1	0	410	1050	4800	3.8	1.8	1.3	0.6	475	3500	2
S2	69.9	215	360	1515	1.7	1	0.7	0.3	185	3400	2.3
S5	357	225	62.5	2.1	0.17	0.22	0.07	0.1	25	3300	3.9
S6	414	32	5.5	0.5	0.02	0.18	0.05	0.01	1.51	1023	4.3
S3	696	570	215	38	0.45	0.4	0.2	0.1	67	3000	3.4
S7	1200	83	4.2	0.12	0.01	0.1	0.05	0.01	1.13	626	4.3

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):

$$S{O}_4{}^{2-}$$E15

where,$$S{O}_4{}^{2-}$$ is the adsorbed concentration, C is the equilibrium concentration, Q₀ denotes the maximum absorption capacity and $$S{O}_4{}^{2-}$$is 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=\frac{Q_0{K}_{l}C}{1+K_{l}C}$$E16

where,

$$q$$= Cartesian coordinates$$K_l$$;

$$\varphi \frac{\partial C}{\partial t}+{\rho }_b\frac{\partial q}{\partial t}=D\frac{{\partial }^{2}C}{{\partial }^2{x}_j}-v_j\frac{\partial C}{\partial x_j}\pm S$$= dispersion coefficient$$x_j{}_{}$$;

$$(m)$$= porosity;

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

$$(m^2/s)$$= time (s);

$$\varphi$$= bulk density (kg/m3);

$$v_j$$= 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:

$$t$$E17

Substituting Equation 17 into Equation 16 yields the following equation:

$${\rho }_b$$E18

where,

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

$$\frac{\partial q}{\partial t}=\frac{\partial q}{\partial C}\frac{\partial C}{\partial t}$$E19

Substituting Equation 19 into Equation 18 yields the following equation:

$$\varphi \frac{\partial C_{}}{\partial t}+{\rho }_b\left(\frac{\partial q}{\partial C_{}}\frac{\partial C}{\partial t}\right)=D\frac{{\partial }^{2}C}{\partial x_j{}^2}-v_j\frac{\partial C_{}}{\partial x_j}\pm S$$E20

Equation 20 can now be rewritten as:

$$C$$E21

The rearrangement of Equation 21 yields:

$$\frac{\partial q}{\partial C}=\frac{Q_0{K}_L}{{\left(1+K_{L}C\right)}^2}$$E22

where, terms in $$\varphi \frac{\partial C}{\partial t}+{\rho }_b\frac{Q_0{K}_L}{{\left(1+K_{L}C\right)}^2}\frac{\partial C}{\partial t}=D\frac{{\partial }^{2}C}{\partial x_j{}^2}-v_j\frac{\partial C}{\partial x_j}\pm 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:

$$\varphi \frac{\partial C}{\partial t}=D\frac{{\partial }^{2}C}{\partial x_j{}^2}-v_j\frac{\partial C}{\partial x_j}-{\rho }_b\frac{Q_0{K}_L}{{\left(1+K_{L}C\right)}^2}\frac{\partial C}{\partial t}\pm S$$E23

In Equation 23, $$\frac{\partial C}{\partial t}=D\frac{{\partial }^{2}C}{\partial x_j{}^2}-v_j\frac{\partial C}{\partial x_j}+\left(1-\varphi \right)\frac{\partial C}{\partial t}-{\rho }_b\frac{Q_0{K}_L}{{\left(1+K_{L}C\right)}^2}\frac{\partial C}{\partial t}\pm S$$and $$\partial C/\partial 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):

$$\frac{\partial C}{\partial t}=D\frac{{\partial }^{2}C}{\partial x_j{}^2}-v_j\frac{\partial C}{\partial x_j}-{\rho }_b\frac{Q_0{K}_L}{{\left(1+K_{L}C\right)}^2}\frac{\partial C}{\partial t}$$E24

where,

$${\rho }_b=1$$ = any of the dependent variable;

$$\varphi =100\%$$E25

= time;

$$\frac{\partial }{\partial t}\left(\rho \psi \right)+\frac{\partial }{\partial x_j}\left(\rho u_j\psi -{\Gamma }_{\psi }\frac{\partial \psi }{\partial x_j}\right)=S_{\psi }$$= PHOENICS-term for density;

$$\psi$$= velocity component in the $$t$$ direction;

$$\rho$$= diffusive exchange coefficient for$$u_j$$;

$$x_j$$= source rate of$${\Gamma }_{\psi }$$.

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 $$\psi$$ 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 Cu²⁺, Mn²⁺, Fe²⁺, Cd²⁺, Pb²⁺, Cr²⁺, Mo²⁺, Zn²⁺ and SO₄^2- 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 pollutant	Kl(L/mg)	Q0(mg/g)	ρb(kg/m3)	References
Cu	2.103	8.598	1	Doulati Ardejani et al., (2011)
Mn	0.347	2.381	1	Doulati Ardejani et al., (2011)
Fe	1.90	7.15	1	Not published data
Cd	2.05	3.65	1	Not published data
Pb	1.45	5.63	1	Not published data
Cr	0.65	3.30	1	Not published data
Mo	2.0	6.20	1	Not published data
Zn	2.25	6.15	1	Not published data
SO4	1.65	8.75	1	Not published data
H+	1.93	5.75	1	Not 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_{\psi }$$produced by the oxidation of ferrous iron may be hydrolysed and precipitated as amorphous ferrihydrite according to Equation 25:

$$\psi$$E26

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 [$$m^2/s$$] and [$$F{e}^{3+}$$] 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 ($$F{e}^{2+}+\frac{1}{4}{O}_2+\frac{5}{2}{H}_{2}O\to Fe{(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 $$F{e}^{2+}$$ have been previously studied and the necessary rate expression was found. Singer & Stumm (1970) proposed the following rate law:

$$H^{+}$$E27

where,

$$P_{O_2}=0.21bar$$rate constant ($$F{e}^{2+}$$);

$$r_{Fe}=\frac{\partial [F{e}^{2+}]}{\partial t}=-k{[O{H}^{-}]}^2{P}_{O_2}[F{e}^{2+}]=\frac{-k{P}_{O_2}{K}_w[F{e}^{2+}]}{{[H^{+}]}^2}$$Concentration of hydroxyl ions;

$$k=$$Concentration of ferrous ions;

$$mo{l}^{-2}{L}^2{\min }^{-2}ba{r}^{-1}$$Concentration of hydrogen ions.

Equation 26 can be rewritten as follow:

$$[OH]=$$E28

$$[F{e}^{2+}]=$$E29

$$[H^{+}]=$$has the unit of inverse time$$r_{Fe}=-\frac{\partial [F{e}^{2+}]}{\partial t}=k_1[F{e}^{2+}]$$ and it varies between$$k_1=\frac{-k{P}_{O_2}{K}_w}{{[H^{+}]}^2}=\frac{k^{*}}{{[H^{+}]}^2}$$ and $$k_1$$ (quoted in Hummel et al., 1985). Substituting Equation 28 into Equation 27 yields the following equation:

$$(1/t)$$E30

Based on the Equation 25, the stoichiometric ratio between protons and ferrous mass formation rate, $$1.6\times {10}^{-13}$$, is:

$$13\times {10}^{-13}$$E31

Considering Equation 30, the rate expression for $$r_{Fe}=-\frac{\partial \left[F{e}^{2+}\right]}{\partial t}=\frac{k^{*}[F{e}^{2+}]}{{[H^{+}]}^2}$$ may be expressed as:

$$K_{Fe}$$E32

Incorporating the kinetics expressions 29 and 31, the following transport equations can then be written to describe reactive transport of $$K_{Fe}=\frac{[H^{+}]}{[F{e}^{2+}]}=3.58\times {10}^{-2}(\frac{mol{H}^{+}}{grF{e}^{2+}})\begin{array}{c}{}^{{}^{}}\\ {}_{{}_{{}_{}}}\end{array}$$ of $$H^{+}$$ in the Shour River.

$$r_{H^{+}}=k_{Fe}\frac{\partial [F{e}^{2+}]}{\partial t}\begin{array}{c}{}^{{}^{}}\\ {}_{}\end{array}$$E33

$$F{e}^{2+}$$E34

The values of $$H^{+}$$ and $$\begin{array}{c}{}^{{}^{{}^{{}^{{}^{{}^{}}}}}}\\ \end{array}\frac{\partial \left[F{e}^{2+}\right]}{\partial t}=D_x\frac{{\partial }^2\left[F{e}^{2+}\right]}{\partial x^2}-v\frac{\partial \left[F{e}^{2+}\right]}{\partial x}-\frac{k}{{[H^{+}]}^2}[F{e}^{2+}]$$ are $$\frac{\partial \left[H^{+}\right]}{\partial t}=D_x\frac{{\partial }^2\left[H^{+}\right]}{\partial x^2}-v\frac{\partial \left[H^{+}\right]}{\partial x}-\frac{k*k_{Fe}}{{[H^{+}]}^2}[F{e}^{2+}]$$and$$k_{Fe}$$, 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 more$$2\times {10}^{-13}$$, $$F{e}^{2+}$$, $$H^{+}$$, and$$F{e}^{3+}$$.

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 $$S{O}_4{}^{2-}$$- 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.

Acknowledgments

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.