Assessing Groundwater Contamination Risk and Detection of Unknown Sources Using a Multi-Component Reactive Transport Model

Haniel Torres

Discipline of Environmental Earth Science, James Cook Univ., Townsville, Australia.

DOI: 10.4236/gep.2020.85009                

Abstract

One of the most serious and important environmental issues related to the mining sector in Central Queensland is the contamination of abandoned mine sites. Representative of this issue is the abandoned Mount Morgan gold mine. The potential dispersal of acid mine drainage (AMD), a product of more than 100 million tons of sulphide-rich waste rock, into the surrounding environment, is the most challenging environmental problem currently facing this abandoned mine site. The abandoned Mount Morgan gold mine has multiple pollutant species that involve complex geochemical processes. The present study simulated the flow and transport processes founded on hydrological and geochemical conditions of the real-life field at the mine site. To assess the groundwater contamination risk and detect unknown pollution sources, few chemical species such as Iron and Sulphur were considered as the contaminants. The flow model was simulated using the computer code MODFLOW, and PHT3D was used for the simulation of advection, dispersion and chemical reactions of constituents dissolved in this groundwater system, and to mimic the reactive chemical transport processes in the polluted groundwater. To improve on results from other studies (Datta et al., 2017; Scotney, 2016; Doyle, 2016), a calibrated model was a main focus for this study. Field concentration measurements were matched with the flow simulation outcomes to calibrate the model. The results obtained showed a great potential to model transport of contaminants in the groundwater system using a real-world situation.

Keywords

Groundwater Contamination, Source Characterization, PHT3D, Chemically Reactive Transport Simulation, Mine Site Contamination

1. Introduction

Detection of groundwater contamination is considered to be of vital importance to manage and protect groundwater from anthropogenic pressures (Soltani et al., 2017). In recent decades, groundwater contamination has become a critical problem for water resources management in many countries (Rahmati et al., 2015). At mining sites, the principal potential environmental polluters are rock waste materials and tailings, these solutes interact with rain and leach to the aquifer. When exposed mining concludes, nature starts to re-establish the basic groundwater and surface water regimes, at which point the mine becomes flooded (Atanacković et al., 2013). The flooding creates pools at lower elevations, leading to degradation of mine water quality (Younger et al., 2004), The design of groundwater models can provide on-site characteristics of the subsurface contaminant source in abandoned mine sites, and help to reduce uncertainties that govern groundwater flows and contaminant transport, as well as the uncertainty about the most likely location and magnitude of the unknown contamination source. The assessment of groundwater contamination risk involves the detection of pollutants’ movements in the subsurface and the characterisation of unknown sources in terms of their magnitude, location and timing, thereby evaluating the groundwater’s vulnerability (Esfahani & Datta, 2016; Chadalavada et al., 2011). The MODFLOW (flow) and PHT3D (reactive transport) simulation codes are widely used to predict spatial and temporal flows, as well as the concentration values in a contaminated aquifer.

The Mount Morgan abandoned mine site, comprising complex and non-linear physical and geochemical groundwater processes, has been the subject of many researches involving the application of developed methodologies for characterisation of contaminant sources (Datta et al., 2017). However, the modelling of reactive chemical species has not been done. Therefore, the aims of the present study were to establish whether the reactive transport simulation model PHT3D, along with the groundwater simulation model MODFLOW, were implemented and validated to achieve reliable prediction of the contaminant transport processes reactive contaminants, such as Iron, Sulphur and Pyrite, in a contaminated aquifer, and whether this model identify pollution sources at the Mount Morgan abandoned mine site.

The first step is to understand the groundwater flow and transport processes (Bear, 1972), considering important parameters such as location, intensity and duration of the contamination activity. The definition of a correct flow and transport model ensures that the correct spatial and temporal distribution of contaminant concentration is well distributed throughout the site. MODFLOW (McDonald & Harbaugh, 1963) is the flow simulation model used in this case. PHT3D is a transport model used to simulate chemical reaction processes in the complex study area.

1.1. Flow Model

According to Zijl et al. (2017), the groundwater flow model is governed by the following two elementary partial-differential equations.

The equations leading groundwater flow through a porous medium are the water balance equation (Equation (1.1)):

$s\frac{\partial h}{\partial t}+\nabla \cdot q=0$ (1.1)

and for a fluid in an aquifer, the flow of groundwater is governed by Darcy’s Laws (Equation (1.2)) of physics and thermodynamics by:

$q=-k\cdot \nabla h$ (1.2)

where s(L−|) is the specific storage, h(L) is the head, q(L T − 1) is the hydraulic conductivity tensor (in short the conductivity), t(T) is time, V(L − 1) is the gradient operator. Vh is the vector with cartesian components (δh/δx, δh/δh, δh/δz), and V(L − 1) is the divergence operator (V∙q is the scalar (δq_x)/δx + (δq_y)/δy + (δq_z)/iz).

The water balance equation describes the physical law of mass conservation (mass cannot be created or destroyed in the groundwater flow field) and the momentum balance equation that describes the physical law of conservation momentum.

Chunmiao and Gordon (2002) provided an equation for steady state conditions applied to the elemental volume as represented in Equation (1.3).

$\frac{\partial \left(\rho V_x\right)}{\partial x}\Delta x+\frac{\partial \left(\rho V_y\right)}{\partial y}\Delta y+\frac{\partial \left(\rho V_z\right)}{\partial z}\Delta z=0$ (1.3)

where ρ stands for the density (kg/m³) and V for the volume (m³). “x”, “y” and “z” are the three Cartesian dimensions.

1.2. Transport Model

Once the solutes are dissolved in the water system, they move in the same direction as the water. According to Schulze-Makuch (2011), the transport of a solute by the bulk movement of groundwater is expressed by the advective transport equation (Equation (1.4)):

$J=v_x\ast C{n}_e$ (1.4)

where J is the mass flux per unit time, V_x is the average linear groundwater velocity in the direction of flow, C refers to the concentration in mass per unit volume of the solution and ne is the effective porosity of the geological medium.

