Open access

Geochemical Modelling of Water Quality and Solutes Transport from Mining Environments

Written By

Bronwyn Camden-Smith, Raymond H. Johnson, Peter Camden- Smith and Hlanganani Tutu

Submitted: April 24th, 2014 Published: September 9th, 2015

DOI: 10.5772/59234

Chapter metrics overview

2,237 Chapter Downloads

View Full Metrics

1. Introduction

Water quality has become more and more important over the years owing to increased demand for water by the growing population. This has created some strain on the supply side as most water resources have become polluted, largely through anthropogenic activities such as agriculture, mining, manufacturing and domestic and industrial effluent among others. In South Africa, for instance, mining accounts for the largest part of water pollution due, largely, to abandoned mine sites that have, over the years, continued to release contaminants into surface and groundwater unabated.

The major challenge that abandoned mine sites, for example gold and coal mines, pose is that of acid mine drainage (AMD), a phenomenon that has been well documented [1-3]. There are four sub-reactions involved in AMD [1,4], the first reaction being the oxidation of sulphide to sulphate and ferrous ion:

2FeS2+ 7/2O2+H2OFe2++ 4SO42-+ 4H+E1

The resulting ferrous iron undergoes oxidation as it enters receiving surface water systems. This reaction is dependent on pH and oxygen and proceeds slowly under acidic conditions (pH< 3):

4Fe2++ O2+ 4H+4Fe3++ 2H2OE2

This reaction is catalysed by the acidophilic Thiobacillus ferrooxidans bacteria. The slow rate of this reaction makes it the overall rate-determining step for the acid-generating process. The emanating acidity (H+) can under neutralising reactions on contact with minerals such as calcite (CaCO3), resulting in elevated pH and sustained concentrations of Fe2+and SO42-.

Ferric ion (Fe3+) produced in reaction 2 precipitates as ferric hydroxide (Fe(OH)3), releasing more acidity:

4Fe3++ 12H2O4Fe (OH)3+ 12H+E3

On losing a water molecule and undergoing structural arrangement, iron hydroxide forms goethite (FeOOH), a common iron(III) oxide. The fourth reaction involves the oxidation of additional pyrite by ferric ion, Fe3+:

FeS2+ 14Fe3++ 8H2O15Fe2++ 2SO42-+ 16H+E4

This is a self-sustaining reaction that perpetuates the pyrite oxidation process. Essentially, once a little portion of pyrite has been oxidised to Fe3+, the oxidation front begins to expand into the unoxidised zones of the ore.

In the Witwatersrand Basin, South Africa, the gold-containing ores also contain a host of other minerals, including a variety of silicates, pyrite (FeS2), sphalerite (ZnS), uraninite (U3O8), gersdorffite (NiAsS), galena (PbS), arsenopyrite (FeAsS) and monazite ((Ce, La, Pr, Nd, Th, Y)PO4) [5]. The incidence of AMD in this basin results in the leaching of sulphates and elements such as Zn, Ni, As, U, Th and rare earths from these ores [6-9]. These contaminants can undergo a number of reactions that can lead to their transformation, adsorption and desorption along the flow path [10].

This chapter is aimed at demonstrating how these processes can be comprehensively assessed using geochemical modelling, a useful tool that can be used to gain insight into and understanding of geochemical processes. It complements analytical techniques (that are based on laboratory experiments and field data), enabling interpretation of complex situations that conventional methods cannot handle as well as the prediction and anticipation of how geochemical systems would evolve. This has drawn interest in the use of geochemical modelling to predict geochemical processes over the long term, design and optimise remediation strategies, manage environmental impact, conduct risk assessment and identify parameters of importance in geochemical systems [11].

Several geochemical modelling codes have been developed over the years and used to model various scenarios. These include, among others: the Geochemist’s Workbench (University of Illinois-Urbana Champaign), Minteq (US EPA), PHREEQC (US Geological Survey), EQ3/6 (Lawrence Livermore National Laboratory), Mineql, Netpath (US Geological Survey), MEDUSA (KTH, Sweden), Aquachem (Schlumberger Water Services) and WHAM (NERC, UK). For this study, the PHREEQC Interactive code was used.

It is common and good practice in modelling to set goals and conceptualise the modelling scenario when conducting simulations. For this work, the goal or aim has been mentioned above. An overview of the various applications of geochemical modelling, using a conceptual model of a polluted abandoned mine site is presented (Figure 1), explaining the key concepts and assumptions behind the simulations using a compartmentalised approach of the site.

Figure 1.

Conceptual model of an abandoned gold mine site (Camden-Smith and Tutu, 2014)

The figure shows a tailings storage facility (TSF) as a source of contaminants owing to AMD generation when the material is contacted with oxygenated rainwater. The leachates released either flow into groundwater or as surface flow, discharging into an adjacent pond. The pond drains into an adjacent wetland due to overflow and groundwater flow and, subsequently, into a natural stream that discharges into an important drainage system to the south of Johannesburg. Some of the important processes along the flow path of released contaminants include changes in elemental speciation; neutralisation; adsorption onto and desorption from hydrous iron oxides (HFOs); and evaporation of leachate solutions, resulting in precipitation of secondary minerals. These processes have been described by different types of models, e.g. speciation-solubility, forward modelling (including neutralisation and adsorption reactions), reaction path models (including evaporation sequences), reactive transport modelling (involving a 1-dimensional flow through system in which contaminants are flushed) [12-15], and inverse modelling (involving finding reactive minerals given an initial and final solution) [13-15].


2. Materials and methods

