Numerical Modeling of a Laboratory-Scale Waste Rock Pile Featuring an Engineered Cover System

: Improved design to reduce contaminant mass loadings from waste rock piles is an increasingly important consideration. In certain cases, an engineered cover system containing a ﬂow control layer (FCL) may be used to mitigate the release of metals out of a pile using capillary barrier e ﬀ ects (CBEs), diverting water away from reactive materials below. In this study, a reactive transport model was calibrated to observational data from a laboratory experiment designed to evaluate a cover system. The results show that the numerical model is capable of capturing ﬂow rates out of multiple drainage ports and matching key e ﬄ uent concentrations by including the spatial distribution of hydraulic parameters and mineral weathering rates. Simulations were also useful for characterizing the internal ﬂow pathways within the laboratory experiment, showing the e ﬀ ectiveness of the cover in diverting the ﬂow away from the reactive waste rock and identifying secondary CBEs between di ﬀ erent rock types. The numerical model proved beneﬁcial in building an improved understanding of the processes controlling the metal release and conceptualizing the system. The model was expanded to investigate the robustness of the cover system as a function of the applied inﬁltration rate, supporting that water diversion will occur with inﬁltration rates representative of ﬁeld conditions.


Introduction
The long-term storage of waste rock-material considered not economically viable, which must be removed from the subsurface to access valuable resources during the mining process-is an increasingly important consideration for mine operations globally. When exposed to atmospheric conditions, sulfide minerals within waste rock oxidize and may cause the release of acid rock drainage (ARD), or if sufficient neutralizing minerals are present, contaminated neutral drainage (CND) [1][2][3][4][5][6]. To mitigate the release of metals and contaminants from ARD and CND, engineered design structures such as cover systems, subaqueous disposal, blending and layering, or co-disposal techniques may be applied as a means of prevention [7][8][9][10][11][12]. In addition to the laboratory-and field-scale tests required to assess these techniques, reactive transport modeling (RTM) may be used as a tool to investigate the processes involved in the release and attenuation of metals from waste rock piles through scenario analyses

Site Background
The Lac Tio Mine is an open-pit mine featuring a massive HI deposit hosted in anorthosite gangue material. The mine site has been in operation since 1950 and is located 43 km north of Havre-Saint-Pierre, Quebec, Canada. Waste rock type at this site varies between anorthosite and highly mineralized material with an HI content >70 wt. % [25]. Nickel, a primary focus of previous studies, comes primarily from Ni-bearing pyrite ((Ni x Fe 1-x )S 2 ), where x is spatially distributed and correlated to Ni abundance, varying between 0 and 10 wt. % throughout the HI deposit. Additionally, Ni is present in the form of millerite (NiS) embedded within pyrite grains [24]. The release of Ni from waste rock is retarded through sorption mechanisms in both plagioclase and ilmenite rock [33]. The drainage from Lac Tio waste rock is maintained at near neutral conditions due to the weathering of silicate minerals, specifically: calcic plagioclase which is mainly close to labradorite in composition (Na 0.4 Ca 0.6 Al 1.6 Si 2.4 O 8 ), enstatite (MgSiO 3 ), and biotite (K(Mg,Fe 2+ ) 3 (Al,Fe 3+ )Si 3 O 10 (OH,F) 2 ) [34]. Silicate dissolution to maintain near-neutral pH in the absence of carbonate minerals has been demonstrated previously [35,36], and is typically incongruent, preferentially releasing cations such as Ca and Na from plagioclase minerals, Mg from pyroxenes [37,38] and Mg, Fe, and K from sheet silicates [39]. Further detailed analyses of the weathering mechanisms and specific dissolution rates have been previously determined for Lac Tio waste rock [34].
In order to improve the drainage quality, an option under investigation at the mine consists of designing an FCL using crushed anorthosite and sand. The objectives of the FCL are to divert water from reactive materials using contrast in grain size to create capillary barrier effects (CBEs) and to reduce the mass loadings of contaminants out of the pile to minimize the water treatment requirements. Previous studies [14,40,41] found that a slope of 5 • was most effective for the FCL to divert precipitation while minimizing erosion. To investigate this design under controlled conditions, a heterogeneous laboratory-scale model was constructed [31] utilizing the relevant waste rock lithologies and compositions. Weekly flush experiments were performed on this laboratory model inclined at various angles between 0 • and 5 • over a period of 24 weeks between October 2017 and April 2018.
Information regarding the results of these transient tests can be found in previous works [31]. In this study, the system covered with the FCL inclined at a 5 • angle was subjected to a constant precipitation rate with leachate collected at different sampling ports along the base of the tank. This experimental study aimed to emulate the flow dynamics for the covered Lac Tio waste rock under high infiltration rates. The objective of the experiment was to quantitatively investigate the performance of the FCL and to help in quantifying hydrogeochemical reaction rates, which can potentially be used for future upscaling purposes.