Equation (1.4) assumes only advection is affecting the transport process which is not an accurate description of the processes. The solutes can interact with the aquifer and microbes along the way; mass transport depends on the behaviour of the solutes. Conservative solutes move with the water (in the same direction and speed) with no interactions along the way. In contrast, if the pollutants are non-conservative or reactive, the prediction of the solutes through the groundwater system is much more complex due to interactions along the way (sorption, adsorption, dissolution, precipitation).

After the solutes enter the groundwater system, the solutes will deal with many processes. Dispersion describes the spread of solute from the point of injection ending in a contaminant plume spreading through the system rather than just a line of contaminant. The different flow paths the water particles take in a geological medium cause dispersion. In porous media, there are faster flow paths because the paths have less friction or the water goes through a direct route or into a slower flow path due to a route closer to the grain boundaries, with more friction causing a mechanical mixture and dilution of the solute within the bulk movement of groundwater.

Advection simply describes the transport due to the natural groundwater flow system (Darcy’s Law), where solutes move with the water at an equal rate to the average linear velocity of groundwater. This is controlled by the hydraulic conductivity and hydraulic gradient of the aquifer. In addition, the average linear groundwater velocity differs from the advective velocity of mass in real scenarios. Indeed, most practical problems assume that the dissolved mass and groundwater will move at the same advection rate.

To resolve all these interactions, the movement of a contaminant, k, in groundwater flow systems is given by the following partial differential equation (Equation (1.5)) in space and time (Rushton & Redshaw, 1979):

$\frac{\partial \left(\theta C^k\right)}{\partial t}=\underset{\text{Dispersion-term}}{\underset{︸}{\frac{\partial }{\partial x_i}\left[\theta D_{ij}\frac{\partial C^k}{\partial x_j}\right]}}-\underset{\begin{array}{c}\text{Advection}\\ \text{Term}\end{array}}{\underset{︸}{\frac{\partial }{\partial x_i}\left(\theta v_i{C}^k\right)}}+\underset{\begin{array}{c}\text{Sink}/\text{Sour}\\ \text{Term}\end{array}}{\underset{︸}{q_3{C}_s^k}}+\underset{\text{Reaction-term}}{\underset{︸}{{{\displaystyle \sum }}^{\text{}}{R}_n}}$ (1.5)

where C^k is the dissolved concentration of the species k (mol/L), θ is the porosity of the subsurface environment, t is the time (s), x is the distance along the Cartesian axis, D_ij is the dispersion coefficient (m²/s). v_i represents the linear pore water velocity (m/S), q_s is the volumetric flow rate per unit volume of aquifer (s − 1) representing fluid sources.

1.3. Reactions

According to Barry (1992), many geochemical reactions, involving sorption, biodegradation, oxidation/reduction and precipitation/dissolution, occur when pollutants enter the groundwater system and mix with the surrounding water. At the Mount Morgan mine site, the main minerals comprising the waste rocks and tailings are pyrite [FeS₂], quartz [SiO₂], jarosite [KFe₃(SO₄)₂(OH)₆], pyrrhotite [Fe(1−x)S], and chalcopyrite [CuFeS₂], with minor amounts of kaolinite [Al₂Si₂O₅(OH)₄] and albite [NaAlSi₃O₈] (Edraki et al., 2005).

Pyrite is considered as primary pollutants in the sources mainly because it is highly reactive in acidic and oxidative environment, as in the Mount Morgan mine site (Gasparon et al., 2017). The following equation is meant to be incorporated into the reactive transport of chemical species when developing the transport model.

Pyrite

${\text{FeS}}_2+2{\text{H}}^{\text{+}}+2{\text{e}}^{-}={\text{Fe}}^{+2}+2{\text{HS}}^{-}$

${\log }_k=-18.479$

${\text{delta}}_h=11.300\text{kcal}$

Groundwater flow models are necessary tools to identify pollutant sources in the groundwater system. In the Mount Morgan case, the inverse problem is distinguished due to not having the specification of all parameters as they occur in reality. This is the case when the majority of the present polluted sites will not have enough significant field measurements to comprehend the pollution process. In effect, hydraulic conductivities, boundary conditions and initial conditions are imposed and assumed. Unknown properties of the studied system are determined using inverse models that allow the user to adjust the properties to check the match between the simulated and observed parameters. In other words, the method is performed by running the simulation process backwards in time.

The fundamental principle to obtain the solution lies in the absence of enough field data, such as heads and concentrations. In addition, some parameters such as the hydraulic conductivity are not unique across the field. The type of modelling used is a calibration of inverse modelling due to the estimation of model parameters; this process is done by adjusting parameters (Zijl et al., 2017) and the numerical models are founded on the discretisation of the continuous equations through a grid.

In the modelling process, the estimation of parameters in a groundwater flow model is generally the most crucial step (Anderson & Woessner, 1992), and the identification of the groundwater pollution source is usually inferred from the field concentration measurements at the study area. Parameter estimation can be done through a series of methods, from manual calibration to complex automatic data assimilation algorithms (Zhou et al., 2014). The methodology used in this study correspond to the indirect inversion method that uses an algorithm to adjust the hydraulic conductivities until the simulated heads coincide with the measured heads. Since it is a versatile and adaptive method, any type of data can be used as input, including soft data (Hunt et al., 2007), with the objective to reduce the difference between the observed and simulated values. When defining the hydraulic conductivity values in the field scale, they must include changes in the groundwater flow that can change with either directly measured quantities (flows in or out of wells, pressures or heads measured in wells and observation wells) or indirectly measured quantities (recharge rates, meteorological data) (Zijl et al., 2017).

2. Methodology

2.1. Study Area

