Investigating the Influence of Structure and Heterogeneity in Waste Rock Piles on Mass Loading Rates—A Reactive Transport Modeling Study

Raymond, Katherine E.; Seigneur, Nicolas; Su, Danyang; Mayer, K. Ulrich

doi:10.3389/frwa.2021.618418

ORIGINAL RESEARCH article

Front. Water, 17 May 2021

Sec. Water and Hydrocomplexity

Volume 3 - 2021 | https://doi.org/10.3389/frwa.2021.618418

This article is part of the Research Topic Modeling-based Approaches for Water Resources Problems View all 7 articles

Investigating the Influence of Structure and Heterogeneity in Waste Rock Piles on Mass Loading Rates—A Reactive Transport Modeling Study

$\nKatherine E. Raymond$ Katherine E. Raymond^*

Nicolas Seigneur^†

Danyang Su

K. Ulrich Mayer

Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, BC, Canada

Placement methods and material availability during waste rock pile (WRP) construction may create significant heterogeneities in physical and geochemical parameters (such as grain size, permeability, mineralogy, and reactivity) and influence the internal pile structure. Due to the enormous scale of WRPs, it is difficult to capture the influence of heterogeneities on mine drainage composition and evolution. Although laboratory- or field-scale experimental studies have provided much insight, it is often challenging to translate these results to full scale WRPs. This study uses a numerical modeling approach to investigate the influence of physical and chemical heterogeneities, structure, and scale on the release of acid rock drainage (ARD) through 2D reactive transport simulations. Specifically, the sensitivity of drainage quality to parameters including grain size distribution, sulfide mineral weathering rates, abundance and distribution of primary minerals, and pile structure as a function of construction methods are investigated. The geochemical model includes sulfide oxidation, pH buffering by calcite dissolution, and ferrihydrite and gypsum as secondary phases. Simulation results indicate that the implications of heterogeneity and construction method are scale-dependent; when grain size distribution trends observed in a pile's core are applied to the entirety of a pile, results between push- and end-dumping methods vary substantially—however, predicted drainage for different construction methods become more similar when features such as traffic surfaces, structural variation, and multiple benches are also considered. For all scales and construction methods investigated, simulated results demonstrate that pile heterogeneity and structure decrease peak mass loading rates 2 to 3-fold, but cause prolonged ARD release compared to the homogeneous case. These findings have implications for the economics of planning water treatment facilities for life of mine and closure operations.

Introduction

Numerous industries globally depend on a reliable supply of commodities produced from mining activities. During the mining process, large amounts of materials that are not economically viable (called waste rock) must first be removed to access the valuable commodities below. After drilling, blasting, and removing the waste rock, it must be stored long-term—often in large piles above ground, where materials are exposed to atmospheric weathering conditions. A waste rock pile (WRP) may be composed of more than 500 million m³ of material, with a wide range of grain sizes (McCarter, 1990). Due to the enormous volumes of materials, heavy equipment such as haul trucks and bull dozers are required for pile construction. Materials deposited from these large machines result in random, heterogeneous waste rock placement, though on a larger-scale, materials may be strategically mixed or confined to prescribed zones to meet closure plan strategies (e.g., MEND, 1998; Miller et al., 2006). Construction methods and type of material thus create significant variations in the distribution of grain size and mineralogy (with implications for sulfide content, metal concentrations, and pH-buffering capacity), which affects drainage quality (Smith et al., 1995; Munroe et al., 1999; Fala et al., 2003; Tran et al., 2003; Herasymuik et al., 2006; Stockwell et al., 2006; Lahmira et al., 2017). While the relationships between these characteristics in a pile and the potential to generate acid rock drainage (ARD), as well as related metal release and attenuation have been studied extensively (e.g., Sherlock et al., 1995; Stumm and Morgan, 1996; Nordstrom and Alpers, 1999; Akcil and Koldas, 2006; Amos et al., 2015; Dold, 2017), the influence of heterogeneity and spatial distribution of these characteristics on ARD has received less attention to date. This is due to the large scale of WRPs and typically insufficient sampling resolution to fully capture a pile's heterogeneity. Furthermore, typical laboratory-based methods (such as humidity-cells) or field-cell experiments do not capture additional structural features such as fining upwards or traffic surfaces (characterized by smaller grain size and lower porosity due to compaction from heavy equipment driving along exposed weathering surfaces during the construction process). To this end, numerical models have proven suitable tools for investigating the influence of physical and chemical heterogeneities (Molson et al., 2005; Fala et al., 2013; Lahmira et al., 2017; Pedretti et al., 2017; Appels et al., 2018; Muniruzzaman and Pedretti, 2020), as well as scale (Wilson et al., 2018; Vriens et al., 2020b) on ARD generation and attenuation processes in WRPs. However, at the time of writing a sensitivity analysis that examines the roles of physical and chemical parameter distribution, with a particular emphasis on the effect of pile construction methods, structural features, and pile scale, has not been completed.

This paper investigates the impacts of physical and chemical heterogeneities, scale and internal pile structure on drainage quality from a WRP using the reactive transport code MIN3P-HPC (Su et al., 2020). First, an overview of parameters controlling the development of ARD are discussed. MIN3P-HPC model outputs from simulated 2D cross sections of WRPs are then analyzed to explore the combined impacts of parameter heterogeneity and internal pile structures on basal drainage quality and mass loadings, as a function of scale. Correlations between spatial distributions of controlling parameters due to grain size (including hydraulic conductivity, porosity, mineral abundance, and weathering rates) as well as the effects of structural features (i.e., grain size distribution and segregation due to construction methods and traffic surfaces) are examined. The influence of pile scale is then explored by comparing results from a simulated single bench WRP to a pile with multiple benches and a larger volume of material. Heterogeneous pile simulations are compared to a reference case with average parameter values assumed uniform throughout the pile (i.e., a homogeneous case). Results from these simulations are discussed to evaluate effects on temporal and spatial development of ARD due to physical and chemical heterogeneities, structure, pile scale, and their combined effects on peak mass loadings and duration of ARD release. While mass loadings provide a measure for comparison between the homogeneous case and simulations with varying degrees of heterogeneity and structure, the purpose is to illustrate differences in solute release behavior in order to provide a suitable foundation for future site-specific numerical modeling studies rather than provide predictive modeling of release rates under defined conditions.

Physical Heterogeneity

