Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Raymond, Katherine Elizabeth
Master of Science - MSc
University of British Columbia
Design strategies for waste rock piles (WRPs) to mitigate or reduce the potential for poor drainage quality which have associated environmental and economic liabilities are an important consideration for mining operations. Due to their size and complexity, conceptualizing internal processes and forecasting drainage quality from WRPs has proved challenging. Reactive transport modeling (RTM) is a tool that can help improve understanding of complex interactions of physical and geochemical processes within WRPs and compare the effects of pile design on drainage water quality. In this thesis, RTM is used to investigate the effect of physical and chemical heterogeneities, pile construction methods and scale on mass loadings in drainage from full-scale WRPs through a series of stylistic models. RTM is also applied to a laboratory-scale model to simulate an engineered cover system as part of the long-term design plan for the Tio Mine, Quebec, Canada. Results from the stylistic models indicate that heterogeneity and structural features in WRPs can cause significant variations in effluent quality, when compared to a homogeneous case. Simulated peak loadings are reduced due to the presence of heterogeneities and internal structure, at the same time prolonging the release of poor-quality drainage. The results also indicate that pile construction methods (e.g. end- and push-dumping) affect drainage release in different ways; however, these differences vanish for larger multi-lift piles. Simulations suggest that multiple lifts and the presence of buried traffic surface homogenize mass loadings, independent of construction methods. Simulations also indicate that peak mass loadings from waste rock piles are reduced by a factor of 2-3 when heterogeneities and structure are considered, which has implications for the design of water treatment systems. These results were put to practice in simulations of a laboratory-scale model, where heterogeneities proved essential for capturing effluent quality from different material zones, and for conceptualizing internal flow paths occurring within the dry cover system in place. In addition, results from this research demonstrate the capabilities of the recently developed unstructured grid code MIN3P-HPC for capturing complex shapes and geometries, which may prove useful for future investigations of capillary barrier effects and other design strategies in WRPs.
APPLICATIONS OF REACTIVE TRANSPORT MODELING FOR ASSESSING THE EFFECTS OF HETEROGENEITIES AND INTERNAL STRUCTURE ON MASS LOADINGS FROM WASTE ROCK PILES by Katherine Elizabeth Raymond B.Sc., The University of Waterloo, 2015 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Geological Sciences) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2020 © Katherine Elizabeth Raymond, 2020 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis entitled: Applications of Reactive Transport Modeling for Assessing the Effects of Heterogeneities and Internal Structure on Mass Loadings from Waste Rock Piles Submitted by Katherine Raymond, in partial fulfillment of the requirements for the degree of Master of Science in Geological Sciences. Examining Committee: Prof. K. Ulrich Mayer, The University of British Columbia, Geological Sciences Supervisor Prof. Roger D. Beckie, The University of British Columbia, Geological Sciences Supervisory Committee Member Prof. Bruno Bussière, Université du Québec en Abitibi-Témiscamingue, Mines and the Environment Supervisory Committee Member Prof. Bern Klein, The University of British Columbia, NBK Mining Engineering Additional Examiner iii Abstract Design strategies for waste rock piles (WRPs) to mitigate or reduce the potential for poor drainage quality which have associated environmental and economic liabilities are an important consideration for mining operations. Due to their size and complexity, conceptualizing internal processes and forecasting drainage quality from WRPs has proved challenging. Reactive transport modeling (RTM) is a tool that can help improve understanding of complex interactions of physical and geochemical processes within WRPs and compare the effects of pile design on drainage water quality. In this thesis, RTM is used to investigate the effect of physical and chemical heterogeneities, pile construction methods and scale on mass loadings in drainage from full-scale WRPs through a series of stylistic models. RTM is also applied to a laboratory-scale model to simulate an engineered cover system as part of the long-term design plan for the Tio Mine, Quebec, Canada. Results from the stylistic models indicate that heterogeneity and structural features in WRPs can cause significant variations in effluent quality, when compared to a homogeneous case. Simulated peak loadings are reduced due to the presence of heterogeneities and internal structure, at the same time prolonging the release of poor-quality drainage. The results also indicate that pile construction methods (e.g. end- and push-dumping) affect drainage release in different ways; however, these differences vanish for larger multi-lift piles. Simulations suggest that multiple lifts and the presence of buried traffic surface homogenize mass loadings, independent of construction methods. Simulations also indicate that peak mass loadings from waste rock piles are reduced by a factor of 2-3 when heterogeneities and structure are considered, which has implications for the design of water treatment systems. These results were put to practice in simulations of a laboratory-scale model, where heterogeneities proved essential for capturing effluent quality from different iv material zones, and for conceptualizing internal flow paths occurring within the dry cover system in place. In addition, results from this research demonstrate the capabilities of the recently developed unstructured grid code MIN3P-HPC for capturing complex shapes and geometries, which may prove useful for future investigations of capillary barrier effects and other design strategies in WRPs. v Lay Summary In order to extract valuable resources such as metals and gems, mines must remove massive amounts of non-valuable rock (referred to as “waste rock”) to reach them. These piles are often stored above ground and may be up to 500,000,000 m3 of rock (a volume which could fill 200,000 Olympic-sized swimming pools). Similar to the mechanism of rust forming on iron, when removed from the ground and exposed to oxygen in the air, waste rock may break down and release potentially harmful metals, which require water treatment. This research studies storage design strategies for waste rock piles in order to minimize contaminant release and water treatment requirements. In this thesis, computer models are used to conceptualize what is happening within these piles to improve long-term waste rock storage design strategies. vi Preface This thesis consists of four chapters with intent to submit Chapters 2 and 3 to internationally peer-reviewed journals at the time of writing. Supplementary information for each chapter is found in the Appendices of this thesis. I was responsible for developing the scope of this thesis together with Mayer KU. I set up and carried out all computer simulations with guidance from Su D and Seigneur N, with some additional help initially from Modzelewski H regarding high performance computing on the UBC EOAS Clusters. I wrote all chapters in this document with edits provided by Seigneur N and Mayer KU. Chapter 2 is a parameter sensitivity analysis using the reactive transport modeling code MIN3P-HPC to investigate the influence of physical and chemical heterogeneities combined with structural trends and scale on the long-term behaviour of mine waste stockpiles built using different construction methods (specifically, push- and end-dumping). Stylistic models of waste rock piles were run for 200 years and compared to simulate differences in long-term effluent quality behaviour under multiple scenarios. As mentioned, I set up and executed simulations with assistance from Su D and Seigneur N. Additionally, Seigneur N provided MATLAB codes which were used in data presentation. I interpreted results with assistance from Mayer KU and Seigneur N. Chapter 3 involved collaboration between myself, Poaty B, Plante B, Mayer KU, and Bussière B to numerically model a laboratory experiment conducted by Poaty (2020) at the Université du Québec en Abitibi-Témiscamingue. The work focused on the assessment of a proposed engineered cover system for a full-scale waste rock pile under controlled conditions. For this chapter, lab vii measurements used in the modeling were completed by Poaty B, with additional information provided from previous characterization of the Tio waste rock by Plante B and Bussière B. An additional steady-state flow experiment was completed specifically for this project by Poaty B, which was used to calibrate the numerical model. I set up and calibrated the numerical model with assistance from Su D and Seigneur N. An additional feature for the MIN3P-HPC code to observe effluent chemistry at multiple drainage ports was developed for this study by Su D. Using this approach, I was able to calibrate the numerical model to match the outflow and chemistry results provided from the laboratory model. Interpretation of simulated results and discussion were led by me with support from Mayer KU, Seigneur N, Plante B, Poaty B, and Bussière B. viii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ....................................................................................................................... viii List of Tables .............................................................................................................................. xiii List of Figures ............................................................................................................................. xiv List of Abbreviations and Symbols ............................................................................................xv Acknowledgements .................................................................................................................... xvi Dedication ................................................................................................................................. xviii Chapter 1: Introduction ................................................................................................................1 1.1 Extraction and Storage of Waste Rock ....................................................................... 2 1.2 Reactive Transport Modeling ..................................................................................... 4 1.2.1 MIN3P-HPC ........................................................................................................... 5 1.2.2 Using RTM to Investigate Mine Waste Problems .................................................. 5 1.3 Thesis Objectives and Organization ........................................................................... 6 Chapter 2: Numerical Modeling of Mass Loadings from Waste Rock Piles Affected by Heterogeneity and Structure .........................................................................................................8 Abstract ....................................................................................................................................... 8 2.1 Introduction ................................................................................................................. 9 2.2 Pile Heterogeneities, Structure, and Scale Overview ............................................... 10 2.2.1 Physical Heterogeneity ......................................................................................... 11 ix 2.2.2 Chemical Heterogeneity ........................................................................................ 12 2.2.3 Pile Structure ......................................................................................................... 12 2.2.3.1 End-dumping ................................................................................................ 13 2.2.3.2 Push-dumping ............................................................................................... 14 2.2.4 Pile Scale ............................................................................................................... 15 2.3 Methods..................................................................................................................... 16 2.3.1 MIN3P Model ....................................................................................................... 16 2.3.2 Domain Geometry, Boundary and Initial Conditions ........................................... 17 2.3.3 Parameter Heterogeneity ....................................................................................... 18 2.3.3.1 Physical Heterogeneity ................................................................................. 19 2.3.3.1.1 Hydraulic Conductivity ......................................................................... 19 2.3.3.1.2 Porosity .................................................................................................. 20 2.3.3.2 Chemical Heterogeneity ............................................................................... 21 2.3.3.2.1 Reactive minerals .................................................................................. 21 2.3.3.2.2 Heterogeneity in NP:AP Ratios ............................................................. 21 2.3.3.2.3 Pyrite Oxidation Rate ............................................................................ 22 2.3.4 Pile Structure ......................................................................................................... 26 2.3.5 Pile Scale ............................................................................................................... 26 2.3.6 Assumptions and Limitations ............................................................................... 27 2.4 Simulation Results .................................................................................................... 30 2.4.1 Sensitivity to Parameter Heterogeneity ................................................................ 30 2.4.1.1 Physical Heterogeneity ................................................................................. 33 2.4.1.2 Chemical Heterogeneity ............................................................................... 33 x 2.4.2 Pile Structure ......................................................................................................... 33 2.4.2.1 Random vs. Fining Upwards Trends ............................................................ 33 2.4.2.2 Traffic Surfaces and Increased Structural Complexity ................................ 35 2.4.3 Pile Scale ............................................................................................................... 39 2.5 Discussion ................................................................................................................. 40 2.5.1 Parameter Sensitivity ............................................................................................ 40 2.5.2 Pile Structure and Scale ........................................................................................ 41 2.5.3 Heterogeneity and Sampling Procedures .............................................................. 42 2.6 Conclusions ............................................................................................................... 43 Chapter 3: Numerical Modeling of a Laboratory-Scale Waste Rock Pile Featuring an Engineered Cover System ...........................................................................................................46 Abstract ..................................................................................................................................... 46 3.1 Introduction ............................................................................................................... 47 3.1.1 Site Background .................................................................................................... 48 3.2 Methods..................................................................................................................... 50 3.2.1 Laboratory Methods .............................................................................................. 50 3.2.1.1 System Overview and Construction ............................................................. 50 3.2.1.2 Flow Measurements ..................................................................................... 51 3.2.1.3 Geochemical Monitoring .............................................................................. 52 3.2.2 Numerical Methods ............................................................................................... 52 3.2.2.1 Reactive Transport Code MIN3P-HPC ........................................................ 52 3.2.2.2 Flow Model and Calibration ........................................................................ 53 3.2.2.3 Reactive Transport Model and Calibration for Leachate Concentrations .... 56 xi 3.2.3 Assumptions and Limitations ............................................................................... 59 3.3 Results ....................................................................................................................... 60 3.3.1 Drainage Rates and Internal Fluid Flow Pathways ............................................... 60 3.3.2 Solute Concentrations in Drainage ....................................................................... 64 3.4 Discussion ................................................................................................................. 68 3.4.1 Implications of Results from Flow Modeling ....................................................... 68 3.4.2 Implications of Results from Reactive Transport Modeling ................................. 69 3.5 Conclusions ............................................................................................................... 72 Chapter 4: Conclusion .................................................................................................................74 4.1 Limitations of RTM .................................................................................................. 76 4.2 Recommended Future Work ..................................................................................... 78 4.3 Implications for Future Use of RTM in Mine Waste Management .......................... 78 Bibliography .................................................................................................................................81 Appendices ....................................................................................................................................98 Appendix A Supplementary Material for Chapter 2 ..................................................................... 98 A.1 Domain Geometry, Boundary and Initial Conditions ........................................... 98 A.2 Parameter Heterogeneity ....................................................................................... 98 A.3 Heterogeneity in NP:AP Ratios .......................................................................... 102 A.4 Pile Structure ....................................................................................................... 102 A.5 Soil Hydraulic Function Parameters ................................................................... 103 A.6 Total vs. Local Drainage Quality ........................................................................ 106 Appendix B Supplementary Material for Chapter 3 ................................................................... 111 B.1 System Overview and Construction.................................................................... 111 xii B.2 Steady-State Flow Experiment ........................................................................... 112 B.3 Grain Size Distribution Curves ........................................................................... 115 B.4 Numerical Model Inputs ..................................................................................... 121 xiii List of Tables Table 2.1 Structural features for end- and push-dumping construction methods ......................... 14 Table 2.2. Mineral rate expressions for stylistic models .............................................................. 25 Table 2.3. Variation in maximum mass loading rates and reduced sulfate release times with parameter heterogeneity .......................................................................................................... 32 Table 2.4. Variation in mass loading rates and duration of release with traffic surfaces and structure ................................................................................................................................................. 39 Table 3.1. Physical parameters for variably saturated flow separated by zone ............................ 55 Table 3.2. Composition of infiltrating water applied to the model ............................................... 56 Table 3.3. Initial mineral volume fractions considered in numerical model ................................ 57 Table 3.4. Minerals and their dissolution-precipitation reactions considered in numerical model. ................................................................................................................................................. 58 Table 3.5. Results for measured versus simulated drainage rates ................................................. 62 Table 3.6. Calibrated effective rate coefficients for mineral weathering reactions ...................... 68 xiv List of Figures Figure 1.1. Construction of an experimental pile with observable fining upwards along sloped surfaces at Lac Tio Mine, QC, Canada. .................................................................................... 4 Figure 2.1. Illustrated conceptual models of push- and end-dumping construction methods in a WRP ........................................................................................................................................ 15 Figure 2.2. Conceptual model for single bench, 2D WRP and triangular mesh ........................... 18 Figure 2.3. Cross-sections showing the initial conditions for the single bench models looking with spatial distributions of hydraulic conductivity, pyrite volume fraction, and effective rate coefficients for pyrite within the pile for various scenarios ................................................... 24 Figure 2.4. Hydraulic conductivity distribution for multiple bench models in cross-section ....... 29 Figure 2.5. Comparison of the effects of heterogeneity in basal sulfate mass loading rates, influenced by physical and chemical heterogeneities ............................................................. 31 Figure 2.6. Sulfate basal drainage evolutions shown for varied pile structure and scale ............. 37 Figure 2.7 pH cross sections at t = 3 years comparing heterogeneous versus structure models .. 38 Figure 3.1.Schematic of laboratory-scale waste rock pile adapted from Poaty (2020). ............... 51 Figure 3.2. Conceptual model with labeled boundary conditions and material property zones ... 55 Figure 3.3. Measured vs. simulated cumulative drainage ............................................................. 61 Figure 3.4. Simulated results of calibrated flow model ................................................................ 63 Figure 3.5. Results of simulated vs. observed effluent concentrations in mol L-1, color-code per reservoir .................................................................................................................................. 67 xv List of Abbreviations and Symbols ARD Acid rock drainage CBE Capillary barrier effect CND Contaminated neutral drainage FCL Fluid control layer FD Frequency domain HI Hematite-ilmenite ICP-MS Inductively coupled plasma mass spectrometry ML Metal leaching RTM Reactive transport modeling WRP Waste rock pile wt. % Weight percent 𝑑𝑑𝑚𝑚 Mean grain size diameter IAP Ion activity product K Hydraulic conductivity K� Mean hydraulic conductivity 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒 Effective rate constant K𝑚𝑚/K25 Equilibrium constant 𝑆𝑆𝑤𝑤𝑤𝑤 Residual saturation Ψ Pressure head α𝐾𝐾 Scale factor for hydraulic conductivity αφ Scale factor for volume fraction θ Porosity φ Volume fraction α van Genuchten air entry suction term 𝑛𝑛 van Genuchten pore size distribution term xvi Acknowledgements First and foremost, thank you to Uli Mayer for being an incredible supervisor and mentor. I gained so much more than I could have expected from this degree because of your considerate approach and dedication to improving society through research—whether it’s in our weekly meetings involving hand calculations, detailed edits for abstracts or full theses, advising roles within the Department, or the projects selected. Thank you for pushing me to work harder when I needed to, and for understanding when I couldn’t. I could not imagine a better advisor. Thank you infinitely to Danyang Su, who always had patience and bountiful knowledge for the embarrassing volumes of emails I sent asking for help with MIN3P. Even more impressive was the mastery of your additions to the code which were essential for attaining many of the results presented in this work, and for the speed at which the models ran. The accomplishments shown in this work could not have been done without you. Another trillion thank-you’s (and sangria grandes) to Nicolas Seigneur, for the Matlab and math wizardry bestowed, and all the advice and edits on my first, second, third, and fourth drafts. I would likely still be plotting results rather than writing this acknowledgement without your help. Thank you to several other professors and researchers who made this work possible: to my committee member Bruno Bussière and also to the UQAT team, Benoît Plante and Bissé Poaty, for providing me with the information and guidance with the work completed in Chapter 3 and for the hospitality (and Gary) you gave me during my exchange to Rouyn-Noranda. Thank you also to my other committee member, Roger Beckie, for providing guidance and challenging me to improve my understanding of hydrogeology through course and thesis work. xvii This research was funded by several grants from the Natural Sciences and Engineering Research Council of Canada (NSERC) including the Toward Environmentally Responsible Resource Extraction Network (TERRE-NET) and the Alexander Graham Bell Canada Graduate Scholarship-Master’s (CGS M). Thank you to especially to the organizers of TERRE-NET: David Blowes, Carol Ptacek, and Steve Holland for funding this research and providing a connected network of students, with whom I have made treasured friends and memories with. The opportunity to get to know other students interested in mine wastes while learning industry skills and networking with current and future professionals was invaluable to my experience. The grant also funded several opportunities for me to attend conferences and present my research, which was an opportunity especially beneficial for progressing this work and my career. A special thank-you to Carol, for first convincing me to switch into hydrogeology as an undergraduate student, and to Dave for advising through my undergraduate career and beyond. The mentorship and opportunities you have both given me (from co-op jobs to this thesis project) have truly influenced the course of my life. Thank-you also to David Wilson, for being a kindred spirit and model modeler. Our conversations about MIN3P and waste rock were always thought provoking and gave new perspectives and excitement for this research. I value all the professional and personal support; I hope we can build some models together soon. Last but never least, thank you forever to my friends and family who have been there to support and encourage me always. Many of you knew I would be doing a Master’s even before I did. To acknowledge all the support I’ve received over the last few years on a personal front would be far longer than this thesis, but was nonetheless essential. I am so lucky to have so many special people in my life; the joy from these shared moments is a cherished and constant source of motivation to keep going. xviii Dedication This thesis is dedicated to my parents, Frank and Mary, who have shown nothing but patience and love to me through all my endeavors—even when they change course ten minutes later. It is through your unconditional support that I have been allowed to be fully creative and curious. Through your love, I am driven to share it with the world. Love you more—ok, same. 1 Chapter 1: Introduction The release of metals from mine waste materials through acid rock drainage (ARD) or contaminated neutral drainage (CND) is the greatest environmental liability facing the mining industry today (Tremblay and Hogan 2001; Canadian Mining Innovation Council 2013), with financial implications for water treatment. While negative environmental effects of acidic drainage from reactive sulfides have been observed as early as during Roman times (INAP 2012), research in strategies for prevention and remediation of metal leaching (ML) from mine wastes is relatively recent. Although significant research over the last several decades has advanced understanding of the important relationships between the production of ARD/CND and factors including sulfide content, reactive surface area, mineral weathering rate evolution, pH-buffering capacity, and mobilization/retardation mechanisms of metals (Sherlock et al. 1995; Blowes et al. 2003; Sapsford et al. 2009; Langman et al. 2014; Erguler and Kalyoncu Erguler 2015; Nordstrom et al. 2015), there are still significant uncertainties when forecasting long-term drainage behaviour, and preventing or remediating poor drainage quality. Part of these uncertainties can be attributed to the complexity and uniqueness that each mine site has with respect to factors such as regional geology, climate, and scale of operations, combined with the often limited data points available to characterize mine waste and operational costs associated with best practice long-term closure strategies. This thesis focuses on the use of reactive transport modeling (RTM) as a tool to help characterize the processes controlling ARD release from waste rock piles (WRP) and to compare the effect of different mine waste storage techniques on their mass loadings. RTM allows the modeler to utilize data obtained from standard measurements typically collected during waste rock material characterization into a numerical model to help investigate internal processes that may be occurring 2 or identify areas where additional data collection may be required (in particular for cases when discrepancies between expected behaviour and model outputs exist). As mentioned in the Preface, this thesis contains two research chapters in preparation for individual peer-reviewed journal submissions. As such, each chapter contains its own introduction and relevant background information tailored to the material presented. However, a brief introduction to waste rock, ML/ARD/CND processes, and presently available tools in RTM used in this thesis to understand and mitigate potential risks from these processes are provided below. 1.1 Extraction and Storage of Waste Rock During the mining process, large volumes of material that are considered not economically valuable (referred to as “waste rock”) must be removed in order to extract the valuable commodities. Once this material has been blasted and removed from the subsurface, it is often stored at ground surface in large stockpiles. WRPs may be composed of more than 500 million cubic metres of material, and range in grain size from clay-sized particles to automobile-sized boulders (McCarter 1990). Piles may also be quite variable in mineralogical composition depending on site-specific material. The overall heterogeneity developed during WRP construction includes physical and chemical properties such as grain size, porosity, hydraulic conductivity, mineralogical composition, and reactive surface area (Fala et al. 2003; Tran et al. 2003; Akcil and Koldas 2006; Stockwell et al. 2006; Amos et al. 2015; Wilson et al. 2018a). In addition, the deposition and construction methods of a pile will influence the spatial distribution of these variable physical and chemical heterogeneities (e.g. Figure 1.1). RTM has previously been used as a tool to investigate the influence of these heterogeneities on the long-term hydrogeological behaviour of unsaturated flow and ARD generation and ML within a WRP (e.g. Fala et al. 2003, 2005, 2006; Blackmore et al. 2014; Lahmira et al. 2016, 2017; Appels et al. 2018; Wilson et al. 2018a,b). A more complete investigation of the different physical and chemical heterogeneities and their influence on ARD/ML is provided in Chapter 2. 3 When removed from the subsurface and exposed to atmospheric conditions, sulfide minerals within the waste rock may oxidize and generate acidity, causing additional mineral weathering and leaching of metals from waste rock (Nicholson et al. 1988; Blowes et al. 2003). Therefore, many prevention or reclamation strategies for mine wastes aim to reduce the exposure of sulfides and reactive minerals to atmospheric oxygen (which fuels sulfide oxidation) or water stemming from precipitation (which transports metals and acidity through and out of a waste rock pile) (Price and Errington 1998; Tremblay and Hogan 2001; Haug 2002; MEND 2004a; Johnson and Hallberg 2005; INAP 2009). Effluent quality produced from a pile depends on multiple factors such as sulfide content, oxidation rates, presence of pH-buffering minerals, and fluid flow pathways; however, many of these processes that occur within a WRP are difficult to measure and interpret, given the scale and complexity of WRPs. To this end, RTM allows for an integrated assessment of these complex physical and hydrogeochemical processes, which may facilitate the development of an improved understanding of the governing processes that control metal leaching. In addition, RTM can be used to validate or disprove hypotheses about dominant processes, as well as characterize and quantify processes occurring within WRPs (Lefebvre et al. 2001; Steefel et al. 2005; Fala et al. 2006). 4 Figure 1.1. Construction of an experimental pile with observable fining upwards along sloped surfaces at Lac Tio Mine, QC, Canada. Photo provided by Bruno Bussière, 2018. 1.2 Reactive Transport Modeling RTM combines transport equations and chemical processes to model a variety of complex hydrogeochemical mechanisms in the subsurface. These RTM codes may be used to build and contrast conceptual models of a system, explore “what-if” scenarios, aid in data interpretation, and compare different storage and management techniques for environmental problems (Mayer et al. 2019). When calibrated using laboratory and field data, RTM may also be used as a tool for forecasting and bracket expected behaviour for a system—although characterization of the model and interpretation of results must be approached with prudence. RTM codes are tested through model verification and validation. Model verification involves ensuring governing equations used in the RTM code are solved correctly, which is evaluated by comparing simulated results from different codes and by comparing results to those of established benchmark studies (Mayer et al. 2015; Steefel et al. 2015b). Model validation involves ensuring that 5 the conceptual model is representative of processes occurring in the field or laboratory by comparing simulation outputs to measured results (e.g. Fala et al. 2013; Blackmore et al. 2014, 2018a; Lahmira et al. 2017; Pedretti et al. 2017; Wilson et al. 2018a). This thesis used the RTM code MIN3P-HPC (Su et al. submitted) for investigations of mine waste storage with a specific focus on aspects of heterogeneity, construction methods, scale, and the design of an engineered cover system. A brief summary of the code capabilities is provided below. 1.2.1 MIN3P-HPC MIN3P-HPC (Su et al. submitted) is a multi-component RTM code for variably saturated porous media. Initially developed as MIN3P by Mayer et al. (2002), several features and upgrades have been made over the last two decades. MIN3P-HPC is an upgrade to the parallelized version of MIN3P-THCm (Bea et al. 2016; Su et al. 2017), which includes the ability to use structured or unstructured grids in one to three spatial dimensions, allowing for better representation of complex geometries (such as the sloped surfaces and angles within WRPs, as presented in this thesis). The code couples transport and complex geochemical reactions using a global implicit solution approach, enforcing global mass balance between species in the gas, water and solid phases (Su et al. 2017). In addition, governing equations are discretized using a conservative finite volume method. Details regarding numerical methods, instructions for use, and applications of MIN3P-HPC can be found in Xie et al. (2015) and Su et al. (2018). 1.2.2 Using RTM to Investigate Mine Waste Problems In order to characterize mine waste and assess long-term mass loadings out of the pile, often laboratory- or field-scale experiments including humidity cell tests, column tests, or field barrel tests are completed, with results then often scaled up to full scale WRPs. However, many times these scale-up procedures rely on empirical scale factors, which have been found insufficient to accurately predict long term drainage quality (Parbhakar-Fox and Lottermoser 2015; Amos et al. 2015; Wilson 6 et al. 2018b). Furthermore, these smaller-scale tests are often unable to capture some of the structural features of a WRP, such as fining upwards trends in grain size or traffic surfaces (discussed further in Chapter 2). Supplemented by laboratory- and field-scale measurements, RTM can be used to take a “bigger picture” approach, integrating aspects of scale, spatial trends and structure (Steefel et al. 2005; Šimůnek et al. 2008). Since the exact structure and composition at each point in the pile cannot be known, well-considered simplifications and assumptions must be made to capture the dominant processes occurring within a WRP (Steefel et al. 2015a). These assumptions manifest themselves as parametrizations of the hydrogeochemical processes considered, dimensionality, and methods for solving flow (including uniform flow, multi-domain, or discrete approaches) and transport. Similar to laboratory- and field-based tests, it is important to understand assumptions and limitations of numerical modeling tools. For example, the idea of a “non-unique” solution is important when interpreting numerical model results, referring to a solution that can be reached through multiple different combined processes, though not all are considered valid based on pre-existing knowledge and system understanding. Provided a considerate approach is taken and supported by real world measurements, RTM is a powerful tool to assist professionals in developing an understanding of the system and for scale-up and structural considerations for mine waste storage design. 1.3 Thesis Objectives and Organization This thesis is divided into two research chapters which address WRP storage design strategies using the RTM code MIN3P-HPC. Work was driven by the following research questions: i. How does physical and chemical parameter heterogeneity in a WRP control the overall effluent quality from a pile and its evolution over time, when compared to a homogeneous case? 7 ii. What differences are observed for the geochemical evolution in a WRP and effluent quality as a function of construction method (specifically, push-dumping vs. end-dumping), internal structure (e.g. traffic surfaces, fining upwards trends), and scale (e.g. single bench vs. multi bench WRPs)? iii. What influences do pile height and shape have on the overall drainage quality? iv. Can RTM methods reproduce drainage and effluent quality observed in a laboratory model WRP, designed with an inclined cover system? If so, can the model be used to shed light on processes occurring with the WRP, which cannot be easily measured? Finally, can the model be used to assess how successful the cover is at diverting water away from potentially metal leaching materials? While investigating these research questions, this work also aimed at testing the capabilities of the unstructured grid in the recently developed MIN3P-HPC code. Scenario analyses in Chapter 2 address questions i – iii using stylistic models of single bench and multi-bench WRPs with spatially varied physical and chemical parameters, arranged with spatial trends of random and/or fining upwards grain size distributions and traffic surfaces to capture the influence of construction methods. In Chapter 3, a numerical model was calibrated to a laboratory-scale model constructed with an engineered cover system as per the design plan for the full-scale WRP at the Tio mine. This chapter additionally explores model validation by reproducing observations from a real-world example to results simulated using MIN3P-HPC. While investigating these research questions, this thesis also demonstrates capabilities of RTM as a tool to conceptualize these complex systems and quantitatively interpret observational data. 8 Chapter 2: Numerical Modeling of Mass Loadings from Waste Rock Piles Affected by Heterogeneity and Structure Abstract Placement methods and material availability during waste rock pile (WRP) construction may create significant heterogeneities in physical and geochemical parameters (such as grain size, permeability, mineralogy and reactivity) and influence the internal pile structure. Due to their enormous scales, it is difficult to capture the influence of WRP heterogeneities using laboratory- or field-scale studies. This study uses a numerical modeling approach to investigate the influence of physical and chemical heterogeneities, structure, and scale on acid rock drainage (ARD) and attenuation through 2D reactive transport simulations. Specifically, the sensitivity of drainage quality to parameters including grain size distribution, sulfide mineral weathering rates, abundance and distribution of primary minerals, and pile structure as a function of construction methods are investigated. The geochemical model includes sulfide oxidation, pH buffering by calcite dissolution, and ferrihydrite and gypsum as secondary phases. Simulation results indicate that the implications of heterogeneity and construction method are scale-dependent; when grain size distribution trends observed in a pile’s core are applied to the entirety of a pile, results between push- and end-dumping methods vary substantially—however, results between different construction methods become more similar when features such as traffic surfaces, structural variation, and multiple benches are also considered. For all scales and construction methods considered, simulated results demonstrate that pile heterogeneity and structure decrease peak mass loading rates 2 to 3-fold but cause prolonged ARD release compared to the homogeneous case. These findings have implications for the economics of planning water treatment facilities for life of mine and closure operations. 9 2.1 Introduction Numerous industries globally depend on a reliable supply of commodities produced from mining activities. During the mining process, large amounts of materials that are not economically viable (called waste rock) must first be removed to access the valuable commodities below. After drilling, blasting, and removing the waste rock, it must be stored long-term—often in large piles above ground, where materials are exposed to atmospheric weathering conditions. A waste rock pile (WRP) may be composed of more than 500 million m3 of material, with a wide range of grain sizes (McCarter 1990). Due to the enormous volumes of materials that constitute them, heavy equipment such as haul trucks and bull dozers are required for pile construction. Materials deposited from these large machines result in random, heterogeneous waste rock placement, though on a larger-scale may be strategically mixed or confined to prescribed zones, based on rock type for closure plan strategies (e.g. MEND 1998; Miller et al. 2006). Construction methods and source of material thus create significant variations in the distribution of grain size and mineralogy (with implications for sulfide content, metal concentrations, and pH-buffering capacity), which affects drainage quality (Smith et al. 1995; Munroe et al. 1999; Fala et al. 2003; Tran et al. 2003; Herasymuik et al. 2006; Stockwell et al. 2006; Lahmira et al. 2017). While the relationships between these characteristics in a pile and the potential to generate acid rock drainage (ARD), related metal release, and attenuation have been studied extensively (e.g. Sherlock et al. 1995; Stumm and Morgan 1996; Nordstrom and Alpers 1999; Akcil and Koldas 2006; Amos et al. 2015; Dold 2017), the influence of heterogeneity and spatial distribution of these characteristics on ARD has received less attention to date. This is heavily due to the large scales of WRPs and typically insufficient sampling resolution to fully capture a pile’s heterogeneity. Furthermore, typical laboratory-based methods (such as humidity-cell or field-cell experiments) do not capture additional structural features such as fining upwards or traffic surfaces. To this end, numerical models have proven suitable tools for investigating the influence of physical 10 and chemical heterogeneities (Molson et al. 2005; Fala et al. 2013; Lahmira et al. 2017; Pedretti et al. 2017; Appels et al. 2018), as well as scale (Wilson et al. 2018a) on ARD generation and attenuation processes in WRPs. However, a sensitivity analysis that examines the roles of physical and chemical parameter distribution, with a particular emphasis on the effect of pile construction methods, structural features, and pile scale, has not been completed. This paper investigates the impacts of physical and chemical heterogeneities, scale and pile-internal structure on drainage quality from a WRP using the reactive transport code MIN3P-HPC (Su et al. submitted). First, an overview of parameters controlling the development of ARD are discussed. MIN3P-HPC model outputs from simulated 2D cross sections of WRPs are then analyzed to explore the combined impacts of parameter heterogeneity and internal pile structures on basal drainage quality, as a function of scale. Correlations between spatial distributions of controlling parameters due to grain size (including hydraulic conductivity, porosity, mineral abundance, and weathering rate) as well as the effects of structural features (i.e. grain size distribution and segregation due to construction methods and traffic surfaces) are examined. The influence of pile scale is then explored by comparing results from a simulated single bench WRP to a pile with multiple benches and a larger volume of material. Simulations of heterogeneous WRPs are compared to a reference case with average parameter values assumed uniform throughout the pile. Results from these simulations are discussed to evaluate effects on temporal and spatial development of ARD due to physical and chemical heterogeneities, structure, and pile scale, and their combined effects on peak mass loadings and duration of ARD release. 2.2 Pile Heterogeneities, Structure, and Scale Overview WRPs exhibit large ranges in grain size, from clay-sized particles to automobile-sized boulders (McCarter 1990; Morin et al. 1991; Smith et al. 1995; Stockwell et al. 2006; Anterrieu et al. 2010). While significant differences in physical and chemical parameters may further be characterized by a 11 pile’s source rock, the rock type itself may be quite variable, and its defining attributes are often site-specific. To begin conceptualizing and quantifying the influence of parameter heterogeneity on the evolution of ARD and metal leaching within a pile, this study will use mean grain size diameter as a proxy. Firstly, the relationships between grain size and the parameters investigated (including hydraulic conductivity, porosity, mineralogical abundance and weathering rate) are described. Following these defined relationships, observations made regarding trends in grain size and their spatial distributions due to construction methods (specifically, end- versus push-dumping) are also explored. This investigation considers the presence of traffic surfaces (which are characterized by smaller grain size and lower porosity due to compaction from heavy equipment driving along exposed weathering surfaces during the construction process) and local patterns in grain size distributions such as random or fining upwards trends, which vary between construction methods. In addition, WRPs are often built in multiple benches, repeating these structural distributions within each bench as they grow in height. Pile scale is investigated in this study by comparing results from simulations in single bench models to larger piles with multiple benches, changing the total amount of material and pile shape. The relationships between grain size distribution and heterogeneous physical and chemical parameters, structural features, and scale provide the foundation for the stylistic models developed for this research. 2.2.1 Physical Heterogeneity In this study, the physical heterogeneities influenced by grain size include hydraulic conductivity and porosity. Hydraulic conductivity may vary by several orders of magnitude in a pile and typically follows a lognormal distribution (Smith et al. 1995; Tietje and Hennings 1996). Variations in hydraulic conductivity values throughout the pile are expected to affect contaminant residence times and migration pathways (Jung and Navarre-Sitchler 2018). Porosities throughout WRPs may also be extremely variable (Smith et al. 1995; Lefebvre et al. 2001; Azam et al. 2007; Zhang 2017; Dimech 12 et al. 2018). In this study, porosity is assumed to have a log-linear relationship with hydraulic conductivity, following the work of Smith et al. (1995). During pile construction, infilling with fines and compaction along traffic surfaces may reduce both hydraulic conductivity and porosity, as observed in previous studies such as Bellehumeur (2001) and Cash (2014), effects of which are also considered in our study. 2.2.2 Chemical Heterogeneity WRPs may exhibit highly heterogeneous mineral abundances depending on the geology of the deposit and waste rock source material. Even if a pile is considered relatively homogeneous in terms of rock type, naturally-occurring variations in mineralogy can lead to the development of pockets of higher reactivity and acidity (e.g. Atherton 2017; Pedretti et al. 2017). To investigate the influence of chemical heterogeneities on ARD-related processes, we consider here local variations in NP:AP ratios (where NP is the neutralization potential from pH-buffering minerals, and AP is the acid production potential from sulfide minerals) and sulfide weathering rates. The availability of sulfide minerals and associated oxidation rates have been shown to depend on particle size. For example, several waste rock studies have found that sulfide mineral availability is higher in the more fine-grained materials (Aranda 2010; Erguler and Kalyoncu Erguler 2015; Elghali et al. 2019). Furthermore, the reactivity of sulfide minerals is controlled by mineral properties and reactive surface areas, with an inverse relationship between pyrite oxidation rate and particle diameter (Nicholson et al. 1988; Komnitsas et al. 1995; Cornell and Schwertmann 2003; Blowes et al. 2003). 2.2.3 Pile Structure Rock placement and construction of a WRP controls the internal structure, influencing and overprinting the physical and chemical heterogeneities mentioned previously. While a combination of construction methods may be used within a single pile, for the purpose of this study, differences between classically described features of push- and end-dumping construction methods were 13 investigated. Features of these construction methods are summarized in Table 2.1 and illustrated in Figure 2.1. The investigation of engineered waste rock storage techniques (such as covers, capillary barriers or the placement of more acid generating rock within non-acid generating cover material) to mitigate ARD generation (Yanful and Woyshner 1993; O’Kane et al. 1998; Haug 2002; MEND 2004a; Fala et al. 2005; Amos et al. 2009) was beyond the scope of this study. 2.2.3.1 End-dumping In end-dumping methods, rock is typically dumped over a sloped pile face via haul truck. This technique results in observable effects of gravity sorting: larger particles and boulders are driven by momentum towards the toe of the pile, whereas significantly more of the finer materials remain close to the surface (Nichols 1986; Morin et al. 1991; Smith et al. 1995; Smith and Beckie 2003; Stockwell et al. 2006; Wilson 2011, 2017; Amos et al. 2015). Nichols (1986) described three distinct horizons in an end-dumped pile: (1) the upper horizon, which contains a high concentration of fine materials; (2) an evenly distributed/evenly graded horizon of fine and coarse particles along the remainder of the slope up to the toe; and (3) a bottom horizon beyond the slope’s toe containing a high concentration of coarse fractioned materials and boulders. The upper horizon contains a higher portion of fines due to gravity sorting, which is further compacted at the surface due to heavy equipment traffic as the pile grows laterally. Relative proportions of upwards fining trends throughout the pile are consistent regardless of the height of the pile bench; however, the degree of particle size segregation depends on a pile’s slope length and origin of material (since rock type will affect influencing factors such as density, strength, and shape). These structural trends (illustrated in Figure 2.1) have been observed in several WRPs (Herasymuik et al. 2006; Stockwell et al. 2006; Azam et al. 2007; McLemore et al. 2009; Barsi 2017; Zhang 2017; Vriens et al. 2019). 14 2.2.3.2 Push-dumping In push-dumping methods, rock is dumped onto the surface of a pile, followed by pushing and levelling out materials using tractors and shovels. As a result, fine-grained and sub-horizontal compacted layers occur in the internal structure due to operating equipment along traffic surfaces as the pile grows (Lefebvre et al. 2001; Fala et al. 2003; Anterrieu et al. 2010). Coarse fractions are often embedded in a collection of fine aggregates at the slope crest, while being pushed along the pile surface. As a result, a random distribution (without observable fining upwards trends) within the core occurs, overlain by traffic surfaces (Anterrieu et al. 2010). Materials pushed over the slope crest may experience the effects of gravity sorting, which may cause fining upwards along the batters similar to end-dumping, but spatially more restricted. These features characteristic to push-dumping construction methods have been observed by Fines (2006), Poisson et al. (2009) and Mele et al. (2013). Table 2.1 Structural features for end- and push-dumping construction methods Construction Method Description/Features End-dumping  Internal structure of inclined layers along angle of repose (Wilson 2017)  Single traffic surface at the top of the pile (bench)  Greater segregation due to gravity sorting throughout pile  Upper 10-15% contains highest concentration of fines (Nichols 1986) Push-dumping  Higher angled slopes  Internal sub-horizontal layers from multiple traffic surfaces  Less segregation due to gravity sorting and limited to batters; larger materials are deposited randomly throughout pile as they get lodged during pushing and resurfacing (Nichols 1986) 15 Figure 2.1. Illustrated conceptual models of push- and end-dumping construction methods in a WRP. 2.2.4 Pile Scale Internal structures defined in the previous section were described in terms of construction of a single bench. However, piles are often built in multiple terraced benches, where internal features are repeated, separated by traffic surfaces as shown in Figure 2.1 (e.g. Piteau Associates Engineering 16 Ltd. 1991; Knight Piésold Ltd. 2012; Dominion Diamond Mines 2014; StrataGold Corporation 2014). Typically, the process of scale-up is not quite as simple with respect to quantifying the differences in drainage quality. Often, laboratory- or field-based experiments such as humidity-cells, column tests, or field-barrels are utilized to characterize mass loadings from mineral leaching rates using empirical scale factors for a full-scale pile, which may be insufficient for accurate water quality prediction (Parbhakar-Fox and Lottermoser 2015; Amos et al. 2015; Wilson et al. 2018a; Vriens et al. 2019). While these issues regarding scale-up have been addressed and improved upon in several studies (e.g. Malmström et al. 2000; Kempton 2012; Plante et al. 2014), these laboratory- and field-based experiments do not consider structural influences such as fining upwards, traffic surfaces, or differences between the core and batters of a pile, and do not fully address the contribution to effluent quality from pile structure and scale. Supported by recent advancements to reactive transport codes and computing capabilities, numerical modeling provides a suitable tool for investigating the influence of pile structure and scale-up related issues and is used in this study. 2.3 Methods 2.3.1 MIN3P Model MIN3P (Mayer et al. 2002) is a three-dimensional multicomponent reactive transport model designed for applications in variably saturated porous media. The unstructured grid version used in this study, MIN3P-HPC (Su et al. submitted), was developed based on the parallel version of the MIN3P-THCm code (Bea et al. 2016; Su et al. 2017). The code uses a polygonal (2D, e.g., triangular, quadrilateral) or polyhedral (3D, e.g., tetrahedral, prismatic, hexahedral) mesh, although a standard regular gridded mesh available in previous versions may also be applied. The versatility of the triangular mesh allows structural features such as sloped batters and deposition following the angle of repose within a WRP to be captured, while previous numerical models of piles using MIN3P were constrained to only the core of a pile due to limitations related the structured grid. 17 MIN3P-HPC utilizes a node-centered finite volume approach. Values are assigned to each node and then interpolated for the interfaces of adjacent cells, which defines a representative elemental volume; the degree to which these values are upscaled is therefore controlled by cell size and mesh discretization. Richards’ equation is used to solve flow in the unsaturated/saturated subsurface, while a global implicit solution approach is used to solve the reactive transport equations (Mayer et al. 2002). Additional details regarding the governing equations and numerical methods of MIN3P-HPC are provided in Mayer et al. (2015) and Su et al. (2018), respectively. 2.3.2 Domain Geometry, Boundary and Initial Conditions To determine the influence of parameter heterogeneity on ARD in single bench WRPs, triangular meshes with an average resolution of 1.0 m were generated (Figure 2.2). Each triangular cell represents an average area of approximately 0.4 m2. In addition, models including traffic surfaces and internal features representative of those found within push- and end-dumped piles were also investigated for structure and scale. Additional mesh configurations for single bench models (20 m height and 2D cross-sectional area of 1,456 m2) and multiple bench models (60 m height and 2D cross-sectional area of 10,032 m2) were completed for the models with increased structural complexity. Examples of these meshes are included for the single-bench models accounting for pile structure in Appendix A. All simulations assumed a recharge rate of 300 mm/year applied to the pile surface and slopes, comparable to recharge rates observed in semi-arid northern climates (e.g. Lefebvre et al. 2001; Neuner et al. 2013). A constant head boundary of -0.5 m was assigned to the base of the pile to allow free drainage out of the pile into a basal collection system (not simulated). Gases considered include O2 and CO2, which were held constant at values of 21% and 0.5%, respectively. These values were chosen based on observations at the Diavik Diamond Mine and thought to be representative of a real-world, well aerated pile with low sulfide and carbonate content (Amos et al. 2018). All single bench 18 simulations were run for 300 years under steady-state flow conditions to monitor long-term pile behavior and effluent response to heterogeneities. Figure 2.2. Conceptual model for single bench, 2D WRP and triangular mesh. Input parameters shown here are consistent in all models investigated for the single bench scenario. For multi-bench models, the same infiltration rate was applied along the top and slopes of exposed surfaces, with the same constant head boundary assigned at the base. Simulations assumed O2 and CO2 remained constant at values of 21% and 0.5%, respectively. The triangular mesh was refined to include traffic surfaces and internal structural features for investigations of structure and traffic surfaces and are shown in Appendix A. 2.3.3 Parameter Heterogeneity Heterogeneity was accounted for by relating physical and chemical parameters to a mean grain size diameter assigned to each node. In this way, parameter distributions of hydraulic conductivity, porosity, mineralogical assemblage and sulfide mineral oxidation rates were determined to mimic a range of parameter values observed in real-world piles. For all simulations, an average mean grain size diameter (𝑑𝑑𝑚𝑚) of 5 mm was considered, with local mean grain size diameter (𝑑𝑑𝑚𝑚,𝑖𝑖 , where the subscript 𝑖𝑖 denotes the node) values ranging between 0.05 mm and 170 mm, following a lognormal distribution as observed in previous WRP case studies (e.g. Smith et al. 1995; Tietje and Hennings 1996; Erguler and Kalyoncu Erguler 2015; Zhang 2017). Although other representative grain size diameters (such as 𝑑𝑑5, 𝑑𝑑10, 𝑑𝑑80 etc.) may be used to derive values for parameters investigated in this study (Chapuis 2012; Zhang 2017), mean grain size diameter was 19 used for consistency with the geometrical estimation method for pyrite abundance (described below). Furthermore, these relationships were used to develop a range of representative parameter values for the stylistic models, based on information from multiple mine sites, rather than to calculate absolute values in the context of a site-specific calibration. To identify the influence of varying spatial distributions of mean grain size and the related parameters, five random, lognormal distributions of 𝑑𝑑𝑚𝑚,𝑖𝑖 were compared to five random, lognormal distributions with fining upwards trends. Going forward, these two configurations are referred to as “random” and “fining upwards”. These 𝑑𝑑𝑚𝑚,𝑖𝑖 values were produced using the rnorm() function in Base R (R Documentation n.d.) to generate random variates for the random distributions, or random variates scaled by a elevation to account for the effects of gravity sorting in fining upwards cases. Additional details on these rnorm() functions are included in Appendix A. These models with random and fining upwards trends represent heterogeneous distributions displaying features common to the core of a pile from push- and end-dumping construction methods respectively. Both scenarios were investigated for parameter sensitivity. Additional information such as parameter distributions for mean grain size diameter, hydraulic conductivity, porosity, pyrite volume fraction and rates used in the simulations are also included in Appendix A. 2.3.3.1 Physical Heterogeneity 2.3.3.1.1 Hydraulic Conductivity The relationship between hydraulic conductivity (K) and 𝑑𝑑𝑚𝑚 was developed given the relationship between K and a representative grain size diameter related to the Hazen and Kozeny- Carman equations (Freeze and Cherry 1979) (Equation 2.1). Considering the large spread in grain size diameters typical for waste rock, the representative grain size diameter was assumed to be 𝑑𝑑50 (or 𝑑𝑑𝑚𝑚) using an approach by Freeze and Cherry (1979) for materials with a higher mean grain size 20 than is typical for soils. In this context, the mean hydraulic conductivity was adjusted using a scale factor α𝐾𝐾 to yield a mean hydraulic conductivity of 10-5 m s-1, with local hydraulic conductivity values Ki ranging between 10-9 and 10-1 m s-1. These values were consistent with observations made from WRPs (e.g. Morin et al. 1991; Smith et al. 1995; Azam et al. 2007; McLemore et al. 2009; Fala et al. 2013; Neuner et al. 2013). An example of the initial hydraulic conductivity distributions for random distribution and fining upwards trends are given in Figure 2.3. K𝑖𝑖 = α𝐾𝐾 �𝑑𝑑𝑚𝑚,𝑖𝑖𝑑𝑑𝑚𝑚 �2 (2.1) 2.3.3.1.2 Porosity Porosity (𝜙𝜙𝑖𝑖) throughout the pile was assumed to be related to K𝑖𝑖 for all simulations. For simplicity, a log-linear relationship between hydraulic conductivity and porosity (Equation 2.2) was chosen based on field observations by Smith et al. (1995), which typically demonstrated larger porosities in coarser grained waste rock. A range in porosity values between 0.1 to 0.5 with an average porosity of 0.3 was used based on field observations from WRPs (Morin et al. 1991; Lefebvre et al. 2001; Fala et al. 2003; Neuner et al. 2013; Zhang 2017; Dimech et al. 2018). 𝜙𝜙𝑖𝑖 = 𝜙𝜙𝑚𝑚𝑖𝑖𝑚𝑚 + � 𝜙𝜙𝑚𝑚𝑚𝑚𝑚𝑚 − 𝜙𝜙𝑚𝑚𝑖𝑖𝑚𝑚log10Kmax − log10Kmin� (log10K𝑖𝑖 − log10K𝑚𝑚𝑖𝑖𝑚𝑚 ) (2.2) Models assuming a homogeneous porosity were initially run and compared to results of models with spatially distributed porosity values defined by the relationship in Equation 2.2. Differences in results between cases with a constant porosity and spatially distributed porosities as a function of K were minimal. Since a heterogeneous distribution of porosity influenced by grain size distribution and hydraulic conductivity is considered more realistic of real-world conditions, only results from the heterogeneous porosity cases are presented here. 21 2.3.3.2 Chemical Heterogeneity 2.3.3.2.1 Reactive minerals While a variety of site-specific minerals must be considered when evaluating the long-term behaviour of drainage quality for a WRP, a simplified representation was used in our conceptual analysis to assess the impact of chemical heterogeneity on ARD generation and attenuation. Four minerals were considered in this study: pyrite, calcite, gypsum, and ferrihydrite. This simplified mineralogy allows for conceptualization of the effects of spatially varied NP:AP ratios, sulfide weathering rates, as well as the role of secondary mineral formation. Average initial pyrite and calcite contents of 3.4 and 0.9 wt. % respectively were assumed to achieve potentially acid generating (PAG) conditions. These abundances fall in range with values observed at mine sites (e.g. Lawrence and Wang 1996; Azam et al. 2007; McLemore et al. 2009; Aranda 2010). Gypsum and ferrihydrite were permitted to precipitate to account for the effects of secondary mineral formation on iron and sulfate release. 2.3.3.2.2 Heterogeneity in NP:AP Ratios To investigate the effects of spatial distributions in NP:AP ratios (here, the NP comes from calcite, and the AP comes from pyrite), cases were evaluated against models that varied the spatial distribution of pyrite volume fractions, while calcite remained homogeneously distributed. This variation in NP:AP assessed the implications of local pockets of elevated AP on the overall basal drainage quality. Initial pyrite volume fractions (φ𝑖𝑖0) were considered inversely proportional to 𝑑𝑑𝑚𝑚 (Equation 2.3) to account for an increased sulfide mineral availability in the fines, as observed by Aranda (2010). This relationship was then scaled appropriately using a scale factor (αφ) to maintain an average 3.4 wt. % pyrite in all cases. The lognormal distribution of pyrite volume fractions assigned to each node is given in Appendix A. 22 φ𝑖𝑖0 = αφ 1𝑑𝑑𝑚𝑚,𝑖𝑖 (2.3) 2.3.3.2.3 Pyrite Oxidation Rate In this simplified ARD conceptual model, pyrite dissolution was considered the rate-determining step controlling drainage quality. Pyrite oxidation kinetics may be influenced by multiple factors such as microbial activity, gas or heat transfer, surface area, etc. (Nordstrom 1982; Nicholson et al. 1988; Blowes et al. 2003; Amos et al. 2015). In these numerical simulations, factors influencing oxidative pyrite dissolution rates are captured in the effective rate coefficient 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒. To investigate the influence of weathering rate variability, pyrite oxidation rates were spatially varied, while calcite weathering rates and precipitation/dissolution rates for gypsum and ferrihydrite were considered as quasi-equilibrium reactions (Table 2.2). Local effective rate constants (𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑖𝑖0 ) [mol L-1 bulk s-1] were assumed to be inversely related to grain size using the geometrical estimation method (Aranda 2010) as shown in Equation 2.4. This method considers the relationship between grain size diameter and reactive surface area, which is accounted for the in 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 term. 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑖𝑖0 = α𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒 𝑑𝑑𝑚𝑚𝑑𝑑𝑚𝑚,𝑖𝑖 (2.4) 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑖𝑖0 s were scaled by α𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒 to maintain an average initial rate constant (𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 ) of 1x10-9 mol pyrite L-1 bulk s-1, producing pyrite oxidation rates comparable to those reported in a literature review by Lorca et al. (2016). Total pyrite dissolution rate evolves over time following a two-third power law, modelling a surface-controlled reaction with a progressively decreasing particle radius (Lichtner 1996; Seigneur et al. 2019). Reactions pathways and rate expressions for minerals considered in all simulations are summarized in Table 2.2. Initial distributions of rate coefficients for pyrite oxidation are shown in Figure 2.3. 23 24 Figure 2.3. Cross-sections showing the initial conditions for the single bench models looking with spatial distributions of hydraulic conductivity, pyrite volume fraction, and effective rate coefficients for pyrite within the pile for various scenarios. For the heterogeneous case (random and fining upwards), one of the five configurations that was used to test parameter sensitivity is shown. Porosity is also spatially varied following the relationship to hydraulic conductivity shown in Equation 2.2. Histograms showing the nodal distribution for these parameters are included in Appendix A 25 Table 2.2. Mineral rate expressions for stylistic models. K𝑚𝑚 is the equilibrium constant [-] for mineral m; 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 is the average effective rate constant [mol L-1 bulk s-1]; φi0 is the initial volume fraction [-] and φi is volume fraction updated through time as per the two-thirds dissolution method; IAP𝑚𝑚 is the ion activity product [-] . Subscript 𝑖𝑖 is used to identify values specified for the ith node. Mineral Reaction logK𝑚𝑚 Rate Expression 𝜑𝜑0 [-] 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 [mol L-1 bulk s-1] pyrite FeS2 + H2O + 7/2 O2 → Fe2+ + 2SO42− + 2H+ --- −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 � 𝑑𝑑𝑚𝑚𝑑𝑑𝑚𝑚,𝑖𝑖� �φiφi0�23 0.0127 1×10-9 calcite CaCO3 → Ca2+ + CO32− 8.475 −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 ⎭⎪⎬⎪⎫ × �1 − IAP𝑚𝑚K𝑚𝑚 � 0.00587 5×10-7 gypsum Ca2+ + SO42− + 2H2O → CaSO4 ∙ 2H2O 4.58 −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 0 1×10-7 ferrihydrite Fe(OH)3 + 3H+ → Fe3+ + 3H2O -4.891 −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒0 0 1×10-7 26 2.3.4 Pile Structure The influence of pile structure was explored through the addition of structural complexities including traffic surfaces and internal grain size trends for push- and end-dumping construction methods, as shown in Figure 2.3. In modeled simulations, traffic surfaces were assumed to occur within 0.5 m thick surfaces with a compaction factor that was scaled based on the 𝑑𝑑𝑚𝑚,𝑖𝑖. This scaling follows a non-linear relationship given in Appendix A for hydraulic conductivity, porosity, and sulfide oxidation rates due to the compaction and crushing of materials from heavy equipment (Nichols 1986; Herasymuik et al. 2006; Stockwell et al. 2006; Wilson 2017). However, it was assumed that pyrite volume fractions are dependent on source material during deposition, and therefore were not correlated to smaller 𝑑𝑑𝑚𝑚,𝑖𝑖 values in traffic surfaces. For the “heterogeneous” cases, one of the five random and fining upwards trending piles were chosen for comparison to complex structural models, since they help to describe general grain size trends observed in the core of push- and end-dumped piles, respectively. However, as described in Section 2.2.3, WRPs typically have a combination of these trends and include local traffic surfaces, dependent on the construction method. Therefore, the “structure” models shown in Figure 2.3 are thought to be more realistic based on real-world pile observations. For the push-dumped structure model, multiple traffic surfaces were included at 5 m height increments throughout the core, with random grain size distributions below these compacted surfaces. Fining upwards trends were included along the slopes at the end of each traffic surface, representative of gravity-sorting along the batters. In the end-dumped structure model, multiple zones defined along the angle of repose represent different materials deposited as a pile is constructed (where each zone has a different average 𝑑𝑑𝑚𝑚, though the total average parameter values for the entire pile are maintained), with a single traffic surface at the top of the pile. In comparing the heterogeneous and structure simulations for push- and end-dumping stylistic models, 27 implications of grain size trends, traffic surfaces, and construction methods are conceptualized and discussed. 2.3.5 Pile Scale To identify the influence of increasing pile size and height, the “structure” models featured in Figure 2.3 were expanded upon to the “multi-bench” models in Figure 2.4. These simulations include the same structural features as described in the push- and end-dumping structure models, but with three stacked benches separated by traffic surfaces as described in Section 2.2.4. These models were used to investigate whether the implications of structure and traffic surfaces from single bench model results hold up under multiple benches and over a larger pile scale. In this conceptual investigation, scale is simplified compared to traditional construction method features previously discussed. For example, all modeled WRPs share average parameter values between single-bench and multi-bench models to compare to the homogeneous model. This allows for isolated observations to be made with respect to pile size and dimensions, as well as the structural implications which may be missed under traditional laboratory-based testing in site-specific studies. 2.3.6 Assumptions and Limitations As mentioned, infiltration was considered constant and flow was assumed to be at steady-state for all simulations. In reality, the wet-up phase would begin during the initial pile construction period which may take several years to complete, and flow may fluctuate due to seasonal and climatic variations that are site-dependent (Smith et al. 1995). For the purposes of this study, the assumption of steady state flow was justified to assess the effects of physical heterogeneity and structure on drainage evolution, considering the time frame for simulations of several hundred years. To limit the complexity of the simulations, residual saturation (𝑆𝑆𝑤𝑤𝑤𝑤 = 0.22) and van Genuchten parameters (𝛼𝛼 = 5.0 m-1; 𝑛𝑛 = 2.7) were kept constant in all simulations, determined based on values observed at other 28 mine sites (Lefebvre et al. 2001; Fala et al. 2003; Azam et al. 2007; Neuner et al. 2013; Peterson 2014). Although these parameters deviate with grain size, saturated hydraulic conductivity values provide a dominant control on the calculated effective saturations and residence times within the model; therefore, homogeneous residual saturation and van Genuchten parameters were considered acceptable for the purposes of this investigation. Two cases that considered heterogeneity for these soil hydraulic function parameters were tested to validate this assumption for one random and fining upwards model and are included in Appendix A. The effect of varying these parameters was limited. In addition, isothermal conditions were assumed, and modeled piles were considered well aerated with gas concentrations held uniform through time, ensuring that pyrite oxidation was not limited by oxygen content. The simulation of highly reactive piles, leading to oxygen depletion, internal heating, and/or the accumulation of CO2 was beyond the scope of this study. It is possible that zones of low hydraulic conductivity could limit oxygen availability in these zones due to variations in saturation (e.g. Parbhakar-Fox et al. 2013), which could further contribute to spatial variation in pyrite oxidation rates under different model configurations. The mineralogy chosen for simulations was simplified and considers pyrite, calcite, gypsum and ferrihydrite, excluding processes such as weathering of aluminosilicate minerals and additional secondary mineral precipitation which would act to additionally buffer low pH values observed from pyrite oxidation. However, the purpose of this study was to investigate the ways in which physical and chemical heterogeneities (controlled by structure and spatial distributions in grain size) and scale influence ARD and overall drainage quality. While the simplifications made here may be important contributing factors for site-specific studies, they are considered sufficient for the purposes of this study to isolate the effects of heterogeneity, structure, and scale compared to the homogeneous model. 29 Figure 2.4. Hydraulic conductivity distribution for multiple bench models in cross-section. These cases expand upon the structure single bench cases from Figure 2.3 to include three stacked benches, with a homogeneous multi-bench model for comparison. 30 2.4 Simulation Results 2.4.1 Sensitivity to Parameter Heterogeneity The influence of spatially distributed parameters controlled by 𝑑𝑑𝑚𝑚 on simulated mass loading rates in basal drainage over time is depicted in Figure 2.5 and summarized in Table 2.3. Results between the five different random and fining upwards distributions did not vary significantly, suggesting that exact spatial location of heterogeneities is not significant for total drainage quality. In the following sections, references to “maximum mass loading rates” are frequently made, with respect to sulfate, which is used as a marker in this study for ARD-related release. These mass loadings assume complete mixing of basal effluent waters and are normalized per kilogram of waste rock represented in the pile. The addition of each parameter and its influence on total mass loading rates are explained using results from stylistic models below. In Table 2.3, a column for ‘Years Until Reduced Mass Release’ is presented, defined at mass loading rates <1 mg kg-1 week-1 (considerably smaller than maximum mass loading peaks, which ranged from 35 to 62 mg kg-1 week-1). While this does not suggest the end of ARD-related processes, this selected value was used to compare tail duration between homogeneous and heterogeneous model results. Additional information regarding the presence of calcite and depletion through time within simulated piles is found in Appendix A. 31 Figure 2.5. Comparison of the effects of heterogeneity in basal sulfate mass loading rates, influenced by physical and chemical heterogeneities. Plots on the left show results from heterogeneous models with random distributions; models that include fining upwards trends are shown on the right. From the legend: K is hydraulic conductivity, 𝛗𝛗𝒑𝒑𝒑𝒑𝒑𝒑 is the pyrite volume fraction (controlling differences in NP:AP ratios), and keff is the effective rate of pyrite oxidation, “c” refers to constant (or homogeneous spatial distribution for that parameter). Curves with the same colours represent simulations with similar trends and heterogeneities considered, but a different randomly generated spatial distribution. Specifically, Red: reference/homogeneous case. Yellow: K is spatially varied, 𝛗𝛗𝒑𝒑𝒑𝒑𝒑𝒑 and keff remain constant. Green: K and keff spatially varied, 𝛗𝛗𝒑𝒑𝒑𝒑𝒑𝒑 is constant. Blue: K and 𝛗𝛗𝒑𝒑𝒑𝒑𝒑𝒑 are spatially varied, keff is constant. Purple: K, 𝛗𝛗𝒑𝒑𝒑𝒑𝒑𝒑, and keff are spatially varied. In all cases where K is not constant, 𝝓𝝓 is also spatially varied. 32 Table 2.3. Variation in maximum mass loading rates and reduced sulfate release times with parameter heterogeneity Hydraulic Conductivity Volume Fraction Pyrite Rate Max Mass Loading Rate [mg1 kg-1 week-1] Decrease from Homogeneous Case Standard Deviation Years Until Reduced Mass Release a 1 2 3 4 5 Mean Homogeneous Constant Constant 61.8 --- --- 53 Random Constant Constant 52.2 52.0 51.9 52.0 51.4 51.9 16% 0.31 55 𝑓𝑓�𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒� 39.4 39.0 39.3 38.6 38.8 39.0 37% 0.32 139 𝑓𝑓 �1𝑑𝑑𝑚𝑚� Constant 37.6 37.7 38.0 36.8 36.9 37.4 39% 0.52 141 𝑓𝑓(𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒 ∙ 𝜑𝜑𝑖𝑖) 36.1 39.2 40.1 40.4 40.0 39.1 37% 1.77 138 Fining Upwards Constant Constant 54.2 54.0 54.5 54.3 54.0 54.2 12% 0.18 54 𝑓𝑓�𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒� 49.1 48.2 49.4 49.2 48.4 48.9 21% 0.53 118 𝑓𝑓 �1𝑑𝑑𝑚𝑚� Constant 40.0 41.1 41.8 41.4 41.5 41.1 33% 0.69 126 𝑓𝑓(𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒 ∙ 𝜑𝜑𝑖𝑖) 51.1 50.0 51.9 50.7 47.8 50.3 19% 1.56 121 a Average time duration [years] when average mass loading rate for SO42- returns to <1 mg kg-1 week-1 for mass outflux in modeled outputs. For multiple configurations, the average number of years between all 5 cases was used. 33 2.4.1.1 Physical Heterogeneity For cases of homogeneous chemical parameters, Figure 2.5 suggests that heterogeneity in hydraulic conductivity and porosity may influence residence times, varying the release of pyrite oxidation products. This results in decreased maximum mass loading rates and solute release rates compared to the homogeneous reference case, observed in Table 2.3. While the transport of pyrite oxidation products may cause minor changes in ARD evolution, given the assumptions listed in Section 2.3.6 simulated results suggest that variations in hydraulic conductivity and porosity alone do not significantly influence the duration of sulfate release; in Figure 2.5 the sulfate release tails (yellow) only extends marginally longer than the homogeneous reference case (red). 2.4.1.2 Chemical Heterogeneity In addition to the physical heterogeneities discussed in the previous section, the influence of chemical heterogeneities was assessed in Figure 2.5 by isolating tested heterogeneities in three scenarios: (1) NP:AP ratio heterogeneity, (2) pyrite oxidation rate heterogeneity, and (3) a combination of the two. Results from all scenarios where chemical heterogeneities are considered produce fairly similar results amongst the random and fining upwards cases, with the exception of scenario (1) (particularly for the fining upwards models). Cases (2) and (3) produce the most comparable results with respect to maximum mass loading rates and duration of release, which suggests that the spatial distribution of sulfide oxidation rates has a greater impact than local variations of mineral abundance in modeled simulations. 2.4.2 Pile Structure 2.4.2.1 Random vs. Fining Upwards Trends Though parameter distributions followed the same statistical trends for both the random and fining upwards models, differences in simulation results suggest that the evolution of the overall 34 basal effluent quality is subject to a both heterogeneity and structure. In all cases, the maximum mass loading rates for sulfate were greater for fining upwards trends compared to a random distribution. With respect to physical heterogeneity, a higher abundance of fines in the uppermost section retains higher moisture contents and restricts flow through lower hydraulic conductivity, causing increased residence times at the top of the pile. When chemical heterogeneities are considered, lower NP:AP ratios and higher pyrite oxidation rates are abundant in the fines due to the inverse relationship to 𝑑𝑑𝑚𝑚 used in the numerical models. This causes a build-up of poor quality water near the pile surface, which then must travel through the entire pile before exiting at a basal drainage point. As a result, a time-lagged but rapid-growing peak was observed for mass loading rates in fining upwards heterogeneous models. In the random heterogeneous models, the spatial variations in hydraulic conductivity, NP:AP ratios, and pyrite oxidation rates are more evenly distributed throughout the pile; as a result, pyrite oxidation products are observed at an earlier time and continue to be released for a longer duration, but have lower maximum mass loading rates compared to the fining upwards models. In cases where all investigated heterogeneities are considered (including hydraulic conductivity, NP:AP ratio, and sulfide oxidation), decreases of 13% and 32% in maximum mass loading rates were observed for the fining upwards and random models, when compared to the homogeneous model respectively. With respect to the homogeneous model, the heterogeneous models led to an extended duration until reduced drainage release was reached at 57 and 71 years, respectively, for the fining upwards and random models. While the heterogeneous models isolate for the influences of physical and chemical heterogeneities for fining upwards vs. random spatial distributions, WRPs often include a combination of these distributions combined with compacted traffic surfaces and multiple benches. These contributions from traffic surfaces and complex pile structure are addressed in the next section. 35 2.4.2.2 Traffic Surfaces and Increased Structural Complexity In order to investigate the influence of traffic surfaces and increased structural complexity on simulation results, Figure 2.6 shows the total sulfate mass loading rates, which were compared for the structure models for push- and end-dumping construction methods from Figure 2.3. The results can further be conceptualized through pH distributions in cross-section shown in Figure 2.7. While the differences in mass loading rates between the heterogeneous models (i.e. random versus fining upwards trends) are significant and suggest that spatial distributions of grain size may play an important role on drainage quality, results from the structure models (which include traffic surfaces and increased structural complexity) become very similar for push- and end-dumping construction methods. Furthermore, the addition of heterogeneities and structure resulted in a reduction in peak mass loading rates by one third from the homogeneous reference case for both scenarios considered. More detailed differences in behaviour for each construction method are described below. 4.2.2.1 Push-Dumping Models The single bench “structured” case push-dumping model (Figure 2.3) includes four traffic surfaces with decreased grain size, porosity, hydraulic conductivity, and increased pyrite oxidation rate overlying randomly distributed coarser material, with fining upwards trends along the slopes. As discussed previously, the same average parameter values were maintained in all simulations to ensure a fair comparison to the homogeneous reference case. However, due to structure, maximum sulfate mass loading rates are reduced, and release occurs over a longer duration. Differences in mass loading rates between random and structured heterogeneous models for push-dumping are less pronounced compared to the end-dumping models. The random spatial distribution in grain size in the heterogeneous models results in spatial variations in ARD-related processes laterally 36 and with depth, averaging out the release of oxidation products from sulfide oxidation. The addition of multiple compacted traffic surfaces and structure in the model acts to further average out these trends, resulting in a 36% decrease in peak mass loading rates compared to the homogeneous model. In addition, fining upwards along the pile slopes resulted in a slight time delay in mass loading, although less significant than the fining upwards and end-dumping construction models. 4.2.2.2 End-Dumping Models The single bench “structured” case end-dumping model (Figure 2.3) maintains the fining upwards grain size trends from the heterogeneous model but assumes zones of different average parameter values following the angle of repose, overlain by a single traffic surface. While the addition of a traffic surface influences the parameter values similarly through the mechanisms described for the push-dumping model, it has less influence on the results since there is only one at the pile surface—which already hosts a larger portion of fines due to the effects of gravity sorting. Instead, the presence of different zones (representative of different materials during pile construction) contributes most meaningfully to sulfate mass loading rates, as it acts to even out parameter distributions laterally through the pile rather than only by depth (as in the heterogeneous model). As a result, the maximum sulfate mass loading rate is greatly reduced compared to the homogeneous model (38% in the model including structure compared to 12% in the random heterogeneous model), with an increased duration of poor-quality drainage release. Although simulated sulfate release is evened out by the presence of material zoning, the time-lagged peak characteristic of fining upwards trends is still observed in the end-dumping model including structure. 37 Figure 2.6. Sulfate basal drainage evolutions shown for varied pile structure and scale. Results for total sulfate mass loading rates are normalized per kilogram of waste rock are shown for different scenarios illustrated in Figure 2.4. As the pile becomes more complex through the addition of traffic surfaces and internal structural trends, differences between push- and end-dumping construction are dampened. 38 Figure 2.7 pH cross sections at t = 3 years comparing heterogeneous versus structure models. Cross sections showing spatial distributions of pH throughout single bench models from heterogeneous models (compared to the homogeneous reference case) to conceptualize the influence of parameter heterogeneities. 39 Table 2.4. Variation in mass loading rates and duration of release with traffic surfaces and structure Case Max Mass Loading Rate [mg1 kg-1 week-1] Decrease from Homogeneous Case Difference Between Push vs. End Dump Methods Years Until Reduced Mass Releasea Single Bench Models Homogeneous 61.8 - - 53 Heterogeneous Push Dump 43.6 29% 25% 125 End Dump 54.6 12% 114 Structure Push Dump 39.8 36% 4% 127 End Dump 38.4 38% 139 Multiple Bench Models Homogeneous 36.2 - - 60 Heterogeneous Push Dump 20.5 43% <0.1% 139 End Dump 20.6 43% 132 a Average time duration [years] when average mass loading rate for SO42- returns to <1 mg kg-1 week-1 for mass outflux in modeled outputs. 2.4.3 Pile Scale To account for the change in size and shape, a multi-bench homogeneous case was used as reference for the multi-bench heterogeneous models. Comparing the two homogeneous models in Figure 2.6, the homogeneous multi-bench scenario showed reduced mass loading rates compared to the single bench model. Since reported mass loading rates are normalized per kilogram of waste rock in models, the total mass of sulfate released would be higher in the multi-bench scenario due to a larger volume of material. However, it is important to note that this 40% decrease in normalized mass loadings rate is significant, influenced only by size and shape when comparing the homogeneous models. Specifically, the addition of pile height over multiple benches increases duration of ARD-related product release and flushing of materials. Since the multi-bench models were adapted from the single-bench models including structure over a three-bench scenario, the processes influencing sulfate mass loading and ARD-generating 40 processes described in the previous section (i.e. influence of traffic surfaces and internal structures) also apply to the multi-bench models. However, as shown in Figure 2.6 and Table 2.3, not only is there a reduction in maximum mass loading rates for sulfate in the multi-bench models due to this influence of size and shape, but the peak release rates become essentially equal for both push- and end-dumping models. Generally, the effects of gravity sorting on a time-lagged release are observed in both multi-bench models, although much more pronounced in the case of the end-dumping scenario, with multiple delays observed in the sulfate mass loading curve. Most importantly, results suggest that the addition of heterogeneity, pile structure, and multiple benches significantly influence predicted effluent quality. Compared to the single bench homogeneous model, stylistic multi-bench models reduced sulfate maximum mass outflux rates by a factor of 3—an amount that would likely have significant implications for long-term treatment and closure planning. 2.5 Discussion 2.5.1 Parameter Sensitivity Modeled results from the previous section suggest that physical and chemical heterogeneities in waste rock piles likely affect drainage water chemistry significantly. If numerical models are used to make comparisons of different WRP storage strategies, incorporating these physical and chemical heterogeneities (particularly, mineral weathering rates and their relation to grain size distribution) causes substantial variation in peak mass loading rates and duration of contaminant release. While parameter heterogeneity is an important consideration, findings suggest that its influence on results is ultimately dependent on the spatial distribution, internal structures, and pile scale. 41 2.5.2 Pile Structure and Scale Results from overall basal drainage discussed in Section 2.4.2 suggest that differences in structure due to construction methods in a well-aerated pile matters to a certain scale, but does not strongly influence the outflow chemistry over the entirety of the pile once more complex internal structures and traffic surfaces are considered. In the random and fining upwards heterogeneous models, these spatial trends which are found in the cores of WRPs constructed using push- and end-dumping construction methods respectively produced quite different results in sulfate mass loading rates. However, the mass loading rate curves for different construction methods become much more comparable for the single bench structure model, and almost identical for the multi-bench models. This is an important finding; while spatial trends may influence drainage evolution through the mechanisms described above and may be dependent on the construction method (as observed in the heterogeneous models), the addition of traffic surfaces and complex internal structures which are common in WRPs (i.e. structure and multi-bench models) suggest that the choice of construction method may in fact not have a significant effect on drainage water chemistry and evolution. Simulated results are most significantly affected by the inclusion of structure, related parameter heterogeneities, and scale. In the multi-bench model, peak sulfate mass loading rates decrease by almost a factor of 2 from the multi-bench homogeneous case, and a factor of 3 from the single-bench homogeneous case. This is a significant decline in overall mass loading rates and has important implications for long term mine planning and environmental risk. On one hand, the homogeneous model arguably provides a higher factor of safety, since it predicts higher maximum mass loading rates. However, this scenario will fail to adequately predict long term mass loading rates and duration and may result in higher initial capital costs for building water treatment facilities than needed (on top of additional factors of safety implemented through regulatory 42 requirements). For the mass loading rates discussed, it is implied that basal drainage is entirely contained and perfectly mixed, which may not be the case in real world waste rock piles. The implications of the investigated conceptual models’ heterogeneity regarding monitoring and sampling techniques are discussed below. 2.5.3 Heterogeneity and Sampling Procedures In Section 2.4.1, it was noted that that there were no significant variations in overall mass sulfate loading between the five different spatial distributions in models with random heterogeneity and models considering fining upwards trends. Material characterization such as mean and standard deviation of parameters will influence the numerical outputs, but it would be nearly impossible (and certainly not cost effective) to determine the exact spatial distributions within a WRP given current measurement and sampling techniques. However, examining the pH distributions throughout the piles in Figure 2.7 and additional effluent water quality curves presented in Appendix A, the influence of heterogeneity on local drainage quality may have further implications for field sampling procedures. For example, seepage samples may be collected to monitor the drainage quality from a WRP following freshet or large rainfall events, but it is important to consider what these measurements represent. The consequences of this depend on the usage of these results: for example, is seepage chemistry used as a check for guidelines for water entering the receiving environment only, or being used as an input to validate models and make predictions? Is drainage assumed to be well mixed using a liner below the base of the pile? Furthermore, if a waste rock pile is not lined and discharge seeps directly into the ground below, it becomes impossible to know what the real evolution of the drainage quality is. In other sampling methods such as the collection of pore water samples using porous cups installed within the pile, Figure 2.7 and additional figures in Appendix A suggests that variation due to local acidic or 43 neutral pockets may cause a skewed perception of drainage quality that will come out of the pile. The results of these simulations strengthen the argument for solid, representative sample collection practices, and consider what implications exist if samples are not, especially in the case of highly heterogeneous material such as waste rock. 2.6 Conclusions The purpose of this study was to investigate the influence of spatially varied distributions of physical and chemical properties throughout WRPs and their internal structures on ARD-generating processes. Stylistic models simulated using MIN3P-HPC suggest that heterogeneity and pile structure are important to the overall basal drainage evolution observed from a pile. In the single bench heterogeneous models, it was suggested that push-dumping construction methods may be beneficial to reduce maximum solute concentrations entering the receiving environment over end-dumping methods—which caused a buildup of poor drainage quality at the top of the pile and resulted in a delayed spike in high concentrations of sulfate and iron concentrations from pyrite oxidation. However, this idea was challenged by the addition of structure (i.e. inclusion of zones, varied spatial trends, and traffic surface) and multiple benches which reduced solute release up to 3-fold for both construction types when compared to the single-bench, homogeneous model, proposing little difference in construction methods on maximum effluent concentrations observed. This study also suggests that oxygen limitation is not the only explanation for lower mass loading rates, as these results were achieved assuming well-aerated conditions for all simulated piles (i.e. O2 was not limited in pyrite oxidation). Due to these variations and their influence on simulated results, care must be taken when considering the effects of heterogeneity and structure on long-term closure planning. Furthermore, field data collected to constrain and develop conceptual 44 models should be representative and should allow conceptualization of localized heterogeneities and internal structure. While capturing the sensitivities of model outputs to parameter heterogeneity, this study also aimed to demonstrate the capabilities of the recently developed MIN3P-HPC code to simulate large-scale mine waste dumps. With model upgrades introduced in this version, it possible to simulate complex models with irregular geometry within a reasonable time, particularly if utilizing the parallelized version of the code (in this study, most models were simulated for 300 years, requiring only a few hours of CPU time). Further advancements to MIN3P and other reactive transport codes along with technological progressions will improve certain shortcomings and better capture heterogeneity, though will also come with challenges associated with the opportunity to add more information to RTMs. As demonstrated in this study, modeled results may be very sensitive to the inclusion of heterogeneity and structure. Further investigations for these future models and the influence of incorporating heterogeneity will need to continue to be updated. Finally, this work further supports the use for numerical modeling as a conceptualization tool in understanding potential WRP behaviour. The idea that WRPs are heterogeneous is not a new concept, but the most effective ways to capture the effects of heterogeneity on the overall pile evolution and ARD-related processes from field or lab-based experiments still have many uncertainties (as discussed at the beginning of this paper). Numerical models cannot provide a substitute for necessary field data, supplementary techniques, and lab-based experiments which have shown to support and constrain site-specific material behaviour and results. However, RTMs can be used as an effective tool to conceptualize different waste rock storage strategies, which can be explored at an initial stage before further investment and narrow down the most effective long-term storage methods and develop more robust closure strategies. However, simulated model 45 outputs can vary greatly depending on pile size, internal structures, and parameter heterogeneities considered as shown in this study, and are ultimately the result of the empirical formulas used in the code and the input parameters. It is the responsibility of the user to carefully consider and constrain models with added complexities such as heterogeneity and structure with solid, representative field measurements. 46 Chapter 3: Numerical Modeling of a Laboratory-Scale Waste Rock Pile Featuring an Engineered Cover System Abstract Improved design of stockpiles and storage techniques to reduce contaminant mass loadings from waste rock piles (WRP) is an increasingly important consideration for mining operations globally. In certain cases, the implementation of an engineered cover system may be a suitable technique to mitigate the release of metals through acid rock drainage (ARD) or contaminated neutral drainage (CND) out of WRPs. This method often relies on using differences in soil hydraulic parameters to create a capillary barrier effect (CBE) to divert water away from reactive materials below. In this study, a numerical model using the reactive transport code MIN3P-HPC was calibrated to observational data from a laboratory-scale model constructed by Poaty (2020) to investigate the effectiveness of a proposed engineered cover system for a full-scale WRP design plan. Results showed the numerical model is capable of capturing flow rates out of multiple drainage ports and matching key effluent concentrations and loadings by including heterogeneities in soil hydraulic parameters and mineral weathering rates. These results would not have been achievable using a homogeneous model. In addition, simulations were useful for characterizing internal flow pathways within the laboratory model, showing the effectiveness of the cover and also identifying a secondary CBE between different rock types within the laboratory model. The simulated results were useful in quantifying reaction rates and understanding the processes controlling metal release and attenuation. The numerical model proved beneficial in conceptualizing the system and could be applied to other WRP designs and testing of cover systems in the future. 47 3.1 Introduction The long-term storage of waste rock—material considered not economically viable, which must be removed from the subsurface to access valuable resources during the mining process—is an increasingly important consideration for mine operations globally. When exposed to atmospheric conditions, available sulfide minerals within waste rock may oxidize and cause the release of acid rock drainage (ARD), or if sufficient neutralizing minerals are present, contaminated neutral drainage (CND) (Sherlock et al. 1995; Lawrence and Wang 1996; Blowes et al. 2003; MEND 2004b; Price 2009; Vriens et al. 2018). To mitigate the release of metals and contaminants from ARD and CND, engineered design structures such as cover systems, subaqueous disposal, blending and layering, or co-disposal techniques may be applied as a means of prevention (e.g. Price and Errington 1998; MEND 2001, 2004a; Haug 2002; Johnson and Hallberg 2005; INAP 2009). In addition to the laboratory- and field-scale tests required to assess these techniques and characterize materials, reactive transport modeling (RTM) may be used as a tool to investigate processes involved in the release and attenuation of metals from waste rock piles through scenario analyses and/or to quantitatively analyze and interpret data from laboratory experiments and field observations. For example, RTM has been used in the context of mine wastes to compare outcomes under multiple scenarios, investigate complex features such as heterogeneity, contrast different storage techniques, and bracket expected behaviour for forecasting in a number of previous studies (e.g. Linklater et al. 2005; Molson et al. 2005; Fala et al. 2013; Mayer et al. 2015, 2019; Amos et al. 2015; Shokri et al. 2016; Lahmira et al. 2017; Blackmore et al. 2018b; Wilson et al. 2018b). In this study, the reactive transport code MIN3P-HPC (Su et al. submitted) was used to create a numerical model to supplement the ongoing research for an engineered cover system applied at 48 the Tio Mine site (Rio Tinto Iron Titanium). At this site, near-neutral pH drainage may at times release Ni concentrations above local regulatory limits (0.50 mg/L) (i.e. Quebec Province, Canada) (Plante et al. 2011a, 2014), requiring water treatment. To control Ni and release of other metals, a fluid control layer (FCL) composed of unreactive crushed anorthosite available at the site and sand was designed to divert water from the metal-containing waste material, which is dominated by below-grade hematite-ilmenite (HI) ore. The characterization of Lac Tio waste rock and the effectiveness of the FCL at diverting water while minimizing erosion has been well-studied in the past (Aubertin et al. 2009; Dubuc et al. 2013; Plante et al. 2014; Dimech et al. 2018, 2019; Martin et al. 2017; Poaty 2020). However, a quantitative two-dimensional analysis of flow and reactive transport processes has not been conducted to date. To this end, a numerical model was constructed and calibrated for the quantitative analysis of the laboratory experiment conducted by Poaty (2020), which was designed to evaluate the implementation of an FCL at the full-scale waste rock pile at the Tio site. Specifically, RTM was used to characterize the internal flow regime within the laboratory model constrained by observational data and measured material properties and to improve understanding of the influences of physical and chemical heterogeneities (such soil hydraulic parameters and mineral precipitation-dissolution kinetics) between different material zones controlling the effluent quality, which varied spatially throughout the experiment. 3.1.1 Site Background The Tio Mine is an open-pit mine featuring a massive HI deposit hosted in anorthosite gangue material. The mine site has been in operation since 1950 and is located 43 km north of Havre-Saint-Pierre, Quebec, Canada. Waste rock at this site varies between anorthosite and highly mineralized material with HI content >70% (Plante et al. 2014). Nickel, a primary focus of previous studies, comes primarily from Ni-bearing pyrite ((NixFe1-x)S2), with x varying spatially 49 between 0 – 10% throughout the HI deposit. Additionally, Ni is present in the form of millerite (NiS) embedded within pyrite grains (Plante et al. 2011a). Release of Ni from waste rock is retarded through sorption mechanisms in both plagioclase and ilmenite rock (Plante et al. 2010). The drainage from Tio waste rock is maintained at near neutral conditions due to the weathering of silicate minerals, specifically: calcic plagioclase which is mainly close to labradorite in composition (Na0.4Ca0.6Al1.6Si2.4O8), enstatite (MgSiO3), and biotite (K(Mg,Fe2+)3 (Al,Fe3+)Si3O10(OH,F)2) (Plante et al. 2011b). Silicate buffering to maintain near-neutral pH in the absence of carbonate minerals has been demonstrated previously by Li (2000) and Jambor et al. (2007). Silicate dissolution is typically incongruent, preferentially releasing cations such as Ca, Na, and Al from plagioclase minerals, Mg from pyroxenes (Oelkers and Schott 2001; Zakaznova-Herzog et al. 2008) and Mg, Fe, and K from sheet silicates (Stoch and Sikora 1976). Further detailed analyses of the weathering mechanisms and specific dissolution rates previously determined for Tio waste rock can be found in Plante et al. (2011b). In order to improve drainage quality, an option under investigation at the mine consists of designing a FCL using crushed anorthosite and sand. Objectives of the FCL are to divert water from reactive materials using contrast in grain size to create a capillary barrier effect (CBE) and to reduce the mass loadings of contaminants out of the pile to minimize water treatment requirements. Previous studies by Fala et al. (2003, 2005); Bréard-Lanoix (2017) found that a slope of 5˚ was most effective for the FCL to divert surface precipitation while minimizing erosion. To investigate this design under controlled conditions, a heterogeneous laboratory-scale model was constructed by Poaty (2020) utilizing the relevant waste rock lithologies and compositions. Weekly flush experiments were performed on this laboratory model inclined at various angles between 0 – 5˚ over a period of 24 weeks between October 2017 and April 2018. Information regarding 50 results of these transient tests can be found in Poaty (2020). In this study, the system covered with the FCL inclined at a 5˚ angle was subjected to a constant application rate with leachate collected at different sampling ports along the base of the tank. This experimental study aimed to emulate the flow dynamics for covered Tio waste rock under high infiltration rates. The objective of the experiment was to quantitatively investigate the performance of the FCL and to help quantifying hydrogeochemical reaction rates, which can potentially be used for future upscaling purposes. 3.2 Methods 3.2.1 Laboratory Methods 3.2.1.1 System Overview and Construction The laboratory-scale model was constructed by Poaty (2020) by depositing crushed waste rock material from the Tio mine site using a hand-held scoop in 4 sections: 3 zones of HI material (labeled S1 – S3) and 1 anorthosite zone at the end, each 80 cm in height, in a 60 cm wide tank (Figure 3.1). Overlying these layers, a 10 cm thick sand layer and 8 cm thick crushed anorthosite layer were placed to create the FCL. Anorthosite material was crushed to <10 mm for the FCL using a crusher at the Tio mine site and compacted using a Proctor hammer with constant mass to approach target porosities and to create a CBE, diverting water from the HI material below (Poaty, personal communication, 2020). Grain size distribution curves for the materials in each of the sections can be found in Appendix B . Along the base, 8 drainage ports (1 cm diameter, referred to as reservoirs 1 – 8) were installed equipped with tubing, allowing water to drain into collection buckets below (tubing lengths for each port are provided in the Appendix B ). Effluents collected from each reservoir were weighed daily to determine flow rates; subsequently samples were collected for geochemical analysis. With the tank exposed to atmospheric conditions along the top, initial conditions throughout the pile were assumed to be well-aerated with respect to O2 (21%) 51 and CO2 (0.0412%) (US Department of Commerce 2019). Additional details on the laboratory model design and construction can be found in Poaty (2020). 3.2.1.2 Flow Measurements The experiment was performed under steady-state flow rates on the laboratory model over a 17-day period in February 2019 with a sprinkler system applying a constant and uniform water application across the entire surface with the laboratory model angled at 5˚. Drainage was collected in buckets below each reservoir and weighed to determine daily flow rates from each port. In addition, frequency domain (FD) sensors placed throughout the pile were used to measure volumetric water content using dielectric constants (Dimech et al. 2019). While FD probes were not calibrated specifically for this 17-day experiment, variations in measured apparent dielectric constant values were used to support bucket weight measurements and to determine when steady-state flow conditions in the laboratory model were reached. While the flow rates obtained from the bucket weight measurements provided information on how much water is leaving at each reservoir, they cannot be used alone to characterize the internal flow pathways within the experiment. Figure 3.1.Schematic of laboratory-scale waste rock pile adapted from Poaty (2020). 52 3.2.1.3 Geochemical Monitoring Leachates were collected in buckets below reservoirs 1 – 8, in which pH and electrical conductivity were monitored daily. Once these parameters had stabilized, water samples were collected from the buckets and filtered with 0.45 μm mesh filters into to pre-acidified (with 2% HNO3) 125 mL PVC bottles for full-suite geochemical analyses using ICP-MS. Samples were collected and analyzed over a period of 10 days after steady-state flow conditions were achieved, which was considered sufficient to represent stabilized effluent quality out of the laboratory model within the scope of this study. In between the transient flush experiments performed by Poaty (2020) and the steady-state experiment completed for this study, the laboratory model was left inactive. During this period of inactivity, materials were left to freely oxidize; therefore, samples were not collected between days 1 through 7 due to budgetary constraints and since these samples were not considered representative for quasi-steady state drainage conditions. Early time drainage was expected to represent release of weathering products that had built-up during the several months of inactivity between transient and steady-state flow tests. For additional details on the breakthrough curves or transient behaviour of the laboratory model, see Poaty (2020). 3.2.2 Numerical Methods 3.2.2.1 Reactive Transport Code MIN3P-HPC MIN3P-HPC (Su et al. submitted) is a multicomponent reactive transport code which builds on MIN3P-THCm (Bea et al. 2016; Su et al. 2017), coupling multicomponent transport and geochemical reactions in variably saturated porous media using a structured or unstructured grid (i.e. a polygonal (2D) or polyhedral (3D) mesh). The code uses Richards’ equation to solve flow in the unsaturated/saturated subsurface, while a global implicit solution approach is used to solve 53 the reactive transport equations (Mayer et al. 2002). Additional details about the code may be found in Mayer et al. (2002), Xie et al. (2015), and Su et al. (2017). In this study, a 2D numerical model solving for steady-state flow and transient geochemical conditions was designed to quantitatively interpret data collected from the laboratory model presented above. The conceptual model of the experiment including the boundary conditions is depicted in Figure 3.2, while the mesh configuration used in provided in Appendix B . The versatility of the unstructured grid used in this study allowed for direct inclusion of the angled surfaces (such as the 5˚ rotation of the FCL and sloped edges following the angle of repose in the HI zones (Figure 3.2), which could not have been captured adequately using a standard rectangular-based mesh. 3.2.2.2 Flow Model and Calibration Previous works by Poaty (2020) and Dimech et al. (2019) utilized suction and moisture content measurements (collected using Watermark® Granular Matrix sensors and ECH2O EC-5 soil moisture sensors, respectively) from transient flow experiments to better understand the internal flow regimes. In this study, measured flow rates from each reservoir were used to calibrate the numerical model, which was constrained by measured material parameters and the system geometry. Figure 3.2 illustrates the material property zones and boundary conditions used in the numerical model. Although a water application of 180 L day-1 was intended from the sprinkler system, the total outflow collected from the reservoirs averaged only 158 L day-1 under steady-state conditions. This difference in mass in versus mass out is attributed to a combination of factors including evaporation and difficulties in maintaining pump pressure and water delivery (additional details are provided in Appendix B ). For consistency, a steady-state water application of 158 L day-1 was applied in the model to match the total measured outflow rate under steady-state 54 conditions. This water application was uniformly distributed across the surface of the experiment, corresponding to an infiltration rate of 108 mm day-1. Constant pressure heads were specified to calibrate steady-state drainage rates at each reservoir. Simulated drainage rates proved to be very sensitive to specified pressure heads at each of the sampling ports; pressure head changes on the order of a tenth of a mm influenced drainage rates by over 1 L day-1 in some instances. Boundaries along all remaining sides of the domain were assumed to be no flow (confined by the tank walls of the physical model). Since no flow was measured out of reservoirs 6 and 7 in the physical model, these locations were also treated as no flow boundaries in the numerical model (Figure 3.2). The variably saturated flow model was further constrained by measured physical parameters as summarized in Table 3.1. Porosity measurements were taken directly from the laboratory model, while saturated hydraulic conductivity and unsaturated hydraulic soil parameters were obtained from previous column experiments and characterization of the waste rock and FCL (Intissar 2009; Bréard-Lanoix 2017; Dimech et al. 2018; Poaty 2020). Residual volumetric moisture contents reported by Bréard-Lanoix (2017) and Peregoedova (2012) for sand and waste rock materials, respectively, were used to determine residual saturations for these materials. These parameters were used directly and were not varied during model calibration. 55 Table 3.1. Physical parameters for variably saturated flow separated by zone Zone Saturated Hydraulic Conductivity [ms-1] Porosity [-] Residual Saturation [-] 𝜶𝜶 [m-1] n [-] Air Entry Pressure [m] HI S1 5.5 × 10−4a 0.39c 0.13ac 5.96d 3.61d 0.04a HI S2 5.5 × 10−4a 0.55c 0.09ac 5.96d 3.61d 0.04a HI S3 5.5 × 10−4a 0.58c 0.09ac 5.96d 3.61d 0.04a Anorthosite 5.0 × 10−4b 0.35c 0.14ac 3.60a 1.56a 0.04ab Sand 7.8 × 10−4a 0.37c 0.08ce 1.64a 4.99e 0.30ce Crushed Anorthosite 5.0 × 10−5b 0.29c 0.10ac 20.0b 1.56d 0.04b a Peregoedova 2012; b Martin et al. 2017; c Poaty 2020; d Demers 2013 e Bréard-Lanoix 2017; f Intissar 2009 Figure 3.2. Conceptual model with labeled boundary conditions and material property zones. Along the top surface (blue), an infiltration rate of 108 mm day-1 is applied. Red dashed lines surrounding the other boundaries represent no flow boundaries. Yellow symbols with arrows represent pressure head boundaries where effluent was collected from each reservoir. Pressure head (Ψ) boundaries were calibrated to achieve drainage rates observed from the laboratory experiment under steady-state flow conditions. Reservoirs 6 and 7 were treated as no flow boundaries in the numerical model, since no outflow was observed at either port, however their locations are labeled here. Corresponding hydraulic heads from left to right: -2.2 cm, -2.3 cm, -2.2 cm, -2.5 cm, -2.1 cm, -2.0 cm. 56 3.2.2.3 Reactive Transport Model and Calibration for Leachate Concentrations The reactive transport model considered a free-phase diffusion coefficient of D0 = 2.4 × 10-9 m2 s-1 for all dissolved species, while longitudinal and transverse vertical dispersivities of 0.05 m were used. In the gas phase, a free phase diffusion coefficient of 2.1 x 10-5 m2 s-1 was used for all gases. Effective diffusion coefficients were calculated based on the Millington-Quirk relationship (Mayer et al. 2002). Only the most relevant chemical components were considered in the model including SO42-, Ni+, Fe2+, Fe3+, Ca2+, K+, Mg2+, Na+, Al3+, H4SiO4, CO32-, H+ and O2(aq). Aqueous complexation of these components was also included yielding a total of 26 dissolved species (see Appendix B for included hydrolysis and complexation reactions). Measured composition of the tap water applied to the surface of the physical experiment was specified along the infiltration boundary in the numerical model as specified flux (Cauchy) boundary (Table 3.2). Table 3.2. Composition of infiltrating water applied to the model Species Concentration [mg L-1] SO42- 65 Ni+ 0.007 Fe2+ -- Fe3+ 0.105 Ca2+ 24 K+ 0.859 Mg2+ 3.41 Na+ 5.95 Al3+ 0.02 H4SiO40 7.84 pH 6.87 Alkalinity [as CaCO3] 22 O2(aq) 0.21 (as pO2 in atm) 57 Boundary conditions in the reservoirs with active drainage were specified as free exit boundary conditions, while all other boundaries were considered impermeable, consistent with the presence of the tank walls. While the mineralogy of the various materials is complex, the mineralogy considered for the numerical model has been simplified and only minerals which are considered to effectively control the system were included. Table 3.3 and Table 3.4 summarize the primary and secondary minerals, their initial volume fractions, reaction pathways, and relevant equilibrium constants considered in the numerical model. The majority of these volume fractions were taken from previous characterization of waste rock by (Demers et al. 2013), though volumes of pyrite was calculated using results of C-S analyses completed by Poaty (2020). Table 3.3. Initial mineral volume fractions considered in numerical model Mineral Volume Fraction HI S1 HI S2 HI S3 Anorthosite Sand Crushed Anorthosite Ni-bearing Pyrite 0.1649a 0.1524a 0.1060a 0.0004a --- 0.0096a Labradorite 0.2054b 0.2054b 0.2054b 0.2054b --- 0.2054b Enstatite 0.1158b 0.1158b 0.1158b 0.1158b --- 0.1158b Biotite 0.0086b 0.0086b 0.0086b --- --- --- Ferrihydrite --- --- --- --- --- --- Al(OH)3 --- --- --- --- --- --- SiO2 --- --- --- --- --- --- Kaolinite --- --- --- --- --- --- a Calculated from measured S wt. % (Poaty 2020) b Demers et al. 2013 58 Table 3.4. Minerals and their dissolution-precipitation reactions considered in numerical model. Effective rate constants (𝒌𝒌𝒆𝒆𝒆𝒆𝒆𝒆,𝒎𝒎) were dependent on material zone and are presented in the results section in Table 3.6. Mineral Reaction Rate Expression 𝐥𝐥𝐥𝐥𝐥𝐥𝐊𝐊𝒎𝒎 Ni-bearing Pyrite Ni0.14Fe0.86S2 + H2O + 3.5O2→ 0.86Fe2+ + 0.14Ni+ + 2SO42− + 2H+ −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚[O2] - Labradorite Na0.4Ca0.6Al1.6Si2.4O8 + 1.6H2O + 6.4H+→ 0.4Na+ + 0.6Ca2+ + 1.6Al3+ + 0.6H4SiO4 −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚[H+]0.55a - Enstatite MgSiO3 + 2H+ + H2O → Mg2+ + H4SiO4 −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚[H+]0.25b - Biotite K(Mg2Fe)AlSi3O10(OH)2 + 10H+→ 2Mg2+ + K+ + 2Fe2+ + 3H4SiO4 −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚[H+]0.25c - Ferrihydrite Fe(OH)3 + 3H+ → Fe3+ + 3H2O −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚 �1−IAPK � -4.891d Al(OH)3 (am) Al3+ + 3H2O → Al(OH)3 + 3H+ −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚 �1−IAPK � -10.80d SiO2 (am) H4SiO4 → SiO2 + 2H2O −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚 �1−IAPK � 3.018d Kaolinite 2Al3+ + 2H4SiO4 + H2O → Al2Si2O5(OH)4 + 6H+ −𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑚𝑚 �1−IAPK � -7.435d a Blum and Stillings 1995 b Brantley and Chen 1995 c Nagy 1995 d MINTEQA2 - Allison et al., 1991 Nickel in the HI material is sourced from both pyrite and millerite grains (found in inclusions within the pyrite). However, simulations considered pyrite as the only Ni source, with dissolution rates calibrated using measured SO4 release rates. Although previous work to characterize Tio waste rock identified adsorption as an important Ni-attenuation mechanism (Plante et al. 2010; Plante et al. 2011b; Demers et al. 2013, Plante et al, 2014), this process was not explicitly included in the current numerical model. In this study, the mineral weathering rates in the numerical model were calibrated to match release rates observed in the drainage of the physical model and constitute effective release rates, implicitly accounting for the effect of adsorption. Accounting for sorption explicitly was beyond the scope of the present study. Previous work by Plante et al. (2011b) showed that both S and Ca release rates were stoichiometrically associated with sulfide (primarily in the form of pyrite, although trace amounts of chalcopyrite and millerite have been detected) and plagioclase weathering, respectively, due to 59 the lack of secondary mineral precipitation. To match the effluent chemistry of the system, pyrite and labradorite weathering rates were calibrated to match SO4 and Ca release rates. The presence of Ni impurities in pyrite has been noted to increase pyrite oxidation rates in previous studies (Nicholson et al. 1988; Lehner and Savage 2008). Enstatite and biotite were added to the model to represent additional pH-buffering reactions by silicate minerals (using Mg and K release rates to calibrate, respectively). The secondary minerals ferrihydrite, SiO2 (am), Al(OH)3 (am) and kaolinite have been identified by (Plante et al. 2011b) and were included to control Fe, Al, and Si concentrations in the drainage waters. 3.2.3 Assumptions and Limitations Building a conceptual model for RTM requires making assumptions and simplifications depending on what data is available and the overall purpose of the model. For example, Table 3.1 shows that porosity was determined for the various materials emplaced in the experimental vessel, while other input parameters were taken from other Tio column and field experiments. In reality, values for each material property zone may vary slightly compared to parameters obtained from the companion experiments. Parameters may also vary spatially and show heterogeneities within each material property zone. In addition, the geochemical reaction network represents a simplification of the actual reactions occurring in the physical experiment and mineral reaction rates which were assumed to be constant in time within each material property zone. This assumption is considered justified given the timescale of this quasi steady-state experiment, and is also supported through previous work by Plante et al. (2011a) who determined through humidity cell tests that silicate weathering rates in fresh Tio waste rock are very similar to weathering rates in material that has undergone 25 years of natural weathering. 60 In addition, the numerical model was built in 2D, while the laboratory experiment contains features that deviate from a true 2D-system. Specifically, the drainage ports in the laboratory experiment were located in the center of the base of the experiment in y-direction. In the 2D numerical model, each reservoir had to be taken as a line of drainage rather than a 1 cm diameter central port. Despite these simplifications, it must be emphasized that most model input data was available through measurements from the well constrained experiment and were kept constant during the model calibration process. The only parameters varied during the modeling were the pressure heads at the drainage ports and mineral weathering rates with the goal to reproduce drainage rates and concentrations from the various reservoirs under quasi-steady state conditions. 3.3 Results 3.3.1 Drainage Rates and Internal Fluid Flow Pathways FD probe sensor measurements showed stable volumetric water contents within the laboratory model were reached after 2 days, while measured flow rates reached quasi-steady-state after approximately 8 days. Simulated results shown in Figure 3.3 and Table 3.5 were able to capture drainage rates with less than 1% difference at all points. The drainage rate out of reservoir 3 (approximately 70 L day-1) constitutes 44% of the total water volume leaving the model. Based on examining drainage rates alone, it is unclear whether a higher volume of water leaving reservoir 3 is caused by water leaking through the FCL into the HI material, or whether water is being diverted by the FCL, as intended, and then rerouted to reservoir 3. 61 Figure 3.3. Measured vs. simulated cumulative drainage. Blue points represent measured cumulative drainage from the laboratory experiment, while black lines represent simulated cumulative drainage, calibrated using pressure heads at each reservoir. Dashed lines represent the simulated cumulative drainage but shifted to account for the delay of the development of steady-state flow conditions in the physical experiment. 62 Table 3.5. Results for measured versus simulated drainage rates Reservoir 8 7 6 5 4 3 2 1 Total Flow Rate [L d-1] Measured 7.98 0.00 0.00 21.20 3.88 69.15 29.33 26.50 158.03 Simulated 7.98 0.00 0.00 21.13 3.86 69.38 29.05 26.54 157.94 Percent Difference -0.08% 0.00% 0.00% 0.33% 0.51% -0.33% 0.94% -0.16% 0.06% Percentage of Total Flow Measured 5.05% 0.00% 0.00% 13.42% 2.46% 43.75% 18.56% 16.77% 100.0% Simulated 5.05% 0.00% 0.00% 13.38% 2.45% 43.93% 18.39% 16.80% 100.0% In Figure 3.4, modeled results of Darcy flux vectors are presented by scaled arrows, which are used to conceptualize flow through the laboratory model. Simulations suggest that in the steady-state flow experiment, a water table forms at z ≈ -2 cm, approximately corresponding to the elevation of reservoir 8. The pressure head boundaries used to calibrate the model (Figure 3.2) were positive for reservoirs 1-5, facilitating outflow. Pressure head at reservoir 8 was slightly negative; however, outflow can still occur under these conditions since a drainage tube was attached to the outflow port at the base of the tank. In the HI zones in Figure 3.4, smaller arrows show that there is less water flowing downwards compared to the anorthosite zone. These simulated flow pathways suggest that while there is some water leaking through the FCL into the HI material below, the FCL angled at 5˚ effectively diverts a significant portion of the water under quasi-steady state flow conditions. Furthermore, model results suggest a secondary CBE exists between the HI material and the uncrushed anorthosite layer (Figure 3.4). Based on these conditions, a flat water table forms along the base of the model and causes some of the drainage water to move upwards (to the left) from the uncrushed anorthosite layer back into the HI material, 63 contributing to the larger flow rates observed at reservoirs 2 and 3. These results suggest that leakage through the FCL is limited and that the large outflows at reservoirs 3 are due to internal rerouting of flow paths due to CBEs. Figure 3.4. Simulated results of calibrated flow model. Black arrows represent Darcy flux vectors simulated in the numerical model indicating direction of fluid flow pathways and their relative magnitude, overlain on simulated water saturation throughout the cross section of the experiment. The left bottom corner was set to z = 0 cm; a simulated flat water table at z = -2.2 cm is shown in white. 64 3.3.2 Solute Concentrations in Drainage Results of measured and simulated water quality at each reservoir are shown in Figure 3.5. Observed solute concentrations in drainage suggest that quasi-steady state conditions were achieved after day 12. Therefore, average solute concentrations between day 12 and day 17 were used to calibrate mineral weathering rates within the numerical model (Table 3.6). Observations between day 8 and day 12 show elevated solute concentrations and decreased pH, capturing the tail end of the flushing of weathering products originating from oxidation reactions that took place while the experimental setup was idle. Early time behavior in the model also show variations in concentrations, due to the specified initial conditions. Simulated concentrations in drainage between day 0 and day 10 are affected by displacement and mixing of resident pore water with infiltrating water and weathering reactions. These early time results are not relevant for the model calibration and interpretation, performed for quasi-steady state conditions, and are therefore represented with dashed lines (Figure 3.5). Simulations included the diffusion of O2 through the gas phase and results suggest that O2 was not depleted within the model domain during the experiment duration, and subsequently the oxidation of Ni-bearing pyrite was not inhibited. The lack of O2 depletion, even for fully saturated regions near the base of the domain, can be explained by short residence times and fast water flow velocities. Pyrite oxidation rate coefficients were calibrated to reproduce SO4-release observed in the experiment. Pyrite grains have been observed to typically contain between 0 – 10 wt. % Ni at the Tio site (Plante et al. 2014), and a pyrite composition of Ni0.14Fe0.86S2 (corresponding to 7 wt% Ni) was found to represent the observed Ni release rates most accurately, consistent with the observed Ni contents in drainage. The current approach utilizing effective Ni-release rates, implicitly accounting for sorption, was able to adequately reproduce observed Ni-release, although 65 previous work had shown that sorption processes are significant in fresh rock at the Tio site and important for controlling the release of Ni (Plante et al. 2011a, 2014). It is possible that the effect of adsorption in this experiment was less pronounced than in previous studies, since multiple pore volumes had passed through the material during previous experiments (Poaty 2020), which may have contributed to saturating available adsorption sites. In addition, high volumetric flow rates were applied to the laboratory model for the steady-state flow experiment, which may have reduced typical sorption processes occurring in the Tio waste rock due to rapid flushing. Another possibility is that some sorption-processes may have been accounted for implicitly by the calibrated Ni-content in pyrite used in the numerical model. As previously mentioned, Ni may range from 1 – 10 wt. % content in pyrite in Tio waste rock (Plante et al. 2014), therefore it is possible that the HI material in the laboratory model is greater than the 7 wt. % used here. Additional laboratory testing may be able to deduce actual Ni-release rates and the contribution of adsorption towards effective Ni-release rates in this experiment; however, this is beyond the scope of the present study. As mentioned in the methods section, labradorite, enstatite, and biotite were the primary minerals considered for pH buffering. Volume fractions of minerals presented in Table 3.3 were not depleted during the 17 days simulated by the numerical model and were therefore not considered a limitation for the model reaching a quasi-steady state. Rate coefficients for these phases were calibrated to match the release of major ions for each mineral. Rate coefficients for pyrite and silicates are summarized for each zone in Table 3.6, consistent with values found in the literature. While simulated concentrations in drainage waters reproduce measured values for SO4, Ca, K, Mg, Na and pH fairly well for most reservoirs (Figure 3.5) —Si was above measured values. 66 Simulations suggest that pore water was undersaturated with respect to SiO2 (am) throughout the domain. Although kaolinite precipitation was considered in the simulations, its precipitation was not able to reproduce the observed Si concentrations. Elevated Si values are likely attributed to the incongruent weathering process for silicates mentioned previously, which were not captured in the simulations. Because aluminosilicate weathering occurs congruently in the numerical model, in order to match major ions (Ca, Mg, and K, which are not controlled by secondary mineral precipitation), while buffering pH and Si release are over-predicted. Previous studies (Plante et al. 2010, 2011b) provide evidence for the presence of additional silicate minerals in the HI material such as pigeonite, muscovite, chlorite, orthoclase, and spinel, which were not considered. These minerals may also contribute to the silicate weathering and may at least in part explain discrepancies between simulated and measured results. However, one of the main objectives of this study was reproduce the release of solutes from the different reservoirs with an emphasis on SO4, Ni and major cations, which was successfully accomplished (Figure 3.5). Simulated concentrations for Al and Fe were both below measured values of approximately 3×10-6 mol/L and 5×10-6 mol L-1, respectively. Speciation calculations for observed water chemistry indicates supersaturated conditions with respect to ferrihydrite and kaolinite, implying that precipitation of these secondary phases is kinetically inhibited, partially due to low concentrations of Al and Fe. 67 Figure 3.5. Results of simulated vs. observed effluent concentrations in mol L-1, color-code per reservoir. Effective rate coefficients for mineral dissolution-precipitation rates (𝒌𝒌𝒆𝒆𝒆𝒆𝒆𝒆,𝒎𝒎) in Table 3.6 were calibrated through an iterative process to reproduce observed concentrations in the effluents from the laboratory experiment. 68 Table 3.6. Calibrated effective rate coefficients for mineral weathering reactions Mineral 𝒌𝒌𝒆𝒆𝒆𝒆𝒆𝒆,𝒎𝒎 HI S1 HI S2 HI S3 Anorthosite Sand Crushed Anorthosite Ni-bearing Pyrite 9.6×10-7 3.0×10-7 3.5×10-7 1.0×10-7 --- 8.0×10-7 Labradorite 2.3×10-6 8.0×10-7 8.0×10-7 3.5×10-8 --- 9.0×10-8 Enstatite 7.0×10-11 2.0×10-11 1.0×10-11 1.0×10-11 --- 1.0×10-11 Biotite 3.0×10-11 3.0×10-11 3.0×10-11 3.0×10-11 --- 4.5×10-11 Ferrihydrite 7.0×10-9 7.0×10-9 7.0×10-9 7.0×10-9 7.0×10-9 7.0×10-9 Al(OH)3 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 SiO2 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 Kaolinite 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 1.0×10-8 3.4 Discussion 3.4.1 Implications of Results from Flow Modeling For this study, the daily infiltration rates considered were equivalent to one month of precipitation at the field site. In reality, saturation within the full-scale pile would be lower and undergo transient infiltration processes, leading to differences in flow dynamics. Nevertheless, the most important results from the numerical model suggest that preferred fluid flow pathways exist controlled by capillary barriers, based on measured soil hydraulic parameters from previous Tio waste rock studies. Locally-sourced sand in the FCL was used to achieve a desired thickness for the engineered cover system for the full-scale piles (Martin et al., 2017). The numerical model in Figure 3.4 suggests a CBE between the sand and HI material, which results in the majority of water being diverted through the sand layer towards the clean anorthosite material. A subsequent study by 69 Dimech et al. (2019) investigating the effectiveness of the FCL in a field-scale experimental pile for the Tio site also found that fluid diversion occurred mainly through the sand layer, which is consistent with the results shown in the numerical model here, as well as modeling completed by Broda et al. (2017). These findings lend support to the FCL’s ability to divert most water away from the reactive HI material. Additionally, results from the numerical model suggest the presence of a secondary CBE between the HI material and uncrushed anorthosite material, which had not been previously identified. The development of the secondary CBE may be explained by increased water saturations caused by the high infiltration rates applied to the laboratory model for the steady-state flow experiment. The secondary capillary barrier is beneficial in keeping water away from the metal-rich HI material and retaining it within the “clean” anorthosite zone. This interpretation is supported by the agreement between observed and simulated drainage rates (Figure 3.3 and Table 3.5), and also by the geochemical results, discussed in the section below. However, additional tests are required to determine whether this secondary CBE exists for lower infiltration rates. Numerical modeling facilitated to account for (1) heterogeneities in physical flow parameters between different material zones and (2) the pile geometry, owing to the unstructured grid capabilities of the upgraded code. These factors were particularly important to investigate the influence of CBEs in the engineered cover system, since in this case the process of fluid diversion is controlled by contrast in grain size, associated hydraulic parameters and pile geometry. 3.4.2 Implications of Results from Reactive Transport Modeling Both simulated and observed effluent concentrations presented Figure 3.5 support the simulated internal fluid flow pathways depicted in Figure 3.4. It was noted above that the highest outflow rates (both measured and simulated) were observed at reservoir 3. If the majority of this 70 water was leaking past the FCL and moving through the HI material zones (contrary to the simulated flow model), the effluent quality would be expected to lie closer to observations in reservoirs 5 and 4. However, actual measured concentrations from reservoir 3 fall somewhere in between reservoirs 5 – 4 and 2 – 1, supporting the hypothesis that a fraction of the water released at reservoir 3 moves through the HI material, while the remainder percolates primarily through the anorthosite zone and is then redirected upwards towards reservoir 3 following hydraulic gradients. Furthermore, reservoirs 1 and 2 have very similar drainage concentrations with comparatively improved water quality relative to other reservoirs in both modeled and measured results. While this is expected for reservoir 1 which lies below the anorthosite material, reservoir 2 is located under the HI material. This finding supports results of the flow model, suggesting that water draining at reservoir 2 primarily percolates through the anorthosite material, but then moves upwards and to the left towards reservoir 2, driven by the hydraulic gradient. Additionally, the Darcy flux vectors shown in Figure 3.4 are smaller in magnitude in HI S1 towards reservoir 8 compared to simulated flow towards reservoirs 5 and 4. Longer residence times and contact between infiltrating water and HI material are likely responsible for the increased concentrations observed from reservoir 8. However, simulations suggest that physical flow alone was insufficient to capture differences in effluent concentrations observed between reservoirs 8 through 1; in Table 3.6, differences in rate coefficients for each material property zone are summarized. For pyrite, higher oxidative dissolution rates were required in HI S1 compared to S2 and S3 in order to match SO4 concentrations observed at each reservoir. Consequently, rate coefficients were also higher in zone S1 for labradorite and enstatite, although since rate expressions for these minerals (as well as biotite) were pH dependent, the overall dissolution rates were also higher due to greater acidity produced from faster sulfide oxidation rates. Table 3.1 71 shows that measured porosities from the laboratory model varied quite significantly between HI S1, S2, and S3 with values of 0.39, 0.55, and 0.58, respectively. These variations in porosity suggest heterogeneity (or potentially differences in compaction) between each zone, which may have influenced the reactive surface area and weathering rates. Porosity values were lower in zones from left to right following deposition history. It is possible that materials in S1 underwent greater compaction than the other HI material zones due to the location, influencing amount (due to more overall rock), reactivity or availability of sulfide grains in this zone, although it is speculative whether this would have a greater influence over naturally occurring chemical heterogeneities between zones. While variations between laboratory- and field-scale rates were beyond the scope of this study, the rates used in the numerical model were consistent with values from the literature and may be used in the future to improve upon the existing knowledge of reaction rates observed in Tio waste rock. Overall, simulated effluent quality was considered to capture measured effluent chemistries from each reservoir quite well. Similar to the physical heterogeneity being essential for capturing the influences of CBEs, differences in mineral weathering rates between each zone were essential for capturing differences between each reservoir. Results suggest that for infiltration that is not diverted from the FCL, the degree of compaction and grain size may influence the residence time and resulting effluent quality released from these zones. Simulated results would not have been captured with a homogeneous model; this suggest the importance of considering physical and chemical heterogeneity in numerical models going forward, and also the implications of compaction and grain size and their influence on effluent quality released in physical studies. 72 3.5 Conclusions RTM can be a useful tool for quantitative interpretation of experimental data, identifying fluid migration pathways, bracketing reaction kinetics, or forecasting expected solute release rates and mass loadings. However, with many experiments aimed to explain complex hydrogeochemical interactions, a thoughtful approach to understanding the dominant processes in the system and justifiable assumptions and limitations must be considered. In this study, simplifications were carefully made to ensure that the project objectives could be met and to capture the rate-controlling processes to better understand the functioning and effectiveness of the FCL. The assumptions were considered justified by reproducing laboratory data including spatially distributed drainage rates and solute concentrations, while keeping most material properties constant based on measured parameters. Although some simplifications were made, physical and chemical heterogeneities between different material zones were essential for capturing measured results in simulations. The main benefit of the numerical model was its ability to characterize internal flow processes (which are often unknown in many systems given limited data points) in the laboratory model, which improved understanding of the functioning of the proposed design for the WRP. Additionally, simulations identified a secondary CBE between the HI and anorthosite waste rock, which had not been considered previously but may prove useful in diverting water from the HI material. Simulated results from this study can also be used to better quantify mineral weathering rates of the fresh HI material used in the laboratory model for pyrite and the aluminosilicate minerals identified from previous studies. Finally, results demonstrate the capabilities of MIN3P-HPC as a tool for simulating flow and reactive transport processes in complex hydrogeochemical systems such as WRPs. The recent 73 implementation of unstructured grid capabilities into the model allows the user to incorporate more complex features, such as the influence of angled layers like the FCL or material zones following the angle of repose, which may contribute to improved water quality leaving the WRP in this study. Future advancements in computing power and improvements made to complex multi-phase interactions within RTM codes will further improve capabilities of the models to characterize other mine waste stockpiles and may be useful to contrast and improve different techniques for the design planning and storage of mine wastes. 74 Chapter 4: Conclusion The use of site-specific design strategies to minimize mass loadings from waste rock and reduce water treatment requirements is important for mining operations to reduce potential leaching of metals from drainage due to ARD/CND processes. The overarching objective of this thesis was to investigate the use of RTM (specifically using the code MIN3P-HPC) to evaluate the effect of different construction methods and design strategies for WRPs on drainage water quality and mass loadings. Although study-specific conclusions have been presented in Chapter 2 and Chapter 3, the findings from both chapters are summarized here and related to one another. In addition, perspectives on some broader implications for the future use of RTM in long-term mine waste management are provided below. In Chapter 2, scenario analyses were performed through a series of stylistic models and modeling results were compared to each other to investigate (1) the influence of spatial distributions of physical and geochemical parameters within WRPs; (2) the effect of spatial trends and features associated with end-dumping vs. push-dumping construction methods; and (3) the effects of pile scale, with a specific focus on single-lift versus multi-lift piles. These simulations provided a comparison of different waste rock storage techniques and assessed how heterogeneities, structural features and scale influence maximum mass loading rates and duration of poor-quality drainage release from WRPs. Results from this study suggest not only that the spatial distributions of physical and geochemical parameters, including hydraulic conductivity, abundance of sulfide minerals and oxidation rates play an important role in shaping peak and long-term mass loadings from WRPs, but also indicate that heterogeneities and structural features should be included in RTM to improve long-term forecasts of pile behaviour in the real-world. These study results are consistent with previous findings from WRP studies by Atherton (2017) 75 and Pedretti et al. (2017), who identified localized zones of increased acidity generated due to the presence of heterogeneities; however, the current study further highlights the role of structure associated with construction methods and pile size on mass loadings from WRPs. Results from simulations demonstrate that the heterogeneity actually reduces normalized peak mass loading rates in acid-generating conditions; if materials were to weather at the same rate throughout the pile, the effects compound and result in higher maximum concentrations in drainage, as can be seen from results of the homogeneous model. This result is similar to well-known concepts of the effect of heterogeneity in groundwater flow systems regarding the effects of dispersion and diffusion on breakthrough curves. The impact of these heterogeneities will depend on their spatial distributions and magnitude given site-specific considerations, but results enforce the need for representative sampling procedures—an important consideration for long-term mine waste management practices with direct implications on peak water treatment requirements. Additionally, the ability to capture different pile scales using MIN3P-HPC models was demonstrated in Chapter 2 through simulations of single bench vs. multi-bench piles. Results from this comparison demonstrated that scale, spatial trends (such as upward fining) and features such as traffic surfaces are also important to consider. When investigating the influence of internal structure at the single bench scale, differences were observed between push- and end-dumping methods. However, when combined over multiple benches, differences in peak mass loading rates were negligible between construction methods. While this finding emphasizes the importance of capturing scale and structural features within a WRP, the actual impact is contextual and it is important to consider the scope and objectives of the numerical model produced. In some instances, a simplistic 1D homogeneous model may be sufficient, while in other cases, 2D or 3D 76 models with complex structural features and multi-phase hydrogeochemical interactions may be considered essential to provide meaningful forecasts. If the site-specific impact of these internal features is unknown during the design phase of a project, RTM may help guide mine site planning to better understand how aspects of the system fit together. Scenario analyses can be beneficial by providing guidance and insights, when evaluating which parameters have the most significant impact on ARD release, its evolution and longevity. In Chapter 3, MIN3P-HPC was used to conceptualize the effectiveness of the FCL engineered cover system in the laboratory-scale model built by Poaty (2020). In this work, the numerical model was calibrated to reproduce real-world experimental results, specifically basal drainage rates and effluent geochemistry. Modeling results and observations agreed fairly well, constrained by measured material property parameters and weathering rates that were consistent to literature values. In doing so, the numerical model built with MIN3P strengthened and enhanced hypotheses developed from results of the laboratory model—the ability to achieve a good agreement between modeled and measured results reinforced the conceptual model that grew from simulations regarding the ability of the FCL to divert significant amounts of water away from the HI material effectively. Furthermore, these results strengthen the argument that considering heterogeneity and internal structure is important. Although zones were homogeneous in the numerical model used in Chapter 3, differences in weathering rates in the HI zones S1 – S3 proved important to matching effluent geochemistry. 4.1 Limitations of RTM Although this thesis demonstrated that the use of RTMs can be beneficial and provide insights, it must be emphasized that a thoughtful approach must be taken when using numerical models in order for the results to be meaningful. As with solving many problems in hydrogeology, it is 77 important to start with a conceptual model of the system and define boundary conditions and the different controlling components while keeping in mind how these align with the objectives for the problem in question. When creating a numerical model to forecast mass loading rates or bracket expected behaviour, it is important to consider long-term processes that may be occurring and influence release rates such as sorption and precipitation-dissolution reactions when calibrating the model. In Chapter 2, SO4 mass loading rates were used as a proxy to compare forecasted release rates from ARD in different construction methods. This approach was useful as a comparative tool for oxidation rates influenced by heterogeneity and pile structure between different structural configurations; however, for the purposes on this study, oxygen was assumed to remain at atmospheric conditions. While beyond the scope of work and objectives presented here, this assumption may not be considered reasonable for highly reactive waste rock piles, in which oxygen depletion occurs. However, these results demonstrate the role of heterogeneity, structure and to some degree scale, which was the objective of this work. In addition to these considerations, it should be noted that some RTM codes are more capable of capturing certain processes than others. In the case of this thesis, the recently developed capabilities provided with MIN3P-HPC to capture sloped features could be tested in both research chapters. Before implementation of unstructured grid capabilities, simulations completed with MIN3P were limited to the pile core (Wilson et al. 2018a) and restricted to rectangular geometries. It is important to define which processes may be important or required during the planning/conceptualization stage when choosing which RTM code may be most applicable for the problem in question. It can be expected that future upgrades and technological advancements will 78 see improved features available in RTMs and provide more powerful capabilities for the user to solve complex subsurface problems. 4.2 Recommended Future Work In Chapter 2, the stylistic models provided insight into influences of heterogeneity, scale, and structures in RTM. In the future, construction of laboratory-scale models (or ideally, experimental to full-scale piles) built using the same waste rock material designed with the features described for push- and end-dumping methods could be compared to improve knowledge on the effect of construction methods on controlling ARD. Results could also be used to validate numerical models and provide additional insights into which parameters are most important to consider when assessing the performance of waste rock pile systems using RTM methods in the future. For Chapter 3, additional testing with reduced infiltration rates applied to the laboratory-scale model may provide insights into the influence that water saturation plays in CBEs and fluid flow paths within the experimental model. Furthermore, additional simulations where the FCL is composed entirely of sand or crushed anorthosite under various infiltration rates may help clarify discrepancies between the interpretations related to the functioning of the FCL, as presented by Dimech et al. (2019) and shown in the current simulation results (i.e. whether flow primarily moves through the crushed anorthosite or sand layer). 4.3 Implications for Future Use of RTM in Mine Waste Management Beyond the contributions provided individually through Chapter 2 and Chapter 3, the overall results of this thesis demonstrate that RTMs are now able to deal with complex geometries, structure and heterogeneities and can be used for the investigation of mine waste weathering in a number of cases. During the planning stage for pile design or reclamation activities, numerical models could be employed similarly to what was completed here to compare proposed site-specific 79 design strategies. Such simulations could be useful to help determine which methods may have the greatest benefits in reduction/mitigation of ARD/CND generation and release of poor-quality drainage waters. Furthermore, many of the parameters required to constrain numerical models are often already required for standard laboratory/field tests and waste rock characterization. In such cases, models can be considered as a supplementary tool to guide in the planning process for long-term waste rock storage with limited additional sampling requirements. Chapter 3 provides a relevant example in this context; the parameters utilized in the numerical model were taken from previous characterization studies done for the Tio mine waste rock. In this instance, an additional steady-state experiment was performed re-using a pre-existing laboratory-scale model to meet the objectives of the numerical study. However, the design of the numerical study using RTM may in some cases be tailored to match data available or with minimum additional requirements to match the objectives of the investigation. Similar to laboratory or field studies, which aim to classify and understand waste rock systems, the samples collected, testing required, and assumptions made for numerical models will be site-specific in order to match the objectives of the modeling. In addition, this thesis has demonstrated that recent improvements in computing power and parallel computing capabilities, as available in MIN3P-HPC, have increased the utility of RTMs for mine waste management. Models presented here were fairly complex; in the case of Chapter 2, multiple configurations simulated over several hundred years were run in only a few hours of wall clock time. Although the geochemical systems were simplified, these runtimes are a significant improvement in comparison to previous efforts and demonstrate that advancements made in recent years, such as those made possible by the parallelized version of MIN3P-HPC, allow simulations to be run on a workstation or supercomputer to achieve accelerated runtimes using multiple processors. As RTM becomes more accessible and more widely used for scenario 80 analyses or as a data interpretation tool, optimization of required parameters to successfully characterize the system will become clearer and more integrated with current standard laboratory/field practices. The work presented in this thesis has demonstrated some applications of numerical modeling tools for conceptualization and planning of mine waste rock stockpiles. As technological advancements continue to improve, RTM will be available as a flexible tool to help specialists in developing improved understanding of complex hydrogeochemical systems and to improve the design and performance of long-term storage facilities, by reducing environmental liabilities and water treatment requirements associated with long-term mine waste storage. 81 Bibliography Akcil A, Koldas S (2006) Acid Mine Drainage (AMD): causes, treatment and case studies. J Clean Prod 14:1139–1145 Allison, J. D., W. Brown and K.J. Novo-Gradac (1991). MINTEQA2/PRODEFA2, A Geochemical Assessment Model for Environmental Systems: Version 3.0 User’s Manual. Environmental research laboratory, U.S Enviromentla protection Agency, EPA/600/3-91/021 Amos RT, Blowes DW, Bailey BL, Sego DC, Smith LJD, Ritchie IAM (2015) Waste-rock hydrogeology and geochemistry. Appl. Geochemistry 57:140–156. https://doi.org/10.1016/j.apgeochem.2014.06.020 Amos RT, Blowes DW, Smith L, Sego DC (2009) Measurement of Wind-Induced Pressure Gradients in a Waste Rock Pile. Vadose Zo J 8:953–962. https://doi.org/10.2136/vzj2009.0002 Amos RT, Wilson D, Bailey BL, Blowes DW, Sego DC, Smith L (2018) Conceptual Overview of Mechanistic Scale-up in Waste Rock: Lessons Learned from Diavik and Potential for Application to Other Mine Sites. In: The 25th British Columbia MEND Metal Leaching/Acid Rock Drainage Workshop. Vancouver, Canada, 28-29 November 2018 Anterrieu O, Chouteau M, Aubertin M (2010) Geophysical characterization of the large-scale internal structure of a waste rock pile from a hard rock mine. Bull Eng Geol Environ 69:533–548. https://doi.org/10.1007/s10064-010-0264-4 Appels WM, Ireson AM, Barbour SL (2018) Impact of bimodal textural heterogeneity and connectivity on flow and transport through unsaturated mine waste rock. Adv Water Resour 112:254–265. https://doi.org/10.1016/j.advwatres.2017.12.008 82 Aranda CA (2010) Assessment of waste rock weathering characteristics at the Antamina mine based on field cells experiment. Master's Thesis. University of British Columbia, Vancouver, BC Atherton C (2017) An Investigation of Heterogeneity and the Impact of Acidic Regions on Bulk Effluent from a Deconstructed Low Sulfide Waste-Rock Pile. Master's Thesis. University of Waterloo, Waterloo, ON Aubertin M, Cifuentes E, Apithy SA, Bussière B, Molson J, Chapuis RP (2009) Analyses of water diversion along inclined covers with capillary barrier effects. Can Geotech J 46:1146–1164. https://doi.org/10.1139/T09-050 Azam S, Wilson WG, Herasymuik G, Nichol C, Barbour LS (2007) Hydrogeological behaviour of an unsaturated waste rock pile: A case study at the Golden Sunlight Mine, Montana, USA. Bull Eng Geol Environ 66:259–268. https://doi.org/10.1007/s10064-006-0077-7 Barsi D (2017) Spatial Variability of Particles in Waste Rock Piles. Master's Thesis. University of Alberta, Edmonton, AB Bea SA, Mayer KU, MacQuarrie KTB (2016) Reactive transport and thermo-hydro-mechanical coupling in deep sedimentary basins affected by glaciation cycles: model development, verification, and illustrative example. Geofluids 16:279–300. https://doi.org/10.1111/gfl.12148 Bellehumeur TM (2001) Mechanisms and Spatial Variability of Rainfall Infiltration on the Claude Waste Rock Pile. Master's Thesis. University of British Columbia, Vancouver, BC Blackmore S, Pedretti D, Mayer KU, Smith L, Beckie RD (2018a) Evaluation of single- and dual-porosity models for reproducing the release of external and internal tracers from heterogeneous waste-rock piles. J Contam Hydrol 214:65–74. https://doi.org/10.1016/j.jconhyd.2018.05.007 83 Blackmore S, Smith L, Ulrich Mayer K, Beckie RD (2014) Comparison of unsaturated flow and solute transport through waste rock at two experimental scales using temporal moments and numerical modeling. J Contam Hydrol 171:49–65. https://doi.org/10.1016/j.jconhyd.2014.10.009 Blackmore S, Vriens B, Sorensen M, Power IM, Smith L, Hallam SJ, Mayer KU, Beckie RD (2018b) Microbial and geochemical controls on waste rock weathering and drainage quality. Sci Total Environ 640–641:1004–1014. https://doi.org/10.1016/j.scitotenv.2018.05.374 Blowes DW, Ptacek CJ, Jambor JL, Weisener CG (2003) The Geochemistry of Acid Mine Drainage. In: Mineralogy of Ore Deposits. Elsevier Inc., pp 149–204 Blum AE, Stillings LL (1995) Chemical weathering rates of silicate minerals. Rev Mineral 31:291–351 Brantley SL, Chen Y (1995) Chemical Weathering Rates of Pyroxenes and Amphiboles. Rev Mineral 31:119–172 Bréard-Lanoix ML (2017) Characterization of hydrogeological properties of the flow control layer constructed on the top of the experimental waste rock pile in the Lac Tio mine (in French). Master's Thesis. Polytechnique Montréal, Montréal, QC Canadian Mining Innovation Council (2013) Environmental Analysis of the Mining Industry in Canada. http://www.cmic-ccim.org/wp-content/uploads/2013/07/HatchScopingReport.pdf. Accessed March 26, 2020 Cash A (2014) Structural and Hydrologic Characterization of Two Historic Waste Rock Piles. Master's Thesis. University of Alberta, Edmonton, AB Chapuis RP (2012) Predicting the saturated hydraulic conductivity of soils: a review. Bull Eng Geol Env 71:401–434. https://doi.org/10.1007/s10064-012-0418-7 84 Cornell RM, Schwertmann U (2003) The Iron Oxides, 2nd edn. Wiley-VCH. ISBN: 3-527-30274-3 Demers I, Molson JW, Bussière B, Laflamme D (2013) Numerical modeling of contaminated neutral drainage from a waste-rock field test cell. Appl Geochemistry 33:346–356. https://doi.org/10.1016/j.apgeochem.2013.02.025 Dimech A, Chouteau M, Aubertin M, Bussière B, Martin V, Plante B (2019) Three‐Dimensional Time‐Lapse Geoelectrical Monitoring of Water Infiltration in an Experimental Mine Waste Rock Pile. Vadose Zo J 18:1–19. https://doi.org/10.2136/vzj2018.05.0098 Dimech A, Chouteau M, Aubertin M, Bussière B, Martin V, Plante B (2018) 3D Time-lapse geoelectrical monitoring of moisture content in an experimental waste rock pile: Validation using hydrogeological data. Proc Symp Appl Geophyics to Eng Environ Probl SAGEEP 2018-March: https://doi.org/10.4133/sageep.31-009 Dimech A, Chouteau M, Chou ÉTK, Aubertin M, Martin V, Bussière B, Plante B (2017) Monitoring water infiltration in an experimental waste rock pile with time-lapse ERT and multi-parameter data collection. Proc Symp Appl Geophyics to Eng Environ Probl SAGEEP. https://doi.org/10.4133/sageep.30-009 Dold B (2017) Acid rock drainage prediction: A critical review. J Geochemical Explor 172:120–132. https://doi.org/10.1016/j.gexplo.2016.09.014 Dominion Diamond Mines (2014) Ekati Diamond Mine Waste Rock and Ore Storage Management Plan Version 4.1, May 2014 Dubuc J, Pabst T, Aubertin M (2013) An Assessment of the Hydrogeological Response of the Flow Control Layer Installed on the Experimental Waste Rock Pile At the Lac Tio Mine. 2–7 85 Elghali A, Benzaazoua M, Bussière B, Bouzahzah H (2019) Determination of the available acid-generating potential of waste rock, part II: Waste management involvement. Appl Geochemistry 100:316–325. https://doi.org/10.1016/j.apgeochem.2018.12.010 Engineering ToolBox (2004) Evaporation from Water Surface. https://www.engineeringtoolbox.com/evaporation-water-surface-d_690.html. Accessed 13 Feb 2020 Engineering ToolBox (2003) Mollier Diagram. https://www.engineeringtoolbox.com/psychrometric-chart-mollier-d_27.html. Accessed 13 Feb 2020 Erguler ZA, Kalyoncu Erguler G (2015) The effect of particle size on acid mine drainage generation: Kinetic column tests. Miner Eng 76:154–167. https://doi.org/10.1016/j.mineng.2014.10.002 Fala O, Aubertin M, Molson JW, Bussière B, Wilson GW, Chapuis R, Martin V (2003) Numerical Modelling of Unsaturated Flow in Uniform and Heterogeneous Waste Rock Piles. 6th International Conference of Acid Rock Drainage 12–18 Fala O, Molson JW, Aubertin M, Dawood I, Bussière B, Chapuis RP (2013) A numerical modelling approach to assess long-term unsaturated flow and geochemical transport in a waste rock pile. Int J Mining, Reclam Environ 27:38–55. https://doi.org/10.1080/17480930.2011.644473 Fala O, Molson JW, Aubertin M, Bussière B, Chapuis RP (2006) Numerical simulations of long term unsaturated flow and acid mine drainage at waste rock piles. 7th International Conference of Acid Rock Drainage 1:582–597. https://doi.org/10.21000/jasmr06020582 86 Fala O, Molson JW, Aubertin M, Bussière B (2005) Numerical modelling of flow and capillary barrier effects in unsaturated waste rock piles. Mine Water Environ 24:172–185. https://doi.org/10.1007/s10230-005-0093-z Fines PE (2006) Hydrologic Characterization of Two Full-Scale Waste Rock Piles. Master's Thesis. University of British Columbia, Vancouver, BC Freeze A, Cherry J (1979) Groundwater. Prentice-Hall Haug MD (2002) A Review of of Non-traditional Dry Cover (MEND Report 2.21.3b). 142 Herasymuik GM, Azam S, Wilson GW, Barbour LS, Nichol C (2006) Hydrological characterization of an unsaturated waste rock dump. In: Proceedings of the 59th Canadian geotechnical conference, Vancouver, BC, Canada INAP (2012) Chapter 2: The Acid Rock Drainage Process (GARD Guide). http://www.gardguide.com/index.php/Chapter_2 INAP (2009) GARD Guide. http://www.gardguide.com/index.php?title=Main_Page. Accessed 11 Feb 2020 Intissar R (2009) Apport des méthodes électriques à la détermination de la structure interne d`une halde à stériles. Master's Thesis. Polytechnique de Montréal, Montréal, QC Jambor JL, Dutrizac JE, Raudsepp M (2007) Measured and computed neutralization potentials from static tests of diverse rock types. Environ Geol 52:1173–1185. https://doi.org/10.1007/s00254-006-0555-z Jodeiri Shokri B, Doulati Ardejani F, Ramazi H, Moradzadeh A (2016) Predicting pyrite oxidation and multi-component reactive transport processes from an abandoned coal waste pile by comparing 2D numerical modeling and 3D geo-electrical inversion. Int J Coal Geol 164:13–24. https://doi.org/10.1016/j.coal.2016.03.014 87 Johnson DB, Hallberg KB (2005) Acid mine drainage remediation options: A review. Sci Total Environ 338:3–14. https://doi.org/10.1016/j.scitotenv.2004.09.002 Jung H, Navarre-Sitchler A (2018) Physical heterogeneity control on effective mineral dissolution rates. Geochim Cosmochim Acta 227:246–263. https://doi.org/10.1016/j.gca.2018.02.028 Kempton H (2012) A review of scale factors for estimating waste rock weathering from laboratory tests. In: Proceedings from the 9th International Conference on Acid Rock Drainage. Ottawa, Canada, 20-26 May 2012 Knight Piésold Ltd. (2012) Waste Dumps and Stockpiles Preliminary Design. Taseko Mines Limited New Prosperity Gold-Copper Project Komnitsas K, Xenidis A, Adam K (1995) Oxidation of pyrite and arsenopyrite in sulphidic spoils in Lavrion. Miner Eng 8:1443–1454. https://doi.org/10.1016/0892-6875(95)00109-3 Kool JB, Parker JC, van Genuchten MT (1985) Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: I. Theory and numerical studies. Soil Sci Soc Am J 49:1348–1354 Lahmira B, Lefebvre R, Aubertin M, Bussière B (2016) Effect of heterogeneity and anisotropy related to the construction method on transfer processes in waste rock piles. J Contam Hydrol 184:35–49. https://doi.org/10.1016/j.jconhyd.2015.12.002 Lahmira B, Lefebvre R, Aubertin M, Bussière B (2017) Effect of material variability and compacted layers on transfer processes in heterogeneous waste rock piles. J Contam Hydrol 204:66–78. https://doi.org/10.1016/j.jconhyd.2017.07.004 Langman J, Moore M, Ptacek C, Smith L, Sego DC, Blowes DW (2014) Diavik Waste Rock Project: Evolution of Mineral Weathering, Element Release, and Acid Generation and 88 Neutralization during a Five-Year Humidity Cell Experiment. Minerals 4:257–278. https://doi.org/10.3390/min4020257 Lawrence RW, Wang Y (1996) Determination of neutral potential for acid rock drainage prediction. MEND Proj 1.16.3 1-89 Lefebvre R, Hockley D, Smolensky J, Geélinas P (2001) Multiphase transfer processes in waste rock piles producing acid mine drainage. 1: Conceptual model and system characterization. J Contam Hydrol 52:137–164. https://doi.org/10.1016/S0169-7722(01)00156-5 Lehner S, Savage K (2008) The effect of As, Co, and Ni impurities on pyrite oxidation kinetics: Batch and flow-through reactor experiments with synthetic pyrite. Geochim Cosmochim Acta 72:1788–1800. https://doi.org/10.1016/j.gca.2008.02.003 Li MG (2000) Acid rock drainage prediction for low-sulfide, lowneutralization potential mine wastes. In: Proceedings of the 5th International Conference on Acid Rock Drainage. Denver, CO, USA, pp 567–580 Lichtner PC (1996) Continuum formulation of multicomponent–multiphase reactive transport. In: Lichtner PC, Steefel CI, Oelkers EH (eds) Reactive Transport in Porous Media, Reviews in Mineralogy. pp 1–81 Linklater CM, Sinclair DJ, Brown PL (2005) Coupled chemistry and transport modelling of sulphidic waste rock dumps at the Aitik mine site, Sweden. Appl Geochemistry 20:275–293. https://doi.org/10.1016/j.apgeochem.2004.08.003 Lorca ME, Mayer KU, Pedretti D, Smith L, Beckie RD (2016) Spatial and Temporal Fluctuations of Pore-Gas Composition in Sulfidic Mine Waste Rock. Vadose Zo J 15:vzj2016.05.0039. https://doi.org/10.2136/vzj2016.05.0039 89 Malmström ME, Destouni G, Banwart SA, Strömberg BHE (2000) Resolving the scale-dependence of mineral weathering rates. Environ Sci Technol 34:1375–1378. https://doi.org/10.1021/es990682u Martin V, Bussière B, Plante B, Pabst B, Aubertin M, Medina F, Bréard-Lanoix ML, Dimech A, Dubuc J, Poaty B, Wu R, Mackenzie J, Broda S, Chen D (2017) Controlling water infiltration in waste rock piles: Design, construction, and monitoring of a large-scale in-situ pilot test pile. Proc 70th Can Geotech Soc Conf Mayer KU, Alt-Epping P, Jacques D, Arora B, Steefel CI (2015) Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage. Comput Geosci 19:599–611. https://doi.org/10.1007/s10596-015-9476-9 Mayer KU, Amos RT, Blowes DW, Wilson D, Bea SA, Dipple GM, Romano C, Su D, Raymond K (2019) Reactive Transport Modeling in Mine Waste Management – Part 1: Examination of its Role Based on Past and Present Applications. In: Proceedings of the 26th Annual BC/MEND Metal Leaching/Acid Rock Drainage Workshop. Vancouver, BC Mayer KU, Frind EO, Blowes DW (2002) Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water Resour Res 38:13-1-13–21. https://doi.org/10.1029/2001wr000862 McCarter MK (1990) Design and operating considerations for mine waste embankments. In: Surface Mining, 2nd edn. Littleton, CO, pp 890–899 McLemore VT, van Zyl D, Ayakwah GF, Anim K, Boakye K, Ennin F, Felli P, Fredlund DG, Gutierrez LAF, Nunoo S, Tachie-Menson S, Viterbo VC (2009) Literature review of other rock piles: characterization, weathering, and stability. Questa Rock Pile Weathering Stability Project. New Mex Bur Geol Miner Resour OF-Report 517 90 Mele M, Servida D, Lupis D (2013) Characterisation of sulphide-bearing waste-rock dumps using electrical resistivity imaging: The case study of the Rio Marina mining district (Elba Island, Italy). Environ Monit Assess 185:5891–5907. https://doi.org/10.1007/s10661-012-2993-2 MEND (2004a) Design, Construction, and Performance Monitoring of Cover Systems for Waste Rock and Tailings Volume 1. MEND Report 2.21.4a MEND (1998) Blending and Layering Waste Rock to Delay, Mitigate, or Prevent Acid Rock Drainage and Metal Leaching: A case study review. Natural Resources Canada, Ottawa, Canada MEND (2004b) Review of Water Quality Issues in Neutral pH Drainage : Examples and Emerging Priorities for the Mining Industry in Canada. 58. MEND Report 10.1 Miller S, Rusdinar Y, Smart R, Andrina J, Richards D (2006) Design and Construction of Limestone Blended Waste Rock Dumps – Lessons Learned from a 10-Year Study at Grasberg. In: R.I. Barnhisel (Ed.), Proceedings of 7th International Conference on Acid Rock Drainage, March 26-30,. St. Louis, MO, American Society of Mining and Reclamation, Lexington, KY. Molson JW, Fala O, Aubertin M, Bussière B (2005) Numerical simulations of pyrite oxidation and acid mine drainage in unsaturated waste rock piles. J Contam Hydrol 78:343–371. https://doi.org/10.1016/j.jconhyd.2005.06.005 Morin KA, Gerencher E, Jones CE, Konasewich DE (1991) Critical Literature Review of Acid Drainage from Waste Rock. MEND Rep 1111 175. https://doi.org/http://dx.doi.org/10.1016/j.jct.2009.07.005 Munroe EA, Mclemore VT, Kyle P (1999) Waste rock pile characterization, heterogeneity, and geochemical anomalies in the Hillsboro Mining District, Sierra County, New Mexico. J Geochemical Explor 67:391–405 91 Nagy KL (1995) Dissolution and Precipitation Kinetics of Sheet Silicates. Rev Mineral 31:173–233 Neuner M, Smith L, Blowes DW, Sego DC, Smith, LJD, Fretz N, Gupton M (2013) The Diavik waste rock project: Water flow through mine waste rock in a permafrost terrain. Appl geochemistry 36:222–233. https://doi.org/10.1016/j.apgeochem.2012.03.011 Nichols RS (1986) Rock Segregation in Waste Dumps. In: Proceedings of the International Symposium on Flow-through Rock Drains 105–120 Nicholson RV, Gillham RW, Reardon EJ (1988) Pyrite oxidation in carbonate-buffered solution: 1. Experimental kinetics. Geochim Cosmochim Acta 52:1077–1085. https://doi.org/10.1016/0016-7037(88)90262-1 Nordstrom DK (1982) Aqueous pyrite oxidation and the consequent formation of secondary iron minerals. Acid Sulfate Weather 37–56. https://doi.org/10.2136/sssaspecpub10.c3 Nordstrom DK, Alpers CN (1999) Geochemistry of Acid Mine Waters. In: Plumlee GS, Logsdon MJ, Filipek LF (eds) The Environmental Geochemistry of Mineral Deposits, Part A: Processes, Techniques and Health Issues. Reviews in Economic Geology Vol. 6A, Society of Economic Geologists, Inc. pp 133–160 Nordstrom DK, Blowes DW, Ptacek CJ (2015) Hydrogeochemistry and microbiology of mine drainage: An update. Appl Geochemistry 57:3–16. https://doi.org/10.1016/j.apgeochem.2015.02.008 O’Kane M, Wilson GW, Barbour SL (1998) Instrumentation and monitoring of an engineered soil cover system for mine waste rock. Can Geotech J 35:828–846 Oelkers EH, Schott J (2001) An experimental study of enstatite dissolution rates as funcion of pH, temperature, and aqueous Mg and Si concentration, and the mechanism of 92 pyroxene/pyroxenoid dissolution. Geochim Cosmochim Acta 65:1219–1231. https://doi.org/10.1016/S0016-7037(00)00564-0 Parbhakar-Fox A, Lottermoser B, Bradshaw D (2013) Evaluating waste rock mineralogy and microtexture during kinetic testing for improved acid rock drainage prediction. Miner Eng 52:111–124. https://doi.org/10.1016/j.mineng.2013.04.022 Parbhakar-Fox A, Lottermoser BG (2015) A critical review of acid rock drainage prediction methods and practices. Miner Eng 82:107–124. https://doi.org/10.1016/j.mineng.2015.03.015 Pedretti D, Mayer KU, Beckie RD (2017) Stochastic multicomponent reactive transport analysis of low quality drainage release from waste rock piles: Controls of the spatial distribution of acid generating and neutralizing minerals. J Contam Hydrol 201:30–38. https://doi.org/10.1016/j.jconhyd.2017.04.004 Peregoedova A (2012) Étude expérimentale des propriétés hydrogéologiques des roches stériles à une échelle intermédiaire de laboratoire. Master's Thesis. Polytechnique Montréal, Montréal, QC Peterson HE (2014) Unsaturated Hydrology, Evaporation, and Geochemistry of Neutral and Acid Rock Drainage in Highly Heterogeneous Mine Waste Rock at the Antamina Mine, Peru. Doctoral Thesis. University of British Columbia, Vancouver, BC Piteau Associates Engineering Ltd. (1991) Investigation and Design of Mine Dumps. British Columbia Mine Waste Rock Pile Research Committee Plante B, Benzaazoua M, Bussière B (2011a) Kinetic Testing and Sorption Studies by Modified Weathering Cells to Characterize the Potential to Generate Contaminated Neutral Drainage. Mine Water Environ 30:22–37. https://doi.org/10.1007/s10230-010-0131-3 93 Plante B, Benzaazoua M, Bussière B, Biesinger MC, Pratt AR (2010) Study of Ni sorption onto Tio mine waste rock surfaces. Appl Geochemistry 25:1830–1844. https://doi.org/10.1016/j.apgeochem.2010.09.010 Plante B, Benzaazoua M, Bussière B (2011b) Predicting Geochemical Behaviour of Waste Rock with Low Acid Generating Potential Using Laboratory Kinetic Tests. Mine Water Environ 30:2–21. https://doi.org/10.1007/s10230-010-0127-z Plante B, Bussière B, Benzaazoua M (2014) Lab to field scale effects on contaminated neutral drainage prediction from the Tio mine waste rocks. J Geochemical Explor 137:37–47. https://doi.org/10.1016/j.gexplo.2013.11.004 Poaty B (2020) Comportement hydrogéochimique des stériles d’une halde expérimentale visant à limiter la production du drainage neutre contaminé par contrôle des écoulements. Doctoral Thesis. Université du Québec en Abitibi-Témiscamingue, Rouyn-Noranda, QC Poisson J, Chouteau M, Aubertin M, Campos D (2009) Geophysical experiments to image the shallow internal structure and the moisture distribution of a mine waste rock pile. J Appl Geophys 67:179–192. https://doi.org/10.1016/j.jappgeo.2008.10.011 Price WA (2009) Prediction manual for drainage chemistry from sulphidic geologic materials. CANMET Mining and Mineral Sciences Laboratories. Smithers, BC 1-579 Price WA, Errington JC (1998) Guidelines for Metal Leaching and Acid Rock Drainage At Minesites in British Columbia. BC Ministry of Energy Mines. Smithers, BC R Documentation rnorm function. https://www.rdocumentation.org/packages/compositions/versions/1.40-3/topics/rnorm. Accessed 21 Mar 2020 94 Sapsford DJ, Bowell RJ, Dey M, Williams KP (2009) Humidity cell tests for the prediction of acid rock drainage. Miner Eng 22:25–36. https://doi.org/10.1016/j.mineng.2008.03.008 Seigneur N, Mayer KU, Steefel CI (2019) Reactive Transport in Evolving Porous Media. Rev Mineral Geochemistry 85:197–238. https://doi.org/10.2138/rmg.2019.85.7 Sherlock EJ, Lawrence RW, Poulin R (1995) On the neutralization of acid rock drainage by carbonate and silicate minerals. Environ Geol 25:43–54. https://doi.org/10.1007/BF01061829 Šimůnek J, van Genuchten MT (2008) Modeling Nonequilibrium Flow and Transport Processes Using HYDRUS. Vadose Zone J 7:782–797. https://doi.org/10.2136/vzj2007.0074 Šimůnek J, van Genuchten MT, Šejna M (2008) Development and Applications of the HYDRUS and STANMOD Software Packages and Related Codes. Vadose Zone J 7:587–600. https://doi.org/10.2136/vzj2007.0077 Smith LJ, Beckie RD (2003) Hydrologic and geochemical transport processes in mine waste rock. Environ Asp Mine Wastes Short Course Ser Miner Assoc 51–72 Smith LJ, Lopez DL, Beckie RD, Morin KA, Dawson RF, Price WA (1995) Hydrogeology of Waste Rock Dumps. MEND Assoc Proj PA-1 Sobek AA, Schuller WA, Freeman JR, Smith RM (1978) Field and laboratory methods applicable to overburdens and minesoils. EPA-600/2-78-054. Cincinnati, Ohio Steefel CI, Appelo CAJ, Arora B, Jacques D, Kalbacher T, Kolditz O, Lagneau V, Lichtner PC, Mayer KU, Meeussen JCL, Molins S, Moulton D, Shao H, Šimůnek J, Spycher N, Yabusaki SB, Yeh GT (2015a) Reactive transport codes for subsurface environmental simulation. Comput Geosci 19, 445–478. https://doi.org/10.1007/s10596-014-9443-x 95 Steefel CI, DePaolo DJ, Lichtner PC (2005) Reactive transport modeling: An essential tool and a new research approach for the Earth sciences. Earth Planet Sci Lett 240:539–558. https://doi.org/10.1016/j.epsl.2005.09.017 Steefel CI, Yabusaki SB, Mayer KU (2015b) Reactive transport benchmarks for subsurface environmental simulation. Comput Geosci 19:439–443. https://doi.org/10.1007/s10596-015-9499-2 Stoch L, Sikora W (1976) Transformations of micas in the process of kaolinitization of granites and gneisses. Clays Clay Miner 24:. https://doi.org/10.1346/ccmn.1976.0240402 Stockwell J, Smith L, Jambor JL, Beckie RD (2006) The relationship between fluid flow and mineral weathering in heterogeneous unsaturated porous media: A physical and geochemical characterization of a waste-rock pile. Appl Geochemistry 21:1347–1361 StrataGold Corporation (July, 2014) Eagle Gold Project: Waste Rock and Overburden Facility Management Plan, Version 2014-01. Stumm W, Morgan JJ (1996) Aquatic chemistry: chemical equilibria and rates in natural waters, 3rd edn. John Wiley and Sons, New York Su D, Mayer KU, MacQuarrie KTB (submitted) MIN3P-HPC: a high performance unstructured grid code for subsurface flow and reactive transport simulation. Submitted to Math Geosci Su D, Xie M, Mayer KU, Kerry TB (2017) MIN3P-THCm User Manual for Parallel Configuration (Draft). User Guide. Su D, Xie M, Mayer KU, Macquarrie KTB (2018) MIN3P-THCm-USG User Manual (Draft) Tietje O, Hennings V (1996) Accuracy of the saturated hydraulic conductivity prediction by pedo-transfer functions compared to the variability within FAO textural classes. Geoderma 69:71–84. https://doi.org/10.1016/0016-7061(95)00050-X 96 Tran B, Miller S, Williams DJ, Fines P, Wilson GW (2003) Geochemical and Mineralogical Characterisation of Two Contrasting Waste Rock Dumps — The INAP Waste Rock Dump Characterisation Project. 6th Int Conf Acid Rock Drain 939–945 Tremblay GA, Hogan CM (2001) Mine Environment Neutral Drainage (MEND) Manual 5.4.2d: Prevention and Control, Report 5.4.2. 352 US Department of Commerce, NOAA (2019) Earth System Research Laboratory Global Monitoring Division - Global Greenhouse Gas Reference Network. Average CO2 for the Month of February 2019 from Mauna Loa, Hawaii. https://www.esrl.noaa.gov/gmd/dv/data/index.php?category=Greenhouse%2BGases&frequency=Monthly%2BAverages&parameter_name=Carbon%2BDioxide. Accessed July 27 2019 Vriens B, Arnault M St., Laurenzi L, Smith L, Mayer KU, Beckie RD (2018) Localized Sulfide Oxidation Limited by Oxygen Supply in a Full-Scale Waste-Rock Pile. Vadose Zo J 17:180119. https://doi.org/10.2136/vzj2018.06.0119 Vriens B, Peterson H, Laurenzi L, Smith L, Aranda C, Mayer KU, Beckie, RD (2019) Long-term monitoring of waste-rock weathering at the Antamina mine, Peru. Chemosphere 215:858–869. https://doi.org/10.1016/j.chemosphere.2018.10.105 Wilson D, Amos RT, Blowes DW, Langman JB, Smith L, Sego DC (2018a) Diavik Waste Rock Project: Scale-up of a reactive transport model for temperature and sulfide-content dependent geochemical evolution of waste rock. Appl Geochemistry 96:177–190. https://doi.org/10.1016/j.apgeochem.2018.07.001 Wilson D, Amos RT, Blowes DW, Langman JB, Ptacek CJ, Smith L, Sego DC (2018b) Diavik waste rock project: A conceptual model for temperature and sulfide-content dependent 97 geochemical evolution of waste rock – Laboratory scale. Appl Geochemistry 89:160–172. https://doi.org/10.1016/j.apgeochem.2017.12.007 Wilson GW (2011) Rock Dump Hydrology: An overview of full‐scale excavations and scale‐up experiments conducted during the last two decades. In: Bell LC, Braddock B (eds) Proceedings of the Seventh Australian Workshop on Acid and Metalliferous Drainage. 307–322 Wilson GW (2017) Management of Acid Rock Drainage. In: Hawley PM, Cunning J (eds) Guidelines for Mine Waste Dump and Stockpile Design. CSIRO Publishing, Australia, 259–268 Xie M, Su D, Mayer KU, Macquarrie K (2015) A Three-dimensional Numerical Model for Multicomponent Reactive Transport in Variably Saturated Porous Media: Verification and emonstration reports (Draft). 1–165 Yanful EK, Woyshner MR (1993) Design of a composite soil cover for an experimental waste rock pile near Newcastle, New Brunswick, Canada. Can Geotech J 30:578–587 Zakaznova-Herzog VP, Nesbitt HW, Bancroft GM, Tse JS (2008) Characterization of leached layers on olivine and pyroxenes using high-resolution XPS and density functional calculations. Geochim Cosmochim Acta 72:69–86. https://doi.org/10.1016/j.gca.2007.09.031 Zhang S (2017) Relationship between Particle Size Distribution and Porosity in Dump Leaching. Master's Thesis. University of British Columbia, Vancouver, BC 98 Appendices Appendix A Supplementary Material for Chapter 2 A.1 Domain Geometry, Boundary and Initial Conditions Figure A.1. Single-lift structure model mesh configurations with traffic surfaces and zones. With the increased complexity added to structure models, the average cell area is reduced to 0.28 m2 and 0.32 m3 for push- and end-dumping models, respectively. A similar mesh sizing was used for the multibench models, however given the larger scale the mesh becomes unclear at this scale. Therefore figures of the multibench model mesh configurations have been omitted. A.2 Parameter Heterogeneity The rnorm() function was applied to generate a random distribution of grain size diameters, following format in Equation 1 below. Using rnorm(), a list of random grain size values are generated when a mean and variance are provided (5 mm and 3.87 mm, respectively). An example node distribution as a function of pile elevation (z) is shown below. 𝑟𝑟𝑛𝑛𝑟𝑟𝑟𝑟𝑟𝑟(𝑥𝑥,𝑟𝑟𝑚𝑚𝑚𝑚𝑛𝑛, 𝑣𝑣𝑚𝑚𝑟𝑟) (𝐴𝐴. 1) 99 Figure A.2. Nodal values for mean grain size diameter distributed for a random WRP simulation case For the fining upwards case, the following were used to generate a random distribution influenced by pile height. An example node distribution for a single lift model used for the random, fining upwards trends is shown below. A similar approaching using rnorm() was taken, however to account for the effects of fining upwards, the slope was used in the scaling factor as a function of elevation calculated using the standard slope equation below, resulting in the adjusted rnorm() function shown in Equation A.3. 𝑟𝑟 = (𝑑𝑑𝑚𝑚 𝑚𝑚𝑖𝑖𝑚𝑚 − 𝑑𝑑𝑚𝑚 𝑚𝑚𝑚𝑚𝑚𝑚)(𝑧𝑧𝑚𝑚𝑖𝑖𝑚𝑚 − 𝑧𝑧𝑚𝑚𝑚𝑚𝑚𝑚) (𝐴𝐴. 2) 𝑟𝑟𝑛𝑛𝑟𝑟𝑟𝑟𝑟𝑟(𝑥𝑥,𝑟𝑟𝑚𝑚𝑚𝑚𝑛𝑛, 𝑣𝑣𝑚𝑚𝑟𝑟) + (𝑟𝑟 × 𝑧𝑧 × 0.5) (𝐴𝐴. 3) This slope was summed to the generation of random grain size distributions to generate a nodal distribution like the following example: 100 Figure A.3. Nodal values for mean grain size diameter distributed for a fining upwards WRP simulation case These distributions were then shifted to confirm that 𝑑𝑑𝑚𝑚 was 5 mm before being used to calculate other parameters. Histograms showing the lognormal distributions of physical and chemical parameters are given below in Figure A.4. 101 Figure A.4. Initial condition nodal distributions for conceptual models. Lognormal distribution assigned to nodes for hydraulic conductivity, volume fraction, and pyrite rate, and a normal distribution for porosity. Five cases were evaluated for each, which are shown in the overlapping histograms. The left-side histograms show parameter distributions for random trend piles, while the right-side are the random, fining upwards trends (i.e. the heterogeneous structure models). 102 A.3 Heterogeneity in NP:AP Ratios Additional details regarding the AP:NP ratios according to the Sobek (1978) method (standard ABA method) have been included here for relevance to industry. The calculation for acid potential (AP) and neutralizing potential (NP) for the values used in the models in this study are given below (Lawrence and Wang 1996): 𝐴𝐴𝐴𝐴 = 𝑆𝑆𝑆𝑆𝑆𝑆𝑓𝑓𝑆𝑆𝑟𝑟 𝑐𝑐𝑟𝑟𝑛𝑛𝑐𝑐𝑚𝑚𝑛𝑛𝑐𝑐 (%) × 1000 𝑘𝑘𝑘𝑘100 × 𝑟𝑟𝑟𝑟𝑆𝑆𝑚𝑚𝑐𝑐𝑆𝑆𝑆𝑆𝑚𝑚𝑟𝑟 𝑤𝑤𝑚𝑚𝑖𝑖𝑘𝑘ℎ𝑐𝑐 𝑟𝑟𝑓𝑓 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3𝑚𝑚𝑐𝑐𝑟𝑟𝑟𝑟𝑖𝑖𝑐𝑐 𝑤𝑤𝑚𝑚𝑖𝑖𝑘𝑘ℎ𝑐𝑐 𝑟𝑟𝑓𝑓 𝑠𝑠𝑆𝑆𝑆𝑆𝑓𝑓𝑆𝑆𝑟𝑟 Assuming the value of 4 wt. % pyrite: 3.4 𝑤𝑤𝑐𝑐 % 𝐹𝐹𝑚𝑚𝑆𝑆2 = 1.8 𝑤𝑤𝑐𝑐 % 𝑆𝑆 𝐴𝐴𝐴𝐴 = 1.8 𝑤𝑤𝑐𝑐 % 𝑆𝑆 × 1000 𝑘𝑘𝑘𝑘100 × 100 𝑘𝑘/𝑟𝑟𝑟𝑟𝑆𝑆32 𝑘𝑘/𝑟𝑟𝑟𝑟𝑆𝑆 𝐴𝐴𝐴𝐴 = 1.8 𝑤𝑤𝑐𝑐. % × 31.25 𝑘𝑘𝑘𝑘 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑐𝑐𝑟𝑟𝑛𝑛𝑛𝑛𝑚𝑚−1 𝐴𝐴𝐴𝐴 = 56.3 𝑘𝑘𝑘𝑘 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑐𝑐𝑟𝑟𝑛𝑛𝑛𝑛𝑚𝑚−1 For calcite, an NP contribution of 1000 kg CaCO3 per tonne is assumed. 𝑁𝑁𝐴𝐴 = 100 × 𝑤𝑤𝑐𝑐. % 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑘𝑘𝑘𝑘 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑐𝑐𝑟𝑟𝑛𝑛𝑛𝑛𝑚𝑚−1 𝑁𝑁𝐴𝐴 = 100 × 0.9 𝑤𝑤𝑐𝑐. % 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑘𝑘𝑘𝑘 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑐𝑐𝑟𝑟𝑛𝑛𝑛𝑛𝑚𝑚−1 𝑁𝑁𝐴𝐴 = 9 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑘𝑘𝑘𝑘 𝐶𝐶𝑚𝑚𝐶𝐶𝑂𝑂3 𝑐𝑐𝑟𝑟𝑛𝑛𝑛𝑛𝑚𝑚−1 Therefore, 𝑁𝑁𝐴𝐴:𝐴𝐴𝐴𝐴 = 1: 6.3 A.4 Pile Structure In order to account for the reduction in hydraulic conductivity influenced by traffic surfaces, a scale factor was applied depending on the original 𝑑𝑑𝑚𝑚 value following a non-linear relationship 103 shown in Table A.1. Lower 𝑑𝑑𝑚𝑚 values were assigned a higher scale factor as they were less influenced by compaction due to already having smaller grain size and closer packing (and vice-versa for larger hydraulic conductivity values). This is applied to Equations 1 and 4 as seen in Equations A.4 and A.5, respectively. Resulting hydraulic conductivity values were comparable to those observed by Lorca et al. (2016). The calculated values for Ki were used to generate porosity using Equation 2. K𝑖𝑖 = α𝑡𝑡α𝐾𝐾K� �𝑑𝑑𝑚𝑚,𝑖𝑖𝑑𝑑m �2 (A. 4) 𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒,𝑖𝑖0 = α𝑡𝑡α𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒 𝑑𝑑𝑚𝑚𝑑𝑑𝑚𝑚,𝑖𝑖 (A. 5) Table A.1 Traffic surface scale factor based on hydraulic conductivity value Original 𝒅𝒅𝒎𝒎 Value Scale Factor (𝛂𝛂𝒕𝒕) Applied to K Calculation >5 mm 0.01 1.58 to 5 mm 0.05 0.5 to 1.58 mm 0.1 0.158 to 0.5 mm 0.5 A.5 Soil Hydraulic Function Parameters Variability in soil hydraulic function parameters in the models include residual water saturation (𝑆𝑆𝑤𝑤𝑤𝑤) and van Genuchten fitting parameters 𝛼𝛼 and 𝑛𝑛. 𝑆𝑆𝑤𝑤𝑤𝑤 is the residual saturation remaining in pore spaces when pores are drained and saturated hydraulic conductivity approaches 0. Parameters 𝛼𝛼 and 𝑛𝑛 are empirical parameters inversely related to air-entry tension and the pore size width, respectively (Kool et al. 1985). Below the results of two test cases with spatial variability in these parameters were compared to the original case where these values are assumed homogeneous. These parameters followed a log-linear relationship with hydraulic conductivity (K) to maintain average values used in the homogeneous cases (i.e. 𝑆𝑆𝑤𝑤𝑤𝑤 = 0.22, 𝛼𝛼 = 5.0 m-1, and 𝑛𝑛 = 2.7) (Lefebvre et al. 2001; Fala et al. 2003; Azam et al. 2007; Šimůnek and van 104 Genuchten 2008; Neuner et al. 2013), with distributions shown in Figure A.5. 𝑆𝑆𝑤𝑤𝑤𝑤 and K are considered inversely proportional since smaller grains with lower K are associated with a lower porosity and smaller pore throat diameters and retain more water. Here 𝛼𝛼 and 𝑛𝑛 are considered directly proportional to K. Results in Figure A.6. demonstrate that variable spatial distributions of these parameters did not have a significant effect on the overall basal drainage quality, and these parameters were therefore kept homogeneous for subsequent simulations. 105 Figure A.5. Nodal distributions for soil hydraulic function parameters. This range follows a normal distribution, where 𝑺𝑺𝒘𝒘𝒑𝒑 was assumed to follow an inverse log-linear relationship to 𝐊𝐊, while van Genuchten parameters 𝜶𝜶 and 𝒏𝒏 followed a direct log-linear relationship to 𝐊𝐊. 106 Figure A.6. Mass loadings from homogeneous vs. heterogeneous soil hydraulic function parameters (distributions given in Figure A.5). Results from these investigations suggest soil hydraulic function parameters did not have a significant influence on modeled pile total basal outfluxes in this study. A.6 Total vs. Local Drainage Quality In Figure A.5, five simulations were considered for each scenario investigated in this sensitivity analysis (represented through overlapping curves of the same colour). Although the acid generating pockets vary in locality within these simulations, total basal mass release rates follow a very similar trend. While the spatial locations of these heterogeneities may not largely influence overall mass loading rates from basal drainage, it does create significant variations in 107 local effluent chemistry. Figures A.7. – A.9. demonstrate this by showing effluent solute concentrations at every 1 m along the base of the pile. In the homogeneous case, drainage quality is consistent in all locations beneath the pile core. Along the slopes, decreases in mass loading rates occur towards the toe, since each parcel of water entering the pile has less contact time with weathering pyrite before exiting the base. However, in the heterogeneous cases, we see considerable variation in solute release rates and pH of drainage quality throughout the pile. 108 Figure A.7. Variations in drainage quality due to location under pile. Individual basal drainage locations have been plotted for one simulation of each case considering all heterogeneities investigated. Red represents basal drainage between the 0 – 55 m; purple is between 55 – 65 m, blue is between 65 – 75 m, and green is between 75 – 86 m. 109 Figure A.8. Calcite volume fractions along pile base through time. These cases correspond to the model configurations shown in Figure 3 of the main paper. In the homogeneous model, calcite is expended at the same rate throughout the pile. For push- and end-dumping models however, though initial concentrations and weathering rates for calcite are homogeneous, depletion occurs at different rates throughout the pile due to spatial variations in hydraulic conductivity, pyrite volume fraction and weathering rate. 110 Figure A.9. Pyrite volume fraction along pile base through time. Pyrite volume fractions at points along the pile base are shown here, corresponding to the model configurations shown in Figure 3 of the main paper. 111 Appendix B Supplementary Material for Chapter 3 B.1 System Overview and Construction Figure B.1 Photos of laboratory model. Top left: Front view of the laboratory scale model. HI rock is the dark grey material on the left which was deposited in three segregated zones; anorthosite is the light grey material (to the right of the HI); sand is orange-brown in colour along the top with crushed anorthosite (light grey) above to form the FCL. The tank was pivoted at an angle of 5˚ for the duration of the experiments with steady-state flow rates. Top right: Drainage tubes 1 cm in diameter that allow effluent to freely drain into collection buckets below. Bottom: Set up of tubing with drainage holes installed along the top of the lab model to distribute water applied at 180 L day-1. Table B.1. Tube Length at Reservoir Ports Reservoir Tube Length [cm] 1 72 2 93 3 120 4 96 5 128 6 52 7 85 8 70 112 B.2 Steady-State Flow Experiment For the steady-state flow experiment, a constant rate of 180 L day-1 of tap water was applied to the model. However, an average outflow rate of 158 L day-1 was measured from collection buckets in reservoirs 1 through 8. In this case, 22 L day-1 is unaccounted for. The following equation was used to obtain a rough estimate of evaporation rates from an open tank. gh = θA(xs − x) (Engineering ToolBox 2004) Where, 𝑘𝑘ℎ = amount of evaporated water (kg/hr) θ = (25 + 19v) evaporation coefficient (kg/m2 hr) A = area xs = max humidity ratio of saturated air (kg/kg) (temperature dependent) x = humidity ratio of air (kg/kg) Values of xs were obtained from the Mollier plot, below: 113 Figure B.2. Mollier diagram for obtaining maximum humidity ratio. Obtained from Engineering ToolBox 2003. 114 Table B.2. Estimated evaporation rate range from tank model Temperature RH (%) x (kg/kg) (Mollier) xs (kg/kg) Evap co. (kg m-2 h-1) Area (m2) gh [L day-1] Average 22.5 29.8 0.0049 0.0150 25 1.482 8.98 Min 21.7 15.0 0.0024 0.0155 11.65 Max 23.1 42.9 0.0073 0.0160 7.74 From these calculations, 35 – 53% attributed to evaporation—however in reality it is likely lower. Difficulties were encountered while attempting to maintain a constant infiltration rate from the sprinkler with the pump available for this experiment (Plante, personal communication, 2019). While the flow rates at each port were fairly consistent, it is likely that a lower flow rate was actually applied to the laboratory model from the sprinkler system shown in Figure B.1, therefore discrepancies in volumetric mass balance are likely significantly attributed to this. 115 B.3 Grain Size Distribution Curves Figure B.3. Grain size distribution of hemo-ilmenite section 1 01020304050607080901001 10 100 1000 10000 100000% PassingGrain Size Diameter (µm)HI S1 116 Figure B.4. Grain size distribution of hemo-ilmenite section 2 01020304050607080901001 10 100 1000 10000 100000% PassingGrain Size Diameter (µm)HI S2 117 Figure B.5. Grain size distribution of hemo-ilmenite section 3 01020304050607080901001 10 100 1000 10000 100000% PassingGrain Size Diameter (µm)HI S3 118 Figure B.6. Grain size distribution of anorthosite material 01020304050607080901001 10 100 1000 10000 100000% PassingGrain Size Diameter (µm)Anorthosite 119 Figure B.7. Grain size distribution of sand material (FCL) 01020304050607080901001 10 100 1000 10000% PassingGrain Size Diameter (µm)Sand 120 Figure B.8. Grain size distribution of crushed anorthosite material (FCL) 01020304050607080901001 10 100 1000 10000 100000% PassingGrain Size Diameter (µm)Crushed Anorthosite 121 B.4 Numerical Model Inputs Figure B.9. Triangular mesh with zone outlines used for numerical model. The mesh file used in numerical simulations consisting of 16602 cells at a sloped at a 5 ̊ angle. Each cell has an average area and representative elemental volume of 1.4 cm3. This unstructured grid allows for the influence of a 5° FCL slope to be investigated, as well as the effects of internal sloped boundaries following the angle of repose. 122 Table B.3. Aqueous complexation and hydrolysis reactions used in model simulations Aqueous Complexation Reaction 𝐥𝐥𝐥𝐥𝐥𝐥𝟏𝟏𝟏𝟏𝐊𝐊𝟐𝟐𝟐𝟐 OH− + H+ ↔ H2O (𝑆𝑆) -13.9980 O2 (𝑘𝑘) ↔ O2 (𝑚𝑚𝑎𝑎) 2.8980 HSO4− ↔ H+ + SO42− 1.9870 H2S (𝑚𝑚𝑎𝑎) ↔ H+ + HS− 6.9940 S2− + H+ ↔ HS− -12.9180 FeCO3+ ↔ Fe3+ + CO32− 9.7200 FeOH+ + H+ ↔ Fe2+ + H2O (𝑆𝑆) 13.2000 FeOH2+ + H+ ↔ Fe3+ + H2O (𝑆𝑆) -2.1900 FeOH2+ + 2H+ ↔ Fe3+ + 2H2O (𝑆𝑆) -5.6700 FeOH3 (𝑚𝑚𝑎𝑎) + 3H+ ↔ Fe3+ + 3H2O (𝑆𝑆) -12.5600 FeSO4+ ↔ Fe3+ + SO42− 4.0400 HCO3− ↔ H+ + CO32− 10.3300 H2CO3 (𝑚𝑚𝑎𝑎) ↔ 2H+ + CO32− 16.6810
Thesis/Dissertation
eng
Geological Sciences
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Applications of reactive transport modeling for assessing the effects of heterogeneities and internal structure on mass loadings from waste rock piles
Text