The Mount Morgan abandoned mine site is located in Central Queensland, approximately 35 km southwest of Rockhampton and 500 km north of Brisbane, Australia (see Figure 1). This mine was operative for more than a century, from 1882 to 1991. During this period, it produced a total of 247 t of gold, 37 t of silver, 387,000 t of copper and over 100 Mt of sulphide-rich waste rock (Edraki et al., 2005). The mine site includes an open pit consisting of two lakes containing acidic water, with a surface area of around 34.5 ha and a maximal depth of 300 m. The main acid mine drainage (AMD) storage area is the open pit, a product of the oxidation of sulphide minerals, such as pyrite, pyrrhotite, chalcopyrite and sphalerite (Ulrich et al., 2003; Wels et al., 2006), and has a capacity of 11,555 ML (Taylor et al., 2002). Therefore, the open pit is considered to be the main source of the contamination of the groundwater system. The potential of this pit to cause pollution to the groundwater network was recognised by Wels et al. (2006) who defined a seepage of 8.0 L/s generated from the backfilled and flooded open cut pit. The seepage interception system (SIS) collects an estimated 80% of all seepage, while the remaining 20% (or ~3 L/s) enters the Dee River and reaches the underlying aquifer. The distribution of heads makes the water flow from the open pit to the Dee River through the aquifer, as shown in Figure 2. Indeed, the open cut seepage is the biggest issue for the river.

Figure 1. Plan view of historic mine site (map images: map data © 2018 Google Earth Pro.

Figure 2. Conceptual model of groundwater flow at Mount Morgan (Wels et al., 2004).

The argument to be developed is based on many assumptions that give rise to uncertainties. One uncertainty involves the definition of the source of contamination. The actual contaminant source comes from a range of different sources and is very site specific. In several studies, as presented in this document, the open cut pit is defined as the main source of groundwater contamination, yet, there are other sources that are not included in the simulation. In addition, hydrogeological uncertainties can be caused by spatial distribution and variability of hydraulic parameters, especially hydraulic conductivity in the field. For example, it is assumed that the aquifer has a rectangular box geometry with potential monitoring wells, defined as columns of model grid cells in which the hydraulic conductivity is spatially distributed, while it can be assumed that the recharge rates are steady throughout the model and the leak is equally probable at any location with no significant constraints. Therefore, to predict the contaminant transport, initially we assume that the concentration of the source and the hydrological conditions are constants to know to where the contaminant is moving and how quickly it is moving.

A large restriction involves field concentration measurements, in some cases they are not well defined or have not even been tested. This situation becomes more complex with the incorporation of reactive transport of chemical species. However, it is likely field measurements only are used at monitoring sites and interpolated them to each node using interpolation techniques.

To test these assumptions, a large number of iterations with several plausible scenarios, with respect to the values of hydraulic conductivity and other parameters, as well as concentration measurements, are necessary to achieve a suitable calibrated model.

However, the computational viability performance of the methodology represents a typical obstacle when simulations are intended to be developed.

2.2. Groundwater Flow and Contaminant Transport Simulation Model

The proposed methodology is based on a performance evaluation conducted to validate the process. Two steps are implemented: 1) the simulation of a groundwater regime in the study area, and 2) the calibration of flow and transport model.

The Mount Morgan abandoned mine site is subject to simulation. The open pit is known to be considered the contaminant source polluting the groundwater system. The simulation is done by considering a specific period. As a result, concentration measurements of each chemical species are generated at different chosen locations, which are considered observed values in the field due to the impossibility of obtaining these values in real scenarios.

The present study has focused on determining the groundwater flow and transport contaminant of two chemical species: Iron (Fe) and Sulphur (S) across the abandoned mine site for a total of five years. The reason for choosing these two species is because the AMD is one of the biggest environmental problems facing the mining industry globally (Akabzaa et al., 2007), and pyrite oxidation may contribute significant amounts of Iron to the drainage that leads to the contamination of groundwater systems (Equeenuddin et al., 2010). Therefore, defining the contaminant transport of these chemicals and the generation of pyrite from such sources is a major challenge to control AMD production.

2.3. Flow Model

The present study used a developed flow model by Jha and Datta (2012), and recently adapted by Datta et al. (2017), for the same study area. The flow model was previously established in GMS using MODFLOW with a total dimension of 4000 m² and variable depths. The grid approach was used to construct the MODFLOW simulation in GMS. The grid approach involves working directly with the 3D grid to apply sources/sinks, a model parameter, on a cell-by-cell basis (Aquaveo, 2013a). With the grid approach, values are manually assigned to the grid (Aquaveo, 2013c). The grid in the study area consisted of 100 rows and 100 columns, each cell measuring 40 m by 40 m in plain view. The aquifer was simulated in a model with a four layers’ flow regime and conceptualised by making use of available data in the computational grid. The layering of the subsurface was found to be due to the difference in the velocity of groundwater produced by contrast in hydraulic conductivity among various layers (Domenico & Schwartz, 1998). The accuracy and precision of the simulated model depends on the sizing of the grid. An area with high grid resolution means significantly more data describing the model, and its accuracy will also be higher. However, this requires significant computational effort to construct the model and update the solution in each cell. For this reason, the grid approach was defined to give a balance in data accuracy and computational efficiency.

The parameters considered to build the flow model in terms of the aquifer’s physical and hydrogeological characteristics can be found in Table 1.

An accurate simulation of the groundwater system in the abandoned mine site is possible due to the large amount of hydrogeological available data for this area, making it an ideal and solid base for this study. The study area described in Section 2.1 denotes specific geologic and hydrogeologic conditions that allow the flow movement due to advection, dispersion and chemical reactions such as sorption. The main constituent of the lithology in the area is composed of volcanic and intrusive unconsolidated material with very low permeability. However, due to constant mining and excavation in this area, the permeability might be increased in the shallow bedrocks in which the groundwater can flow; this aspect is reflected in the model. According to Jha and Datta (2012), most of the groundwater flow in the site occurs in the permeable waste material form mining activities and in shallow bedrocks.

The area consists of four layers with changing depths, variable hydrogeological parameters and different boundary conditions, creating a complex flow model.

Table 1. Flow Model parameters.

For this reason, the Layer-Property Flow (LFP) Package was used in MODFLOW. The LFP is an internal flow package that computes conductance components and the rate of water movement into and out of storage. This package was used to determine horizontal hydraulic conductivity, vertical and horizontal anisotropy at each layer. It was assumed that the last mentioned parameters remained constant for each layer, except for layer 2, which was divided into eight zones with different horizontal hydraulic conductivity values, as shown in Figure 3. According to Table 2, zone 8 in layer 2 has the highest hydraulic conductivity mainly due to its location and geology properties. The deepest layer was considered to have the lowest hydraulic conductivity due to its proximity to the confined bedrock.