WRPs exhibit large ranges in grain size, from clay-sized particles to automobile-sized boulders (McCarter, 1990; Morin et al., 1991; Smith et al., 1995; Stockwell et al., 2006; Anterrieu et al., 2010). In this study, the physical heterogeneities that are influenced by grain size which are considered include hydraulic conductivity and porosity. Hydraulic conductivity may vary by several orders of magnitude in a pile and typically follows a lognormal distribution (Smith et al., 1995; Tietje and Hennings, 1996). Variations in hydraulic conductivity values throughout a pile are expected to affect contaminant residence times and migration pathways (Jung and Navarre-Sitchler, 2018; Vriens et al., 2020a). Porosities throughout WRPs may also be extremely variable (Smith et al., 1995; Lefebvre et al., 2001; Azam et al., 2007; Zhang, 2017; Dimech et al., 2018). In this study, porosity is assumed to have a log-linear relationship with hydraulic conductivity, following the work of Smith et al. (1995). During pile construction, infilling with fines and compaction along traffic surfaces may reduce both hydraulic conductivity and porosity, as observed in previous studies such as Bellehumeur (2001) and Cash (2014)—effects of which are also considered in our study.

Chemical Heterogeneity

WRPs may exhibit highly heterogeneous mineral abundances depending on the geology of the deposit and waste rock source material. Even if a pile is considered relatively homogeneous in terms of rock type, naturally-occurring variations in mineralogy can lead to the development of pockets of higher reactivity and acidity (e.g., Atherton, 2017; Pedretti et al., 2017). To investigate the influence of chemical heterogeneities on ARD-related processes, we consider here local variations in NP:AP ratios (where NP is the neutralization potential from pH-buffering minerals, and AP is the acid production potential from sulfide minerals) and sulfide weathering rates. The availability of sulfide minerals and associated oxidation rates have been shown to depend on particle size. For example, several waste rock studies have found that sulfide mineral availability is higher in the more fine-grained materials (Aranda, 2010; Erguler and Kalyoncu Erguler, 2015; Elghali et al., 2019). Furthermore, the reactivity of sulfide minerals is controlled by mineral properties and reactive surface areas, with an inverse relationship between sulfide mineral oxidation rate and particle diameter (Nicholson et al., 1988; Komnitsas et al., 1995; Blowes et al., 2003; Cornell and Schwertmann, 2003; Vriens et al., 2020b).

Pile Structure