Solid material was sampled from the oxidised zone of the TSF (a depth of up to 1 m). The distinct orange colour characterises this material, an indication of oxidation of pyrite which is dark greyish in colour. The material clods were pulverised and sieved to < 2mm prior to further treatment.

Due to the unavailability of rainwater during sampling (conducted during the dry season), an artificial rainwater solution was prepared according to the procedure described by Anderson et al. [16], but with modification to suit composition of local rainwater in the study area [17]. This was used to leach the tailings at an optimised solid:liquid ratio of 1 g:20 mL (the actual leaching was done on 20 g of solid using 400 mL of artificial rainwater) for 24 h. The leachate was then filtered using a 0.45 µm cellulose nitrate filter paper and geochemical parameters (pH, electrical conductivity, redox potential and temperature) and metals analysed using appropriate probes and inductively coupled plasma optical emission spectroscopy (ICP-OES), respectively.

Three water samples were collected from different points in the pond and mixed to make a composite. Geochemical parameters were determined immediately as previously outlined. The composite sample was then filtered using a 0.45 µm cellulose nitrate filter paper and separated into portions for determination of anions and metals by ion chromatography and ICP-OES, respectively.

Hydrous ferric oxides (HFOs) were extracted and determined according to Breit et al. [18]. From the HFO concentrations, the amounts of weak and strong adsorptions sites, Hfo_w and Hfo_s, were determined according to Dzombak and Morel [19].


3. Data treatment and basis for modelling

The above two samples, namely leachate from oxidised tailings and pond water composite will be discussed and utilised in the case studies related to different geochemical models. The analytical data used for modelling are presented in Table 1. For simplicity, only major metals and a selection of minor metals have been presented. The advantages and disadvantages of this will be discussed in the section on speciation modelling.

As mentioned earlier, the PHREEQC geochemical modelling code was used for modelling [13]. This code, developed at the US Geological Survey, is available for free online. It has wide modelling capabilities, with the possibility of using a number of databases, or incorporating new information into its existing databases. The following is worth mentioning regarding databases.

A geochemical model is said to be as good as the quality of the database used to construct it. A database would contain aqueous, mineral and gas phase constituents and reactions. This is usually referred to as the basis. If these constituents, consisting of components and their related species, are not present in the database, they have to be defined in the input file. Their definition includes thermodynamic parameters such as equilibrium constants, entropy and enthalpy values. Thus, the more comprehensive a database, the better the modelling results.

The detailed descriptions related to the different geochemical models are presented for each case study.

Parameter Pond Water Rainwater Leachate
Geochemical parameters
pH 3.57 3.64
pe 6.4 7.4
Temperature (°C) 18.7 19.8
Conductivity (µS cm-1) 1790 700
Chemical Analyses (mg L-1)
Cl- 47 <20
SO4 2- 1245 257
Al 33.1 29.1
Ca 195.0 2.4
Co 2.1 0.6
Cu 0.5 0.5
Fe 2.3 0.0
K 10.1 6.7
Mg 96.6 16.3
Mn 27.3 0.8
Na 43.8 0.5
Ni 3.5 1.2
U 0.2 0.2
Zn 4.2 0.2

Table 1.

Analysed pond water and artificial rainwater leachate of oxidised tailings


4. Results and discussion

This section presents and discusses the results for the different models obtained from simulating the evolution of the above water samples. Further descriptions of particular models have been included together with the construction protocols of those models.