The recharge inputs to the system only occur in the first layer, which was divided into seven zones (Figure 4). The main source for groundwater recharge is known to come from precipitation. The study area is located in a region with an average rainfall of roughly 815.5 mm/year (MLA, 2008); this value needs to be converted in m/day, which is required in MODFLOW. Due to runoff generation, evapotranspiration losses and soil properties, the amount of water infiltrating the soil is reduced after each rainfall event. For this study, seven different zones in layer 1 as shown in Table 3 and Figure 4 were considered. In some areas, 80% of average yearly rainfall is considered input due to the permeability and topography of the area, such as mining waste dumps that are buried and act as a pathway for the flow of water. Physical properties played a main role in defining the recharge zones.

Another source for groundwater recharge is the infiltrating water coming from the pit. As the pit is flooded, it was considered as a constant head boundary matching the water level in the pit. The pit is currently formed by two lakes that are located in layers 2 and 3 respectively. Considering only one contaminant source as the open pit, the contaminant plume is not altered by other sources such as dumps, and it is practical to obtain actual field measurements. As a result, the accuracy of this process is improved.

The Dee River was also modelled as part of the study area. This river interacts with groundwater, receiving water from the aquifer due to the slope of the water table. Specific head is applied to the river with varying head stages from 194.45 m to 180.49 m. The water table of the river was updated to match real conditions.

Previous researches in this area considered the Dee River as a constant head boundary condition with a standing water level (Scotney, 2016).

Within any one study area there are a number of different groundwater recharges but this study only considered precipitation and the open pit (recharge lake).

Figure 3. Hydraulic conductivity zones for this study in Layer 2.

Figure 4. Recharge zones in Layer 1 according to Doyle (2016).

Table 2. Hydraulic conductivity at each zone according to Doyle (2016) in Layer 2.

Table 3. Recharge rates at each zone according to Doyle (2016).

2.4. Model Calibration

The flow model parameters were calibrated under steady-state flow conditions. This study used the model provided by Datta et al. (2017) in which errors oscillated between 3.20 m to 24.38 m. To improve the error range and reduce the difference between simulated and observed hydraulic heads, some alterations have been applied only to the horizontal hydraulic conductivity and the recharge zones in layer 2. Similar to Doyle’s (2016) and Scotney’s (2016) realistic values, which considered 30% of total annual rainfall as a maximum infiltration rate in each zone, the current study applied the highest recharge value of 5.75e−04 to zone 4, representing 30% of total annual rainfall (1.92e−03) (see Figure 4 and Table 3). In addition, hydraulic conductivity values at each zone were also replaced to find the best match between simulated and observed values, and eventually impacted positively in the reduction of the residual head. These values were also obtained from Doyle’s (2016) and Scotney’s (2016) studies, as shown in b and Table 2. An increase or decrease in hydraulic conductivity values before and after each observation well was an essential part in the calibration process.

2.5. Transport Model

Mass transport processes such as dispersion spread the contaminants in groundwater beyond the region they normally would occupy due to advection alone. For this reason, the present study not only considered a model based on advection and dispersion, it also considered chemical reactions. The dispersion occurs in a porous medium because of two main processes, diffusion and mixing. Due to velocity variation or mechanical dispersion, diffusion and mixing are included in different packages when using PHT3D, and the processes are added into the model.

The contaminant species included in this simulation were defined to be iron and sulfur due to their importance in generating AMD. The AMD tends to occur naturally as part of the rock weathering process but is exacerbated by large-scale earth disturbance characteristics of mining (Dowding & Mills, 2000). After being exposed to air and water, oxidation of metal sulfides (pyrite) within the surrounding rock and overburden generate acidity (Ferguson & Morin, 1991).

The processes by which the contaminants move through a porous media are complex. The processes are expressed mathematically in the model, and the field data, included in Table 4, are necessary for application to the differential equations. The concentration for iron and sulphur were taken from groundwater quality observed at Mount Morgan. The water quality of the open cut presents an iron concentration of 284 mg/L, sulphur concentration of 13118 mg/L and pH of 2.67—the average of quarterly monitoring between June 2003 and August 2004 (Wels et al., 2006). These values were converted in mol/L as required in the PHT3D simulation.

The Basic Transport Package within PHT3D enables the user to define input parameters, such as equilibrium species and equilibrium mineral phases, along

Table 4. Transport Model parameters.

with the stress period and time steps. These parameters are required by the entire transport model. For this study, iron (Fe) and sulfur (S) were defined as equilibrium species, and pyrite (FeS₂) as the equilibrium mineral phase. In addition, a total of 5 years was considered, with time steps of 36.5 days. The Basic Transport Package allows the user to define the problem, specify the boundary and initial conditions, determine the time step, prepare mass balance information and print the simulation results (Aquaveo, 2013b).

The Advection Package solves the concentration change using the Third Order TVD scheme (ULTIMATE), and the Dispersion Package was modified using values from Table 4. The longitudinal and transverse dispersivity, as well as diffusion coefficient, TRPT (ratio of longitudinal to transverse dispersivity) and TRVT (ratio of vertical to longitudinal dispersivity) were inputted into the Dispersion Package.

The Source/Sinks Mixing Package can define the contaminant source location in the model domain. For this study, the open pit acts as the only source of contamination in the abandoned mine site and is represented by 109 cells in layer 2 and 353 cells in layer 3 (462 cells in total). Each cell had to be defined as a specific storage point (PHT3D: Point SS) with concentration values of 0.0050855 mol/L (284 mg/L) for iron and 0.013655 mol/L (13,118 mg/L) for sulphur. The infiltration of contaminant to the subsoil can also be established as a constant recharge with a concentration of 0.009 mg/L as used by Doyle (2016) in a previous study in the same area. The purpose of assigning the recharge concentration was to simulate the transport of contaminant emitted from the landfill in the model. This recharge represents leachate from the landfill in the form of filtration and was assigned directly to the recharge polygons in the conceptual model.