Design of Experiment and Construction
The laboratory experiment [31] was constructed by depositing crushed waste rock material from the Lac Tio mine site using a hand-held scoop in 4 sections: 3 zones of HI material (labeled S1-S3) and 1 anorthosite zone at the end, all of which were deposited mimicking the angle of repose expected to develop during the deposition of waste rock at the field scale. Each zone was 80 cm in height and was emplaced in the 60 cm wide tank ( Figure 1). Overlying these layers, a 10 cm thick sand layer and 8 cm thick crushed anorthosite layer were placed to create the FCL. Anorthosite material was crushed to <10 mm for the FCL using a crusher at the Lac Tio mine site. The emplaced cover material was compacted using a Proctor hammer with a constant mass to approach target porosities and to create a CBE, with the goal to divert water from the HI material below. Grain size distribution curves for the materials in each of the sections can be found in Appendix A. Along the base, 8 drainage ports (1 cm diameter, referred to as reservoirs [1][2][3][4][5][6][7][8] were installed and equipped with tubing to guide the drainage water into collection buckets below ( Figure 1). Effluent volumes were determined for each reservoir; subsequently, the samples were collected for geochemical analysis. With the tank exposed to atmospheric conditions along the top, initial conditions throughout the experiment were assumed to be well aerated with respect to O 2 (21%) and CO 2 (0.0412%) [42]. Additional details on the laboratory model design and construction can be found in Poaty [31].

Flow Measurements
The experiment was performed to create steady-state flow conditions, applying precipitation with a sprinkler system uniformly across the surface at a constant rate over a 17-day period in February 2019 (photos of the model setup are included in Appendix A). A rate of 180 L day −1 (corresponding to an application rate of 108 mm day −1 ) was selected. This application rate represents the equivalent of one month of precipitation entering the pile under field conditions. The high-water application rate was selected in order to meet the experimental time constraints and to achieve steady state conditions quickly, while providing the benefit of testing the engineered cover system under extreme conditions. Drainage was collected in the buckets below each reservoir and weighed to determine the daily flow rates from each port. In addition, the frequency domain (FD) sensors placed throughout the experiment were used to measure the volumetric water content using dielectric constants [29]. While the FD probes were not calibrated specifically for this 17-day experiment, variations in the measured apparent dielectric constant values were used to support the bucket weight measurements and to determine when the steady-state flow conditions in the laboratory experiment were reached. While the flow rates obtained from the bucket weight measurements provided information on how much water was leaving at each reservoir, they could not be used alone to characterize the internal flow pathways within the experiment.

Geochemical Monitoring
Leachates were collected in buckets below reservoirs 1-8, in which pH and electrical conductivity were monitored daily. Once these parameters stabilized, water samples were collected from the buckets and filtered with 0.45 µm mesh filters into pre-acidified (with 2% HNO 3 ) 125 mL polyvinyl chloride (PVC) bottles for full-suite geochemical analyses using ICP-MS. Samples were collected and analyzed over a period of 10 days after the steady-state flow conditions were achieved, which was considered sufficient to represent the stabilized effluent water quality within the scope of this study. After the completion of transient flush experiments [31] and prior to the experiment completed for this study, the laboratory set-up was left inactive. During this period of inactivity, the materials were left to oxidize freely; therefore, the samples were not collected between days 1 through 7, since these samples were not considered representative for geochemically stable drainage conditions. Early time drainage was expected to represent the release of weathering products that had built up during several months of inactivity between transient flush experiments and the experiment conducted for this study. Results on the transient behaviour of the laboratory experiment have been published elsewhere [31].

Reactive Transport Code MIN3P-HPC
MIN3P-HPC [23] is built on the parallel version of MIN3P-THCm [43,44], with implementation of unstructured grid capabilities. The key features of the code include coupled multicomponent transport and biogeochemical reactions in variably saturated porous media [44][45][46][47]. The code uses Richards' equation to solve the flow in the unsaturated/saturated subsurface. A global implicit method and direct substitution approach was implemented for the solution of the multicomponent advection-dispersion equation and biogeochemical reactions [44,48]. Additional details about the code may be found in [44,45,48].
In this study, a 2D numerical model solving for the steady-state flow and transient geochemical conditions was set up to quantitatively interpret the data collected from the laboratory experiment presented above. The conceptual model of the experiment including the boundary conditions is depicted in Figure 2, while the mesh configuration used is provided in Appendix B. The versatility of the unstructured grid available in MIN3P-HPC allowed for the direct inclusion of the angled surfaces (such as the 5 • rotation of the FCL and sloped edges following the angle of repose in the HI zones, as depicted in Figure 2), which could not have been captured adequately using a standard rectangular-based mesh. Details regarding the calibration process for flow and geochemistry are provided in the following subsections.

Flow Model and Calibration
Previous works [27,29] utilized suction and moisture content measurements (collected using Watermark ® Granular Matrix sensors and ECH2O EC-5 soil moisture sensors, respectively) from transient flow experiments to better understand the internal flow regimes. In this study, the measured flow rates from each reservoir were used to calibrate the numerical model. Figure 2 illustrates the material property zones and boundary conditions used. Although a water application of 180 L day −1 evenly distributed across the surface was intended from the sprinkler system, the total outflow collected from the reservoirs only averaged 158 L day −1 under steady-state flow conditions. This difference in mass-in versus mass-out is attributed to a combination of factors including evaporation and difficulties in maintaining pump pressure and water delivery. A steady and uniform water infiltration (158 L day −1 , corresponding to an infiltration rate of 108 mm day −1 ), was applied in the model to match the total measured outflow under the steady-state conditions. Constant pressure heads were specified to calibrate the steady-state drainage rates at each reservoir. The simulated drainage rates proved to be sensitive to specified pressure heads at each of the sampling ports. Since no flow was measured out of reservoirs 6 and 7 during the experiment, these locations were also treated as no flow boundaries in the numerical model ( Figure 2). The variably saturated flow model was further constrained by measured physical parameters as summarized in Table 1. Porosity measurements were taken directly from the materials used in the laboratory experiment, while saturated hydraulic conductivity and unsaturated hydraulic soil parameters were obtained from previous column experiments and characterization of the waste rock and FCL [28,31,41,49]. Residual volumetric moisture contents reported in previous studies [41,50] were used to determine residual saturations for sand and waste rock materials. The van Genuchten-model was used to describe functional relationships between the pressure head, saturation and relative

Flow Model and Calibration
Previous works [27,29] utilized suction and moisture content measurements (collected using Watermark ® Granular Matrix sensors and ECH 2 O EC-5 soil moisture sensors, respectively) from transient flow experiments to better understand the internal flow regimes. In this study, the measured flow rates from each reservoir were used to calibrate the numerical model. Figure 2 illustrates the material property zones and boundary conditions used. Although a water application of 180 L day −1 evenly distributed across the surface was intended from the sprinkler system, the total outflow collected from the reservoirs only averaged 158 L day −1 under steady-state flow conditions. This difference in mass-in versus mass-out is attributed to a combination of factors including evaporation and difficulties in maintaining pump pressure and water delivery. A steady and uniform water infiltration (158 L day −1 , corresponding to an infiltration rate of 108 mm day −1 ), was applied in the model to match the total measured outflow under the steady-state conditions. Constant pressure heads were specified to calibrate the steady-state drainage rates at each reservoir. The simulated drainage rates proved to be sensitive to specified pressure heads at each of the sampling ports. Since no flow was measured out of reservoirs 6 and 7 during the experiment, these locations were also treated as no flow boundaries in the numerical model ( Figure 2). The variably saturated flow model was further constrained by measured physical parameters as summarized in Table 1. Porosity measurements were taken directly from the materials used in the laboratory experiment, while saturated hydraulic conductivity and unsaturated hydraulic soil parameters were obtained from previous column experiments and characterization of the waste rock and FCL [28,31,41,49]. Residual volumetric moisture contents reported in previous studies [41,50] were used to determine Minerals 2020, 10, 652 6 of 24 residual saturations for sand and waste rock materials. The van Genuchten-model was used to describe functional relationships between the pressure head, saturation and relative permeability [51], based on the parameters derived in previous studies [30,41,49,50]. For the HI material, zone-specific soil hydraulic function parameters were not available, and material-specific data from previous studies [30,50] were used uniformly for all zones. These parameters were used directly and were not varied during the model calibration.  [50] 0.04 [30] The values from the literature (rather than direct measurements such as porosity) were considered homogeneous within the individual material zones. In reality, the hydraulic conductivity and additional soil hydraulic parameters within these zones are likely heterogeneous; influenced by the natural heterogeneity of the material and potentially also affected by compaction during the construction of the laboratory model. However, since direct measurements of these parameters in the laboratory model were not available, the outflow was calibrated using pressure heads-rather than inferring the hydraulic conductivity distribution-was considered justified. In particular, this approach was chosen because the boundary conditions at the drainage ports are not well defined, since tubing to guide drainage was attached to the ports and since it was unclear whether all drainage ports performed equally. While this approach does not allow to infer the exact hydraulic conductivity distribution throughout the experimental vessel, the hydraulic parameters were adequately constrained by the related literature data and measurements, providing a suitable approach for investigating the diversion of water as a result of CBEs, which was the main purpose of this study.

Reactive Transport Model and Calibration for Leachate Concentrations
The reactive transport model considered a free-phase diffusion coefficient of D 0 = 2.4 × 10 −9 m 2 s −1 for all dissolved species, while longitudinal and transverse vertical dispersivities of 10 −3 m and 5 × 10 −4 m were used. In the gas phase, a free phase diffusion coefficient of 2.1 × 10 −5 m 2 s −1 was used for all gases. Effective diffusion coefficients were calculated based on the Millington-Quirk relationship [48]. Only the most relevant chemical components were considered in the model including SO 4 2− , Ni + , Fe 2+ , Fe 3+ , Ca 2+ , K + , Mg 2+ , Na + , Al 3+ , H 4 SiO 4 , CO 3 2− , H + and O 2 (aq). The measured composition of the tap water applied to the surface of the physical experiment was specified along the infiltration boundary in the numerical model as a specified flux boundary ( Table 2). Concentrations in Table 2 are provided as the total component concentrations, except for the components H + , CO 3 2− and O 2 (aq), which are specified in terms of pH, alkalinity and pO 2 , respectively. Aqueous speciation amongst these components is computed based on the hydrolysis and complexation reactions defined in the model database, which is derived from the USGS wateq4f database [52]. Chloride is present in the tap water as a result of chlorination, although it was not measured since Cl − has very limited reactivity and its omission has no significant effect on the simulation results. Although rainwater is more dilute than tap water, both tap water and rain water are substantially more dilute compared to pore water within the waste rock and drainage water; therefore, no significant impacts on the weathering reactions are anticipated by using tap water in the experiment. Boundary conditions in the reservoirs with active drainage were specified as free exit boundary conditions. While the mineralogy of the various materials is complex, the mineralogy considered for the numerical model has been simplified to include only the minerals which are considered to effectively control drainage chemistry. Table 3 summarizes the volume fractions of the primary mineral phases included. The majority of these volume fractions were taken from the previous characterization of waste rock by [53]; however, pyrite volume fractions were calculated based on the results from S analyses [31]. In addition to these primary phases, the secondary phases, ferrihydrite, Al(OH) 3 (am), SiO 2 (am) and kaolinite were allowed to precipitate. Table 4 summarizes the reaction pathways, rate expressions and relevant equilibrium constants of the primary and secondary minerals considered in the numerical model. A simple first order rate expression was used to describe the oxidation of the Ni-bearing pyrite. More complex models were not deemed necessary, in particular considering that the published rate expressions were determined under highly controlled laboratory conditions and are not necessarily directly applicable to field situations [54].   Table 4. Minerals and their dissolution-precipitation reactions considered in the numerical model. Effective rate constants (k eff,m ) were calibrated depending on material zones and are presented in the results section in Table 6.

Mineral Reaction Rate Expression logK m
Ni-bearing pyrite Nickel in the HI material is sourced from both pyrite and millerite grains (found in inclusions within the pyrite). However, the simulations considered pyrite as the only Ni source, with the dissolution rates calibrated using measured SO 4 release rates. Although previous work to characterize Lac Tio waste rock identified sorption as an important Ni-attenuation mechanism [25,33,34,53], this process was not explicitly included in the current numerical model. In this study, the mineral weathering rates in the numerical model were calibrated to match the release rates observed in the drainage of the physical model and constitute effective release rates relevant to the duration of the experiment, implicitly accounting for the effect of sorption. Accounting for sorption explicitly was beyond the scope of the present study.
Previous work has shown that both S and Ca release rates were stoichiometrically associated with sulfide (primarily in the form of pyrite, although trace amounts of chalcopyrite and millerite have also been detected) and plagioclase weathering, respectively, due to limited secondary mineral precipitation [34]. To reproduce the observed effluent chemistry of the system, rate coefficients for mineral dissolution-precipitation reactions ( Table 4) were calibrated using a stepwise approach. Initially, the oxidative dissolution rate of Ni-bearing pyrite was calibrated constrained by the measured concentrations of SO 4 and Ni from each reservoir. During this step, the Ni abundance in pyrite was also adjusted within the 0-10 wt. % range observed at the Lac Tio site. The presence of Ni impurities in pyrite has been noted to increase pyrite oxidation rates in previous studies [59,60]. Then, the most reactive silicates with pH-buffering capacity were considered (labradorite, enstatite, and biotite) to reproduce pH and major ion concentrations. Labradorite weathering rates were calibrated to reproduce Ca release rates and to achieve an initial calibration for pH. Enstatite and biotite were subsequently introduced to represent additional pH-buffering reactions by silicate minerals (using observed Mg and K release rates to calibrate, respectively). The secondary minerals ferrihydrite, SiO 2 (am), Al(OH) 3 (am) and kaolinite have been identified [34] and were included to control the Fe, Al, and Si concentrations in the drainage waters.

Reduced Infiltration Scenario Analysis
Expanding upon the model calibrated for the experiment under the steady-state flow conditions, the influence of reduced infiltration on the internal flow patterns and effectiveness of the FCL was investigated through additional simulations with a 10× lower infiltration rate (10.8 mm day −1 ). Since the outflow rates are unknown for these conditions, three different scenarios for a hypothetical water table were explored: (1) the same water table height as observed in the model for high infiltration rates; (2) a reduced water table height, which is expected for the reduced infiltration conditions; and (3) a scenario with no water table forming within the experimental vessel. This analysis was only performed to investigate the flow processes with a specific focus on the performance of the FCL; an analysis on the effect of drainage chemistry was not performed due to a lack of observational data.

Assumptions and Limitations
Building a conceptual model for RTM requires making assumptions and simplifications depending on what data are available and the overall purpose of the model. For example, Table 1 shows that porosity was determined for the various materials emplaced in the experimental vessel, while other input parameters were taken from other Lac Tio column and field experiments. In reality, the values for each material property zone may vary slightly compared to the parameters obtained from the companion experiments. Parameters may also vary spatially and show heterogeneities within each material property zone. In addition, the geochemical reaction network represents a simplification of the actual reactions occurring in the physical experiment and mineral reaction rates, which were assumed to be constant in time within each material property zone. This assumption is considered justified given the timescale of this experiment, and is also supported through previous work by [24] who determined through humidity cell tests that silicate weathering rates in fresh Lac Tio waste rock are very similar to weathering rates in the material that has undergone 25 years of natural weathering.
In addition, the numerical model was built in 2D, while the laboratory experiment contains features that deviate from a true 2D system. Specifically, the drainage ports in the laboratory experiment were located in the center of the base of the experiment in the y-direction. In the 2D numerical model, each reservoir had to be taken as a line of drainage rather than a 1 cm diameter central port.
Despite these simplifications, it must be emphasized that most model input data were available through measurements from the well constrained experiment and were kept constant during the model calibration process. The only parameters that varied during the modeling were the pressure heads at the drainage ports and the mineral weathering rates with the goal to reproduce drainage rates and concentrations from the various reservoirs under quasi-steady-state flow conditions.

Drainage Rates and Internal Fluid Flow Pathways
FD probe sensor measurements showed that stable volumetric water contents were reached within the laboratory experiment after 2 days, while the measured flow rates reached a quasi-steady state after approximately 8 days. Simulated results shown in Figure 3 and Table 5 capture the drainage rates with less than 1% difference at all points. The drainage rate out of reservoir 3 (approximately 70 L day −1 ) constitutes 44% of the total water volume leaving the model. However, based on examining the drainage rates alone, it is unclear whether a higher volume of water leaving reservoir 3 is caused by water leaking through the FCL into the HI material, or whether water is being diverted by the FCL, as intended, and then rerouted to reservoir 3. In Figure 4, the modeled results of Darcy flux vectors are presented by scaled arrows, which are used to conceptualize flow through the laboratory model. Simulations suggest that in the experiment, a near flat water table forms at z ≈ −2 cm, approximately corresponding to the elevation of reservoir 8. The pressure head boundaries used to calibrate the model (Figure 2) were positive for reservoirs 1-5, facilitating the outflow. The pressure head at reservoir 8 was slightly negative; however, outflow can still occur under these conditions since a drainage tube was attached to the outflow port at the base of the tank. In the HI zones in Figure 4, the smaller arrows show that there is less water flowing downwards compared to the anorthosite zone. These simulated flow pathways suggest that while there is some water leaking through the FCL into the HI material below, the FCL angled at 5 • effectively diverts a significant portion of the water under quasi-steady-state flow conditions. Furthermore, the model results suggest the development of secondary CBEs between the HI material and the uncrushed anorthosite layer (Figure 4). Based on these conditions, a nearly flat water table forms along the base of the model and causes some of the drainage water to move updip (to the left) from the uncrushed anorthosite layer back into the HI material, contributing to the larger flow rates observed at reservoirs 2 and 3. These results suggest that leakage through the FCL is limited and that the large outflows at reservoirs 3 are due to the internal rerouting of flow paths due to CBEs. Minerals 2020, 10, x FOR PEER REVIEW 10 of 24

Solute Concentrations in Drainage
Results of the measured and simulated water quality at each reservoir are shown in Figure 5. The observed solute concentrations in drainage suggest that geochemically stable conditions were achieved after 12 days. Therefore, average solute concentrations between day 12 and day 17 were used to calibrate the mineral weathering rates within the numerical model (Table 6). Table A1 in Appendix C provides statistical information regarding the achievement of stable geochemical conditions for the major elements of concern (S, Ca, K, Mg, Na and Si), justifying this approach. Observations between day 8 and day 12 show the elevated solute concentrations and decreased pH, interpreted as capturing the tail end of the flushing of weathering products originating from oxidation reactions that took place while the experimental setup was idle. Early time behavior in the model also shows variations in concentrations, due to the specified initial conditions. Simulated concentrations in drainage between day 0 and day 10 are affected by the displacement and mixing of resident pore water with infiltrating water and weathering reactions. These early time results are not relevant for the model calibration and interpretation, performed for geochemically stable conditions, and are therefore represented with dashed lines ( Figure 5).
The simulations included the diffusion of O2 through the gas phase. Results suggest that O2 was not depleted within the model domain over the duration of the experiment, and subsequently the oxidation of Ni-bearing pyrite was not inhibited. Pyrite oxidation rate coefficients were calibrated to reproduce the SO4-release observed in the experiment. Pyrite grains have been observed to typically contain between 0-10 wt. % Ni at the Lac Tio site [25], and a pyrite composition of Ni0.14Fe0.86S2 (corresponding to 7 wt. % Ni) was found to represent the observed Ni release rates most accurately, consistent with the observed Ni contents in drainage. The current approach utilizing effective Nirelease rates, implicitly accounting for sorption, was able to adequately reproduce the observed Nirelease, although previous work has shown that sorption processes are significant in fresh rock at the Lac Tio site and important for controlling the release of Ni [24,25]. It is possible that the effect of sorption in this experiment was less pronounced than in previous studies, since multiple pore

Solute Concentrations in Drainage
Results of the measured and simulated water quality at each reservoir are shown in Figure 5. The observed solute concentrations in drainage suggest that geochemically stable conditions were achieved after 12 days. Therefore, average solute concentrations between day 12 and day 17 were used to calibrate the mineral weathering rates within the numerical model (Table 6). Table A1 in Appendix C provides statistical information regarding the achievement of stable geochemical conditions for the major elements of concern (S, Ca, K, Mg, Na and Si), justifying this approach. Observations between day 8 and day 12 show the elevated solute concentrations and decreased pH, interpreted as capturing the tail end of the flushing of weathering products originating from oxidation reactions that took place while the experimental setup was idle. Early time behavior in the model also shows variations in concentrations, due to the specified initial conditions. Simulated concentrations in drainage between day 0 and day 10 are affected by the displacement and mixing of resident pore water with infiltrating water and weathering reactions. These early time results are not relevant for the model calibration and interpretation, performed for geochemically stable conditions, and are therefore represented with dashed lines (Figure 5).
The simulations included the diffusion of O 2 through the gas phase. Results suggest that O 2 was not depleted within the model domain over the duration of the experiment, and subsequently the oxidation of Ni-bearing pyrite was not inhibited. Pyrite oxidation rate coefficients were calibrated to reproduce the SO 4 -release observed in the experiment. Pyrite grains have been observed to typically contain between 0-10 wt. % Ni at the Lac Tio site [25], and a pyrite composition of Ni 0.14 Fe 0.86 S 2 (corresponding to 7 wt. % Ni) was found to represent the observed Ni release rates most accurately, consistent with the observed Ni contents in drainage. The current approach utilizing effective Ni-release rates, implicitly accounting for sorption, was able to adequately reproduce the observed Ni-release, although previous work has shown that sorption processes are significant in fresh rock at the Lac Tio site and important for controlling the release of Ni [24,25]. It is possible that the effect of sorption in this experiment was less pronounced than in previous studies, since multiple pore volumes had passed through the material during previous experiments [31], which may have contributed to saturating available sorption sites. In addition, high volumetric flow rates were applied to the laboratory model for this study, which may have reduced the sorption processes typically occurring in the Lac Tio waste rock due to rapid flushing. Another possibility is that some sorption-processes may have been accounted for implicitly by the calibrated Ni-content in pyrite used in the numerical model. As previously mentioned, Ni may range from 0-10 wt. % content in pyrite in Lac Tio waste rock [25], therefore it is possible that the HI material in the laboratory model is greater than the 7 wt. % used here. Additional laboratory testing may be able to deduce the actual Ni-release rates and the contribution of sorption towards effective Ni-release rates in this experiment; however, this is beyond the scope of the present study.
As mentioned in the methods section, labradorite, enstatite, and biotite were the primary minerals considered for pH buffering. The volume fractions of minerals presented in Table 3 were not depleted during the 17 days simulated by the numerical model and therefore did not affect the progress towards geochemically stable conditions. The rate coefficients for these phases were calibrated to match the release of major ions for each mineral (specifically: the pyrite was calibrated using SO 4 release; labradorite to Ca, biotite to K, and the remaining Mg to enstatite, while considering these influences on effluent pH). These rate coefficients for pyrite and silicates are summarized for each zone in Table 6.
While simulated concentrations in drainage waters reproduce the measured values for SO 4 , Ca, K, Mg, Na and pH fairly well for most reservoirs (Figure 5), the Si concentrations were above the measured values. The simulations suggested that the pore water was undersaturated with respect to SiO 2 (am) throughout the domain. Although the precipitation of kaolinite was considered in the simulations, its precipitation was not able to control Si to the observed concentrations. Previous studies [33,34] provide evidence for the presence of additional silicate minerals in the HI material such as pigeonite, muscovite, chlorite, orthoclase, and spinel. These minerals, with different Si contents and incongruent weathering processes, may also contribute to the silicate weathering and may at least in part explain the discrepancies between the simulated and measured results for Si. However, considering these additional mineral phases would result in increased model complexity and non-uniqueness; these phases were omitted from the analysis. In addition, simulated concentrations for Al and Fe were both below the measured values of approximately 3 × 10 −6 mol L −1 and 5 × 10 −6 mol L −1 , respectively. In the simulations, the precipitation of Al(OH) 3 (am) and ferrihydrite were considered as quasi-equilibrium reactions, while the speciation calculations for observed water chemistry indicate the supersaturated conditions with respect to these phases, implying that their precipitation is kinetically inhibited. Despite these simplifications, the model was able to reproduce the release of solutes of key interest (SO 4 and major cations) from the different reservoirs, meeting the study objectives ( Figure 5). Table 6. Calibrated effective rate coefficients for the mineral weathering reactions (units in (mol dm −3 bulk s −1 ), except for pyrite, which is provided in units of (L H 2 O dm −3 bulk s −1 ); these units account for the scaling from pore volume to bulk volume for the first-order rate oxidation rate).  Table 6 were calibrated through an iterative process to reproduce the observed concentrations in the effluents from the laboratory experiment. Figure 5. Results of the simulated versus observed effluent concentrations in mol L −1 , color coded per reservoir. Effective rate coefficients for the mineral dissolution-precipitation rates (k eff,m ) in Table 6 were calibrated through an iterative process to reproduce the observed concentrations in the effluents from the laboratory experiment.
Although high water application rates dilute the drainage chemistry, mass loading rates are expected to be similar to those under natural precipitation conditions. As mentioned previously, simulated oxygen levels did not decline throughout the domain for the conditions considered; implying that sulfide oxidation rates are not inhibited by high water application rates. However, some effect on the precipitation of secondary phases and metal attenuation is expected, since the formation of secondary minerals is more limited under high infiltration conditions. The impact is expected to most significantly affect trace metal release, which was not the focus of this study.

Reduced Infiltration Scenario Analysis
Simulation results of the flow model for a reduced infiltration rate, better aligned with real world conditions, are depicted in Figure 6. In this scenario, the infiltration rate was reduced by a factor of 10 in comparison to the conditions in the laboratory experiment. Since measured outflow rates were not available for these conditions, the boundary condition at the base of the experimental vessel had to be estimated. It was assumed that similar to the experiment with higher flow rates, a flat water table would develop, however, at a lower elevation. Simulation results indicate that the fluid flow pathways and functioning of the FCL remain very similar to the pattern observed for the higher flow rates (Figure 6), with the FCL effectively capturing most of the infiltration. Although the material zones show decreased water saturation above the water table compared to the simulations calibrated to the observational data from the laboratory experiment under higher flow conditions (compare Figure 6 to Figure 4), the simulated results suggest that water is still diverted through the sand material in the FCL ( Figure 6). Furthermore, a secondary CBE between the HI material and anorthosite is also present under these reduced infiltration scenarios. Simulation results show that the infiltration is diverted through the FCL towards the anorthosite material and that secondary capillary barrier effects (CBEs) between the anorthosite and HI material persist, similar to conditions in the laboratory experiment for higher infiltration rates. Darcy flux vectors are scaled the same as those shown in Figure 4, which are smaller in this case due to the reduced infiltration causing smaller fluxes. A closer look at the flow pathways for water movement from the sand into the anorthosite is included due to the reduced arrow sizes.

Implications of Results from Flow Modeling
Locally sourced sand in the FCL was used to achieve a desired thickness for the engineered cover system for the full-scale piles [30]. The numerical model confirmed the development of a capillary barrier between the sand and the HI material, which results in the majority of water being diverted through the sand layer towards the uncrushed anorthosite material (Figure 4). A related study [29] using a water truck to complete a wetting test on an experimental pile for the Lac Tio site to investigate the effectiveness of the FCL at the field-scale also found that fluid diversion occurred mainly through the sand layer, which is consistent with the results forecasted by the numerical model here, as well as the previously completed modeling [61]. These findings lend support to the FCL's ability to divert most water away from the reactive HI material.
Additionally, the results from all numerical models presented here indicate the development of secondary CBEs between the HI material and uncrushed anorthosite material, which had not been Figure 6. Simulated results of the flow model with a reduced infiltration rate and lowered water table (z ≈ −13 cm). Simulation results show that the infiltration is diverted through the FCL towards the anorthosite material and that secondary capillary barrier effects (CBEs) between the anorthosite and HI material persist, similar to conditions in the laboratory experiment for higher infiltration rates. Darcy flux vectors are scaled the same as those shown in Figure 4, which are smaller in this case due to the reduced infiltration causing smaller fluxes. A closer look at the flow pathways for water movement from the sand into the anorthosite is included due to the reduced arrow sizes.
To evaluate the sensitivity of the results with respect to the assumed basal boundary condition, two additional simulations were conducted with the water table located at the same elevation as seen in the laboratory experiments for high inflow rates and for the conditions without a water table within the experimental vessel (suction at all reservoirs) (for results see Appendix D). These results indicate that the functioning of the FCL and internal flow pattern are not substantially impacted by the water table location and the basal boundary condition, although the drainage volume at the various reservoirs differ in response to these changes.
The low infiltration simulation analysis was only performed to investigate the flow processes with a specific focus on the performance of the FCL; an analysis on the effect of drainage chemistry was not performed due to a lack of observational data.

Implications of Results from Flow Modeling
Locally sourced sand in the FCL was used to achieve a desired thickness for the engineered cover system for the full-scale piles [30]. The numerical model confirmed the development of a capillary barrier between the sand and the HI material, which results in the majority of water being diverted through the sand layer towards the uncrushed anorthosite material (Figure 4). A related study [29] using a water truck to complete a wetting test on an experimental pile for the Lac Tio site to investigate the effectiveness of the FCL at the field-scale also found that fluid diversion occurred mainly through the sand layer, which is consistent with the results forecasted by the numerical model here, as well as the previously completed modeling [61]. These findings lend support to the FCL's ability to divert most water away from the reactive HI material.
Additionally, the results from all numerical models presented here indicate the development of secondary CBEs between the HI material and uncrushed anorthosite material, which had not been previously identified. The existence of these secondary CBEs is dependent on the unsaturated soil hydraulic parameters, which were constrained by the experimental data and field observations ( Table 1). The secondary CBEs are beneficial in diverting water away from the metal-rich HI material and retaining it within the uncrushed anorthosite zone, and therefore may warrant additional future investigation. The presence of secondary CBEs is also supported by the geochemical results (discussed in greater detail below).
For this study, the daily infiltration rates in the experiment to calibrate the numerical model were equivalent to one month of precipitation at the field site. In reality, water saturations within the full-scale pile would be lower and infiltration processes would be transient, possibly leading to differences in flow dynamics. The effect of a reduced infiltration rate was considered through an additional scenario with an infiltration rate of 10.8 mm day −1 (a reduction by a factor of 10). Model results suggest that fluid flow pathways controlled by the FCL and CBEs are very similar for both infiltration rate scenarios, indicating that the FCL will also be functioning under field conditions. Numerical modeling facilitated to account for (1) heterogeneities in the physical flow parameters between different material zones and (2) the pile geometry, owing to the unstructured grid capabilities of the MIN3P-HPC code. These factors were particularly important to investigate the influence of CBEs in the engineered cover system, since in this case the process of fluid diversion is controlled by contrast in grain size, associated hydraulic parameters and pile geometry.

Implications of Results from Reactive Transport Modeling
Both the simulated and observed effluent concentrations presented in Figure 5 support the simulated internal fluid flow pathways depicted in Figure 4. It was noted above that the highest outflow rates (both measured and simulated) were observed at reservoir 3. If the majority of this water was leaking past the FCL and percolating through the HI material zones (contrary to the simulated flow model), the effluent quality would be expected to lie closer to the observations in reservoirs 4 and 5. However, the actual measured concentrations from reservoir 3 fall somewhere in between reservoirs 4-5 and 1-2, supporting the hypothesis that only a portion of the water released at reservoir 3 moves through the HI material, while a substantial fraction percolates through the anorthosite zone and is then redirected updip towards reservoir 3. Furthermore, reservoirs 1 and 2 have very similar drainage concentrations with improved water quality relative to the other reservoirs in both the modeled and measured results. While this is expected for reservoir 1 which lies below the anorthosite material, reservoir 2 is located under the HI material. This finding further supports the results of the flow model, suggesting that water draining at reservoir 2 primarily percolates through the anorthosite material, but then moves updip and to the left towards reservoir 2 ( Figure 4).
Additionally, the Darcy flux vectors shown in Figure 4 are smaller in magnitude in HI S1 towards reservoir 8 compared to the simulated flow towards reservoirs 4 and 5. Longer residence times and contact between infiltrating water and HI material are likely responsible for the increased concentrations observed in the drainage collected from reservoir 8. However, the simulations suggest that physical flow alone was insufficient to capture the differences in effluent concentrations observed between reservoirs 1 and 8. In Table 6, the differences in the rate coefficients for each material property zone are summarized. For pyrite, higher oxidative dissolution rates were required in HI S1 compared to S2 and S3 in order to reproduce the SO 4 concentrations observed at each reservoir. Consequently, the rate coefficients were also higher in zone S1 for labradorite and enstatite, although since the rate expressions for these minerals (as well as biotite) were pH dependent, the overall dissolution rates were also higher due to the greater acidity produced from faster sulfide oxidation rates. Table 1 shows that the measured porosities from the laboratory model varied quite significantly between HI S1, S2, and S3 with values of 0.39, 0.55, and 0.58, respectively. These variations in porosity suggest heterogeneity (or potentially differences in compaction) between each zone, which may have influenced the reactive surface area and the related weathering rates. Porosity values were lower in zones from left to right following the deposition history. It is possible that materials in S1 underwent greater compaction than the other HI material zones due to the location, influencing amount (due to more overall rock), reactivity or availability of sulfide grains in this zone, although it is speculative whether this would have a greater influence over the naturally occurring chemical heterogeneities between zones. Furthermore, the grain size distribution curves (Appendix A) show that S1 displayed a smaller grain size compared to S2 and S3, which may explain the faster weathering rates due to the increased surface area and sulfide abundance [2,60,[62][63][64][65].
Overall, the simulated effluent quality reproduced the measured effluent chemistries from each reservoir quite well. Similar to the physical heterogeneity being essential for capturing the influences of CBEs on flow patterns, differences in mineral weathering rates between the zones were essential for reproducing the variations in chemical loadings collected in drainage from the different reservoirs. Results suggest that for infiltration that is not diverted from the FCL, the degree of compaction and grain size may influence the residence time and the resulting effluent quality released from these zones. Simulated results would not have been able to reproduce the observational data with a homogeneous model, pointing towards the importance of considering physical and chemical heterogeneity in the numerical models going forward. In addition, the results indicate that compaction and grain size should be considered and assessed with respect to their influence on effluent quality.

Conclusions
RTM can be a useful tool for the quantitative interpretation of experimental data, identifying fluid migration pathways, bracketing reaction kinetics, or forecasting expected solute release rates and mass loadings. However, with many experiments aimed at explaining complex hydrogeochemical interactions, a thoughtful approach to understanding the dominant processes in the system and justifiable assumptions and limitations must be considered. In this study, RTM was useful to quantitatively assess the functioning of an engineered FCL, in the process revealing the likely development of beneficial secondary CBEs and providing insights into internal flow migration patterns. Although simplifying assumptions had to be made, the resulting model was able to reproduce the observational laboratory data including the spatially distributed drainage rates and solute concentrations, while keeping most material properties constant based on the measured parameters. The simulations were also useful to estimate the mineral weathering rates for pyrite and the aluminosilicate minerals present in the fresh HI material used in the laboratory experiment. For the simulations conducted, the consideration of physical and chemical heterogeneities between different material zones was essential to reproduce the observations made in the laboratory. The model was also useful to evaluate the effectiveness of the FCL for the conditions of reduced infiltration rates, suggesting that the FCL will function as designed under infiltration rates more representative of field conditions. Finally, the results demonstrate the capabilities of MIN3P-HPC as a tool for simulating the flow and reactive transport processes in complex hydrogeochemical systems such as waste rock piles. The recent implementation of unstructured grid capabilities into the model allows the user to incorporate more complex features, such as the influence of angled layers like the FCL or material zones following the angle of repose. Future advancements in computing power and improvements made to complex multi-phase interactions within RTM codes will further improve the capabilities of the models to characterize mine waste stockpiles and may be useful to contrast and improve different techniques for the design planning and storage of mine wastes.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. model. HI rock is the dark grey material on the left which was deposited in three segregated zones; anorthosite is the light grey material (to the right of the HI); sand is orange-brown in colour along the top with crushed anorthosite (light grey) above to form the FCL. The tank was pivoted at an angle of 5° for the duration of the experiments with steady-state flow rates. Top right: Drainage tubes 1 cm in diameter that allow effluent to freely drain into the collection buckets below. Bottom: Set up of tubing with drainage holes installed along the top of the lab model to distribute the water applied at 180 L day −1 . This unstructured grid allows for the influence of a 5° FCL slope to be investigated, as well as the effects of internal sloped boundaries following the angle of repose.  Figure A3. Triangular mesh with zone outlines used for the numerical model. The mesh file used in numerical simulations consisting of 16,602 cells sloped at a 5 • angle. The cells have an average representative elemental volume of 1.4 cm 3 . This unstructured grid allows for the influence of a 5 • FCL slope to be investigated, as well as the effects of internal sloped boundaries following the angle of repose.