The input files or scripts have not been included, but information on how these are written or generated can be obtained in the PHREEQC manual and other relevant text [20] that can be downloaded from the USGS website (

4.1. Speciation-solubility models

The chemical analysis of a sample will generally only provide total concentrations of the elements in the solution. How the elements are related to one another is not easily determined using analytical techniques. Basically, the aim of a speciation model is to predict whether the metal ions are “free” (associated with water molecules only) or has anions in the solution replaced water molecules in the hydration sphere of the metal ion.

Speciation-solubility models are used to define the concentration and activity of stable species in a system and to determine the saturation states of minerals within a system [14]. Most models assume local equilibria in order to solve the equations involved in determining the species spectrum. This assumption is often valid as the reactions that affect the bulk chemistry of natural water bodies occur relatively fast with respect to the time scale of interest [21]. Chemical equilibrium can be solved in two ways, namely by either minimising the Gibbs’ free energy of a system or by using mass action equations combined with equilibrium constants. The second method is used by the majority of programming software because of a lack of consistent Gibbs’ free energy data [11]. Speciation-solubility models are the simplest geochemical models to define and are the starting models for geochemical reaction models and reactive transport models [14].

For this model, a SOLUTION_SPREAD function was defined in PHREEQC Interactive and the Wateq4f.dat database was used. The output was a long file containing a large amount of information e.g. solution composition following adjustments for charge and electrical balance as well as speciation of elements in solution. Only the results for the latter in pond water will be discussed.

From the input, an output containing the calculated distribution of species in solution was obtained. The results for iron speciation are presented in Table 2.

At a redox potential (pe) value of 6.4 depicted for the solution, iron is expected to exist predominantly as the ferrous valence state. The majority of this iron is in a “free” form, associating through its hydration sphere. The effect of this on activity is quite notable as the activity was found to be significantly lower than molality. The same effect is not observed for the second most abundant species, FeSO4. This neutral species is not expected to have as extensive a hydration sphere as the Fe2+species and as such its molality and activity are approximately equal. The ferric ion is mainly associated with sulphates and hydroxides. The free ferric iron has a much lower activity than molality because of its small hydration sphere and higher charge. There is a large amount of information that can be gained from inspecting the speciation of elements in a solution. A simple speciation-solubility model can prove useful in explaining laboratory experiments. For example, at these pH and pe conditions, a membrane or polymer designed in the hope of capturing ferric ions might be ineffective as only a small percentage of iron exists in the ferric state and secondly, most of the ferric ions exist as a lower charged sulphate species. This could be of particular relevance in the case of uranium capture. In this example, uranium exists predominantly as an uncharged UO2SO4 species. An important consideration is the existence of equilibrium between redox elements. In the initial calculation, the valence states of each redox species are dispersed using the experimental pe. The iron speciation presented here is not necessarily in equilibrium with other redox elements. If the defined solutions are to be used for further modelling, thermodynamic pe as calculated by the software will be used. The redox elements will equilibrate with one another and the dispersal of valence states will be different.

Valence Species Molality Activity
Fe+2 3.04E-05 1.62E-05
FeSO4 1.16E-05 1.17E-05
FeCl+ 2.92E-08 2.48E-08
Fe(HSO4)+ 2.42E-08 2.04E-08
Fe(OH)+ 1.38E-11 1.17E-11
Fe(OH)2 2.10E-19 2.12E-19
Fe(OH)3 - 3.24E-26 2.75E-26
Fe(HS)2 0.00E+00 0.00E+00
Fe(HS)3 - 0.00E+00 0.00E+00
Total Fe(II) Species 4.21E-05
FeSO4 + 1.53E-10 1.30E-10
Fe(OH)2+ 9.47E-11 4.89E-11
Fe(OH)2 + 5.54E-11 4.71E-11
Fe(SO4)2 - 1.51E-11 1.27E-11
Fe+3 1.04E-11 2.98E-12
Fe(HSO4)2+ 1.86E-13 9.46E-14
FeCl+2 1.58E-13 8.16E-14
Fe(OH)3 1.69E-14 1.70E-14
FeCl2+ 5.88E-16 4.97E-16
Fe2(OH)2 4+ 1.27E-18 8.40E-20
FeCl3 5.48E-20 5.52E-20
Fe(OH)4 5.24E-20 4.45E-20
Fe3(OH)4 5+ 1.04E-25 1.50E-27
Total Fe (III) species: 3.29E-10

Table 2.

Species distribution of iron in pond water

The final section of the output lists the saturation indices (SI) of the minerals that could be possibly be associated with the solution. The SI of the mineral is the log of the ratio between the ion activity product (IAP) and the solubility constant (Ksp) [22]. An SI value of zero indicates that the mineral is in equilibrium with the solution while a negative SI value indicates that the system is undersaturated with the mineral. If the mineral was present, then there would be a possibility that the mineral would dissolve. A positive SI value indicates that the system is supersaturated with respect to that mineral and as such, the solution could precipitate that mineral. There are several reasons as to why a mineral can be supersaturated in solution but not precipitate. These include: kinetic constraints (SI values give no indication of the rate of precipitation reactions); calculated SI values are for pure end members of a mineral series and not reflective of the solid solution composition that would precipitate; errors in analysis and unreliable thermodynamic data [14].

A selection of saturation indices for the pond water is presented in Table 3. According to the indices, alunite, goethite and gypsum are close to equilibrium with the pond water. Pyrite is undersaturated as the activity of the sulphide ion is very low as all of the sulphur is speciated as sulphate or bisulphate ion. Gaseous oxygen is undersaturated. For a gas, the SI is equal to the logarithm of the fugacity. Oxygen was not defined as a species in the initial model as it was not analysed. If the solution was in equilibrium with atmospheric oxygen, then the saturation index would be -0.7. In the next section, the solution will be equilibrated with atmospheric gases.

Name in database SI Formula
Al(OH)3(a) -4.76 Al(OH)3
Alunite 0.93 KAl3(SO4)2(OH)6
Anhydrite -0.75 CaSO4
Boehmite -2.57 AlOOH
Epsomite -3 MgSO4:7H2O
Fe(OH)3(a) -5.71 Fe(OH)3
Goethite -0.05 α-FeOOH
Gypsum -0.51 CaSO4:2H2O
Jarosite-K -12.8 KFe3(SO4)2(OH)6
Jarosite-Na -15.93 NaFe3(SO4)2(OH)6
O2(g) -45.33 O2
Pyrite -68.77 FeS2

Table 3.

Selected mineral saturation indices for the pond water sample

Modelled results modelling often reveal more information than analytical results. As indicated earlier, in order to gain a better model of the pond water, atmospheric gases were added as equilibrium phases. Since the pond water samples were obtained from the upper 30 cm of the pond, it is reasonable to assume that they are in equilibrium with oxygen, carbon dioxide and nitrogen. As nitrogen is an inert gas, it was not included as it would not affect the speciation or solubility of the solution. The concentration of gases in the solution was determined by taking into account their atmospheric abundance and using Henry’s Law. In order to model the equilibration of the solution at 1 atm with atmospheric gases, the log(fugacity) or saturation index of oxygen was set to -0.7 and carbon dioxide to -3.5. For the pond water, the model predicted that 0.2940 mmol of O2 and 0.0128 mmol of CO2 were added to the solution in order to equilibrate the solution with the atmosphere. The addition of oxygen caused the pe to increase to 17.58. This is the thermodynamic pe of the solution. The redox species were then defined according to this pe. At higher pe values, ferric species were the major components. If ferrous and ferric concentrations had been quantified, then it would have been possible to decouple ferrous and ferric iron and they would not speciate according to the defined or calculated pe. It is also possible to define the iron concentration as consisting solely of the ferrous species. This could be a reasonable assumption in cases where the pH of the solution is greater than 4 since, at this pH, precipitation of ferric hydroxide would have removed any ferric species from solution. However, it should be noted that for further reactions, the programme will cause the solution to equilibrate and this will result in either the consumption of oxygen (if defined as an element in the initial solution) or the oxidation of ferrous iron in the presence of a constant oxygen supply (if oxygen is defined as an equilibrium phase). The effect on iron speciation of the addition of atmospheric gases on the initial solution is summarised in Figure 2 and Figure 3. At the experimental pe, iron is predominantly a ferrous species as discussed previously. After equilibration with oxygen and carbon dioxide, the iron was found to be predominantly complexed ferric species. The amount of oxygen gas added to the solution to reach equilibrium was 4.54 mg L-1, while the experimental value obtained on site using a dissolved oxygen probe was 4.6 mg L-1, showing that surface water was indeed equilibrated with atmospheric oxygen. The true nature of iron speciation probably lies between the two models. It should be noted that redox disequilibrium is common in natural waters. As the water equilibrates, both within itself and with the atmosphere, the second model becomes more representative of the system.

Figure 2.

Speciation of iron using experimental pe

Figure 3.

Speciation of iron after equilibration of the solution with atmospheric gases

4.2. Forward models

In forward modelling, the final composition of a solution after a reaction or equilibration is calculated [11]. Thermodynamic models utilise equilibrium constraints and are assumed to be valid for homogeneous reactions in which there is sufficient residence time for equilibrium to be established. Heterogeneous reactions such as the dissolution and precipitation of minerals as well as adsorption are often kinetically controlled and models should be adjusted accordingly. Reaction path models are types of forward models that can track reaction progress in small incremental steps. At each step, the species distribution and saturation indices are calculated and equilibrium is maintained by the dissolution or precipitation of defined minerals [11]. These models are used to describe processes such as the addition of minerals to a water system, the mixing of one solution with another (mixing models) and the removal of minerals or components from a system e.g. evaporation sequences [14].

To demonstrate how forward modelling works, the studied pond water was taken through two remediation processes, namely reaction with calcite (CaCO3) and contact with iron oxide surfaces (i.e. adsorption). These two processes were then modelled simultaneously, simulating reactions that are likely to occur in an environment where the two processes co-exist.

Modelling of the addition of calcite to the pond water was undertaken as follows:

  • The pond water was equilibrated with atmospheric gases (oxygen and carbon dioxide)

  • Calcite was added as an equilibrium phase with a saturation index of 0.0 (this was to determine the maximum amount of calcite that will dissolve into the pond water)

  • The step-wise addition of calcite was carried out, with the solution assessed after each step.

  • Minerals that became oversaturated during the addition of calcite were allowed to precipitate.

The final model for the addition of calcite to pond water is presented in Figures 4-10. There is buffering of pH with initial additions of calcite up to just below 250 mg, beyond which the pH increases significantly. The corresponding drop in redox potential (pe) at higher pH values is expected and demonstrated clearly in pe-pH diagrams. The addition of calcite resulted in the precipitation of five minerals, namely diaspore (AlOOH); iron hydroxide (Fe(OH)3); manganite (MnOOH); malachite (Cu2(OH)2CO3) and nickel hydroxide (Ni(OH)2). Diaspore was the first mineral to precipitate. It continued to precipitate until all the aluminium in solution was depleted. A similar path was observed when defining Al(OH)3 as the mineral, although diaspore remains supersaturated during the precipitation. Amorphous iron hydroxide precipitated early on in the profile. The iron concentration also decreased to almost complete depletion. Manganite exhibited a similar profile. As the carbonate content and pH increased in the solution, malachite (Cu2(OH)2CO3) was predicted to precipitate very late in the profile (just before the saturation of calcite was accomplished). Nickel hydroxide also precipitated during this step. The total concentration of the other defined elements remained unaffected. However, the speciation of elements changed with increasing pH. For example, initially uranium existed predominantly as UO2SO4 and UO22+species. After the addition of calcite, the carbonate species UO2(CO3)22-and UO2(CO3)34-were found to predominate.

Figure 4.

Modelled pH profile during addition of calcite to 1 L pond water

Figure 5.

Modelled pe profile during addition of calcite to 1 L pond water

Figure 6.

Modelled aluminium profile during addition of calcite to 1 L pond water

Figure 7.

Modelled iron profile during addition of calcite to 1 L pond water

Figure 8.

Modelled manganese profile during addition of calcite to 1 L pond water

Figure 9.

Modelled copper profile during addition of calcite to 1 L pond water

Figure 10.

Modelled nickel profile during addition of calcite to 1 L pond water

The addition of ferrihydrite (hydrous ferric oxide (HFO)) to the pond water was modelled. Ions in solution adsorb onto the surface of HFO through strong and weak sites [19, 23]. Essentially, strong sites and weak sites depict sites of high and low energy, respectively. Strong complexation of ions to HFO would occur at strong sites while there would be weak complexation at weak sites. An example of adsorption of Zn2+onto HFO is given below:

FeO-+ Zn2+FeOZn+ E5

For every 1 mmol of HFO, a total of 0.2 mmol weak sites and 0.005 mmol strong sites were allocated in the model according to Dzombak and Morel [19] and Appelo et al. [23]. The addition of 1 g of HFO to 1 L of pond water equilibrated with oxygen and carbon dioxide was modelled. The following were input parameters for the surface: 0.0112 mol of weak sites, 5.6 x 10-5 mol of strong sites and a surface area of 600 m2 g-1.

The addition of 1 g of HFO increased the pH of the solution to 5.83, largely due to the removal of H+ ions from solution through adsorption onto HFO. This caused Al(OH)3, diaspore, Fe(OH)3 and manganite to become supersaturated. By setting them as equilibrium phases and allowing them to precipitate, the precipitation of these minerals upon the addition of HFO can be modelled.

The site occupancy of the HFO surface is summarised in Table 4. Approximately, half of the weak sites are occupied by water molecules. Sulphate and hydroxide ions take up the bulk of the remainder of weak sites. Less than 1% are occupied by bivalent oxycations, namely zinc, copper, manganese, nickel. The strong sites are occupied by these oxycations, with the remainder being filled by calcium hydroxide species, hydroxide and uranyl oxide. The sorbed fractions of the metals are summarised in Table 5. This is the percentage reduction in total soluble content of the metal after the addition of 1 g of HFO. Nearly all of the copper and uranium in the pond water has been removed from solution after the addition of HFO. Uranium is sorbed onto the strong sites and copper is split between weak and strong sites. The HFO had no effect on the total concentration of aluminium, chloride, iron, potassium, magnesium and sodium.

Surface complex Number of sites occupied (mmol) Percentage of sites
Weak sites
Hfo_wOH2 + 1.149 51.3
Hfo_wSO4 - 0.589 26.3
Hfo_wOHSO4 2- 0.247 11
Hfo_wOH 0.242 10.8
Hfo_wOZn+ 0.004 0.2
Hfo_wOCu+ 0.004 0.2
Hfo_wOMn+ 0.002 0.1
Hfo_wONi+ 0.002 0.1
Hfo_wO- 0.001 0.1
Strong sites
Hfo_sOZn+ 0.020 35.8
Hfo_sOMn+ 0.011 19.9
Hfo_sONi+ 0.006 11.6
Hfo_sOH2 + 0.006 11.4
Hfo_sOHCa2+ 0.006 10.8
Hfo_sOCu+ 0.004 7.2
Hfo_sOH 0.001 2.4
Hfo_sOUO2 + 0.000 0.9

Table 4.

Surface complexes formed after addition of 1 g HFO to 1 L of tailings pond water

Element Percentage sorbed
Ca 0.1
Cu 98.8
Mn 2.6
Ni 13.5
S 6.4
U 99.9
Zn 37.0

Table 5.

Percentage of total elements adsorbed after addition of 1 g HFO to 1 L of pond water

By allowing minerals that become supersaturated during addition of HFO to precipitate, a different model is observed. The pH only increases to 3.94 (due to the removal of hydroxide ions during precipitation). As such, the speciation of metal ions is different to that discussed above. Less than 0.1% of weak sites are occupied by metal ions. Water molecules occupy 52.5% and sulphate ions occupy 47% of weak sites. Water occupies 97% of the strong sites with uranium, calcium and copper occupying approximately 1% each. The models presented are based on thermodynamic principles only, no kinetics have been included. The relative rates of adsorption onto the HFO surface and of precipitation of minerals will determine the model that better represents the chemistry of the system.

A model combining the addition of calcite with the addition of HFO can be built. These two reactions are important in situations where water remediation strategies are being considered, for example, the use of a permeable reactive barrier. For this combined model, pond water was equilibrated with oxygen and carbon dioxide and reacted with calcite until saturation of calcite was achieved. Minerals that became supersaturated were allowed to precipitate and the resulting solution brought into contact with 1 g HFO. The final predicted water quality for this model is summarised in Table 6. The percentage reduction shown is with respect to the original solution. The pH and pe of the final modelled solution were 7.32 and 13.8, respectively. The addition of calcite increased the calcium concentration in the final solution, hence a negative value. There was a major reduction in the concentrations of aluminium, copper, iron, manganese, nickel, uranium and zinc. Chloride, potassium, magnesium, sodium and sulphate values were not significantly reduced. The strong sites of the HFO had sorbed zinc (41%), calcium (35%), nickel (21%) and copper (1%). Weak sites were dominated by hydroxides and sulphate with 12% occupation for magnesium and smaller amounts of zinc (2%), calcium (1%) and nickel (1%).

Parameter Total concentrations (mg/L) Percentage reduction
Al 0.0 100
Ca 310 -58.9
Cl 47 0.0
Cu 0.0 99.97
Fe 0.0 99.97
K 10.1 0.0
Mg 90.4 6.6
Mn 0.0 100.0
Na 43.9 0.0
Ni 0.3 90.4
SO4 -2 1164 2.5
U 0.0 99.5
Zn 0.1 97.6

Table 6.

Final quality of pond water after addition of 1 g HFO to 1 L initial pond water

Types of reaction path models include: titration, buffering, flush and kinetic reaction path models [11]. The previous example of calcite dissolution is a typical case of a titration or buffering/neutralisation model. Evaporation sequences are some of the most interesting models as they tend to reveal a wide range of possible minerals, thus complementing analytical techniques such as X-ray diffraction (XRD) that may come short in determining some minerals.

An evaporation model of the pond water is presented as Figure 11. The model was prepared assuming equilibration with the atmospheric gases, O2 and CO2. As such, the iron is speciated from the beginning of the calculation as ferric iron. The ferric species is insoluble and was allowed to precipitate at the beginning of the reaction as jarosite (KFe3+3(OH)6(SO4)2). Goethite (α-FeOOH) could also have selected as an initial phase. The model was created by removing water incrementally using the REACTION function and allowing minerals to precipitate by setting them as EQUILIBRIUM_PHASES. The first model allowed for the dissolution of precipitated minerals as the pH changed (Figure 11). This would be relevant for a closed system, for example an evaporation experiment occurring in a beaker (Figure 12). The water would remain in contact with the precipitated minerals and as the pH decreases during evaporation, the minerals become unsaturated in solution (saturation indices decrease to below zero) and the precipitated minerals dissolve. For a more environmentally relevant evaporation profile in an open system, as in the case of the pond water, the subsequent dissolution would not be allowed because the minerals that precipitate on the edge of the pond will no longer be in contact with the remaining water (Figure 12). Therefore, even though the solution becomes unsaturated with respect to the mineral, no resolubilisation occurs because the minerals have left the system. In both cases, alunite (KAl3(SO4)2(OH)6) and jarosite were supersaturated at the start of the evaporation profile. This is not ideal for an evaporation model as we are assuming that the analysed solution is in thermodynamic equilibrium, meaning that a solution in equilibrium would have already precipitated saturated minerals. Minerals can be supersaturated in natural waters as their precipitation could be kinetically hindered. In general, it is preferable that minerals be saturated during evaporation as in the case of gypsum (CaSO4.2H2O) as shown in the profiles in Figures 11 and 12. Notwithstanding, the models are still useful as an indication of the type of minerals that can be precipitated from the water. The minerals presented are pure end members and in nature would precipitate as solid solution minerals, thereby removing some of the transition metals from solution.

Figure 11.

Closed system evaporation model of water (a) Precipitation of gypsum and (b) zoomed in precipitation of jarosite and precipitation and redissolution of alunite.

Figure 12.

Open system evaporation model of water (a) Precipitation of gypsum and (b) zoomed in precipitation of jarosite and precipitation and redissolution of alunite.

4.3. Inverse models

Inverse modelling is also known as mass balance modelling. Given the composition of two water systems along the same flow path and the mineralogy of the rock through which the water flows, inverse modelling can be used to provide a set of possible reactions that occurred to transform the first solution into the second [24]. Mass balance equations are utilised and neither thermodynamic properties nor kinetic restraints are considered [14].

Inverse modelling of the rainwater leachate is presented in Table 7. The initial solution was the blank experiment that was subjected to the same leaching and analytical methods as the leachate but did not contain any solid material. Table 7 summarises the phase transfers from the initial solution into and from the leachate. Phases were defined based on metal fractionation information gained from sequential extractions of the tailings material. It revealed that most of the metals discussed in this chapter exist as readily soluble phases in the oxidised tailings except for uranium which was associated with a reducible phase. The reducible phases are regarded as iron and manganese oxyhydroxides (FeOOH and MnOOH). It is possible that uranium does not exist in the material as a well-defined uranyl hydroxide, but is strongly adsorbed onto the surface of the iron hydroxide. For the sake of simplicity, however, uranium has been allocated as its own phase. The defined phases were chalcanthite (CuSO4.5H2O), epsomite (MgSO4), MnSO4, morenosite (NiSO4.7H2O), β-UO2(OH)2, gypsum, Al2(SO4)3:6H2O, ZnSO4:H2O, Al(OH)3(a), Fe(OH)3(a) and melanterite (FeSO4.7H2O). The phase data for Al2(SO4)3:6H2O was not available in the Wateq4f.dat database and so had to be imported from the llnl.dat database. For simplicity, chloride (Cl-) and nitrate (NO3-) were excluded from the model. Including these species required the definition of more phases, leading to a large number of potential models. No phase including potassium was defined. Model 1 (Table 7) is the simplest model predicted by the program. It contains the minimum number of phases. In the input file, the uncertainty parameter for each solution must be defined. In order to allow for the convergence of the initial solution, a large uncertainty was allocated for the concentrations. This is often required when attempting to apply inverse models to dilute solutions. Model 1 was generated by allowing for the adjustment of calcium, sodium and sulphate values. As such, there were no calcium or sodium containing phases in the model. Model 2 (Table 7) had more phases. Only the sulphate value was adjusted by the program in order to generate the model. The sulphate values in the initial rainwater and leachate solutions were adjusted by 3% and 5%, respectively. As illustrated in this model, the software will use minerals as defined by the user. The iron concentration in both solutions was below the detection limit, but, despite this, melanterite and amorphous iron hydroxide were defined as phases. Model 2 predicts the dissolution of the ferrous sulphate mineral followed by the precipitation of ferric hydroxide. The precipitation of amorphous aluminium hydroxide in Model 1 (and to an extent in Model 2) assists in lowering the pH of the initial solution to equate to that in the final solution. Forward modelling could be undertaken in order to look at the validity of each model. This would involve dissolving the predicted minerals in rainwater to see if this would yield the leachate solution.

Phase Amount (moles) Mass (mg) Action Formula
Model 1
Chalcanthite 7.70 x 10-6 1.9 Dissolve CuSO4.5H2O
Epsomite 6.65 x 10-4 163.9 Dissolve MgSO4.7H2O
MnSO4 1.46 x 10-5 2.2 Dissolve MnSO4
Morenosite 2.05 x 10-5 5.7 Dissolve NiSO4.7H2O
β -UO2(OH)2 8.41 x 10-7 0.3 Dissolve UO2(OH)2
Al2(SO4)3.6H2O 5.75 x 10-4 258.9 Dissolve Al2(SO4)3.6H2O
ZnSO4.H2O 3.06 x 10-6 0.5 Dissolve ZnSO4.H2O
Al(OH)3(a) 7.23 x 10-5 5.6 Precipitate Al(OH)3
Model 2
Chalcanthite 7.70 x 10-6 1.9 Dissolve CuSO4.5H2O
Epsomite 6.65 x 10-4 163.9 Dissolve MgSO4.7H2O
MnSO4 1.46 x 10-5 2.2 Dissolve MnSO4
Morenosite 2.05 x 10-5 5.7 Dissolve NiSO4.7H2O
β-UO2(OH)2 8.41 x 10-7 0.3 Dissolve UO2(OH)2
Gypsum 3.39 x 10-5 5.8 Dissolve CaSO4.2H2O
Al2(SO4)3.6H2O 5.75 x 10-4 258.9 Dissolve Al2(SO4)3.6H2O
ZnSO4.H2O 3.06 x 10-6 0.5 Dissolve ZnSO4.H2O
Fe(OH)3(a) 1.01 x 10-10 1.1 x 10-5 Precipitate Fe(OH)3
Melanterite 1.01 x 10-10 2.8 x 10-5 Dissolve FeSO4.7H2O
Thenardite 3.87 x 10-6 0.5 Dissolve Na2SO4
Al(OH)3(a) 7.23 x 10-5 5.6 Precipitate Al(OH)3

Table 7.

Inverse modelling of rainwater leachate

4.4. Reactive transport models

Modelling the effect of calcite and HFO on pond water chemistry as it leaches through a medium containing these best exemplifies a reactive transport model. For this purpose, a permeable reactive barrier (PRB) containing calcite and HFO was considered (conceptualized in Figure 13), with initial conditions defined as follows:

  • Division of different areas into compartments as follows: zone 1 containing pond water; zone 2 containing 9 cells that had the same amount of calcite and HFO and zone 3 containing 3 cells that have only rainwater in them.

  • Zone 2 above is the PRB barrier and was “filled” with rainwater that was already equilibrated with calcite.

  • Zone 3 did not contain any calcite nor HFO.

The simulation was such that a band of pond water was allowed to flow through the PRB, with the final predicted water quality represented in zone 3. The assumption was that pond water flowed through the zones followed by rainwater.

Figure 13.

Conceptual model and initial conditions of reactive barrier

Simulation of the transport of iron through the zones gave the results presented in Figure 14. The “steps” are time intervals that help to track the changes along the flow path. Iron was found to be transported conservatively, with concentrations remaining just above 4 x 10-5 mol L-1.This is as a result of no precipitation minerals being defined. The speciation of iron changed with change in pH and pe, but the total concentration remained constant.

Figure 14.

Predicted transportation of iron through the barrier

The results for transportation of manganese are presented in Figure 15. Manganese was found to bind to the surface of HFO. Initially, there is a large decrease in manganese concentration with most of it being captured by HFO sites at the front of the barrier (as shown by step 4 in Figure 16). However, as rainwater flowed through after the pond water, it was shown to leach off manganese from the adsorptive sites, resulting in increased concentrations in subsequent solutions (as shown by step 8 in Figure 15). Thereafter, the manganese concentration is spread out in the barrier as it adsorbs onto HFO and is then slowly leached off (step 16 and step 20 in Figure 15). The effects of strong and weak HFO sites on the binding of manganese are shown in Figures 16 and 17, respectively. Weak sites had higher initial manganese (a magnitude higher than for strong sites) adsorbed onto them and this decreased to lower than for strong sites following flushing with rainwater. Strong sites showed a high resistance to desorption of manganese adsorbed onto them as shown by the almost constant adsorbed concentration (Figure 16) while weak sites, on the other hand, showed a low resistance to desorption of manganese by rainwater (Figure 17). The desorption steps showed that the adsorbed concentrations decreased over time, evidence that manganese was lost from solid surface to the water column. This increased manganese content in water is the sustained content observed in Figure 15 as explained before. The observed decrease in manganese concentration on weak sites could be attributed to competition for these sites by calcium and magnesium contained in rainwater (owing to contact with calcite).

Figure 15.

Predicted transport of manganese through barrier

Figure 16.

Manganese adsorption on strong sites

Figure 17.

Manganese adsorption on weak sites


5. Limitations to geochemical modelling

Geochemical modelling is an important tool for predicting reactions. However, there are some limitations to its capability that should be considered. For instance, the quality of the data used for modelling influences the quality of the models obtained. While some manipulations can be done on the data during the modelling process, the best practice requires that analytical data be as accurate and relevant as possible. Another limitation is the availability and suitability of databases. This is encountered when modelling of solutions of high ionic strengths, e.g. brines, has to be conducted. The commonly used Debye-Hückel and Davies expressions have shortcomings in handling such solutions as they are suitable for lower ionic strengths [13], and the more robust Pitzer expression [25] lacks reliable literature of most data, particularly for redox sensitive species. The other challenge is that of a general lack of literature containing kinetic data for a number of important mineral reactions.


6. Conclusion

Geochemical modelling has been shown to be a useful tool in making predictions of contaminant behaviour in the aqueous environment. This has been demonstrated by modelling the same water samples using various models. While it may have some limitations, geochemical modelling remains a robust method that gives a better understanding of contaminants and that can be used for decision making such as in proposing appropriate remediation strategies for contaminated water, predicting potential impacts of contaminants on the environment and for general risk assessment. It can be used to complement analytical techniques, revealing important information where the latter tend to fail.



The authors would like to acknowledge the National Research Foundation of South Africa, the National Academy of Science (through the PEER Program) and the University of the Witwatersrand for the financial support towards this work.


  1. 1. Singer, P.C. and Stumm W. (1970). Acid mine drainage: rate determining step, Science 167, 1121-1123.
  2. 2. Stumm, W. and Morgan, J.J. (1981). Aquatic chemistry, An Introduction Emphasising Chemical Equilibria in Natural Waters, 2nd edition, Wiley, New York.
  3. 3. Nordstrom, D.K., Jenne, E.A. and Ball, J.W. (1979). Redox equilibria of iron in acid mine waters. In: Jenne, E.A. (Ed.): Chemical modeling in aqueous systems. Am. Chem. Soc. Symp. Washington, D.C., Series 93, 51-79.
  4. 4. Blowes, D.W., Ptacek, C.J., Benner, S. G., Waybrant, K.R., Bain, J.G. (1998). Permeable reactive barriers for the treatment of mine tailings drainage water. In: Proceedings of the International Conference and Workshop on Uranium Mining and Hydrogeology, Freiberg, Germany, Vol. 2, 113-119.
  5. 5. Feather, C.E. and Koen, G.M. (1975). The mineralogy of the Witwatersrand Reefs, Minerals Science Engineering, 7, 189-224.
  6. 6. Rosner, T. and van Schalkwyk, A. (2000). The environmental impact of gold mine tailings footprints in the Johannesburg region, South Africa, Bull. Eng. Geol. Env., Vol. 122, 137-148.
  7. 7. Naicker, K., Cukrowska, E. and McCarthy, T. S. (2002). Acid mine drainage arising from gold mining activities in Johannesburg, South Africa and environs, Journal of Environmental Pollution 122, 29-40.
  8. 8. Mphephu, N.F. (2004). Geotechnical environmental evaluation of mining impacts on the Central Rand. PhD Thesis submitted to the University of the Witwatersrand, Johannesburg.
  9. 9. Tutu, H., McCarthy, T.S. and Cukrowska, E.M. (2008). The chemical characteristics of acid mine drainage with particular reference to sources, distribution and remediation: the Witwatersrand Basin, South Africa, as a case study, Applied Geochemistry 23: 3666-3684.
  10. 10. Camden-Smith, B.P.C. and Tutu, H. (2014). Geochemical modelling of the evolution and fate of metal pollutants. Water Science and Technology 69 (5): 1108-14.
  11. 11. Crawford, J. (1999). Geochemical Modelling – A Review of Current Capabilities and Future Directions. Swedish Environmental Protection Agency, Report 262.
  12. 12. Bethke, C.M. (1996). Geochemical reaction modelling – Concepts and applications, Oxford University Press, New York.
  13. 13. Parkhurst, D.L. and Appelo, C.A.J. (1999). User’s guide to PHREEQC (Version2)—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. U.S. Geological Survey Water-Resources Investigations Report 99-4259.
  14. 14. Zhu, C. and Anderson, G. (2002). Environmental Applications of Geochemical Modeling. The Press Syndicate of the University of Cambridge, Cambridge, United Kingdom.
  15. 15. Šráček, O., Černík, M. and Vencelides, Z. (2013). Applications of Geochemical and Reactive Transport Modeling in Hydrogeology, 1st edition, Palacký University, Olomouc, Czech Republic. ISBN 978-80-244-3781-1
  16. 16. Anderson, P., Davidson, C.M., Duncan, A.L., Littlejohn, D., Urea, A.M. and Garden, L.M. (2000). Column leaching and sorption experiments to assess the mobility of potentially toxic elements in industrially contaminated land. J. Environ. Monit. 2: 234-239.
  17. 17. SA Weather Bureau (1998). Climate of South Africa. Weather Bureau Publication WB40, Department of Environmental Affairs and Tourism, Pretoria, 475.
  18. 18. Breit, G.N.; Tuttle, M.L.W.; Cozzarelli, I.M.; Berry, C.J.; Christenson, S.C.; Jaeschke, J.B. (2008). Results of the Chemical and Isotopic Analyses of Sediment and Ground Water from Alluvium of the Canadian River Near a Closed Municipal Landfill, Norman, Oklahoma, Part 2. USGS Open-File Report: 2008-1134
  19. 19. Dzombak, D.A. and Morel, F.M.M. (1990). Surface complexation modeling. Hydrous ferric oxide, John Wiley & Sons, New York.
  20. 20. Parkhurst, D.L. and Appelo, C.A.J. (2013). Description of input and examples for PHREEQC version 3--A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations: U.S. Geological Survey Techniques and Methods, book 6, chap. A43, 497 p. Available at
  21. 21. van der Lee, J. and De Windt, L. (2001). Present state and future directions of modeling of geochemistry in hydrogeological systems. Journal of Contaminant Hydrology, 47 (2-4), 265-2852.
  22. 22. Nordstrom, D. and Munoz, J. (1986). Geochemical Thermodynamics. Palo Alto: Blackwell Scientific Publications.
  23. 23. Appelo, C., van der Weiden, M., Tournassat, C., & Charlet, L. (2002). Surface complexation of ferrous iron and carbonate on ferrihydrite and the mobilisation of arsenic. Environmental Science and Technology, 36 (14), 3096-3103.
  24. 24. Nordstrom, D. (2007). 5.02-Modeling Low-Temperature Geochemical Processes. In H. Holland, & K. Turekian, Treatise on Geochemistry (pp. 1-38). Pargamon: Oxford.
  25. 25. Pitzer K.S. (1979). Theory-Ion interaction approach: in R.M. Pytkowitcz, ed., Activity Coefficients in Electrolyte Solutions, vol. 1, CRC Press, Inc., Boca Raton, Florida, pp. 157-208.

Written By

Bronwyn Camden-Smith, Raymond H. Johnson, Peter Camden- Smith and Hlanganani Tutu

Submitted: April 24th, 2014 Published: September 9th, 2015