The Chemical Reaction Package was included in the simulation. The representations of complex solute-solid relations in transport simulation using sorption processes were included in the model to simulate the contaminant transport in PHT3D. Sorption retards the movement of the plume, and decay (due to biodegradation) reduces the concentration. When dealing with these processes in a complex system like in the Mount Morgan abandoned site, the surface complexation condition should be defined in the model, where the charge of the surface is variable and dependent on the amount and kind of ions sorbed. PHT3D simulates sorption through surface complexation.

The input parameters in the Chemical Reaction Package were defined based on the most commonly used relations for describing equilibrium-controlled reversible sorption. Linear isotherm and First order irreversible kinetic reaction were chosen. Bulk density was defined as 1.6 mg/m³ based on a study conducted by Johnson and Skousen (1995), in which physical properties of different abandoned mine land sites were defined with similar characteristics as Mount Morgan.

3. Results and Discussion

The model showed the hydraulic heads at each grid. The groundwater flow direction was from west to east towards the river. Aligned with Datta et al. (2017) and Petit (2016), the contaminants were assumed to travel with the flow by advection and dispersion through the aquifer from the open pit to the Dee River. As the river is located next to the Mount Morgan abandoned mine site, the groundwater is contributing to the river at approximately 3 L/s (Wels et al., 2009).

The model is a reliable representative of the actual field condition. During the calibration process, computed water levels were matched adequately with observed values. Hydraulic conductivities and recharge values were adjusted in sequential model runs to obtain an almost perfect match between simulated and measured heads. The flow model after the calibration showed error values from 0.1851 m to 5.002 m. The error was determined by the difference between the simulated and observed head values, also called the residual head, shown as colour bars next to each observation well (see Figure 5). The colour of the bars indicates the error, where a green bar signifies the calibration is within the target value, yellow indicates a greater than 100% error, and red indicates an error greater than 200% (Doyle, 2016).

The final calibrated model with an interval of 10 m and a confidence level of 85-90% at each observation well is represented in Figure 5. This model achieved a proper calibration target for the following data: hydraulic head data across the area, groundwater-flow direction, hydraulic-head gradient and contaminant concentration. The difference between simulated and actual field conditions was <5% of the variability in the field data across the model domain. Groundwater was assumed to be under steady state conditions.

The present calibrated model is an improvement on various studies conducted in the same area (Datta et al., 2017; Doyle, 2016; Scotney, 2016) with a maximum reduction in residual head values. Table 5 shows location and residual heads at

Figure 5. Calibrated flow model.

Table 5. Comparison of residual heads obtained by different studies.

each observation well after the calibration for the present study. Negative numbers represent a simulated head above the observed head, while positive values describe the opposite condition. After these adjustments to the parameters, the model was accepted as being adequately calibrated. Calibration is within the target value.

Table 5 and Figure 6 show the difference between different studies in terms of residual heads at each observation well in comparison with other studies.

The unique contaminant source was assumed to be the open pit located at the centre of the model domain. The contaminant species, determined as Iron and Sulphur, were only defined for layer 2 of the model, as represented by time step 36.5 days in Figure 7. Even though the contaminant is present in multiple layers, this study only considered layer 2 as within the largest impacted area in the aquifer. PHT3D displays the concentration values (see the legend to the left of the Figure 7), with concentration values represented by different colours ranging from pink (0 mg/L) to red (maximum concentration value). As the contaminant spreads from the open pit, the colour changes in the border of the pit, showing a decrease in contaminant concentration due to leaching and filtration of contaminant through the soil profile. The balance of the model domain has a concentration value of 0 mg/L.

The reactive transport model developed in this case study was based on several important assumptions: 1) the open pit is the unique contaminant source in the model domain; 2) the two chemical species—iron and sulfur—are assumed to attenuate in concentration through time due to physical, chemical and biogeochemical processes; 3) the reactive solute transport is controlled by sorption interaction resulting in delayed progress of solute; 4) the concentration of the species in solution and the concentration adsorbed to the porous medium are in equilibrium.

Figure 6. Residual heads at different observation well after calibration for different studies.

Figure 7. Transport model for layers 2 and 3 using iron at first time step of 36.5 days.

Only for test purposes, a plume was developed for each species in two scenarios: with and without sorption. When running with sorption, a chemical reaction package and an equilibrium mineral phase have to be defined to account for the effects on reactive solute transport. In this case, pyrite was considered as a product of the reaction between iron and Sulphur (see Figure 8). In contrast, when running the simulation without sorption, no mineral phase is needed, and considering only advection and dispersion processes, two separated plumes for each species are formed, typical of non-reactive solute transport (see Figure 9). Most studies have focused on this type of simulation with non-reactive solute transport; however, this situation does not occur in real world.

The simulation of sorption processes represents a complex solute-solid relation with a reduction of concentration of chemical species. Some solutes, to one degree or another, are removed from the solution and immobilised in the solid matrix, which in this study was pyrite: this process tends to be more common in saturated porous medium. Consequently, simulated concentrations using chemical reactions and sorption for iron and sulphur decrease significantly compared to the simulation without sorption processes. The plume formed by each species is analysed more clearly using Figure 9.

In both cases, with and without sorption, the plume mainly travels in the same direction as the flow, very little dispersion was observed going against the flow.

1) When considering sorption, the exchange between the aqueous and solid phase becomes more pronounced. In other words, high amount of solutes are sorbed onto the solid matrix and left behind as the plume moves away from the source. The plume formed by iron presented concentration values ranging from 1.42 × 10⁻⁶ to 4.28 × 10⁻⁶ mol/L (0.079 to 2.39 mg/L respectively) as shown in Figure 8. In the case of sulphur, concentration values in the plume ranges from 1.18 × 10⁻⁶ to 2.37 × 10⁻⁶ mol/L (0.0378 to 0.0759 mg/L respectively). As the contaminants are not treated as inert compounds, rather they are reactive species, the expected predominant product, in this case pyrite, had a smaller plume than iron and Sulphur.

Figure 8. Plumes for iron (a), sulphur (b) and pyrite (c) at the beginning of the simulation and after five years of simulation with sorption of pyrite.

Figure 9. Plumes for iron (a) and sulphur (b) at the beginning of the simulation and after five years of simulation without sorption of pyrite.

