Science, Faculty of
Non UBC
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
Minerals 11 (6): 565 (2021)
Yi, Xueying
Su, Danyang
Bussière, Bruno
Mayer, K. Ulrich
In order to reduce contaminant mass loadings, thermal cover systems may be incorporated in the design of waste rock piles located in regions of continuous permafrost. In this study, reactive transport modeling was used to improve the understanding of coupled thermo-hydrological and chemical processes controlling the evolution of a covered waste rock pile located in Northern Canada. Material properties from previous field and laboratory tests were incorporated into the model to constrain the simulations. Good agreement between simulated and observational temperature data indicates that the model is capable of capturing the coupled thermo-hydrological processes occurring within the pile. Simulations were also useful for forecasting the pile’s long-term evolution with an emphasis on water flow and heat transport mechanisms, but also including geochemical weathering processes and sulfate mass loadings as an indicator for the release of contaminated drainage. An uncertainty analysis was carried out to address different scenarios of the cover’s performance as a function of the applied infiltration rate, accounting for the impacts of evaporation, runoff, and snow ablation. The model results indicate that the cover performance is insensitive to the magnitude of recharge rates, except for limited changes of the flow regime in the shallow active layer. The model was expanded by performing an additional sensitivity analysis to assess the role of cover thicknesses. The simulated results reveal that a cover design with an appropriate thickness can effectively minimize mass loadings in drainage by maintaining the active layer completely within the cover.
mineralsArticleThermal-Hydrological-Chemical Modeling of a Covered WasteRock Pile in a Permafrost RegionXueying Yi 1,*, Danyang Su 1, Bruno Bussière 2 and K. Ulrich Mayer 1Citation: Yi, X.; Su, D.; Bussière, B.;Mayer, K.U. Thermal-Hydrological-Chemical Modeling of a CoveredWaste Rock Pile in a PermafrostRegion. Minerals 2021, 11, 565.https://doi.org/10.3390/min11060565Academic Editor: Juan AnteloReceived: 7 February 2021Accepted: 22 May 2021Published: 26 May 2021Publisher’s Note: MDPI stays neutralwith regard to jurisdictional claims inpublished maps and institutional affil-iations.Copyright: © 2021 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms andconditions of the Creative CommonsAttribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).1 Department of Earth, Ocean, and Atmospheric Sciences, The University of British Columbia, 2207 Main Mall,Vancouver, BC V6T 1Z4, Canada; dsu@eoas.ubc.ca (D.S.); umayer@eos.ubc.ca (K.U.M.)2 Institut de Recherche en Mines et en Environnement, Université du Québec en Abitibi-Témiscamingue,445 Boulevard de l’Université, Rouyn-Noranda, QC J9X 5E4, Canada; Bruno.Bussiere@uqat.ca* Correspondence: xyi@eoas.ubc.caAbstract: In order to reduce contaminant mass loadings, thermal cover systems may be incorporatedin the design of waste rock piles located in regions of continuous permafrost. In this study, reactivetransport modeling was used to improve the understanding of coupled thermo-hydrological andchemical processes controlling the evolution of a covered waste rock pile located in Northern Canada.Material properties from previous field and laboratory tests were incorporated into the model toconstrain the simulations. Good agreement between simulated and observational temperaturedata indicates that the model is capable of capturing the coupled thermo-hydrological processesoccurring within the pile. Simulations were also useful for forecasting the pile’s long-term evolutionwith an emphasis on water flow and heat transport mechanisms, but also including geochemicalweathering processes and sulfate mass loadings as an indicator for the release of contaminateddrainage. An uncertainty analysis was carried out to address different scenarios of the cover’sperformance as a function of the applied infiltration rate, accounting for the impacts of evaporation,runoff, and snow ablation. The model results indicate that the cover performance is insensitive tothe magnitude of recharge rates, except for limited changes of the flow regime in the shallow activelayer. The model was expanded by performing an additional sensitivity analysis to assess the role ofcover thicknesses. The simulated results reveal that a cover design with an appropriate thickness caneffectively minimize mass loadings in drainage by maintaining the active layer completely withinthe cover.Keywords: reactive transport modeling; mine waste rock; thermal cover; freeze-thaw cycles; permafrost;contaminated drainage1. IntroductionLarge volumes of waste rock (considered non-economically viable materials) are gen-erated and removed from the subsurface during mining operations. This waste rock (whichis usually stored in large and unsaturated stockpiles at ground surface) poses long-termenvironmental risks related to acid rock drainage (ARD) or contaminated neutral drainage(CND), if sufficient pH-buffering minerals are present in the waste rock, such as carbonateminerals. The weathering of potentially acid-generating (PAG) materials, in particular sul-fide minerals, may trigger the release of ARD/CND and dissolved metals [1–3]. Differenttechniques have been proposed to control the release of contaminants and limit adverseimpacts of ARD/CND, such as the placement of covers (including water covers, coverswith capillary barrier effects, elevated water table with monolayer covers, low saturatedhydraulic conductivity covers, store-and-release covers, inorganic or organic coatings, elec-trochemical covers, and thermal covers) [4–8]. Thermal covers (also known as insulationcovers) are a unique option in cold regions where permafrost is typically found due toprevailing cold temperatures, implying that soils remain perennially frozen for at leasttwo consecutive years [5,9–12]. In this scenario, non-potentially acid-generating (NPAG)Minerals 2021, 11, 565. https://doi.org/10.3390/min11060565 https://www.mdpi.com/journal/mineralsMinerals 2021, 11, 565 2 of 26waste rock also extracted from the site can be used as the cover material to encapsulatethe PAG material. The role of the cover is to provide an insulating effect to maintainpermafrost in the underlying waste rock, while the active layer, defined as the top layersubject to seasonal thawing, remains within the cover material. The low temperaturesmaintained in the reactive waste rock control the oxidation of sulfide minerals, whereasfreezing conditions reduce the potential for transport of contaminated water [7]. Therefore,the design for cover thickness will depend on site-specific conditions (e.g., permafrostextent, climate, soil properties) [13,14].In addition to laboratory and field tests, regarded as essential methods to understandwaste rock behavior, reactive transport modeling (RTM) can also provide insights intothe coupled thermo-hydrological-chemical (THC) processes occurring within a wasterock pile (WRP), responsible for controlling conditions in the field. Another practicalapplication of modeling is to assist with improving the design of reclamation methods formine wastes to minimize their environmental impacts. Increasing attention has been givento the capabilities of RTM to address various mine waste related issues [15–23]. Whileseveral codes have been developed to simulate dynamic freezing and thawing processes inpermafrost regions [24–27], only a few studies have applied this functionality to miningenvironments. In this paper, the reactive transport code MIN3P-HPC [28,29] was usedto assess the performance of thermal covers placed on WRPs, based on the example ofthe Meadowbank mine, located in the Kivialliq region in Nunavut, Canada. At this site,a thermal cover system of 4-m thickness was emplaced to avoid potential ARD/CND [30].Laboratory experiments were performed to characterize the Meadowbank materials, andfield measurements were taken to improve understanding of the overall thermal behaviorof the WRP. However, a quantitative two-dimensional analysis of the coupling of THCprocesses has not been conducted to date. The main objectives of this study are to: (1) assessthe long-term THC behavior of a sulfide-bearing WRP in a cold-climate region featuringa thermal cover, (2) assess the sensitivity of the modeling results to uncertainties withrespect to the magnitude of seasonal recharge, and (3) evaluate the functionality andeffectiveness of thermal cover designs on controlling the weathering of waste rock andARD/CND release.2. Conceptual ModelThe simplified conceptual model presented here was developed to capture the cou-pling of the most important processes with respect to freeze-thaw cycling in sulfidic minewaste rock encapsulated by a thermal cover. Under cold-region climate conditions, large-scale WRPs usually allow for the development of permafrost within the waste rock and therelocation of the active layer from the underlying soil to near the pile surface, as shownin Figure 1 [31]. Permafrost established within the waste rock can be conceptualized as anearly impermeable barrier, limiting deep recharge, and in this way altering internal waterrouting and enhancing lateral flow [31]. Unlike in temperate climates, where precipitationis often dominated by rainfall and the duration of ground freezing is relatively short andrestricted to the shallow ground, infiltration into WRPs is more complex in cold regions,due to extended periods of frozen ground during the winter months and the presenceof permafrost at depth. During the winter months, snow can either accumulate at orablate from the pile surface. The accumulation and ablation rates are determined by avariety of parameters such as topographic features, vegetation, climate, and wind [10].In cold regions, spring snowmelt is a particularly significant event for water infiltrationinto piles [32] and has important influences on solute transport and moisture content.Minerals 2021, 11, 565 3 of 26Minerals 2021, 11, x FOR PEER REVIEW 3 of 26 solutes. Oxygen consumption is associated with the oxidation of PAG minerals present in a WRP. However, the cover is designed to minimize the effects of freeze-thaw cycling within the WRP, restricting oxidation processes, oxygen consumption. Therefore, also ox-ygen ingress. As a result, the oxygen content within the pile is expected to remain rela-tively stable at atmospheric levels. Figure 1. Simplified conceptual model of a waste rock pile subjected to a cold-region climate. In addition, heat transfer processes occurring within a WRP are tightly associated with above-mentioned processes. A temperature gradient is present between the frozen core and the seasonally thawed surface layer, which can drive heat transport through con-duction in or out of the pile. The magnitude of the conductive energy flux is controlled by temperature gradients and thermal properties of both NPAG materials in the cover and PAG materials in the core [3]. Another heat transfer mechanism is convection in associa-tion with pore water and pore gas. The limitation on deep infiltration of water due to the presence of permafrost is accompanied by limited convective heat flux by water flow. The oxidation of sulfide minerals is known as strongly exothermic, potentially causing internal heat production. However, the mostly frozen state of the WRP may slow down oxidation processes and associated heat production. In the case of low internal sulfide oxidation rates, convection driven by pore gas is also restricted. 3. Site Description In this study, the Meadowbank Gold Mine is used as an example for quantitatively evaluating cover performance in cold-region climates. The mine is located approximately 300 km from Hudson Bay in the Kivialliq region (65 °N 96 °W), which is subject to typical low Arctic eco-climatic conditions with a mean annual air temperature (MAAT) of −11 °C [35]. The mine site has been in operation by Agnico Eagle Mines Limited since 2010, gen-erating large quantities of both PAG and NPAG waste rock. Previous baseline studies show that the thicknesses of the active layer and the underlying permafrost in the region typically vary from 1.3 to 4 m and 450 to 550 m, respectively [35,36]. Based on these ob-servations, it has been decided to construct a 4-m insulation cover WRP made of NPAG waste rock to control contaminant generation from the WRP. To validate the cover design, two field experimental cells with a cover of 2 m and 4 m were constructed over 1–2 m tailings to aid with closure planning for the reactive tailing storage facilities (TSFs) [12]. Figure 1. Simplified conceptual model of a waste rock pile subjected to a cold-region climate.Cold-region climates also influence chemical reactions and gas transport within thepile. The establishment of permafrost as well as the freezing of porewater can result in asignificant reduction in chemical reaction rates, effectively leading to c ssation of reactionprogress [33,34]. Moreover, the limited thickness of the active lay r and short period ofthawed conditions also limit che ical reactions as well as the release and transport ofsolutes. Oxygen consumption is associated with the oxidation of PAG minerals present in aWRP. However, the cover is designed to minimize the effects of freeze-thaw cycling withinthe WRP, restricting oxidation processes, oxygen consumption. Therefore, also oxygeningress. As a result, the oxygen content within the pile is expected to remain relativelystable at atmospheric levels.In addition, heat transfer processes occurring within a WRP are tightly associatedwith above-mentioned processes. A temperature gradient is present between the frozencore and the season lly thawed surfa e layer, w ich can drive hea t ansport throughconduction in or out of the pile. The magnitude of the conductive energy flux is controlledby temperature gradients and thermal properties of both NPAG materials in the cover andPAG materials in the core [3]. Another heat transfer mechanism is convection in associationwith pore water and pore gas. The limitation on deep infiltration of water due to thepresence of permafrost is accompanied by limited convective heat flux by water flow. Theoxidation of sulfide minerals is known as strongly exothermic, potentially causing internalheat production. However, the mostly frozen state of the WRP may slow down oxidati nprocesses and associated hea produc ion. In the case of low intern l sulfide oxidationrates, convection driven by pore gas is also restricted.3. Site DescriptionIn this study, the Meadowbank Gold Mine is used as an example for quantitativelyevaluating cover performance in cold-region climates. The mine is located approximately300 km from Hudson Bay in the Kivialliq region (65◦ N 96◦ W), which is subject to typ-ical low Arctic eco-climatic conditions with a mean annual air temperature (MAAT) of−11 ◦C [35]. The mine site has been in operation by Agnico Eagle Mines Limited since2010, generating large quantities of both PAG and NPAG waste r ck. Previous baselinestudies show that the thicknesses of the active layer and the underlying permafrost in theMinerals 2021, 11, 565 4 of 26region typically vary from 1.3 to 4 m and 450 to 550 m, respectively [35,36]. Based on theseobservations, it has been decided to construct a 4-m insulation cover WRP made of NPAGwaste rock to control contaminant generation from the WRP. To validate the cover design,two field experimental cells with a cover of 2 m and 4 m were constructed over 1–2 mtailings to aid with closure planning for the reactive tailing storage facilities (TSFs) [12].Since the test results from monitoring the temperatures showed that 4m is effective forthe TSFs, it was assumed that this cover configuration would also work for the WRP. Asshown in Figure 2, the waste rock pile contains five benches, which are covered with 4mNPAG waste rock, resulting in a large pile (118 m in height × 650 m × 1020 m), located ona layer of till. The batters of the pile are sloped at 27◦ (2H:1V) [37].Minerals 2021, 11, x FOR PEER REVIEW 4 of 26 Since the test results from monitoring the temperatures showed that 4m is effective for the TSFs, it was assumed that this cover configuration would also work for the WRP. As shown in Figure 2, the waste rock pile contains five benches, which are covered with 4m NPAG waste rock, resulting in a large pile (118 m in height × 650 m × 1020 m), located on a layer of till. The batters of the pile are sloped at 27° (2H:1V) [37]. A total of fourteen vertically instrumented multipoint thermistor strings were in-stalled into borehole casings within the WRP in 2013 (YSI 44031 from GKM Consultants, Sainte-Julie, QC, Canada) and 2015 (G176-TH002 from Hoskin Scientific Ltd., Burnaby, BC, Canada), to monitor the temperature and thermal behavior of the entire covered sys-tem. Thermistor strings were installed inside borehole casings surrounded by gravel re-jects of ¾ inches. All thermistor sensors had an interchangeable tolerance of ± 0.1 °C. Fig-ure 2 shows the locations of the thermistors distributed across the pile. Detailed infor-mation about the thermistors, including the thermistors’ elevations, depths, string lengths, and the total number of sensors, are available in [37]. Strings that span through all pile benches include RSF-3, RSF-4, RSF-5, and RSF-6, containing more than 26 sensors. Based on the available sensors, cross section AA’ (Figure 2) was selected for the numerical analysis performed in this work. Except for RSF-4, which was installed on a higher bench (and was not measured from 2015 to 2017 due to its loss of functionality), strings were installed on the second or third benches. Daily air temperature and precipitation were also recorded on site at a meteorological station 1km away from the WRP [37]. Additional la-boratory tests were completed at the Université du Québec en Abitibi-Témiscamingue (UQAT, Rouyn, QC, Canada) to characterize the hydrogeological, thermal, physical, and chemical properties of both NPAG and PAG mat rials. Detailed method descriptions and results these la oratory tests are available in [35,37–44]. Figure 2. Plan view of the waste rock pile at the Meadowbank mine site. 4. Methods 4.1. Modeling Approach 4.1.1. Reactive Transport Code MIN3P-HPC MIN3P-HPC [28,29] (a parallel version of the MIN3P-THCm code) is a multicompo-nent reactive transport code designed for simulating the primary thermo-hydro-chemical and mechanical (THCm) processes in variably saturated porous media with unstructured grid capabilities. Key features of the code include density dependent flow and solute Fig re 2. lan view of the waste rock pile at the Meadowbank mi e site.A total of fourteen vertically instrumented multipoint thermistor strings were in-stalled into borehole casings within the WRP in 2013 (YSI 44031 from GKM Consultants,Sainte-Julie, QC, Can da) and 2015 (G176-TH002 from Hoskin Scientific Ltd., Burnaby, BC,Canada), to monitor the temperature and thermal behavior of the entire covered system.Thermistor strings were installed inside borehole casings surrounded by gravel rejects of34 inches. All thermistor sensors had an interchangeable tolerance of ±0.1 ◦C. Figure 2shows the locations of the thermistors distributed across the pile. Detailed informationabout the thermistors, including the thermistors’ elevations, e ths, string lengths, and thetotal num er of ensors, are available in [37]. Strings that span through all pile benchesinclude RSF-3, RSF-4, RSF-5, and RSF-6, containing more than 26 sensors. Based on theavailable sensors, cross section AA’ (Figure 2) was selected for the numerical analysisperformed in this work. Except for RSF-4, which was installed on a higher bench (and wasnot measured from 2015 to 2017 due to its loss of functionality), strings were installed onthe second or third benches. Daily air temperature and precipitation were also recorded onsite at a meteorological station 1km away from the WRP [37]. Additional laboratory testswere completed at the Université du Québec en Abitibi-Témiscamingue (UQAT, Rouyn,QC, Canada) to characterize the hydrogeological, thermal, physical, and chemical proper-ties of both NPAG and PAG materials. Detailed method descriptions and results of theselaboratory tests are available in [35,37–44].Minerals 2021, 11, 565 5 of 264. Methods4.1. Modeling Approach4.1.1. Reactive Transport Code MIN3P-HPCMIN3P-HPC [28,29] (a parallel version of the MIN3P-THCm code) is a multicompo-nent reactive transport code designed for simulating the primary thermo-hydro-chemicaland mechanical (THCm) processes in variably saturated porous media with unstructuredgrid capabilities. Key features of the code include density dependent flow and solutetransport, gas transport, biogeochemical reactions, heat transport, and one-dimensionalhydromechanical coupling [20,28,29,45,46]. In addition, the code has recently been en-hanced to simulate freeze-thaw related processes. The cessation of flow due to freezing isconsidered through lowering the hydraulic conductivity as a surrogate for phase transfor-mation between ice and water phases. Previous studies have shown that the temperaturedependence of reaction rates is adequately described by the Arrhenius equation, even at orslightly below the freezing point [33,47,48]. More recent work has shown that, for sub-zerotemperatures, sulfide mineral reactivity is negligible, implying that the Arrhenius equationscannot represent sulfide oxidation rates for completely frozen conditions [12]. The effectsof freezing on sulfide oxidation can therefore be simplified by defining a threshold temper-ature (site-specific) below which the reaction progress ceases [12]. The Arrhenius equationwas also modified for sulfide mineral oxidation in high temperature environments, toaccount for the decline of microbial activity under these conditions [49]. Furthermore,the energy balance equation in the code was enhanced by including heat generation frompyrite oxidation as a source term. The governing equations for flow, reactive transport,and heat transport, along with details regarding abovementioned relationships can befound in [50]. For this study, the code was used to simulate a base case scenario withspecific focus on coupled thermal-hydrological and geochemical evolution of a coveredWRP, constrained by data from field observations and laboratory experiments, followedby uncertainty analyses. Details regarding the parameterization and assigned boundaryconditions are provided in the following sections.4.1.2. Spatial DiscretizationThe model domain was constrained by the dimensions of the covered waste rockpile at the Meadowbank mine. The 2D-model was set up following cross section A-A’(Figure 2) to facilitate quantitative interpretation of field observations, accounting for pilegeometry. A deep and wide section of natural terrain (991 m × 500 m) underlying the WRPhas been included in the model. The vertical dimension was extended to the base of thecontinuous permafrost where the geothermal heat flux could be reasonably estimated basedon literature [51]. The domain was also extended laterally in order to avoid interferenceof lateral heat flux with the boundaries. Figure 3 shows the model configuration andspatial discretization, along with the locations of the thermistors, which provided the datato constrain the model. The enhanced unstructured grid capabilities of the MIN3P-HPCcode [29] allow for inclusion of multiple benches and sloped surfaces, as well as gridrefinements in areas of interest to simulate the covered pile. The domain was discretizedby using 22,455 nodes. Discretization of the cover and waste rock is more refined than thatof the underlying soil, since the focus on the study is on processes occurring in the WRP.Minerals 2021, 11, 565 6 of 26Minerals 2021, 11, x FOR PEER REVIEW 5 of 26 transport, gas transport, biogeochemical reactions, heat transport, and one-dimensional hydromechanical coupling [20,28,29,45,46]. In addition, the code has recently been en-hanced to simulate freeze-thaw related processes. The cessation of flow due to freezing is considered through lowering the hydraulic conductivity as a surrogate for phase trans-formation between ice and water phases. Previous studies have shown that the tempera-ture dependence of reaction rates is adequately described by the Arrhenius equation, even at or slightly below the freezing point [33,47,48]. More recent work has shown that, for sub-zero temperatures, sulfide mineral reactivity is negligible, implying that the Arrhe-nius equations cannot represent sulfide oxidation rates for completely frozen conditions [12]. The effects of freezing on sulfide oxidation can therefore be simplified by defining a threshold temperature (site-specific) below which the reaction progress ceases [12]. The Arrhenius equation was also modified for sulfide mineral oxidation in high temperature environments, to account for the decline of microbial activity under these conditions [49]. Furthermore, the energy balance equation in the code was enhanced by including heat generation from pyrite oxidation as a source term. The governing equations for flow, re-active transport, and heat transport, along with details regarding abovementioned rela-tionships can be found in [50]. For this study, the code was used to simulate a base case scenario with specific focus on coupled thermal-hydrological and geochemical evolution of a covered WRP, constrained by data from field observations and laboratory experi-ments, followed by uncertainty analyses. Details regarding the parameterization and as-signed boundary conditions are provided in the following sections. 4.1.2. Spatial Discretization The model domain was constrained by the dimensions of the covered waste rock pile at the Meadowbank mine. The 2D-model was set up following cross section A-A’ (Figure 2) to facilitate quantitative interpretation of field observations, accounting for pile geom-etry. A deep and wide section of natural terrain (991 m × 500 m) underlying the WRP has been included in the model. The vertical dimension was extended to the base of the con-tinuous permafrost where the geothermal heat flux could be reasonably estimated based on literature [51]. The domain was also extended laterally in order to avoid interference of lateral heat flux with the boundaries. Figure 3 shows the model configuration and spa-tial discretization, along with the locations of the thermistors, which provided the data to constrain the model. The enhanced unstructured grid capabilities of the MIN3P-HPC code [29] allow for inclusion of multiple benches and sloped surfaces, as well as grid refine-ments in areas of interest to simulate the covered pile. The domain was discretized by using 22,455 nodes. Discretization of the cover and waste rock is more refined than that of the underlying soil, since the focus on the study is on processes occurring in the WRP. Figure 3. Model configuration of 2D cross-section and grid discretization. Figure 3. Model configuration of 2D cross-section and grid discretization.4.1.3. Base Case Specific Input ParametersThe base case was designed to investigate the short- and long-term performance ofthe covered WRP, for comparison with temperature data collected at the site. The modelwas constrained by measured physical, hydrogeological, thermal, and chemical param-eters of both NPAG and PAG materials, summarized in Table 1. The bulk densities,porosities, and saturated hydraulic conductivities were directly derived from laboratoryexperiments [37–39], while thermal conductivities [38] were obtained from a generalizedthermal model reported in [41,43]. The hydraulic conductivity under frozen conditionswas drastically reduced to a value low enough to represent the cessation of flow underpermafrost conditions. The van Genuchten model was used to parameterize the relation-ships between volumetric water content and suction [37] (Figure A1, Appendix A). Theparameters used for the underlying natural soil were obtained from [50].Mineralogical composition in the models were justified based on previous stud-ies [37,38] (Table A1, Appendix A). With specific focus on the generation of acid rockdrainage and generation of heat, only pyrite and carbonates (expressed as equivalentcalcite) were included in the model as primary mineral phases to represent acidificationpotential (AP) and neutralization potential (NP). At this site, pyrite oxidation dominatesthe contribution of mineral dissolution precipitation reactions towards the energy balancedue its strong exothermic nature. The contribution of dolomite towards pH-bufferingreactions was combined with calcite based on total carbonate content to capture the fullNP. The volumetric abundances of both primary minerals were constrained by field ob-servations [37,38] (Table A1, Appendix A). Although abundant, silicate minerals wereneglected, since they are not expected to make significant contributions to pH bufferingand the energy budget of the pile. Secondary mineral phases ferrihydrite and gypsumwere allowed to precipitate to provide solubility controls for iron, carbonate, and sulfate.Therefore, the aqueous components and species considered in the simulations were re-stricted to those involved in these mineral reactions, including SO42−, Fe2+, Fe3+, Ca2+,CO32−, H+, and O2(aq). The initial porewater compositions in both materials were assumedsimilar to dilute rainwater with concentrations of 10−10 mol·L−1 for SO4, Fe(II/III), andCa, and in equilibrium with O2 and CO2 at atmospheric concentrations. Table 2 givesinformation about the reaction stoichiometries, rate expressions, equilibrium constants,and rate constants of included primary and secondary phases. Reaction rates for calcite,ferrihydrite and gypsum reflect quasi-equilibrium conditions, while the reference rateconstant for pyrite oxidation is within the range of values reported in the literature [50].Minerals 2021, 11, 565 7 of 26Previous laboratory tests performed on Meadowbank tailings suggest negligible reactivityof the tailings below 0 ◦C [12]. In line with these findings, it was assumed that mineral reac-tion rates are negligible in the waste rock, if temperatures are below freezing, while abovefreezing, the rate constant for pyrite oxidation was varied with temperature as a functionof the Arrhenius equation and equilibrium constants were temperature-compensated bythe van’t Hoff equation.Table 1. Hydrogeological, thermal and selected geochemical model parameters.Properties Symbols & Units NPAG PAG NaturalTerrainThickness/Height L/m 4 m [37] 118 m [37] 500 m [50]Porous mediumResidual saturation Sr/unitless 0.08 1 [37] 0.1 1 [37] 0.1 [50]Van Genuchten α factor α/m−1 0.22 1 [37] 0.64 1 [37] 4.6 [50]Van Genuchten n factor n/unitless 1.65 1 [37] 1.70 1 [37] 1.35 [50]Saturated hydraulic conductivityat reference temperature (25 ◦C) Ksat/m s−1 7.9 × 10−4 [37] 1.6 × 10−3 [37] 7.4×10−5 [50]Porosity φ/unitless 0.39 [37] 0.31 [37] 0.30 [50]Free-phase aqueous diffusion coefficient D0a /m2s−1 2.1 × 10−9 [15]Heat transportThermal conductivity of water 2 λw/W·m−1·◦C−1 0.6Thermal conductivity of solid λs/W·m−1·◦C−1 4.8 [35] 5.3 [35] 2.0 [50]Heat capacity of water 2 cw/J·kg−1·◦C−1 4180Heat capacity of air ca/J·kg−1·◦C−1 1010Heat capacity of solid cw/J·kg−1·◦C−1 1400 [37] 1600 [37] 1093 [50]Bulk density of solid phase ρs/kg·m–3 1763 [37] 2015 [37] 1750 [50]Threshold temperature for pyrite oxidation Tf/◦C 0 [50]Pyrite oxidationVolume fraction of pyrite Vpyrite /unitless - 0.002 [35] -Standard exothermic heat ∆E r /kJ·mol−1pyrite- 1409 [15] -Arrhenius activation energy ∆Ea/kJ·mol−1pyrite- 64 -Volume fraction of calcite Vcalcite /unitless 0.19 [35] 0.06 [35] -1 Estimated porous-medium parameters based on laboratory data report in [37]. 2 Thermal conductivity and heat capacity of water havebeen applied to both water and ice.Table 2. Primary and secondary minerals with their rate expressions, stoichiometries, equilibrium constants (Ki), andeffective rate constants (keffi ) at 25◦C.Mineral Reaction Rate Expression logKilogkeffi[mol·dm−3·s−1]Primary mineral reactionsPyrite FeS2 + 3.5O2 + H2O→ Fe2+ + SO2−4 + 2H+ R = keffFeS2 - −9.6 [50]Calcite CaCO3 → Ca2+ + CO2−3 R = keffCaCO3[1− IAP10−8.48]−8.48 [52] −7.3 [53]Secondary mineral reactionsgypsum CaSO4·2 H2O→ Ca2+ + SO2−4 + 2H2O R = keffFe(OH)3[1− IAP104.89]4.89 [52] −7.0 [53]ferrihydrite Fe(OH)3 + 3H+ → Fe3+ + 3H2O R = keffCaSO4 H2O[1− IAP10−4.58]−4.58 [52] −7.0 [53]Infiltration into the pile was simulated as seasonally controlled discontinuous events,based on cumulative precipitation data shown in Table 3. Details of the recorded dailyprecipitation rates are available in [37]. In the simulations, recharge was assumed toonly occur when the air temperature was above freezing. Specifically, a constant rate ofMinerals 2021, 11, 565 8 of 26infiltration was applied for average daily air temperatures above freezing, assuming avalue of 0.32 mm day−1 for the base case, derived from cumulative rainfall between 2014and 2017, accounting for evaporation and runoff, which was assumed to comprise 70% ofrainfall. The percentage was selected from the midpoint of the evaporation ranges between50% and 90% in other mining studies located in similar climates [54,55]. In addition,recharge due to snowmelt was considered during the first month after spring thaw. Duringthis period, infiltration from melting snow was set at 0.55 mm day−1, assuming that 50%of snowmelt infiltrates and that 25% of cumulative snowfall remains at the surface whenspring arrives. The percentage of snowmelt represents the contribution of snowmelt tothe infiltration that was derived from the midpoint of the range (10–90%) reported in theliterature [32], and the percentage of accumulated snowfall was selected from the midpointof the model calibration conducted for the Meadowbank mine, which determined that5–45% of snowfall forms the snowpack at the surface at the onset of spring thaw [37].Table 3. Yearly cumulative precipitation (rainfall and snowfall) with corresponding event durations between 2014 and 2017 [37].YearRainfall Snowfall 1Event Duration(Days)CumulativePrecipitation (mm)Event Duration(Days)CumulativePrecipitation (mm)2014 138 116 227 1182015 115 127 250 1782016 122 166 243 1352017 141 137 224 113Average rate (mm day−1) 1.06 0.581 Snowfall reported as rainfall equivalents.The exposure of the covered pile to atmospheric conditions necessitated the incor-poration of seasonal temperature fluctuations. The average daily temperatures (Ta) fromJanuary 2014 to October 2017 are shown in Figure 4, and a sinusoidal function dependenton time (t) captures the seasonal variations adequately [37]:Ta(◦C) = −11.4 ◦C + 23.2 ◦C sin((2pi365t)+ 4.2)(1)Minerals 2021, 11, x FOR PEER REVIEW 8 of 26 The temperature fluctuates around the mean annual air temperature (MAAT = −11.4 °C), and 23.2 °C represents the annual a plitude of temperature variations. The geother-mal heat flux at the model base was ssumed to be 0.087 W m−2 [51], yielding the perma-frost thickness observed in this region. All other boundaries were considered impermea-ble for both water flow and heat flux (Figure A2, Appendix B). The base case simulation was run for 50 years starting from 2014. Figure 4. Observed average daily air temperatures (Avg) for the period from 2014 to 2017 represented by a sinusoidal function (Sin fct). Table 3. Yearly cumulative precipitation (rainfall and snowfall) with corresponding event durations between 2014 and 2017 [37]. Year Rainfall Snowfall 1 Event Duration (days) Cumulative Precipitation (mm) Event Duration (days) Cumulative Precipitation (mm) 2014 138 116 227 118 2015 115 127 250 178 2016 122 166 243 135 2017 141 137 224 113 Average rate (mm day−1) 1.06 0.58 1 Snowfall reported as rainfall equivalents. 4.2. Infiltration Uncertainty Analysis Although infiltration rates used for the base case simulation separately consider con-tributions from rainfall and snowmelt, account for snow ablation, evaporation and runoff, as well as seasonal controls, uncertainty remains due to a lack of site-specific infiltration measurements into the pile. Expanding on the base case study, additional simulations for both increased and reduced infiltration rates were performed to investigate the sensitivity on pile-internal flow patterns, moisture contents, as well as thermal and geochemical evo-lution. The infiltration rates in the uncertainty analysis used the upper limits of the litera-ture ranges for the increased infiltration scenario and the lower limits for the reduced in-filtration scenario, as summarized in Table 4. The model geometry and other parameter inputs remain the same as for the base case. Table 4. Infiltration uncertainty analysis. Case # Rainfall Snowfall Precipitation Infiltration Fraction (% of Rainfall) Rain Infiltration (mm day−1) Snowpack Fraction (% of Snowfall) Infiltration Frac-tion (% of Snowpack) Snowmelt Infiltration (mm day−1) Infiltration 1 (mm day−1) Infiltration 2 (mm day−1) Reduced 10 0.11 5 10 0.02 0.13 0.11 Base 30 0.32 25 50 0.55 0.87 0.32 Increased 50 0.53 45 90 1.78 2.31 0.53 Figure 4. Observed average daily air temperatures (Avg) for the period from 2014 to 2017 represented by a sinusoidalfunction (Sin fct).The temperature fluctuates around the mean annual air temperature (MAAT =−11.4 ◦C),and 23.2 ◦C represents the annual amplitude of temperature variations. The geothermalhe t flux at the model base was ssumed to be 0.087 W m−2 [51], yielding the permafrostthickness observed in this region. All other boundaries were considered impermeable forMinerals 2021, 11, 565 9 of 26both water flow and heat flux (Figure A2, Appendix B). The base case simulation was runfor 50 years starting from 2014.4.2. Infiltration Uncertainty AnalysisAlthough infiltration rates used for the base case simulation separately considercontributions from rainfall and snowmelt, account for snow ablation, evaporation andrunoff, as well as seasonal controls, uncertainty remains due to a lack of site-specificinfiltration measurements into the pile. Expanding on the base case study, additionalsimulations for both increased and reduced infiltration rates were performed to investigatethe sensitivity on pile-internal flow patterns, moisture contents, as well as thermal andgeochemical evolution. The infiltration rates in the uncertainty analysis used the upperlimits of the literature ranges for the increased infiltration scenario and the lower limits forthe reduced infiltration scenario, as summarized in Table 4. The model geometry and otherparameter inputs remain the same as for the base case.Table 4. Infiltration uncertainty analysis.Case # Rainfall Snowfall PrecipitationInfiltrationFraction(% of Rainfall)RainInfiltration(mm day−1)SnowpackFraction(% of Snowfall)Infiltration Fraction(% of Snowpack)SnowmeltInfiltration(mm day−1)Infiltration 1(mm day−1)Infiltration 2(mm day−1)Reduced 10 0.11 5 10 0.02 0.13 0.11Base 30 0.32 25 50 0.55 0.87 0.32Increased 50 0.53 45 90 1.78 2.31 0.531 Infiltration rate during first month of spring thaw. 2 Infiltration rate during other months with above-freezing average daily air temperatures.4.3. Cover Thickness Scenario AnalysisExpanding on previous studies on cover thickness [12], an additional three simulationswere performed with no cover, a 2 m cover, and a 6 m cover. The purpose of thesesimulations was to evaluate and illustrate the effect of the presence and thickness of thecover system on sulfide oxidation, long-term drainage evolution, and mass loadings. TheWRP configuration remained the same as for the base case, except for the presence orthickness of the cover. The material properties of both cover and pile remained unchangedrelative to the base case and all boundary conditions also remained as specified for thebase case.4.4. Assumption and SimplificationsSome simplifications and assumptions had to be made when building the conceptualand numerical models. Infiltration rates were based on seasonal averages rather thanindividual daily rates, and no infiltration was considered when the surficial ground wasfrozen. The magnitude of the infiltration rate was estimated by removing assumed frac-tions of evaporation, runoff, and snow ablation from precipitation data, according to thesuggested ranges reported in the literature. Furthermore, the temperature and infiltrationpatterns assigned to the pile surface were assumed to be repetitive over the 50-year simu-lation period, while all other boundary conditions remain unchanged. Consideration ofthe potential meteorological changes in the long term was beyond the scope of the study.Other input parameters derived from related experiments and the literature were assumedhomogeneous for the various materials. The consideration of the effect of heterogeneitieswas also beyond the scope of the study, and preferential flow was neglected. In addition tothe simplified mineralogy simulated, pyrite oxidation was simplified without considera-tion of dynamic evolution of oxidation processes as a function of porewater geochemistry(complex parallel reactions or multiple oxidation pathways), since geochemical data wereunavailable and geochemical conditions are expected to remain relatively stable due to thepresence of carbonate minerals and the assumption of abundant oxygen [3,15]. A simplezero-order reaction expression was used for pyrite oxidation, since the main focus was ontemperature dependence. In order to highlight and isolate the contribution towards sulfateMinerals 2021, 11, 565 10 of 26mass loadings from the weathering of PAG waste rock, it was also assumed that limitedsulfur content present in the cover material is negligible. In addition, it was assumed thatthe pile is well aerated and that O2 concentrations remain at atmospheric levels throughoutthe whole pile (pO2 was fixed at 0.2 atm), in line with previous observations at WRPslocated in permafrost environments (pore gas including O2 level were found to closelyresemble those of air) [56,57]. This implies that the simulations neglect gas phase diffusion,gas advection, and gas supply through soil gas. The heat convection with pore gases istherefore not included in the simulations. For the conditions at this site, neglecting thisprocess is not expected to have any significant implications on the energy balance, since theheat capacity of air is much lower than that of the waste rock and water. The contributionof water–gas phase transitioning on the energy balance was considered in the model, butthat of ice–water phase transitioning was not included explicitly due to its irreversiblenature which alternates seasonally in the active layer where flow and reactive transportprocesses occur, and its relatively low latent heat of melting when compared to that ofvaporization. The effects of hydrodynamic dispersion and freezing point depression dueto salinity were also neglected.5. Results and Discussion5.1. Base Case ResultsComparing simulated temperature fluctuations from January 2014 to October 2017at RSF-4, RSF-5, RSF-6, RSF-15, and RSF-16 to the measured temperature data at the siteindicates the simulations provide an adequate representation of the thermal regime in thecovered pile (Figure 5). Simulated temperature values from June to November in eachyear show good agreement to the measured values at all stations, especially at the stationswith the most complete data sets (i.e., RSF-4, RSF-5, and RSF-6). Although simulatedtemperatures were lower than the measured values at RSF-15 and RSF-16, the simulationresults adequately capture the two peaks in 2016 and 2017. Overall, simulated resultscapture the penetration of the seasonally induced thermal waves, with a penetration depthlimited to the upper 10–15 m of the WRP. Simulations show that seasonal temperaturefluctuations at greater depth are much lower and lag significantly behind temperaturechanges near the surface due to the effects of heat conduction, consistent with observedpatterns. Both measured and simulated results show that the core of the WRP remainsfrozen and nearly constant in extent, even in the warm seasons for the first 4 years, withonly the cover thawing in summer and early fall. Winter temperature data was limited andonly available over the four-year period at RSF-6, and over the first year at RSF-4 and RSF-5.These field measurements show a multiyear warming trend from 2015 to 2017, most clearlyseen for the winter season, ranging from about 3 to 7 ◦C as a function of depth (as shown atRSF-6), while the simulations yielded identical minimum temperatures over the four-yearperiod. These results indicate the covered pile in the field was warmer than predictedby the model, especially for the shallow depths in the winter months (1.0–2.2 m). Thesedifferences are at least in part caused by insulating effects of the surface snowpack [9],which is only seasonally present and was not considered in the model. Other differences,also visible in the warm seasons, were likely to result from (1) neglecting heterogeneityof material properties in the numerical simulations (such as thermal conductivity), whichinfluence heat transport; (2) wind effects at the batters of the pile (i.e., RSF-15 and RSF-16),leading to gas advection and convective heat transport into or out of the pile; (3) variationsin cover thickness over the pile surface of the constructed WRP; (4) preferential flowthrough large pores in the pile; and (5) the slope orientation. Specifically, the surfacetemperatures along south-facing slopes are expected to be higher than on north-facingslopes. The model fit increases later during the simulations since early results are stronglyaffected by the initial condition relevant to temperature evolution during the processof pile construction, which is difficult to take into consideration. When simulating thepile behavior and implement the conceptual model based on measured data, no attemptwas made to calibrate the model parameters to improve the fit at individual locations toMinerals 2021, 11, 565 11 of 26observed temperature data. Discrepancies between observed and simulated data exist,especially during early time, limiting the predictive capabilities of the model. However,seasonal trends and the dampening of the temperature amplitude with depth are wellrepresented by the model, providing confidence in the modeling approach for waste rockweathering and cover functionality in cold-region climates.Minerals 2021, 11, x FOR PEER REVIEW 11 of 26 Figure 5. Simulated temperatures compared to the observed temperatures at RSF-6, RSF-4, RSF-5, RSF-16, RSF-15 from January 2014 and October 2017. The long-term thermal evolution of the covered pile is shown in Figure 6. Since the seasonally controlled temperature function is assigned for the entire 50-year simulation period, the core temperature of the WRP approaches the average yearly value (i.e., MAAT Figure 5. Simulated temperatures compared to the observed t mperatures at RSF-6, RSF-4, RSF-5,RSF-16, RSF-15 from January 2014 and October 2017.Minerals 2021, 11, 565 12 of 26The long-term thermal evolution of the covered pile is shown in Figure 6. Sincethe seasonally controlled temperature function is assigned for the entire 50-year simula-tion period, the core temperature of the WRP approaches the average yearly value (i.e.,MAAT = −11.4 ◦C) (due to the main controls of climatic condition and geothermal heatflux), surrounded by temperature gradients (due to freeze thaw processes in the activelayer) towards the surface. The simulations indicate that the 0 ◦C isotherm remains in thesurficial 4m NPAG material (above the PAG waste rock) in the long term, implying that theactive layer is maintained within the cover.Minerals 2021, 11, x FOR PEER REVIEW 12 of 26 = −11.4 °C) (due to the main controls of climatic condition a ge thermal heat flux), sur-rounded by temperature gradients (due to freez thaw processes in the active layer) to-wards the surfac . The simulations indicate that th 0 °C isotherm remains in the surficial 4m NPAG material ( bove the PAG waste rock) in the long term, implying that the ctiv layer is maintained within the cover. Figure 6. Simulated contours of temperature, ice saturation, water saturation, pH, dissolved sulfate and calcium in the summer after 50 years (i.e., 2064). The white line and white arrows indicate the interface between waste rock and cover materials. The contours were cropped to more clearly illustrate the differences at the basal drainage location. Temperature fluctuations and seasonal thawing in the active layer allow for the per-colation of water. In Figure 6, the observed ice-free layer near the surface, present in the summer, is another indicator of the active layer extent. The distribution of ice and water Figure 6. Simulated contours of temperature, ice saturation, water saturation, pH, dissolved sulfate and calcium in thesummer after 50 years (i.e., 2064). The white line and white arrows indicate the interface between waste rock and covermaterials. The contours were cropped to more clearly illustrate the differences at the basal drainage location.Minerals 2021, 11, 565 13 of 26Temperature fluctuations and seasonal thawing in the active layer allow for thepercolation of water. In Figure 6, the observed ice-free layer near the surface, present inthe summer, is another indicator of the active layer extent. The distribution of ice andwater in the waste rock pile, and in particular the presence of ice within the cover materialand above the waste rock, confirms that the active layer is restricted to the cover anddoes not extend into the waste rock. The simulations indicate that some limited pondingof water occurs within the cover above an ice layer (especially on the benches near thetoes of the batters), resulting from the blocking effect of permafrost to deeper infiltration,facilitating the accumulation of infiltrating water within the cover, followed by lateraldiversion towards the toes of the pile during the warm seasons. The underlying lowpermeability till blocks infiltration into natural soils at the toe of the WRP, which explainsthe increased water-saturation zones at the base of the slopes.In the absence of geochemical data on pore water chemistry, the model is used toillustrate the expected long-term geochemical evolution of the system (Figure 6). Sulfateand calcium concentrations are good indicators for the occurrence of pyrite oxidationand calcite dissolution, respectively. Very low concentrations of sulfate observed in thesystem indicates a lack of pyrite oxidation and acidity production, except for spatiallyconstrained zones within the steep slopes (Figure A3, Appendix C). These results illustratethe functionality of the cover, showing that sulfide oxidation is very limited, since thePAG waste rock will remain frozen in the longer term. At the top of each tipping face(where the pile bench meets the sloped face, referred to as a “break”), multi-directionalheat conduction enhances the thaw depth (i.e., down from the bench and from the batters)leading to some sulfate release. However, the simulations indicate that basal drainageremains free of sulfate for the entire simulation period. Calcite (which was considered as theprimary pH-buffering mineral) dissolves in the cover mostly in response to slightly acidicinfiltration from rain and snow, establishing equilibrium under atmospheric conditionsand leading to slightly basic pH-values in the near surface pore water (Figure 6).5.2. Uncertainty Analysis to InfiltrationSimulated ice and water saturations for increased and reduced infiltration rates aredepicted in Figure 7 for the summer after 50 years in comparison to the base case results.An ice-free layer of similar thickness is observed near the surface in all infiltration scenarios,indicating that the active layer is maintained within the cover and the waste rock remainsfrozen in the long term, regardless of the infiltration rate. The distributions of watersaturations are directly related to the flow regime and infiltration rates. In the increasedinfiltration scenario, water accumulation within the cover is more pronounced, depictedby the higher average water saturations and more pronounced internal water pondingin Figure 7. The occurrence of the water ponding is related due to the blocking effectof permafrost. Ponding also becomes more obvious at the breaks of each batter, wherethe direction of flow instantly changes to lateral flow. As for the base case scenario, thetoes of the piles become water-saturated due to underlying low-permeability materials,leading to outflow during the summer months. Darcy flux vectors, presented by scaledarrows, are used to visualize flow patterns for the different scenarios. In the scenario witha reduced infiltration rate, flux vectors suggest that water is also laterally diverted throughthe cover during the summer months, although it moves more slowly than for the base case.Conversely, flow is faster in the scenario with an increased infiltration, while flow patternsremain similar. The simulations indicate that, in the long term, basal drainage at the toes ofthe pile is to be expected, with the magnitude of outflow a function of the recharge.Minerals 2021, 11, 565 14 of 26Minerals 2021, 11, x FOR PEER REVIEW 14 of 26 be affected materially, although drainage water volumes will be a function of the infiltra-tion rates in the long term. Figure 7. Simulated contours of ice and water saturations with reduced or increased infiltration rates compared to the base case in the summer after 50 years. The grey arrows represent the Darcy flux vectors. The contours were cropped to more clearly illustrate the differences at the basal drainage location. The white line and white arrows indicate the interface between waste rock and cover materials. 5.3. Sensitivity Analysis to Cover Thickness The simulated results of the internal temperature and ice saturation distributions for the four cover-thickness scenarios are shown in Figure 8. Simulated active layer thickness is similar for the different scenarios in spite of differences of variations of heat capacities between NPAG and PAG materials. This result is not unexpected, considering that heat capacities of the two materials differ by less than 15% (Table 1). However, for a thinner cover, the simulations indicate that the active layer extends into the PAG material. The 0 °C isotherm in the scenarios with no cover and 2 m cover are located within the PAG material, while this isotherm remains near, but above the interface between the cover and the PAG waste rock for the scenarios with 4 m and 6 m covers. In other words, the results suggest that PAG material near the pile surface will be affected by seasonal thaw in the long-term for the scenarios with a 2 m cover or no cover (Figure A4, Appendix C). Alt-hough the cover thickness has little impact on the active layer thickness, it determines whether the PAG rock will experience thaw and allow for the oxidation of sulfides solute transport of reaction products through this material in the warm seasons. The distribution of ice and water saturations are reflective of the flow regime. All four scenarios show a similar active layer extent, indicated by the thickness of the ice-free layer near the surface, but their distributions are different with respect to the cover thickness. Although the waste rock remains mostly frozen for the no-cover case, a shallow active layer forms in the summer, with increased water saturations and pronounced internal ponding, result-ing from deeper infiltration due to the relatively high hydraulic conductivity of the PAG waste rock in combination with the blocking effect of permafrost. For the cases including a cover, the cover with its lower hydraulic conductivity adds more complexity to the ice and water distribution. Zones of water and ice saturations are both observed within the cover material in the 4 m cover and 6 m cover scenarios. The ice layer formed above the interface between the waste rock and cover is thicker in the 6 m cover scenario, whereas Low infiltrationMedium infiltrationHigh infiltrationice saturation [-]0 0.2 0.4 0.6 0.8 1.0water saturation [-]0 0.2 0.4 0.6 0.8 1.0flux vector0 25mAfter 50 yearsFig re 7. Simulated contours of ice and water saturations with reduced or increased infiltration rates compared to thebase case in the summer after 50 years. The grey arrows represent the Darcy flux vectors. The contours were cropped tomore clearly illustrate the differences at the basal drainage location. The white line and white arrows indicate the interfacebetween waste rock and cover materials.The sensitivity analysis for infiltration was designed to assess how uncertaintiesrelated to nfiltration will affect flow processes and th thermal regime with n the p le, withan emphasis on cover performance. Simulatio res lts illustrate that the applied re hargehas an insignificant impact on the thickn ss of th active laye , since flow patterns andlocatio s, as well as the ice distribution , are similar for all infiltrati scenarios. Due toa lack of observati nal dat , furth r analy is of the effects of the recharge n dra agewater c emistry was not performed. However, since the imulation indicate th t activelayer t ickness is not significantly affected, it can be p stul t d th t mass loadings willalso not be affected materially, although drainage water volumes will be a functi n of theinfiltration rates in the long term.5.3. Sensitivity Analysis to Cover ThicknessThe simulated results of the internal temperature and ice saturation distributions forthe four cover-thickness scenarios are sho n in Figure 8. Simulated active layer thicknessis similar for the different scenarios in spite of differences of variations of heat capacitiesbetween NPAG and PAG materials. This result is not unexpected, considering that heatcapacities of the two materials differ by less than 15% (Table 1). However, for a thinnercover, the simulations indicate that the active layer extends into the PAG material. The0 ◦C isotherm in the scenarios with no cover and 2 m cover are located within the PAGmaterial, while this isotherm remains near, but above the interface between the coverand the PAG waste rock for the scenarios with 4 m and 6 m covers. In other words, theresults suggest that PAG material near the pile surface will be affected by seasonal thawin the long-term for the scenarios with a 2 m cover or no cover (Figure A4, Appendix C).Although the cover thickness has little impact on the active layer thickness, it determineswhether the PAG rock will experience thaw and allow for the oxidation of sulfides solutetransport of reaction products through this material in the warm seasons. The distributionof ice and water saturations are reflective of the flow regime. All four scenarios show asimilar active layer extent, indicated by the thickness of the ice-free layer near the surface,but their distributions are different with respect to the cover thickness. Although theMinerals 2021, 11, 565 15 of 26waste rock remains mostly frozen for the no-cover case, a shallow active layer forms inthe summer, with increased water saturations and pronounced internal ponding, resultingfrom deeper infiltration due to the relatively high hydraulic conductivity of the PAG wasterock in combination with the blocking effect of permafrost. For the cases including acover, the cover with its lower hydraulic conductivity adds more complexity to the iceand water distribution. Zones of water and ice saturations are both observed within thecover material in the 4 m cover and 6 m cover scenarios. The ice layer formed above theinterface between the waste rock and cover is thicker in the 6 m cover scenario, whereaselevated water saturations can be seen slightly extending into the PAG waste rock in the2 m cover scenario (i.e., partial thawing of the shallow waste rock). These results confirmthat the active layer is likely to penetrate into PAG material for covers with a thickness ofless than 4 m.Minerals 2021, 11, x FOR PEER REVIEW 15 of 26 elevated water saturations can be seen slightly extending into the PAG waste rock in the 2 m cover scenario (i.e., partial thawing of the shallow waste rock). These results confirm that the active layer is likely to penetrate into PAG material for covers with a thickness of less than 4 m. Figure 8. Simulated contours of temperature, ice saturation, and water distributions for no cover, 2 m cover, and 6 m cover scenarios compared to the base case with 4 m cover in the summer after 50 years. The white line and white arrows indicate the interface between waste rock and cover materials, where applicable. As a function of the presence of the cover and cover thickness, the occurrence of freeze-thaw cycles in the active layer affects geochemical processes within the WRP in different ways. Simulated changes of pyrite volume fractions in the no-cover and 2 m-cover scenarios indicate pyrite oxidation is occurring in the sulfide-bearing rock during the summer months (Figure A5, Appendix C). Pyrite oxidation is more extensive in the no-cover scenario since the active layer is entirely in the PAG material, whereas the active layer in the 2 m cover scenario straddles both NPAG and PAG materials, which leads to a thin zone of pyrite oxidation immediately below the cover. However, with thicker cover designs (4 and 6 m cover scenarios), barely any pyrite is oxidized due to persistent frozen conditions in the PAG waste rock over the 50-year simulation period. The oxidation of pyrite during the warm seasons in the no-cover and 2 m cover scenarios release sulfate, as shown in Figure 9. As expected, sulfate concentrations are the highest in the active layer of the no-cover scenario. In the 2 m cover scenario, the release of sulfate is concentrated in the PAG-material located directly below the cover with some spreading into the cover material, which results from mixing in response to lateral flow as well as solute diffusion. In this case, sulfate release is most visible the uppermost bench (Figure A5, Appendix C), followed by a downslope flushing with infiltrating water. However, the covered pile in the 4 m- and 6 m-cover designs remain sulfate-free, because pyrite oxidation does not occur in the long-term. No cover2m cover4m cover6m coverAfter 50 yearstemperature (℃)0-2 2-4 4-6 6-8 8 10-10-12ice saturation [-]0 0.2 0.4 0.6 0.8 1.0water saturation [-]0 0.2 0.4 0.6 0.8 1.00 25mFigure 8. Simulated contours of temperature, ice saturation, and water distributions for no cover, 2 m cover, and 6 m coverscenarios compared to the base case with 4 m cover in the summer after 50 years. The white line and white ar ows indicatethe interface between waste rock and cover materials, where applicable.As a function of the presence of the cover and cover thicknes , the oc ur ence offre ze-thaw cycles in the active layer affects geochemical proces es within the WRP indif erent ways. Si ulate cha es of pyrite volu e fractio s in the no-cover and 2 m-cover scenarios indicate pyrite oxidation is oc ur ing in the sulfide-bearing rock duringthe summer months (Figure A5, Ap endix C). Pyrite oxidation is more extensive in theno-cover scenario since the active layer is entirely in the PAG material, whereas the activelayer in the 2 m cover scenario strad les both NPAG and PAG materials, whic leads to athin zone of pyrite oxidation im ediately below the cover. However, with thicker coverdesigns (4 and 6 m cover scenarios), barely any pyrite is oxidized due to persistent frozenconditions in the PAG waste rock over the 50-year simulation period. The oxidation ofpyrite during the warm seasons in the no-cover and 2 m cover scenarios release sulfate, asshown in Figure 9. As expected, sulfate concentrations are the highest in the active layerof the no-cover scenario. In the 2 cover scenario, the release of sulfate is concentratedMinerals 2021, 11, 565 16 of 26in the PAG-material located directly below the cover with some spreading into the covermaterial, which results from mixing in response to lateral flow as well as solute diffusion.In this case, sulfate release is most visible the uppermost bench (Figure A5, Appendix C),followed by a downslope flushing with infiltrating water. However, the covered pile in the4 m- and 6 m-cover designs remain sulfate-free, because pyrite oxidation does not occur inthe long-term.Minerals 2021, 11, x FOR PEER REVIEW 16 of 26 Figure 9. Simulated contours of sulfate concentration, calcium concentration, and pH for four cover scenarios in the sum-mer after 50 years. The white line and white arrows indicate the interface between waste rock and cover materials, where applicable. Both the cover material and PAG waste rock contain carbonate materials. The disso-lution of calcite, as the sole carbonate mineral included in the model, was tightly corre-lated with pyrite oxidation and sulfate production as seen in Figure 9. In the no-cover scenario, pyrite oxidation in the active layer promotes calcite dissolution to buffer acidity. Iron, calcium, and sulfate produced from pyrite oxidation along with calcite dissolution drives the secondary precipitation of both ferrihydrite and gypsum (Figure A6 and Figure A7, Appendix C). Although pH is buffered in the no cover and 2 m cover scenarios, pH is lower than in the other scenarios, because pyrite oxidation adds acidity to the system. For the 4m- and 6m-cover scenarios, calcite dissolution leads to mildly basic conditions in the cover. Geochemical reactions occurring within the PAG waste rock control the effluent ge-ochemistry at the toes of the pile (Figure 10). Results for both the 4 m and 6 m cover de-signs indicate that basal drainage remains nearly free of dissolved sulfate over the 50-year simulation period. However, the simulations also show that sulfate mass loadings in basal drainage are higher for the 2 m cover and the no-cover scenarios. In addition, sulfate mass loadings are not equally distributed at each of the toes of the pile due to pile geometry (longer and gentler slope on the right side). Although simplifications have been made regarding the mineralogical composition and corresponding reactions, the sensitivity analysis is helpful at illustrating the effectiveness of thermal covers in a geochemical and mass loading context. The simulations results suggest that for the climate conditions and materials considered, a cover with a thickness of less than 4 m would likely allow the active layer to extend into the PAG material, possibly leading to the release of the contam-inated drainage at the toes of the pile. It is reassuring to note that there were negligible differences in mass loadings between 4 and 6 m cover designs. Figure 9. Simulated contours of sulfate concentration, calcium concentration, and pH for four cover scenarios in thesummer after 50 years. The white line and white arrows indicate the interface between waste rock and cover materials,where applicable.Both the cover material and PAG waste rock contain carbonate materials. The dissolu-tion of calcite, as the sole carbonate mineral included in the model, was tightly correlatedwith pyrite oxidation and sulfate production as seen in Figure 9. In the no-cover scenario,pyrite oxidation in the active layer promotes calcite dissolution to buffer acidity. Iron, cal-cium, and sulfate produced from pyrite oxidation along with calcite dissolution drives thesecondary precipitation of both ferrihydrite and gypsum (Figures A6 and A7, Appendix C).Although pH is buffered in the n cover and 2 m cover scenarios, pH is lower than inthe oth r scenarios, because pyrite oxidation adds acidity to the system. For the 4m- and6m-cover scenarios, cal ite dissolution leads to mildly basic conditions in the cover.Geoche ical reactions occurring within the PAG waste rock control the effluentgeochemistry at the toes of the pile (Figure 10). Results for both the 4 m and 6 m coverdesigns indicate that basal drainage remains nearly free of dissolved sulfate over the50-year simulation period. However, the simulations also show that sulfate mass loadingsin basal drainage are higher for the 2 m cover and the no-cover scenarios. In addition,sulfate mass loadings are not equally distributed at each of the toes of the pile due to pilegeometry (longer and gentler slope on the right side). Although simplifications have beenmade regarding the mineralogical composition and corresponding reactions, the sensitivityanalysis is helpful at illustrating the effectiveness of thermal covers in a geochemical andMinerals 2021, 11, 565 17 of 26mass loading context. The simulations results suggest that for the climate conditions andmaterials considered, a cover with a thickness of less than 4 m would likely allow the activelayer to extend into the PAG material, possibly leading to the release of the contaminateddrainage at the toes of the pile. It is reassuring to note that there were negligible differencesin mass loadings between 4 and 6 m cover designs.Minerals 2021, 11, x FOR PEER REVIEW 17 of 26 Figure 10. Mass loadings of sulfate at two basal drainage locations for four cover design for the 50-year simulation period. 5.4. Implications for Cover Design and Construction The base case simulation results, which are constrained by experimental data and field measurements, provide insights into the coupled processes occurring in sulfide-bear-ing WRPs in permafrost environments. The results forecasted by the model illustrate the influences of freeze-thaw cycles on flow, thermal, and geochemical regimes of the covered pile in the short term, consistent with available observations for a four-year time period. Long-term simulations suggest that, for the given climate conditions and material param-eters, conditions in the WRP will remain stable, subject to seasonally driven dynamic equi-librium conditions. Although some assumptions had to be made regarding infiltration rates and uncertainties regarding the magnitude and temporal distribution of infiltration remains, the simulations indicate that the performance of thermal covers is not very sen-sitive to infiltration. The simulation results suggest that the thickness of the active layer is not significantly affected by uncertainties in infiltration rates, with impacts mostly re-stricted to changes in the shallow active layer flow system. In this context, the results in-dicate that the impacts of infiltration on mass loadings are likely limited, as long as the active layer remains entirely within the thermal cover. Additional sensitivity analysis re-lated to the cover design demonstrates the functions and effectiveness of the thermal cover as a unique reclamation approach commonly used for WRPs in cold-region climates. The simulation results confirm and provide confidence in the cover’s ability to prevent the sulfide-bearing waste rock from being oxidized and to effectively manage the effluent quality in the long-term. Simulated results for mass loadings also imply that the drainage quality is affected by the cover thickness. In general, the overall results from this study indicate that modeling can be used as a guiding tool to help identify and design suitable reclamation approaches. Different scenarios can be simulated in the design planning stage, prior to the costly construction of field experimental studies, and the simulation results can provide insights into the ability of different proposed cover designs to control potential ARD/CND. Reactive transport modeling of thermal-hydrological and geochem-ical processes can be used to address parameter uncertainties, assess different future sce-narios, and illustrate thermal, hydrological, and geochemical processes in WRPs. Figure 10. Mass loadings of sulfate at two basal drainage locations for four cover design for the 50-year simulation period.5.4. Implications for Cover Design and ConstructionThe base case simulation results, which are constrained by experimental data andfield measurements, provide insights into the coupled processes occurring in sulfide-bearing WRPs in permafrost environments. The results forecasted by the model illustratethe influences of freeze-thaw cycles on flow, thermal, and geoche ical regimes of thecovered pile in the short term, consistent with available observations for a four-yeartime period. Long-term simulations suggest that, for the given climate conditions andmaterial parameters, conditions in the WRP will remain stable, subject to seasonally drivendynamic equilibrium conditions. Although some assumptions had to be mad regardinginfiltration rates and uncertainties regarding the magnitude and temporal distribution ofinfiltration remains, the simulations indicate that the performance of thermal covers isnot very sensitive to infiltration. The simulation results suggest that the thickness of theactive layer is not significan ly affected by uncertainties in infiltration rates, with impactsmostly restricted to changes in the shallow active layer flow system. In this context, theresults indicate that the impacts of infiltration on mass loadings are likely limited, aslong as the active layer remains entirely within the thermal cover. Additional sensitivityanalysis related to the cov r d sign demonstrates the functions and effectiven ss of thethermal cover as a unique reclamation approach commonly used for WRPs in cold-regionclimates. The simulation results confirm and provide confidence in the cover’s abilityto prevent the sulfide-bearing waste rock from being oxidized and to effectively managethe effluent quality in the long-term. Simulated results for mass loadings also imply thatthe drainage quality is affected by the cover thickness. In general, the overall resultsfrom this study indicate that modeling can be used as a guiding tool to help identifyand design suitable reclamation approaches. Different scenarios can be simulated in thedesign planni g stage, prior to the costly construction of field experimental studies, and thesimulation results can provide insights into the ability of different proposed cover designsMinerals 2021, 11, 565 18 of 26to control potential ARD/CND. Reactive transport modeling of thermal-hydrological andgeochemical processes can be used to address parameter uncertainties, assess differentfuture scenarios, and illustrate thermal, hydrological, and geochemical processes in WRPs.6. ConclusionsIn addition to experimental or field tests, RTM provides a versatile tool for assessingthe performance of thermal covers in permafrost environments and their impact on limitingsulfide mineral oxidation and mass loadings. MIN3P-HPC, the modeling tool used in thisstudy, has the ability of to simulate thermo-hydrological-geochemical interactions in acovered WRP under permafrost conditions. Simulation results provide insights and alsobuild confidence for using modeling to analyze other similar scenarios through sensitivityand uncertainty analyses.This study demonstrates the capabilities of MIN3P-HPC to simulate the problemsrelated to mine waste in permafrost regions undergoing frequent freeze-thaw cycles. Freeze-thaw processes and their impact on water flow, geochemical reactions, and mass transportare incorporated into the code, facilitating an integrated assessment of these coupledprocesses. The recently added capability of unstructured grids allows to consider complexgeometries with angled surfaces and multiple layers, common at mine sites. Future studieson thermal cover performance in cold-region climates could address additional complexity,by considering physical and chemical heterogeneity, a more detailed geochemical reactionnetwork, and climate warming. In addition, future simulations could consider insulatingeffects of seasonal snow covers, which may aid to further improving the understandingof covered and uncovered WRPs in cold-region climates, and the prognostic value ofthe model.Author Contributions: Conceptualization, X.Y., K.U.M. and D.S.; data curation and methodology,experimental work, laboratory measurements, B.B.; computer simulations by X.Y. with guidancefrom D.S. and K.U.M.; code development by X.Y. with assistance from D.S.; validation, X.Y., K.U.M.and D.S. with site-specific guidance from B.B.; formal analysis and investigation, X.Y. and K.U.M.;resources, B.B., K.U.M., and D.S.; writing—original draft preparation, X.Y.; writing—review andediting, K.U.M. and B.B.; visualization, X.Y. and D.S.; supervision, K.U.M.; funding acquisition,K.U.M. and B.B. 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 ofCanada (NSERC) through the strategic network grant “Toward Environmentally Responsible ResourceExtraction Network (TERRE-NET)” (Program director D.W. Blowes), grant number 479708-2015.Acknowledgments: The Meadowbank mine geotechnical team is greatly appreciated for their in-volvement and collaboration throughout the project. K. Raymond at UBC is also thanked for theassistance with the manuscript.Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the designof the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, orin the decision to publish the results.Minerals 2021, 11, 565 19 of 26Appendix A. Material PropertiesMinerals 2021, 11, x FOR PEER REVIEW 19 of 26 ppendix . aterial Properties Figure A1. Estimated water retention curves with parameters for both NPAG and PAG materials, based on laboratory data reported in [37]. Table A1. Mineralogical composition. Mineral 1 Volume Fraction NPAG PAG Actinolite 0.054 [35] - Calcite 0.052 [35] 0.028 [35] Chlorite 0.178 [35] 0.105 [35] Dolomite 0.117 [35] 0.025 [35] Magnetite - 0.034 [35] Muscovite - 0.052 [35] Quartz 0.047 [35] 0.334 [35] Pyrite - 0.002 [35] Talc 0.198 [35] 0.050 [35] Stilpnomelane - 0.081[35] Other Parameters Ssulfide (%) 0.03 [37] 0.36 [37] Ctotal (%) 1.92 [37] 0.27 [37] AP (kg CaCO3/t) 1.0 [37] 11 [37] NP (kg CaCO3/t) 120 [37] 20 [37] NP/AP 120 [37] 1.8 [37] 1 Calculated volume fractions of mineral phases from weight fraction (wt. %) reported in [35]. i r . sti te ater retention curves with parameters for both NPAG and PAG materials,l boratory data reported in [37].le 1. ineralogical composition.ineral 1Volume FractionNPAG PAGActinolite 0.054 [35] -Calcite 0.052 [35] 0.028 [35]hlorite 0.178 [35] 0.105 [35]olomite 0.117 [35] 0.025 [35]Magnetite - 0.034 [35]Muscovite - 0.052 [35]Quartz 0.047 [35] 0.334 [35]Pyrite - 0.002 [35]Talc 0.198 [35] 0.05 [35]Stilpnomelane - 0.081 [35]Other ParametersSsulfide (%) 0.03 [37] 0.36 [37]Ctotal (%) 1.92 [37] 0.27 [37]AP (kg CaCO3/t) 1.0 [37] 11 [37]NP (kg CaCO3/t) 120 [37] 20 [37]NP/AP 120 [37] 1.8 [37]1 Calculated volume fractions of mineral phases from weight fraction (wt. %) reported in [35].Minerals 2021, 11, 565 20 of 26Appendix B. Initial Conditions and Boundary ConditionsMinerals 2021, 11, x FOR PEER REVIEW 20 of 26 Appendix B. Initial Conditions and Boundary Conditions Figure A2. Initial conditions for ice saturation and temperature. Initial temperature given to the WRP is assumed ho-mogeneous with a temperature of −3 °C, while the temperature of natural soil was simulated with MAAT = −11.4 °C at the surface and geothermal heat flux of −0.087 W m−2 at the base until reaching a steady state. Initial ice saturation of the system was estimated based on the water retention curves for NPAG and PAG materials as well as the underlying soil. Boundary conditions for flow, reactive transport, and heat transport are assigned at each boundary. Transient flow and temperature boundaries were assigned at the surface, while a constant geothermal heat flux boundary was assigned at the base. Other boundaries were assumed as impermeable for flow and transport. Figure A2. Initial conditions for ice saturation and temperature. Initial temperature give t t s edhomogeneous with a t mpe ature of −3 ◦C, while th temperature of ral soil was imulated with MAAT = −11.4 ◦C atthe surf ce and geothermal heat flux of − .087 W 2 at the til reaching a steady state. Initial ice satur ion ofthe system was estimated based on the retention curves for NPAG and PAG materials as well as the underlying soil.ary condit ons for flow, reactive transport, and heat transport are assigned at each boundary. Transient flow a dte rature boundaries were assigned at the surface, while a constant geothermal at flux boundary was assigned at thebase. Other boundaries were assumed as impermeable for flow and transport.Minerals 2021, 11, 565 21 of 26Appendix C. Results of Sensitivity Analysis to Cover ThicknessMinerals 2021, 11, x FOR PEER REVIEW 21 of 26 i . lt f iti it l i t r i Figure A3. Log-scale sulfate concentrations for base case after 50 years (i.e., sulfate concentrations shown in Figure 6 are presented here in log-scale). The purpose of this figure is to highlight the limited sulfate content at the breaks of between some of the benches and slopes. The white arrows indicate the interface between waste rock and cover materials. Figure A4. Comparison of internal temperature fluctuations over time between the base case (4 m cover) and the no cover scenario in the short term. The depth is measured from the surface of PAG waste rock. Figure A3. Log-scale sulfate concentrations for base case after 50 years (i.e., sulfate concentrations shown in Figure 6 arepresented here in log-scale). The pu pose of this figure i to highlight the limited sulfat tent at the breaks of betweensome of the benches and slopes. The white arrows indicate the interface between waste rock and cover materi ls.Minerals 2021, 11, x FOR PEER REVIEW 21 of 26 Appendix C. Results of Sensitivity Analysis to Cover Thickness Figure A3. Log-scale sulfate concentrations for base case after 50 years (i.e., sulfate concentrations shown in Figure 6 are presented here in log-scale). The purpose of this figure is to highlight the limited sulfate content at the breaks of between some of the benches and slopes. The white arrows indicate the interface between waste rock and cover materials. Figure A4. Comparison of internal temperature fluctuations over time between the base case (4 m cover) and the no cover scenario in the short term. The depth is measured from the surface of PAG waste rock. Figure A4. Comparison of internal temperature fluctuations over time between the base case (4 m cover) and the no coversc nario in the short term. The depth is meas red from the surface of PAG waste rock.Minerals 2021, 11, 565 22 of 26Minerals 2021, 11, x FOR PEER REVIEW 22 of 26 Figure A5. Sulfate concentrations for the four cover scenarios. In the 2 m-cover scenario, high sulfate content in the top bench is only observed near the interface between cover and waste rock, indicating that the active layer partially extends into the sulfide-bearing waste rock. The broader spread of sulfate concenatrtions into the cover in the lower benches is due to mixing processes occuring during lateral drainage. The white line and white arrows indicate the interface between waste rock and cover materials. Figure A5. Sulfate concentrations for the four cover scenarios. In the 2 m-cover scenario, high sulfate content in the topbench is only observed near the interface between cover and waste rock, indicating that the active layer partially extendsinto the sulfide-beari g waste rock. The broader spread of sulfate concenatrtions into the cover in the lower benches is dueto mixing processes occuring during lateral drainage. The white line and white arrows i dicate the interface between wasterock and cover materials.Minerals 2021, 11, 565 23 of 26inerals 2021, 11, x FOR PEER REVIEW 23 of 26 Figure A6. Simulated contours of volumetric fraction changes of primary minerals, pyrite and calcite, for no cover, 2 m cover, and 6m cover scenarios compared to the base case with 4 m cover in the summer after 50 years. The negative sign (-) represents mineral dissolution. Only part of the pile is shown for illustrative purposes. The white line and white arrows indicate the interface between waste rock and cover materials. No cover2m cover4m cover6m coverAfter 50 yearsvolume fraction change of pyrite (-) (m3 m-3) volume fraction change of calcite (-) (m3 m-3) 0 0.0002 0.0004 0.0006 0.0008 0 0.0003 0.0006 0.0009 0.0012 0.0015 0.0018 0.0021 0.0024 0.00270 25mFigure A6. Simulated contours of volumetric fraction changes of primary minerals, pyrite and calcite, for no cover, 2 mcover, and 6m cover scenarios compared to the base case with 4 m cover in the summer after 50 years. The negative sign (-)represents mineral dissolution. Only part of the pile is shown for illustrative purposes. The white line and white arrowsindicate the interface between waste rock and cover materials.Minerals 2021, 11, 565 24 of 26Minerals 2021, 1, x FOR PEER REVIEW 24 of 26 Figure A7. Simulated contours of volumetric fraction changes of secondary minerals, ferrihydrite and gypsum, for no cover, 2 m cover, and 6m cover scenarios compared to the base case with a 4 m cover in the summer after 50 years. The positive sign (+) represents mineral precipitation. Only part of the pile is shown for illustrative purposes. The white line and white arrows indicate the interface between waste rock and cover materials. References 1. Amos, R.T.; Blowes, D.W.; Bailey, B.L.; Sego, D.C.; Smith, L.; Ritchie, A.I.M. Waste-rock hydrogeology and geochemistry. Appl. Geochem. 2015, 57, 140–156, doi:10.1016/j.apgeochem.2014.06.020. 2. Vriens, B.; St. Arnault, M.; Laurenzi, L.; Smith, L.; Mayer, K.U.; Beckie, R.D. Localized sulfide oxidation limited by oxygen supply in a full-scale waste-rock pile. Vadose Zone J. 2018, 17, 180119, doi:10.2136/vzj2018.06.0119. 3. Lefebvre, R.; Hockley, D.; Smolensky, J.; Lamontagne, A. Multiphase transfer processes in waste rock piles producing acid mine drainage. 2. Applications of numerical simulation. J. Contam. Hydrol. 2001, 52, 165–186, doi:10.1016/S0169-7722(01)00157-7. 4. Demers, I.; Mbonimpa, M.; Benzaazoua, M.; Bouda, M.; Awoh, S.; Lortie, S.; Gagnon, M. Use of acid mine drainage treatment sludge by combination with a natural soil as an oxygen barrier cover for mine waste reclamation: Laboratory column tests and intermediate scale field tests. Miner. Eng. 2017, 107, 43–52, doi:10.1016/j.mineng.2016.11.017. 5. MEND. Acid mine drainage in permafrost regions: Issues, control strategies and research requirements. In MEND Project 1.61.2; MEND: Smithers, BC, Canada, 1996; pp. 25–35. 6. Bussière, B. Acid mine drainage from abandoned mine sites: Problematic and reclamation approaches. In Advances in Environ-mental Geotechnics; Springer: Berlin, Heidelberg, Germany, 2010; pp. 111–125, doi:10.1007/978-3-642-04460-1. 7. Bussière, B.; Guittonny, M. Hard Rock Mine Reclamation: From Prediction to Management of Acid Mine Drainage; Bussière, B., Guittonny, M., Eds.; CRC Press: Boca Raton, FL, USA, 2020; Chapter 4–9, pp. 1–15, 93–222. No cover2m cover4m cover6m coverAfter 50 yearsvolume fraction change of ferrihydrite (+) (m3 m-3) 0 0.0002 0.0004 0.0006 0.0008volume fraction change of gypsum (+) (m3 m-3) 0 0.0005 0.0015 0.002 0.00370.001 0.0025 0.0030 25mFigure A7. Simulated contours of volumetric fraction changes of secondary minerals, ferrihydrite and gypsum, for no cover,2 m cover, and 6m cover scenarios compared to the base case with a 4 m cover in the summer after 50 years. The positivesign (+) represents mineral precipitation. Only part of the pile is shown for illustrative purposes. The white line and whitearrows indicate the interface between waste rock and cover materials.References1. Amos, R.T.; Blowes, D.W.; Bailey, B.L.; Sego, D.C.; Smith, L.; Ritchie, A.I.M. Waste-rock hydrogeology and geochemistry.Appl. Geochem. 2015, 57, 140–156. [CrossRef]2. Vriens, B.; St. Arnault, M.; Laurenzi, L.; Smith, L.; Mayer, K.U.; Beckie, R.D. Localized sulfide oxidation limited by oxygen supplyin a full-scale waste-rock pile. Vadose Zone J. 2018, 17, 180119. [CrossRef]3. Lefebvre, R.; Hockley, D.; Smolensky, J.; Lamontagne, A. Multiphase transfer processes in waste rock piles producing acid minedrainage. 2. Applications of numerical simulation. J. Contam. Hydrol. 2001, 52, 165–186. [CrossRef]4. Demers, I.; Mbonimpa, M.; Benzaazoua, M.; Bouda, M.; Awoh, S.; Lortie, S.; Gagnon, M. Use of acid mine drainage treatmentsludge by combination with a natural soil as an oxygen barrier cover for mine waste reclamation: Laboratory column tests andintermediate scale field tests. Miner. Eng. 2017, 107, 43–52. [CrossRef]5. MEND. Acid mine drainage in permafrost regions: Issues, control strategies and research requirements. In MEND Project 1.61.2;END: Smithers, BC, Canada, 1996; pp. 25–35.6. Bussière, B. Acid mine drainage from abandoned mine sites: Problematic and reclamation approaches. In Advances in Environmen-tal Geotechnics; Springer: Berlin/Heidelberg, Germany, 2010; pp. 111–125. [CrossRef]7. Bussière, B.; Guittonny, M. Hard Rock Mine Reclamation: From Prediction to Management of Acid Mine Drainage; Bussière, B.,Guittonny, M., Eds.; CRC Press: Boca Raton, FL, USA, 2020; Chapter 4–9; pp. 1–15, 93–222.8. Lin, S.M.; Hayward, G.; Leonard, S.; Gene, S. Electrochemical cover for the prevention of acid mine drainage—A laboratory test.In Proceedings of the Sudbury 2003 Mining and the Environment Conference, Sudbury, ON, Canada, 25–28 May 2003.Minerals 2021, 11, 565 25 of 269. Kyhn, C.; Elberling, B. Frozen cover actions limiting AMD from mine waste deposited on land in Arctic Canada. Cold Reg. Sci.Technol. 2001, 32, 133–142. [CrossRef]10. MEND. Mine waste covers in cold regions. In MEND Report 1.16.5a; MEND: Smithers, BC, Canada, 2009; pp. 35–54.11. Woo, M. Permafrost Hydrology; Springer: Berlin/Heidelberg, Germany, 2012; Volume 110, pp. 5–15. ISBN 9788578110796.12. Boulanger-Martel, V.; Bussière, B.; Côté, J. Thermal behaviour and performance of two field experimental insulation covers tocontrol sulfide oxidation at Meadowbank mine, Nunavut. Can. Geotech. J. 2020, 58, 427–440. [CrossRef]13. Pham, H.N. Heat Transfer in Waste Rock Piles Construted in a Contiuous Permafrost Region. Ph.D. Thesis, University of Alberta,Edmonton, AB, Canada, 2013.14. Stevens, C.W.; Shapka-Fels, T.; Rykaart, M. Thermal Cover Design for Mine Waste Facilities in Cold Regions. 2018. Availableonline: https://www.srk.com/en/publications/thermal-cover956design-for-mine-waste-facilities-in-cold-regions (accessedon 30 September 2018).15. Lefebvre, R.; Hockley, D.; Smolensky, J.; Geélinas, P. Multiphase transfer processes in waste rock piles producing acid minedrainage. 1: Conceptual model and system characterization. J. Contam. Hydrol. 2001, 52, 137–164. [CrossRef]16. Lahmira, B.; Lefebvre, R.; Aubertin, M.; Bussière, B. Effect of heterogeneity and anisotropy related to the construction method ontransfer processes in waste rock piles. J. Contam. Hydrol. 2016, 184, 35–49. [CrossRef]17. Ramasamy, M.; Power, C.; Mkandawire, M. Numerical prediction of the long-term evolution of acid mine drainage at a wasterock pile site remediated with an HDPE-lined cover system. J. Contam. Hydrol. 2018, 216, 10–26. [CrossRef] [PubMed]18. Molson, J.; Aubertin, M.; Bussière, B.; Benzaazoua, M. Geochemical transport modelling of drainage from experimental minetailings cells covered by capillary barriers. Appl. Geochem. 2008, 23, 1–24. [CrossRef]19. Ouangrawa, M.; Molson, J.; Aubertin, M.; Bussière, B.; Zagury, G.J. Reactive transport modelling of mine tailings columns withcapillarity-induced high water saturation for preventing sulfide oxidation. Appl. Geochem. 2009, 24, 1312–1323. [CrossRef]20. Bea, S.A.; Wilson, S.A.; Mayer, K.U.; Dipple, G.M.; Power, I.M.; Gamazo, P. Reactive transport modeling of natural carbonsequestration in ultramafic mine tailings. Vadose Zone J. 2012, 11, vzj2011.0053. [CrossRef]21. Demers, I.; Molson, J.; Bussière, B.; Laflamme, D. Numerical modeling of contaminated neutral drainage from a waste-rock fieldtest cell. Appl. Geochem. 2013, 33, 346–356. [CrossRef]22. Wilson, D.; Amos, R.T.; Blowes, D.W.; Langman, J.B.; Smith, L.; Sego, D.C. Diavik Waste Rock Project: Scale-up of a reactivetransport model for temperature and sulfide-content dependent geochemical evolution of waste rock. Appl. Geochem. 2018, 96,177–190. [CrossRef]23. Raymond, K.E.; Seigneur, N.; Su, D.; Poaty, B.; Plante, B.; Mayer, K.U. Numerical modeling of a laboratory-scale waste rock pilefeaturing an engineered cover system. Minerals 2020, 10, 652. [CrossRef]24. Kurylyk, B.L.; Watanabe, K. The mathematical representation of freezing and thawing processes in variably-saturated, non-deformable soils. Adv. Water Resour. 2013, 60, 160–177. [CrossRef]25. McKenzie, J.M.; Voss, C.I.; Siegel, D.I. Groundwater flow with energy transport and water-ice phase change: Numericalsimulations, benchmarks, and application to freezing in peat bogs. Adv. Water Resour. 2007, 30, 966–983. [CrossRef]26. Zhao, L.; Gray, D.M.; Male, D.H. Numerical analysis of simultaneous heat and mass transfer during Infiltration into FrozenGround. J. Hydrol. 1997, 200, 345–363. [CrossRef]27. Voss, C.I.; Provost, A. Version of September 22, 2010 (SUTRA Version 2.2) Latest Version. USGS Water Resource InvestigationsReport 02-4231 2010. Available online: http://water.usgs.gov/nrp/gwsoftware (accessed on 1 September 2010).28. Mayer, K.U.; Frind, E.O.; Blowes, D.W. Multicomponent reactive transport modeling in variably saturated porous media using ageneralized formulation for kinetically controlled reactions. Water Resour. Res. 2002, 38, 13-1–13-21. [CrossRef]29. Su, D.; Mayer, K.U.; MacQuarrie, K.T.B. Numerical investigation of flow instabilities using fully unstructured discretization forvariably saturated flow problems. Adv. Water Resour. 2020, 143. [CrossRef]30. Awoh, S.A.; Bussière, B.; Batzenschlager, C.; Boulanger-Martel, V.; Lepine, T.; Voyer, E. Design, construction, and preliminaryresults of two insulation covers at the Meadowbank mine. Geo Chicago 2016, 2016, 185–194.31. Collette, L. Cryohydrogeology of a Covered Waste Rock Pile in Permafrost Environment: Large Scale Field Experiment andFreeze-Thaw Numerical Investigations. Master’s Thesis, University of British Columbia, Vancouver, BC, Canada, 2017.32. Mockus, V.; Mullem, J.A.V.; Garen, D.; Woodward, D.E. Chapter 11 Snowmelt. In USDA Part 630 Hydrology National EngineeringHandbook; July 2004; pp. 11–12. Available online: https://directives.sc.egov.usda.gov/OpenNonWebContent.aspx?content=17753.wba (accessed on 25 May 2021).33. Elberling, B. Environmental controls of the seasonal variation ia oxygen uptake in sulfidic tailings deposited in a permafrost-affected area. Water Resour. Res. 2001, 37, 99–107. [CrossRef]34. Meldrum, J.L.; Jamieson, H.E.; Dyke, L.D. Oxidation of mine tailings from Rankin Inlet, Nunavut, at subzero temperatures.Can. Geotech. J. 2001, 38, 957–966. [CrossRef]35. Boulanger-Martel, V.; Bussière, B.; Côté, J.; Gagnon, P. Design, construction, and preliminary performance of an insulation coverwith capillary barrier effects at Meadowbank mine, Nunavut. In Proceedings of the 70th Canadian Geotechnical Conference,Ottawa, ON, Canada, 1–10 October 2017; p. 354.36. Voyer, E.; Robert, S. Waste rock management and closure planning in Northern climate: The Meadowbank Mine, Nunavut. InProceedings of the Symposium 2015 sur l’environnement et les mines, Rouyn-Noranda, Rouyn-Noranda, QC, Canada, 14–17June 2015; Published by the Canadian Institute of Mining, Metallurgy and Petroleum (CIM): Ottawa, ON, Canada, 2015; pp. 1–15.Minerals 2021, 11, 565 26 of 2637. Poirier, A. Étude du Comportement Thermique D’une Halde à Stériles en Milieu Nordique. Master’s Thesis, Université duQuébec en Abitibi-Témiscamingue, Rouyn-Noranda, QC, Canada, 2019.38. Boulanger-martel, V.; Poirier, A.; Abitibi-témiscamingue, Q.; Côté, J.; Bussière, B. Thermal conductivity of Meadowbank’s mine wasterocks and tailings. In Proceedings of the 71th Canadian Geotechnical Conference, Edmonton, AB, Canada, 23–26 September 2018.39. Peregoedova, A. Étude Expérimentale des Propriétés Hydrogéologiques des Roches Stériles à une Échelle Intermédiaire deLaboratoire. Master’s Thesis, École Polytechnique de Montréal, Montréal, QC, Canada, 2012.40. Chapuis, R.P.; Masse, I.; Madinier, B.; Aubertin, M. A drainage column test for determining unsaturated properties of coarsematerials. Geotech. Test. J. 2007, 30, 83–89. [CrossRef]41. Côté, J.; Konrad, J.M. Assessment of structure effects on the thermal conductivity of two-phase porous geomaterials. Int. J. HeatMass Transf. 2009, 52, 796–804. [CrossRef]42. Côté, J.; Konrad, J.M. Indirect methods to assess the solid particle thermal conductivity of Quebec marine clays. Can. Geotech. J.2007, 44, 1117–1127. [CrossRef]43. Côté, J.; Konrad, J.M. A generalized thermal conductivity model for soils and construction materials. Can. Geotech. J. 2005, 42,443–458. [CrossRef]44. Van Genuchten, M.T.; Leij, F.J.; Yates, S.R. The RETC code for quantifying the hydraulic functions of unsaturated soils. Dep. Agric.Agric. Res. Serv. 1991, 65, 1–83.45. Henderson, T.H.; Mayer, K.U.; Parker, B.L.; Al, T.A. Three-dimensional density-dependent flow and multicomponent reactivetransport modeling of chlorinated solvent oxidation by potassium permanganate. J. Contam. Hydrol. 2009, 106, 195–211. [CrossRef][PubMed]46. Su, D.; Ulrich Mayer, K.; MacQuarrie, K.T.B. Parallelization of MIN3P-THCm: A high performance computational framework forsubsurface flow and reactive transport simulation. Environ. Model. Softw. 2017, 95, 271–289. [CrossRef]47. Nicholson, R.V.; Gillham, R.W.; Reardon, E.J. Pyrite oxidation in carbonate-buffered solution: 1. Experimental kinetics. Geochim.Cosmochim. Acta 1988, 52, 1077–1085. [CrossRef]48. Arrhenius, S. Über die Dissociationswärme und den Einfluss der Temperatur auf den Dissociationsgrad der Elektrolyte. Z. Phys.Chem. 1889, 4, 96–116. [CrossRef]49. Xie, X.; Xiao, S.; He, Z.; Liu, J.; Qiu, G. Microbial populations in acid mineral bioleaching systems of Tong Shankou Copper Mine,China. J. Appl. Microbiol. 2007, 103, 1227–1238. [CrossRef] [PubMed]50. Yi, X. Reactive Transport Modeling of Waste Rock Weathering in Permafrost Environments. Master’s Thesis, University of BritishColumbia, Vancouver, BC, Canada, 2021.51. Booshehrian, A.; Wan, R.; Su, X. Hydraulic variations in permafrost due to open-pit mining and climate change: A case study inthe Canadian Arctic. Acta Geotech. 2020, 15, 883–905. [CrossRef]52. Allison, J.D.; Brown, W.; Novo-Gradac, K.J. MINTEQA2/PRODEFA2, A Geochemical Assessment Model for Environmental Systems:Version 3.0 User’s Manual; EPA/600/3-91/021; Environmental Research Laboratory, U.S. Environmental Protectio Agency:Washington, DC, USA, 1991.53. Mayer, K.U.; Alt-Epping, P.; Jacques, D.; Arora, B.; Steefel, C.I. Benchmark problems for reactive transport modeling of thegeneration and attenuation of acid rock drainage. Comput. Geosci. 2015, 19, 599–611. [CrossRef]54. Amos, R.T.; Bailey, B.L.; Pham, H.N.; Fretz, N.; Hannam, S.; Blowes, D.W.; Sego, D.C.; Smith, L. Diavik Waste Rock Project:Northern Aspects and Scaling Predictions. In Proceedings of the 19th British Columbia MEND Metal Leaching/Acid RockDrainage Workshop, Vancouver, BC, Canada, 28–29 November 2012.55. Chevé, N.; Guittonny, M.; Bussière, B. Water budget of field experimental cells with vegetated and non-vegetated soil layersplaced on waste rock. Mine Clos. 2018, 2018, 549–562.56. Pham, N.H.; Sego, D.C.; Arenson, L.U.; Blowes, D.W.; Smith, L. Convective heat transfer in waste rock piles under permafrostenvironment. In Proceedings of the GeoEdmonton 2008: The 61th Canadian Geotechnical Conference, Edmonton, AB, Canada,21–24 September 2008; pp. 940–947.57. Wilson, D.; Amos, R.T.; Blowes, D.W.; Langman, J.B.; Ptacek, C.J.; Smith, L.; Sego, D.C. Diavik waste rock project: A conceptualmodel for temperature and sulfide-content dependent geochemical evolution of waste rock—Laboratory scale. Appl. Geochem.2018, 89, 160–172. [CrossRef]
Article
eng
Reviewed
Vancouver : University of British Columbia Library
Multidisciplinary Digital Publishing Institute
CC BY 4.0
Faculty
Researcher
reactive transport modeling
mine waste rock
thermal cover
freeze-thaw cycles
permafrost
contaminated drainage
Thermal-Hydrological-Chemical Modeling of a Covered Waste Rock Pile in a Permafrost Region
Text