Rock placement and construction of a WRP controls the internal structure, influencing and overprinting the physical and chemical heterogeneities mentioned previously. While a combination of construction methods may be used within a single pile, for the purpose of this study, differences between classically described features of push- and end-dumping construction methods were investigated. Features of these construction methods are summarized in Table 1 and illustrated in Figure 1. The investigation of engineered waste rock storage techniques (such as covers, capillary barriers, or the placement of more acid generating rock within non-acid generating cover material) to mitigate ARD generation (e.g., Yanful and Woyshner, 1993; O'Kane et al., 1998; Haug and Pauls, 2002; MEND, 2004; Fala et al., 2005; Amos et al., 2009; Raymond et al., 2020) was beyond the scope of this study.

TABLE 1

Table 1. Structural features for end- and push-dumping construction methods.

FIGURE 1

Figure 1. Illustrated conceptual models of push- and end-dumping construction methods in a WRP. In the end-dumped pile, different patterns represent distinct materials emplaced along the angle of repose. Darker layers at the top of each bench in the end-dumped pile or within the center of the push-dumped pile represent compacted traffic surfaces.

End-Dumping

In end-dumping methods, rock is typically dumped over a sloped pile face via haul truck. This technique results in observable effects of gravity sorting: larger particles and boulders are driven by momentum toward the toe of the pile, whereas significantly more of the finer materials remain close to the surface (Nichols, 1986; Morin et al., 1991; Smith et al., 1995; Smith and Beckie, 2003; Stockwell et al., 2006; Wilson, 2011, 2017; Amos et al., 2015). Nichols (1986) described three distinct horizons observed due to gravity-sorting in an end-dumped pile: (1) the upper horizon, which contains a high concentration of fine materials (which is further compacted at the surface due to heavy equipment traffic as the pile grows laterally); (2) an evenly distributed/evenly graded horizon of fine and coarse particles along the remainder of the slope up to the toe; and (3) a bottom horizon beyond the slope's toe containing a high concentration of coarse fractioned materials and boulders. Relative proportions of fining upwards trends throughout the pile are consistent regardless of the height of the pile bench; however, the degree of particle size segregation depends on a pile's slope length and origin of material (since rock type will affect influencing factors such as density, strength, and shape). These structural trends are common and have been observed in several WRPs (Herasymuik et al., 2006; Stockwell et al., 2006; Azam et al., 2007; McLemore et al., 2009; Barsi, 2017; Zhang, 2017; Vriens et al., 2019; St-Arnault et al., 2020).

Push-Dumping

In push-dumping methods, rock is dumped onto the surface of a pile, followed by pushing and leveling out materials using tractors and shovels. As a result, fine-grained and sub-horizontal compacted layers occur in the internal structure due to operating equipment along traffic surfaces as the pile grows (Lefebvre et al., 2001; Fala et al., 2003; Anterrieu et al., 2010). Coarse fractions are often embedded in a collection of fine aggregates at the slope crest while being pushed along the pile surface. As a result, a random distribution (without observable fining upwards trends) occurs within the core, overlain by traffic surfaces (Anterrieu et al., 2010). Materials pushed over the slope crest may experience the effects of gravity sorting, which may cause fining upwards along the batters similar to end-dumping, but spatially more restricted. These features characteristic to push-dumping construction have been observed by Fines (2006), Poisson et al. (2009), and Mele et al. (2013).

Pile Scale

Internal structures defined in the previous section were described in terms of construction of a single bench. However, piles are often built in multiple terraced benches, where internal features are repeated, separated by traffic surfaces as shown in Figure 1 (e.g., Anterrieu et al., 2010; Broda et al., 2017; Martin et al., 2017; Vriens et al., 2018; St-Arnault et al., 2020). Typically, the process of scale-up is not quite as simple with respect to quantifying the differences in drainage quality and mass loadings. Often, laboratory- or field-based experiments such as humidity-cells, column tests, or field-barrels are utilized to characterize mass loadings from mineral leaching rates using empirical scale factors for a full-scale pile, which may be insufficient for accurate water quality prediction (Amos et al., 2015; Parbhakar-Fox and Lottermoser, 2015; Wilson et al., 2018; Vriens et al., 2019). While these issues regarding scale-up have been addressed and improved upon in several laboratory studies (e.g., Malmström et al., 2000; Kempton, 2012; Plante et al., 2014), these experiments do not consider structural influences such as fining upwards, traffic surfaces, or differences between the core and batters of a pile, and do not fully address the contribution to effluent quality from pile structure and scale. Supported by recent advancements to reactive transport codes and computing capabilities, numerical modeling provides a suitable tool for investigating the influence of pile structure and scale-up related issues (Wilson et al., 2018; Vriens et al., 2020b).

Methods

Modeling Approach

MIN3P (Mayer et al., 2002) is a three-dimensional multicomponent reactive transport code designed for applications in variably saturated porous media. The unstructured grid version used in this study, MIN3P-HPC (Su et al., 2020), was developed based on the parallel version of the MIN3P-THCm code (Bea et al., 2016; Su et al., 2017). The code uses a polygonal (2D, e.g., triangular, quadrilateral) or polyhedral (3D, e.g., tetrahedral, prismatic, hexahedral) mesh, although a standard regular gridded mesh available in previous versions may also be applied. The versatility of the triangular mesh allows structural features such as sloped batters and deposition following the angle of repose within a WRP to be captured, while previous numerical models of piles using MIN3P were constrained to only the core of a pile due to limitations related to the structured grid.

MIN3P-HPC utilizes a node-centered finite volume approach. Values are assigned to each node and then interpolated for the interfaces of adjacent cells, which defines a representative elemental volume; the degree to which these values are upscaled is therefore controlled by cell size and mesh discretization. A global implicit solution approach is applied to both variably saturated flow and reactive transport equations (Mayer et al., 2002). Additional details regarding the governing equations and numerical methods of MIN3P-HPC are provided in Mayer et al. (2015), Su et al. (2017), and Su et al. (2020).

Domain Geometry, Boundary, and Initial Conditions

In this study, triangular meshes with an average resolution of 1.0 m were used to generate single bench WRPs (Figure 2). Each triangular cell represents an average cross-sectional area of ~0.4 m². Models including traffic surfaces and internal features representative of those found within push- and end-dumped piles were also investigated for structure and scale. Additional mesh configurations for single bench models (20 m height and 2D cross-sectional area of 1,456 m²) and multiple bench models (60 m height and 2D cross-sectional area of 10,032 m²) were completed for models with increased structural complexity. Examples of these meshes are included for the single-bench models accounting for pile structure in Supplementary Figure 1.

FIGURE 2

Figure 2. Conceptual model for single bench, 2D WRP and triangular mesh. Input parameters shown here are consistent in all models investigated for the single bench scenario. For multi-bench models, the same infiltration rate was applied along the top and slopes of exposed surfaces, with the same constant head boundary assigned at the base.

All simulations assumed a recharge rate of 300 mm year⁻¹ applied to the pile surface and slopes, comparable to recharge rates observed in semi-arid northern climates (e.g., Lefebvre et al., 2001; Neuner et al., 2013). A constant hydraulic head boundary of −0.5 m was assigned to the base of the pile to allow free drainage into a basal collection system (not simulated). In all simulations, longitudinal and transverse vertical dispersivity coefficients of 5 × 10⁻² m and 5 × 10⁻⁴ m were chosen, respectively, although the influence of dispersivity was found to be negligible compared the effects of spatially variable hydraulic conductivity and chemical heterogeneity on the scale of grid resolution. Gases considered include O₂ and CO₂, which were held constant at values of 21 and 0.5%, respectively. These values were chosen based on observations at the Diavik Diamond Mine and are thought to be representative of a real-world, well-aerated pile with low sulfide and carbonate content (Amos et al., 2018). All single bench simulations were run for 300 years under steady-state flow conditions to evaluate long-term pile evolution and effluent response to heterogeneities.

Parameter Heterogeneity

Heterogeneity was accounted for by relating physical and chemical parameters to a mean grain size diameter assigned to each node. In this way, parameter distributions of hydraulic conductivity, porosity, mineralogical assemblage and sulfide mineral oxidation rates were determined to mimic a range of parameter values observed in real-world piles. For all simulations, an average mean grain size diameter (d_m) of 5 mm was considered, with local mean grain size diameter (d_m,i, where the subscript i denotes the node) values ranging between 0.05 and 170 mm, following a lognormal distribution as observed in previous WRP case studies (e.g., Smith et al., 1995; Tietje and Hennings, 1996; Erguler and Kalyoncu Erguler, 2015; Zhang, 2017). Although other representative grain size diameters (such as d₅, d₁₀, d₈₀ etc.) may be used to derive values for parameters investigated in this study (Chapuis, 2012; Zhang, 2017), the mean grain size diameter was used for consistency with the geometrical estimation method for pyrite abundance (described below). These relationships were used to develop a range of representative parameter values for the stylistic models, based on information from multiple mine sites, rather than to calculate absolute values in the context of a site-specific calibration.

To identify the influence of varying spatial distributions of mean grain size and the related parameters, five random, lognormal distributions of d_m,i were compared to five random, lognormal distributions with fining upwards trends. Going forward, these two configurations are referred to as “random” and “fining upwards.” The d_m,i values were produced using the rnorm() function in Base R (R Documentation, 2020) to generate random variates for the random distributions, or random variates scaled by elevation to account for the effects of gravity sorting in fining upwards cases. Additional details on these rnorm() functions are included in Section 2 of the Supplementary Material. These models with random and fining upwards trends represent heterogeneous distributions displaying features common to the core of a pile from push- and end-dumping construction methods, respectively. In both situations, material placement is random in nature, which is captured by our approach of assigning grain size and associated parameter distributions randomly. This method does not explicitly account for preferential flow pathways and produces matrix-dominated flow within the pile. This approach is similar to methods previously used by Lahmira et al. (2017) and Wilson et al. (2018) and is in line with recent observations at the Antamina Mine, Peru, where flow through experimental waste rock piles was dominated by matrix flow with relatively small contributions of preferential flow (Blackmore et al., 2014; Pedretti et al., 2017; Seigneur et al., 2021). Although the nature and distribution of mass loadings is affected by preferential flow (especially in the short term), the impact on long-term mass loadings is expected to be more limited. Examples of parameter distributions for random and fining upwards cases are given in Figure 3. Additional information on the parameter distributions used in the simulations are included in the Supplementary Figures 2–4.

FIGURE 3

Figure 3. Cross-sections showing the initial conditions for the single bench models including spatial distributions of hydraulic conductivity, pyrite volume fraction, and effective rate coefficients for pyrite within the pile for various scenarios. For the heterogeneous case (random and fining upwards), one of the five configurations that was used to test parameter sensitivity is shown. Porosity is also spatially varied following the relationship to hydraulic conductivity shown in Equation 2. Histograms showing the nodal distribution for these parameters are shown in Supplementary Figure 4.

Hydraulic Conductivity

The relationship between hydraulic conductivity (K) and d_m was developed given the relationship between K and a representative grain size diameter related to the Hazen equation (Freeze and Cherry, 1979) (Equation 1). Considering the large spread in grain size typical for waste rock, the representative grain size diameter was assumed to be d₅₀ (or d_m) using an approach by Freeze and Cherry (1979) for materials with a higher mean grain size than is typical for soils. In this context, the mean hydraulic conductivity was adjusted using a scale factor α_K to yield a mean hydraulic conductivity of 10⁻⁵ m s⁻¹, with local hydraulic conductivity values K_i ranging between 10⁻⁹ and 10⁻¹ m s⁻¹. These values are consistent with observations made from WRPs (e.g., Morin et al., 1991; Smith et al., 1995; Azam et al., 2007; McLemore et al., 2009; Fala et al., 2013; Neuner et al., 2013).

\begin{array}{l} K_{i} = α_{K} \bar{K} {(\frac{d_{m, i}}{d_{m}})}^{2} & (1) \end{array}

Porosity

Porosity (ϕ_i) throughout the pile was assumed to be related to K_i for all simulations. For simplicity, a log-linear relationship between hydraulic conductivity and porosity (Equation 2) was chosen based on field observations by Smith et al. (1995), which typically demonstrated larger porosities in coarser grained waste rock. A range of porosity values between 0.1 and 0.5 with an average porosity of 0.3 was used based on field observations from WRPs (e.g., Morin et al., 1991; Lefebvre et al., 2001; Fala et al., 2003; Neuner et al., 2013; Zhang, 2017; Dimech et al., 2018).

\begin{array}{r} ϕ_{i} = ϕ_{min} + (\frac{ϕ_{max} - ϕ_{min}}{{log}_{10} K_{max} - {log}_{10} K_{min}}) ({log}_{10} K_{i} - \\ {log}_{10} K_{min}) & (2) \end{array}

Models assuming a homogeneous porosity were initially run and compared to results of models with spatially distributed porosity values defined by the relationship in Equation (2). Differences in these results between cases with a uniform porosity and spatially distributed porosities as a function of K were minimal. Since a heterogeneous distribution of porosity influenced by grain size distribution is considered more realistic of real-world conditions, only results from these cases are presented.

Reactive Minerals and NP:AP Ratios

While a variety of site-specific minerals must be considered when evaluating the long-term behavior of drainage quality for a WRP, a simplified representation was used in this conceptual analysis to assess the impact of chemical heterogeneity on ARD generation and attenuation. Four minerals were considered in this study: pyrite, calcite, gypsum, and ferrihydrite. This intentionally simplified mineralogy allows for conceptualization of the effects of spatially varied NP:AP ratios, sulfide weathering rates, as well as the role of secondary mineral formation, which will affect the drainage quality in more complex systems. Average initial pyrite and calcite contents of 3.4 and 0.9 wt.% respectively were assumed to achieve potentially acid generating (PAG) conditions. These abundances are in line with values observed at mine sites (e.g., Lawrence and Wang, 1996; Azam et al., 2007; McLemore et al., 2009; Aranda, 2010). Gypsum and ferrihydrite were permitted to precipitate to account for the effects of secondary mineral formation on iron and sulfate release.

To investigate the effects of spatial distributions in NP:AP ratios (here, the NP comes from calcite, and the AP comes from pyrite), the spatial distribution of pyrite volume fractions varied while calcite remained homogeneously distributed. This variation in NP:AP assessed the implications of local pockets of elevated AP on the overall basal drainage quality. Initial pyrite volume fractions ( $φ_{i}^{0}$ ) were considered inversely proportional to d_m (Equation 3) to account for an increased sulfide mineral availability in the fines, as observed by Aranda (2010). This relationship was then scaled appropriately using a scale factor (α_φ) to maintain an average 3.4 wt.% pyrite in all cases.

\begin{array}{l} φ_{i}^{0} = α_{φ} \frac{1}{d_{m, i}} & (3) \end{array}

Pyrite Oxidation Rate

In this simplified ARD conceptual model, oxidative pyrite dissolution was considered the rate-determining step controlling drainage quality. Pyrite oxidation kinetics may be influenced by multiple factors such as microbial activity, gas or heat transfer, surface area, etc. (Nordstrom, 1982; Nicholson et al., 1988; Blowes et al., 2003; Amos et al., 2015). In these numerical simulations, factors influencing oxidative pyrite dissolution rates are captured in the effective rate coefficient, k_eff. To investigate the influence of weathering rate variability, pyrite oxidation rates were spatially varied, while calcite weathering rates and precipitation/dissolution rates for gypsum and ferrihydrite were considered as quasi-equilibrium reactions (Table 2) (Blowes et al., 2003; Moncur et al., 2009). Local effective rate constants ( $k_{e f f, i}^{0}$ ) [mol L⁻¹ bulk s⁻¹] were assumed to be inversely related to grain size using the geometrical estimation method (Aranda, 2010) as shown in Equation 4. This method considers the relationship between grain size diameter and reactive surface area, which is accounted for the in k_eff term.

\begin{array}{l} k_{e f f, i}^{0} = α_{k_{e f f}} \frac{d_{m}}{d_{m, i}} & (4) \end{array}

The effective rate coefficients $k_{e f f, i}^{0}$ were scaled by α_{k_eff} to maintain an average initial rate constant ( $k_{e f f}^{0})$ of 1 × 10⁻⁹ mol pyrite L⁻¹ bulk s⁻¹, producing pyrite oxidation rates comparable to those reported in a literature review by Lorca et al. (2016). Total pyrite dissolution rate evolves over time following a two-third power law, representing a surface-controlled reaction with a progressively decreasing particle radius (Lichtner, 1996; Seigneur et al., 2019). Reaction pathways and rate expressions for minerals considered in all simulations are summarized in Table 2.

TABLE 2

Table 2. Mineral rate expressions for stylistic models.

Pile Structure

The influence of pile structure was explored through the addition of structural complexities including traffic surfaces and internal grain size trends for push- and end-dumping methods, as shown in Figure 3. In the simulations, traffic surfaces were assumed to occur within the upper 0.5 m of a lift with a compaction factor scaled based on the d_m,i. This scaling follows a non-linear relationship (Supplementary Table 1) for hydraulic conductivity, porosity, and sulfide oxidation rates due to the compaction and crushing of materials from heavy equipment (Nichols, 1986; Herasymuik et al., 2006; Stockwell et al., 2006; Wilson, 2017). However, it was assumed that pyrite volume fractions are dependent on source material during deposition, and therefore were not correlated to smaller d_m,i values in traffic surfaces.

Referring to Figure 3, one of each of these five random and fining upwards trending piles were chosen and labeled as the “heterogeneous” cases. To investigate the influence of traffic surfaces and the combination of these trends based on the descriptions of each construction method, two single bench “structure” models were developed. For the push-dumped structure model, multiple traffic surfaces were included at 5 m height increments throughout the core, with random grain size distributions below these compacted surfaces. Fining upwards trends were included along the slopes at the end of each traffic surface, representative of gravity-sorting along the batters. In the end-dumped structure model, multiple zones defined along the angle of repose represent different materials deposited as a pile is constructed (where each zone has a different average d_m, though the total average parameter values for the entire pile are maintained), with a single traffic surface at the top of the pile. These different zones are used to represent heterogeneities in source rock material which may be more pronounced based on the mine site and may consist of several different host rock materials (e.g., Azam et al., 2007; McLemore et al., 2009). Comparison between the “heterogeneous” and “structure” cases helps to distinguish the influence between general grain size trends (i.e., random vs. fining upwards) from the combined effects of traffic surfaces and grain size variations, which are conceptualized and discussed below.

Pile Scale

To identify the influence of increasing pile size and height, the “structure” models featured in Figure 3 were expanded upon to the “multi-bench” models in Figure 4. These simulations include the same structural features as described in the push- and end-dumping structure models, but with three stacked benches separated by traffic surfaces. These models were used to investigate whether the implications of structure and traffic surfaces from single bench model results hold up under multiple benches and over a larger pile scale. All modeled WRPs share average parameter values between single-bench and multi-bench models to compare to the homogeneous model. This allows for isolated observations to be made with respect to pile size and dimensions, as well as the structural implications which may be missed under traditional laboratory-based testing in site-specific studies.

FIGURE 4

Figure 4. Hydraulic conductivity distribution for multiple bench models in cross-section. These cases expand upon the structure single bench cases from Figure 3 to include three stacked benches, with a homogeneous multi-bench model for comparison.

Assumptions and Limitations

As mentioned, infiltration was considered constant and flow was assumed to be at steady-state for all simulations. In reality, the wet-up phase would begin during pile construction (which may take several years to complete) and flow may fluctuate due to site-dependent seasonal variations and individual recharge events (Smith et al., 1995). To isolate the impact of heterogeneities and structure on drainage water quality, we chose a steady-state approach for the flow problem, which is justified for the purposes of this numerical experiment and comparisons, given the time frame of simulations over several 100 years. Considering uncertainties associated with the distribution of recharge across the benches and slopes of the pile, we assumed uniform recharge across the entire pile surface. To limit the complexity of the simulations, residual saturation (S_wr = 0.22) and van Genuchten parameters (α = 5.0 m⁻¹; n = 2.7) were kept constant in all simulations, determined based on values observed at other mine sites (e.g., Lefebvre et al., 2001; Fala et al., 2003; Azam et al., 2007; Neuner et al., 2013; Peterson, 2014). Although these parameters deviate with grain size, saturated hydraulic conductivity values provide a dominant control on the calculated effective saturations and residence times within the model; therefore, homogeneous residual saturation and van Genuchten parameters were considered acceptable for the purposes of this investigation. Two cases that considered heterogeneity for these soil hydraulic function parameters were tested to validate this assumption (for one random and fining upwards model) and are included in Section 4 of the Supplementary Figures 5, 6. The effect of varying these parameters was limited.

In addition, isothermal conditions were assumed, and modeled piles were considered well-aerated with gas concentrations held uniform through time, ensuring that pyrite oxidation was not limited by oxygen availability. As a result, diffusion in the gas phase was not simulated. The simulation of highly reactive piles, leading to oxygen depletion, internal heating, and/or the accumulation of high levels of CO₂ was beyond the scope of this study. It is possible that zones of low hydraulic conductivity could limit oxygen availability due to elevated water saturations, limiting gas ingress into a WRP (e.g., Parbhakar-Fox et al., 2013). It is equally possible that the exothermic nature of sulfide mineral oxidation leads to an increase in reaction rates and enhanced oxygen ingress (Lefebvre et al., 2001). The focus here is on aspects of heterogeneity, structure and scale, without the complexities of gas transport limitations. We refer the reader to several studies specifically focusing on influence of gas transport and the exothermic nature of sulfide oxidation in WRPs (Lefebvre et al., 2001; Chi, 2010; Lorca et al., 2016; Steinepreis, 2017).

The mineralogy chosen for simulations was simplified and excludes processes such as weathering of aluminosilicate minerals, ion exchange, adsorption, surface complexation, or additional secondary mineral formation. These latter reactions also play a role in buffering pH and controlling pore water chemistry, and are of particular importance for simulating metal attenuation in more complex systems. However, the purpose of this study was to investigate the ways in which physical and chemical heterogeneities (controlled by structure and spatial distributions in grain size) and scale influence ARD and overall drainage quality. While the simplifications made here may be important contributing factors for site-specific investigations, the conceptual nature of this study is considered suitable for the purposes of isolating the effects of heterogeneity, structure, and scale compared to the homogeneous model.

Simulation Results

Sensitivity to Parameter Heterogeneity

Simulated variations in water saturations and the flow regime for the different cases are shown in Figure 5. While there were some minor horizontal components observed in the heterogeneous cases, the results illustrate that gravity-dominated drainage dominates throughout the WRPs at relatively low water saturations for all cases. This gravity-dominated drainage was also observed in alternative simulations where van Genuchten parameters were varied as a function of grain size. The effect of variable van Genuchten parameters on drainage patterns and rates was found insignificant and was therefore not further considered in our scenario analysis.

FIGURE 5

Figure 5. Saturation profiles with Darcy flux vectors represented as arrows (top right: arrow for scale) which indicate flow direction and their relative magnitude to illustrate influence of physical heterogeneities on flow. In all cases, gravity drainage is the dominant flow control, as flux vectors are predominately downwards in all cases. The boundary condition at the base was set to h = −0.5 m for all cases to mimic discharge into a basal drainage layer.

In order to compare results from different simulations, sulfate mass loadings (as mg ${SO}_{4}^{2 -}$ kg waste rock⁻¹ week⁻¹) (Table 3 and Figures 6, 7) and pH distributions within the pile (Figure 8) are used to compare the influence of heterogeneities on drainage quality from pyrite oxidation and ARD generation to the homogeneous model. Sulfate mass loadings were normalized per kg of waste rock so that volume of drainage water and mass of waste rock in each simulation are inconsequential to comparisons. In addition to the results presented below, information regarding the presence of calcite, pyrite, and differences in spatial depletion through time within simulated piles can be found in the Supplementary Figures 7–9. For visualization of the temporal evolution of pH, animations are also available for download (see animations available for download and Supplementary Table 2).

TABLE 3

Table 3. Variation in maximum mass loading rates and reduced sulfate release times with parameter heterogeneity.

FIGURE 6

Figure 6. Comparison of the effects of heterogeneity on basal sulfate mass loading rates, influenced by physical and chemical heterogeneities. Plots on the left show results from models with “random” material property distributions; models that include “fining upwards” trends are shown on the right. From the legend: K is hydraulic conductivity, φ_pyr is the pyrite volume fraction (controlling differences in NP:AP ratios), and k_eff is the effective rate of pyrite oxidation, and “c” refers to constant (or homogeneous spatial distribution for that parameter). Curves with the same colors represent the various realizations of a case with different randomly generated spatial distribution of parameters. Specifically, Red: reference/homogeneous case. Yellow: K is spatially varied, φ_pyr and k_eff remain constant. Green: K and k_eff spatially varied, φ_pyr is constant. Blue: K and φ_pyr are spatially varied, k_eff is constant. Purple: K, φ_pyr, and k_eff are spatially varied. In all cases where K is not constant, porosity is also spatially varied.

FIGURE 7

Figure 7. Sulfate basal drainage evolutions shown for varied pile structure and scale. Results for total sulfate mass loading rates are normalized per kilogram of waste rock are shown for different scenarios illustrated in Figure 4. As the pile becomes more complex through the addition of traffic surfaces and internal structural trends, differences between push- and end-dumping construction are dampened.

FIGURE 8

Figure 8. pH cross sections at t = 3 years comparing “heterogeneous” vs. “structure” models. To visualize the progression of pH evolution over 50 years, see animations and Supplementary Table 2.

Physical Heterogeneity

For cases with a heterogeneous distribution of hydraulic conductivity and porosity and homogeneous chemical parameters, simulation results show a reduction and delay in peak mass loading rates (Figure 6). Maximum mass loading rates decrease by 6 and 10% compared to the homogeneous reference case (Table 3). However, the simulation results suggest that variations in hydraulic conductivity and porosity alone do not significantly influence the duration of sulfate release; in Figure 6, the tailing of sulfate release only extends marginally longer than the homogeneous reference case. Results between the five different realizations did not vary significantly, suggesting that exact spatial location of physical heterogeneities is not significant for total drainage quality.

Chemical Heterogeneity

In addition to the physical heterogeneities discussed in the previous section, the influences of chemical heterogeneities were assessed by isolating three scenarios: (1) heterogeneous distribution of NP:AP ratio, (2) heterogeneous distribution of pyrite oxidation rates, and (3) a combination of the two. Results from scenarios where chemical heterogeneities were considered in addition to physical heterogeneities show a more pronounced decrease in peak mass loading rates in comparison to the homogeneous case (ranging from 12 to 34%) (Figure 6 and Table 3). In addition, results indicate that chemical heterogeneities lead to a prolonged release of ARD in comparison to the homogeneous case; it takes at least twice as long (118–141 years) until ${SO}_{4}^{2 -}$ mass loadings return to <1 mg kg⁻¹ week⁻¹ than for the homogenous case (53 years) (Table 3). Similar to the test cases restricted to varying physical heterogeneities, the different realizations of random material distributions (i.e., multiple curves of same colors in Figure 6) did not affect the magnitude and timing of mass loading rates significantly.

Pile Structure

Random vs. Fining Upwards Trends

Although parameter distributions followed the same general statistical trends for both the random and fining upwards models, differences in simulation results between these cases suggest that the evolution of the overall basal effluent quality is also affected by internal pile structure. In all cases, the maximum mass loading rates for sulfate were greater for fining upwards trends compared to a purely random distribution. For the fining upwards case, a higher abundance of fines in the uppermost section retains higher moisture contents and restricts flow due to lower hydraulic conductivities, causing increased residence times at the top of the pile. When chemical heterogeneities are considered, lower NP:AP ratios and higher pyrite oxidation rates are abundant in the fines due to the inverse relationship to d_m used in the numerical models. This causes a build-up of poor-quality water near the pile surface (see animations and Supplementary Table 2), which then must travel through the entire pile before exiting at a basal drainage point. As a result, a time-lagged but more rapidly growing peak was observed for mass loading rates in fining upwards models, when compared to the models with random material property distributions (Figure 6 and Table 3). In the random models, the spatial variations in hydraulic conductivity, NP:AP ratios, and pyrite oxidation rates are more evenly distributed throughout the pile (see animations and Supplementary Table 2). As a result, pyrite oxidation products are observed at an earlier time and continue to be released for a longer duration, but have lower maximum mass loading rates compared to the fining upwards models. Specifically, when comparing these “heterogeneous” cases (that varied both pyrite distribution and reactivity) to the homogeneous reference case, fining upwards trends led to an average decrease of peak mass loading rates by 12% compared to a 32% decrease for the random models (Table 3). Additionally, the fining upwards models prolonged sulfate mass loadings exceeding 1 mg kg⁻¹ week⁻¹ for an average of 121 years, while the random cases produced elevated ${SO}_{4}^{2 -}$ mass loadings for a longer time period (averaging 138 years).

While the simulations presented so far account for the influences of physical and chemical heterogeneities and for fining upwards vs. random spatial distributions, WRPs often include a combination of these distributions combined with compacted traffic surfaces and multiple benches. These contributions from traffic surfaces and complex pile structure are addressed in the next section.

Traffic Surfaces and Increased Structural Complexity

In order to assess the influence of traffic surfaces and increased structural complexity on simulation results, Figure 7 shows the total sulfate mass loading rates compared for the “structure” models for push- and end-dumping construction methods (Figure 3). These simulations take into consideration physical (hydraulic conductivity, porosity), as well as chemical (NP:AP ratios and pyrite reactivity) heterogeneities. The results can further be visualized through pH distributions in cross-sections shown in Figure 8, and in animations included in the Supplementary Material. While the differences in mass loading rates between the “heterogeneous” cases (i.e., random versus fining upwards models) are significant and suggest that grain size trends may play an important role on drainage quality (Figure 7), results from the “structure” and “multi-bench” models that consider traffic surfaces and increased structural complexity become very similar for push- and end-dumping construction methods (Figure 7). Furthermore, the consideration of additional structural complexity resulted in a further reduction in peak mass loading rates relative to the homogeneous reference case for both scenarios considered. More detailed differences between the behaviors, as a function of construction method, are described below.

Push-Dumping Models

The single bench “structure” push-dumping model (Figure 3) includes four traffic surfaces with decreased grain size, porosity, hydraulic conductivity, and increased pyrite oxidation rate overlying randomly distributed coarser material within the pile core, with fining upwards trends along the slopes. As discussed previously, the same average parameter values were maintained in all simulations to ensure a fair comparison to the homogeneous reference case (Figure 7). In this model, maximum sulfate mass loading rates are reduced, and release occurs over a longer duration (Table 4). The random spatial distribution in grain size in the heterogeneous models results in spatial variations in ARD-generating processes laterally and with depth, averaging out the release of products from sulfide oxidation. The addition of multiple compacted traffic surfaces and structure in the model acts to further average out these trends, resulting in a 36% decrease in peak mass loading rates compared to the homogeneous model. In addition, fining upwards along the pile slopes resulted in a slight time delay in mass loadings, although less significant than for the fining upwards and end-dumping construction models.

TABLE 4

Table 4. Variation in mass loading rates and duration of release with traffic surfaces and structure.

End-Dumping Models

The single bench “structure” end-dumping model (Figure 3) maintains the fining upwards grain size trends from the heterogeneous model, but assumes zones of different average parameter values following the angle of repose overlain by a single traffic surface. While the addition of a traffic surface influences the parameter values similarly through the mechanisms described for the push-dumping model, it has less influence on the results since there is only one traffic surface present at the pile surface—which already hosts a larger portion of fines due to the effects of gravity sorting. Instead, the presence of different zones (representative of different materials during pile construction) contributes most significantly to sulfate mass loading rates, as it acts to even out parameter distributions laterally through the pile rather than only by depth (as in the model without added structural complexity). As a result, the maximum sulfate mass loading rate is greatly reduced compared to the homogeneous model (38% in the “structure” model compared to 12% in the “heterogeneous” case), with an increased duration of poor-quality drainage release. Although simulated sulfate release is evened out by the presence of material zoning, the time-lagged peak characteristic of fining upwards trends is still observed in the structure end-dumping model. However, it should be noted that these results summarized in Figure 7 and Table 4 represent basal drainage that is entirely contained and perfectly mixed. Therefore, although the addition of zones in the end-dumped structure case presented here resulted in simulated reduced overall maximum mass loading rates, in reality these trends may change due to factors such as seepage discharge within the pile or sampling location, which locally may be more influenced by the factors described in the “heterogeneous” end-dumping (or fining upwards) scenario.

Pile Scale

To account for the change in size and shape, a multi-bench homogeneous case was used as reference for the multi-bench heterogeneous models. Comparing the two homogeneous models in Figure 7, the homogeneous multi-bench scenario showed reduced mass loading rates compared to the single bench model. Since reported mass loading rates are normalized per kilogram of waste rock in models, the total mass of sulfate released would be higher in the multi-bench scenario due to a larger volume of material. However, it is important to note that this 40% decrease in normalized mass loadings rate is significant, influenced only by size and shape when comparing the homogeneous models; the addition of pile height over multiple benches increases the spread and therefore duration of ARD-related product release and flushing of materials.

Since the multi-bench models were adapted from the single-bench models including structure over a three-bench scenario, the influence of traffic surfaces and internal structures on sulfate mass loading and ARD-generating processes described in the previous section also apply to the multi-bench models. However, as shown in Figure 7 and Table 4, not only is there a reduction in maximum mass loading rates for sulfate in the multi-bench models due to this influence of size and shape, but the peak release rates become essentially equal for both push- and end-dumping models. Generally, the effects of gravity sorting on a time-lagged release are observed in both multi-bench models, although much more pronounced in the case of the end-dumping scenario, with multiple delays observed in the sulfate mass loading curve. Most importantly, results suggest that the addition of heterogeneity, pile structure, and multiple benches significantly influence predicted effluent quality. Compared to the single bench homogeneous model, stylistic multi-bench models reduced sulfate maximum mass outflux rates by a factor of 3—an amount that would likely have significant implications for long-term treatment and closure planning.

Discussion

Parameter Sensitivity

Simulation results discussed above suggest that drainage quality is significantly affected by physical and chemical heterogeneities in waste rock piles, while considerations such as dispersivity and van Genuchten parameters had limited influence in these simulated field-scale piles. If numerical models are used to make comparisons of different WRP storage strategies, incorporating these physical and chemical heterogeneities (particularly, mineral weathering rates and their relation to grain size distribution) yields substantial variations in simulated peak mass loading rates and duration of contaminant release. While results summarized in Table 3 demonstrate that parameter heterogeneity is an important consideration, findings suggest that its influence on results is ultimately dependent on the spatial distribution, internal structures, and pile scale.

Pile Structure and Scale

Results from overall basal drainage suggest simulated differences in structure due to construction methods in a well-aerated pile matters to a certain scale, but does not strongly influence the outflow chemistry over the entirety of the pile once more complex internal structures and traffic surfaces are considered. In the random and fining upwards “heterogeneous” models, these spatial trends produced quite different results in sulfate mass loading rates (Table 3). However, the mass loading rate curves for different construction methods become much more comparable for the single bench “structure” model, and almost identical for the “multi-bench” models (Table 4). This is an important finding; while spatial trends may influence drainage evolution through the mechanisms described above and may be dependent on the construction method (as observed in the “heterogeneous” cases), the addition of traffic surfaces and complex internal structures, which are common in WRPs (i.e., “structure” and “multi-bench” models) suggest that the choice of construction method (assuming oxygen is not limited) may in fact not have a significant effect on overall drainage water chemistry and evolution. Simulated results are most significantly affected by the inclusion of structure, related parameter heterogeneities, and scale. In the “multi-bench” model, peak sulfate mass loading rates decrease by almost a factor of 2 from the multi-bench homogeneous case, and a factor of 3 from the single-bench homogeneous case. This is a significant decline in overall maximum mass loading rates and has important implications for long term mine planning and environmental risk. On one hand, the homogeneous model arguably provides a higher factor of safety, since it predicts higher maximum mass loading rates. However, this scenario will fail to adequately predict long-term mass loading rates and duration and may result in higher initial capital costs for building water treatment facilities than needed (on top of additional factors of safety implemented through regulatory requirements). As mentioned previously, for the mass loading rates discussed, it is implied that basal drainage is entirely contained and perfectly mixed, which may not be the case in real world WRPs depending on effluent drainage, use of liners, or other pile design techniques. The implications of the investigated conceptual models' heterogeneity regarding monitoring and sampling techniques are discussed in the next section.

It should be acknowledged that the use of laboratory and/or field measurements is important in site-specific studies to validate model results. While modeled input parameters were constrained by values from measured parameters and their relationships reported in the literature, this study is conceptual in nature. The benefit of this numerical experiment is to provide insight to the role of heterogeneity, structure, and pile dimension in WRP models and to assess the anticipated effects of the parameters investigated here on mass loading rates and their evolution over time. Results from this numerical experiment demonstrate that these factors do play a role in affecting the simulation results and suggest that some investigation of parameter sensitivity may be worthwhile in future site-specific studies.

Heterogeneity and Sampling Procedures

It was previously noted that that there were no significant variations in overall sulfate mass loadings between the five different spatial distributions in models with random heterogeneity and models considering fining upwards trends. Material characterization such as mean and standard deviation of parameters used in model inputs will influence the numerical outputs, but it would be nearly impossible (and certainly not cost effective) to determine the exact spatial distributions within a WRP given current measurement and sampling techniques. However, examining the pH distributions throughout the piles in Figure 8 (as well as animations and additional effluent water quality curves provided in the Supplementary Material), the influence of heterogeneity on local drainage quality may have further implications for field sampling procedures. For example, seepage samples may be collected to monitor the drainage quality from a WRP following freshet or large rainfall events, but it is important to consider what these measurements represent. The consequences of this depend on the usage of these results: for example, is seepage chemistry used as a check for guidelines for water entering the receiving environment only, or is it being used as an input to validate models and make predictions? Is drainage assumed to be well-mixed using a liner below the base of the pile? Furthermore, if a waste rock pile is not lined and discharge seeps directly into the ground below, it becomes impossible to know what the real evolution of the drainage quality is. In other sampling methods such as the collection of pore water samples using soil water solution samplers installed within the pile, Figure 8 and Supplementary Figure 7 suggest that variation due to local acidic or neutral pockets may cause a skewed perception of drainage quality that will be released from the pile.

Conclusions

The purpose of this study was to investigate the influence of spatially varied distributions of physical and chemical properties throughout WRPs and their internal structures on ARD-generating processes. Stylistic models simulated using MIN3P-HPC suggest that heterogeneity and pile structure are important to the overall basal drainage evolution observed from a pile. Results from the single bench “heterogeneous” models suggest that push-dumping construction methods may be beneficial to reduce maximum solute concentrations entering the receiving environment over end-dumping methods—which caused a buildup of poor drainage quality at the top of the pile and resulted in a delayed spike of high sulfate concentrations from sulfide oxidation. However, this idea was challenged by the addition of structure (i.e., inclusion of zones, varied spatial trends, and traffic surface) and “multi-bench” cases, which reduced solute release up to 3-fold for both construction types when compared to the single-bench, homogeneous model, while indicating only small differences between construction methods on maximum sulfate effluent concentrations observed. These results were achieved assuming well-aerated conditions for all simulated piles (i.e., O₂ was not limited in pyrite oxidation). While oxygen availability is an essential control for pyrite oxidation rates, this study demonstrates that oxygen limitation is not the only explanation for lower mass loading rates with increasing scale, demonstrated by significant differences in the results presented here.

While capturing the sensitivities of model outputs to parameter heterogeneity, this study also aimed to demonstrate the capabilities of the recently developed MIN3P-HPC code to simulate large-scale mine waste dumps. With model upgrades introduced in this version, it possible to simulate complex models with irregular geometry within a reasonable time, particularly if utilizing the parallelized version of the code (in this study, most models were simulated for 300 years, requiring only a few hours to complete). Further advancements to MIN3P and other reactive transport codes along with technological progressions will improve abilities to capture heterogeneity, though will also come with challenges associated with the opportunity to add more information to RTMs.

As demonstrated in this study, modeled results may be very sensitive to the inclusion of heterogeneity and structure, and care must be taken regarding their effects on results when used for long-term closure planning. Furthermore, field data collected to constrain and develop the numerical models should be representative and should allow conceptualization of localized heterogeneities and internal structure. Further investigations for these future models and the influence of incorporating heterogeneity will need to continue to be updated.

This work further supports the use for numerical modeling as a conceptualization tool in understanding potential WRP behavior. The idea that WRPs are heterogeneous is not a new concept, but the most effective ways to capture the effects of heterogeneity on the overall pile evolution and ARD-related processes from field or lab-based experiments still have many uncertainties. Numerical models cannot provide a substitute for necessary field data, supplementary techniques, and lab-based experiments which have shown to support and constrain site-specific material behavior and results. However, RTMs can be used as an effective tool to conceptualize different waste rock storage strategies, which can be explored at an initial stage before further investment and narrow down the most effective long-term storage methods and develop more robust closure strategies. Finally, simulated model outputs can vary greatly depending on pile size, internal structures, and parameter heterogeneities considered, and are ultimately the result of the empirical formulas used in the code and the input parameters. It is the responsibility of the user to carefully consider and constrain models with added complexities such as heterogeneity and structure with solid, representative field measurements.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

Conceptualization of this study was completed in collaboration by KR and KM. Computer simulations were completed by KR with guidance and validation from DS, NS, and KM. The original draft was prepared by KR, with review, editing, and supervision from NS and KM. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the strategic network grant Toward Environmentally Responsible Resource Extraction Network (TERRE-NET) led by D. W. Blowes, grant number 479708-2015 and an Alexander Graham Bell Canada Graduate Scholarship-Master's (CGS M), award to KR.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Special thanks to the UBC EOAS IT Department (H. Modzelewski and C. Krzysik) for guidance with high performance computing and set-up and departmental computing resources. This research was completed as part of the requirements for KR's Master's thesis at the University of British Columbia.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2021.618418/full#supplementary-material

References

Akcil, A., and Koldas, S. (2006). Acid mine drainage (AMD): causes, treatment and case studies. J. Clean. Prod. 14, 1139–1145. doi: 10.1016/j.jclepro.2004.09.006