2) When considering without sorption, the plume formed by Iron and Sulphur presented larger concentration than with sorption processes (8.6 × 10⁻⁴ to 5.2 × 10⁻³ mol/L for iron and 1.8 × 10⁻³ to 1.3 × 10⁻² mol/L for Sulphur) (see Figure 9).

The plume shape formed depends mainly on the hydraulic conductivity established in the model domain. To better distinguish the role of hydraulic conductivity in defining the plume shape, Figure 3 and Figure 8 were combined as shown in Figure 10. The plume appears to flow through the aquifer following the high hydraulic conductivity zones. In the case of iron, it is clear that the plume was formed more sharply in zone 9, which has a hydraulic conductivity of 3.456 m/day (see Figure 3), avoiding zones 6 and 7 with lower hydraulic conductivities of 0.001 and 0.216 m/day respectively. A fast flow rate with high dispersion potential was also noticed for Sulphur through zone 9 instead of going through zone 7, as shown in Figure 10.

After 5 years, the contaminant plumes of Iron and Suphur have travelled through different paths in the aquifer, as shown in Figure 10. Therefore, chemical reactions were possible to accurately simulate in this model.

Validation of the Contaminant Transport

The reactive transport was performed to investigate the propagation of iron and sulphur fronts in porous media in the presence of pyrite. Figure 11 shows the measurements in the field of iron concentrations for the current modelling that were obtained from four observation wells close to the open pit by Wels et al. (2006) as part of the Groundwater Quality Program at Mount Morgan. The average of Wels et al.’s (2006) quarterly monitoring between June 2003 and August 2004 is shown in Table 6. The current simulations’ measurements showed the prediction of the plume that is progressively impacting the aquifer. Since the concentration of iron reaching the first monitoring point (MON 01) is lower than the concentration in the open pit and successively lower down gradient in

Figure 10. Contaminant transport for iron (a) and sulphur (b) after five years of simulation showing hydraulic conductivity zones.

Figure 11. Field monitoring points for iron close to the open pit.

Table 6. Concentration values measured in the field compared with simulated concentration values.

MON 03 and MON 04, the contaminant transport is assumed to be accurate. However, Wels et al’s (2006) field measurements in MON 02 were very high in concentration values compared to the first, third and fourth monitoring points. Consequently, a new contaminant source was assumed to be located close to MON 02, the Mundic tailing. According to Wels et al. (2009), one of the Mundic tailings were placed into the historic drainage channel of Mundic Creek (between the open cut and Frog Hollow), with anecdotal evidence suggesting that this tailing was initially deposited in the Mundic drainage without proper containment. In addition, EWL Sciences (2001), through geochemical testing data, concluded that this tailing was highly reactive and can release significant amount of iron.

The comparison between the measured iron concentrations and the outcomes of the reactive transport simulations or selected spatial monitoring points across the study area is shown in Figure 12. The agreement between observed and simulated iron concentrations was very good, particularly five years after beginning the simulation of iron. The transport model was found to be able to reproduce both the shape of the profile and the magnitude of the iron concentrations measured within the heterogeneous porous media. A significant discrepancy

Figure 12. Simulated concentration (Fe) values at different time steps and field measurements in monitoring points close to the open pit.

between the model and the measured field data was seen mainly in the MON 02 which show a field measured iron concentration ten times higher than the simulated value, at 0.035 mol/L and 0.0030 mol/L respectively (see Table 6). This is probably due to the assumption of the existence of a new contaminant source very close to MON 02, the Mundic tailing. For this reason, MON 02 is considering as an outlier for the validation of the contaminant transport model.

Once the flow model was calibrated and the transport model set with the open pit as the source of contaminant with initial concentration for iron and sulphur, monitoring points were chosen. Six observation points were located within a proximity of the source due to the fact that the contaminants’ transport through the aquifer is very slow and concentration values are very small for each specie. Each observation well was systematically selected in order to obtain a change in concentration through time, in other words, the contamination’s location depended on the variation in colour in each cell in the model domain and was expressed in cell IJK identification as shown in Table 7.

According to Figure 13 and Figure 14, Iron concentration has sharply increased during the first 200 days of initialized simulation in the closest point to the source (Obs 1). At this point, the prediction of the concentration of iron fluctuated, showing a decreasing trend after five years with low values at the end of the study time. A decrease in contaminant concentration corroborates the hypothesis that the path taken by the contaminant changes direction. Iron created a plume towards Obs 2 and Obs 3 due to an increase in concentration values at these two points, similar to the Obs 1 at the end of the simulation. The plume has expanded along this way. The plume was least likely to go through the Obs 4 point, as reflected by very low values in concentration after the simulation. In contrast, Obs 5 and Obs 6 are the most likely locations where Iron could make its path.

When referring to sulphur, Obs 1 shows the same pattern as with Iron. Figure 15 shows how the concentration in Obs 1 increased during the first 250 days,

Table 7. Concentration values measured in the field compared with simulated concentration values.

Figure 13. Observation points situated near the contaminant source. (a) Left figure refers to iron plume formation after five years with sorption and (b) Right figure refers to Sulphur plume formation after five years with sorption.

Figure 14. Prediction of iron concentration at each observation point.

to end with a decrease over time. Sulphur tended towards Obs 1 at the beginning of the simulation. However, it changed direction, following Obs 4 which presented high concentration values. Obs 4 was a potential point to contain very high concentrations of sulphur in five years. The concentration values obtained in this study are very small compared to the concentration in the source (0.0050855 mol/L for Iron and 0.013655 mol/L for sulphur), which supports chemical reaction with sorption.

According to Figure 16, the present of pyrite in the aquifer is only reported for Obs 1 which location was located closest to the contaminant source. These

Figure 15. Prediction of sulphur concentration at each observation point.

Figure 16. Prediction of pyrite concentration at each observation point.

detections included values ranging from approximately 1.4 × 10⁻⁸ mol/L to 1.32 × 10⁻⁶ mol/L. Pyrite appeared occasionally at specific time during the first 1460 days of simulation (at the 401.5, 1131.5, 1277.5, 1314 and 1350 days only). After which pyrite concentration values showed a variation with an increasing trend; pyrite started to increase at around 1460 days after the beginning of the simulation.

4. Conclusion

In conclusion, this modelling of contaminants’ transport in the groundwater system using a real-world simulation has been improved compared to earlier attempts to calibrate a simulation model for the complex flow and transport process in an abandoned mine site. The model was calibrated with very good accuracy when recharge rates and hydraulic conductivities were modified and adjusted from past modeling efforts. The changes made did not affect the formation of the plume in the simulation after five years for any of the chemical species but successfully decrease the residual heads at each observation points. Therefore, the flow regime seems to have been calibrated more accurately in this study.

The main contaminant source specified in this study was the open Pit, and it was specified that the contaminants propagate with time, from northwest to southeast following the flow direction as well as the longitudinal and transverse directions. Compared to an advection-dispersion model simulation only, the effect of the same simulation combined with chemical reactions using PHT3D is more realistic. In addition, by adding a mineral phase as a product of iron and sulphur, in this study, pyrite, the concentration of iron and sulphur decreased significantly due to sorption processes. The results of this study have shown a chemical sorption of iron and sulphur in the presence of water, and the formation of a minor amount of pyrite in areas near the source. High amounts of iron and sulphur resulted in their transformation into pyrite, and a very small concentration of these chemical species; therefore; migrated from the pit in this simulation.

Based on this work, more rigorous evaluations can be carried out using PHT3D as a simulation model which incorporates iron and sulphur and other reactive chemical species in complex contaminated aquifer sites, e.g., in abandoned mine sites.

Acknowledgements

The author thanks Dr. Bithin Datta for kindly being my Principal Supervisor of this Postgraduate level minor research project, for his guidance, knowledge and ongoing support throughout the process. I would like to extend my thanks to Neupane Rabindra, a PhD student, for his help with the simulation model and the countless issues that were experienced.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Akabzaa, T., Armah, T., & Banoeng-Yakubo, B. (2007). Prediction of Acid Mine Drainage Generation Potential in Selected Mines in the Ashanti Metallogenic Belt Using Static Geochemical Methods. Environmental Geology, 52, 957-964. https://doi.org/10.1007/s00254-006-0537-1
[2]	Anderson, M. P., & Woessner, W. W. (1992). Applied Groundwater Modeling: Simulation of Flow and Advective Transport. San Diego, CA: Academic Press.
[3]	Atanacković, N., Dragisic, V., Stojkovic, J., Papic, P., & Živanovic, V. (2013). Hydrochemical Characteristics of Mine Waters from Abandoned Mining Sites in Serbia and Their Impact on Surface Water Quality. Environmental Science and Pollution Research, 20, 7615-7626. https://doi.org/10.1007/s11356-013-1959-4
[4]	Aquaveo (2013a). GMS 9.1 Tutorial, PHT3D—Ion Exchange and Surface Complexation. http://gmstutorials-10.4.aquaveo.com/PHT3D-IonExchangeAndSurfaceComplexation.pdf
[5]	Aquaveo (2013b). GMS 9.1 Tutorial, MODFLOW—Conceptual Model Approach I. http://gmstutorials-10.4.aquaveo.com/MODFLOW-ConceptualModelApproach1.pdf
[6]	Aquaveo (2013c). GMS 10.1 Tutorial, MT3DMS—Grid Approach. http://gmstutorials-10.4.aquaveo.com/MODFLOW-GridApproach.pdf
[7]	Barry, D. A. (1992). Modelling Contaminant Transport in the Subsurface: Theory and Computer Programs. In H. Ghadiri, & C. W. Rose (Eds.), Modelling Chemical Transport in Soil: Natural and Applied Contaminants (pp. 105-144). Boca Raton, FL: Lewis, Publishers.
[8]	Bear, J. (1972). Dynamics of Fluids in Porous Media. New York: American Elsevier Pub. Co.
[9]	Chadalavada, S., Datta, B., & Naidu, R. (2011). Optimisation Approach for Pollution Source Identification in Groundwater: An Overview. International Journal of Environment and Waste Management, 8, 40-61. https://doi.org/10.1504/IJEWM.2011.040964
[10]	Chunmiao, Z., & Gordon, D. (2002). Applied Contaminant Transport Modeling (2nd ed.).
[11]	Datta, B., Petit, C., Palliser, M., Esfahani, H., & Prakash, O. (2017). Linking a Simulated Annealing Based Optimization Model with PHT3D Simulation Model for Chemically Reactive Transport Processes to Optimally Characterize Unknown Contaminant Sources in a Former Mine Site in Australia. Journal or Water Resource and Protection, 9, 432-454. http://www.scirp.org/journal/jwarp https://doi.org/10.4236/jwarp.2017.95028
[12]	Domenico, P. A., & Schwartz, F. W. (1998). Physical and Chemical Hydrogeology. New York: Wiley.
[13]	Doyle, J. (2016). Optimal Identification of Unknown Sources of Groundwater Contamination in Abandoned Mine Sites. Undergrad Thesis Submitted to the College of Science & Engineering, James Cook University.
[14]	Dowding, B., & Mills, C. (2000). Natural Acid Rock Drainage and Its Impact upon Background Metal Concentration. http://technology.infomine.com/enviromine/ard/Introduction/Natural.htm
[15]	Edraki, M., Golding, S. D., Baublys, K. A., & Lawrence, M. G. (2005). Hydrochemistry, Mineralogy and Sulfur Isotope Geochemistry of Acid Mine Drainage at the Mt. Morgan Mine Environment, Queensland, Australia. Applied Geochemistry, 20, 789-805. https://doi.org/10.1016/j.apgeochem.2004.11.004
[16]	Equeenuddin, M., Tripathy, S., Sahoo, P., & Panigrahi, M. (2010). Geochemistry of Ochreous Precipitates from Coal Mine Drainage in India. Environmental Earth Sciences, 61, 723-731. https://doi.org/10.1007/s12665-009-0386-9
[17]	Esfahani, H. K., & Datta, B. (2016). Linked Optimal Reactive Contaminant Source Characterization in Contaminated Mine Sites: Case Study. Journal of Water Resources Planning and Management, 142, Article ID: 4016061. https://doi.org/10.1061/(ASCE)WR.1943-5452.0000707
[18]	Ewl Sciences Pty Ltd. (2001). Contaminant Source Study. Mount Morgan. Prepared for Qld Department of Natural Resources and Mines.
[19]	Ferguson, K., & Morin, K. (1991). The Prediction of Acid Rock Drainage—Lessons from the Database. In Proceedings of the Second International Conference on the Abatement of Acid Drainage (pp. 1-4). Montreal, Canada.
[20]	Gasparon, M., Smedley, A., Jong, T., Costagliola, P., & Benvenuti, M. (2017). Acid Mine Drainge at Mount Morgan, Queensland (Australia): Experimental Simulation and Geochemical Modelling of Buffering Reactions. In Water in Mining Environments (pp. 433-436). Cagliari: International Mine Water Association.
[21]	Hunt, R. J., Doherty, J., & Tonkin, M. J. (2007). Are Models Too Simple? Arguments for Increased Parameterization. Ground Water, 45, 254-262. https://doi.org/10.1111/j.1745-6584.2007.00316.x
[22]	Jha, M., & Datta, B. (2012). Simulated Annealing Based Simulation-Optimization Approach for Identification of Unknown Contamination Sources in Groundwater Aquifers. Desalination and Water Treatment, 32, 79-85. https://doi.org/10.5004/dwt.2011.2681
[23]	Johnson, C., & Skousen, G. (1995). Minesoil Properties of 15 Abandoned Mine Land Sites in West Virginia. Journal of Environmental Quality, 24, 635-643. https://doi.org/10.2134/jeq1995.00472425002400040014x
[24]	McDonald, M. G., & Harbaugh, A. W. (1963). MODFLOW, a Modular Three-Dimensional Finite-Difference Ground-Water Flow Model (p. 586). Techniques of Water-Resources Investigations, Book 6, US Geological Survey.
[25]	MLA Meat & Livestock Australia (2008). Mount Morgan Climate History (p. 113).
[26]	Petit, C. (2016). Modelling and Characterization of Groundwater Contamination Sources in Abandoned Mine Sites. Internship Work, Douglas: James Cook University.
[27]	Rahmati, O., Samani, A. N., Mahmoodi, N., & Mahdavi, M. (2015). Assessment of the Contribution of N Fertilizers to Nitrate Pollution of Groundwater in Western Iran (Case Study: Ghorveh-Dehgelan Aquifer). Water Quality, Exposure and Health, 7, 143-151. https://doi.org/10.1007/s12403-014-0135-5
[28]	Rushton, K. R., & Redshaw, S. C. (1979). Seepage and Groundwater Flow: Numerical Analysis by Analogue and Digital Methods. New York, Chichester: Wiley.
[29]	Schulze-Makuch, D. (2011). Groundwater: Advection, Dispersion, Sorption, Degradation, Attenuation Vol. II. In Enciclopedia of Life Support Systems (p. 383). Paris, France.
[30]	Scotney, D. (2016). Optimal Identification of Unknown Sources of Groundwater Contamination in Abandoned Mine Sites. Undergrad Thesis, Douglas: The College of Science & Engineering, James Cook University.
[31]	Soltani, S., Asghari Moghaddam, A., Barzegar, R., Kazemian, N., & Tziritis, E. (2017). Hydrogeochemistry and Water Quality of the Kordkandi-Duzduzan Plain, NW Iran: Application of Multivariate Statistical Analysis and PoS Index. Environmental Monitoring and Assessment, 189, 1-20. https://doi.org/10.1007/s10661-017-6171-4
[32]	Taylor, G., Howse, R., Duivenvoorden, L., & Vicente-Beckett, V. (2002). Downstream Flow Event Sampling of Acid Mine Drainage from the Historic Mt Morgan Mine. Water Science and Technology, 45, 29-34. https://doi.org/10.2166/wst.2002.0376
[33]	Ulrich, T., Golding, S. D., Kamber, B. S., Zaw, K., & Taube, A. (2003). Different Mineralization Styles in a Volcanic-Hosted Ore Deposit: The Fluid and Isotopic Signatures of the Mt Morgan Au-Cu Deposit, Australia. Ore Geology Reviews, 22, 61-90. https://doi.org/10.1016/S0169-1368(02)00109-9
[34]	Wels, C., Findlater, L., Shaw, S., & Laurencont, T. (2004). Mount Morgan Mine—A Case Study of ARD Impacted Groundwater. In Proceedings of Mine Water 2004: Process, Policy and Progress (pp. 235-245). Newcastle upon Tyne, UK.
[35]	Wels, C., Findlater, L., & McCombe, C. (2006). Assessment of Groundwater Impacts at the Historic Mount Morgan Mine Site, Queensland, Australia. In The 7th International Conference on Acid Rock Drainage (pp. 2311-2331). St. Louis, Missouri, USA. https://doi.org/10.21000/JASMR06022311
[36]	Wels, C., Findlater, L., & McCombe, C. (2009). Contaminant Load Balance Study for Mount Morgan Mine, QLD, Australia. In The 8th International Conference on Acid Rock Drainage (pp. 1487-1496). Skelleftea, Sweden.
[37]	Younger, P. L., Wolkersdorfer, C., & Consortium, E. (2004). Mining Impacts on the Fresh Water Environment: Technical and Managerial Guidelines for Catchment Scale Management. Mine Water and the Environment, 23, s2-s80. https://doi.org/10.1007/s10230-004-0028-0
[38]	Zhou, H., Gomez-Hernandez, J. J., & Li, L. (2014). Inverse Methods in Hydrogeology: Evolution and Recent Trends. Advances in Water Resources, 63, 22-37. https://doi.org/10.1016/j.advwatres.2013.10.014
[39]	Zijl, W., De Smedt, F., El-Rawy, M., & Batelaan, O. (2017). The Double Constraint Inversion Methodology: Equations and Applications in Forward and Inverse Modeling of Groundwater Flow. Cham: Springer. https://doi.org/10.1007/978-3-319-71342-7_2