ON THE ROLE OF MULTICOMPONENT DIFFUSION AND ELECTROCHEMICAL MIGRATION FOR REACTIVE TRANSPORT IN POROUS MEDIA by Pejman Rasouli M.Sc., The Eberhard Karls University of Tübingen, 2008 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Geological Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2016 © Pejman Rasouli, 2016 ii Abstract In multicomponent solutions, electrostatic coupling between charged species leads to a process called “electromigration”. Neglecting electromigration results in a charge imbalance and an incomplete and unrealistic description of mass transfer. Although not commonly considered in reactive transport codes, electromigration can strongly affect mass transport processes and can explain unexpected behaviors such as uphill diffusion or isotope fractionation. Including electrostatic coupling in reactive transport codes enables simulation of problems involving mass transport by advection and diffusion, electromigration and geochemical reactions, such as electrokinetic remediation and the geobattery concept associated with buried ore bodies. There are generally two methods for coupling charge and mass continuities. The first method is based on the null-current approach which assumes negligible electric current transmission. The second method considers explicit coupling of mass and electric fluxes. In this study both methods are investigated and their implications for reactive transport are examined. To this end, MIN3P, a fully coupled 3D reactive transport code, was extended by integrating the Nernst-Planck and Gauss-Ampère equations. The implementation of the Nernst-Planck equations was verified by inter-comparison with other existing reactive transport codes based on a set of benchmark problems. At the same time, these benchmark problems illustrate the effect of electric coupling during multicomponent diffusion and electrochemical migration. By explicit coupling of the Nernst-Planck and Gauss-Ampère equations, MIN3P was further enhanced to simulate electrokinetic remediation and the resulting code was tested for desalination problems. In addition, scenario and sensitivity analysis were used to investigate the potential for spontaneous exsolution of gases in response to gas generation at the electrodes of electrokinetic remediation systems. iii Finally, a process-based model linking surface-measureable self-potential signals to electrochemical transport and geochemical reactions associated with buried metallic bodies was developed. The enhanced code provides a reactive transport modeling framework for process-based forward modeling of self-potential signals and associated geochemical signatures of buried ore bodies and allows a quantitative investigation of the “geobattery concept”. The code was tested based on published data from a laboratory experiment involving a buried iron bar and used to evaluate the geobattery concept based on an illustrative example of a buried ore body. iv Preface This thesis has been prepared as a collection of manuscripts either published or in preparation for publication. Although these manuscripts have been co-authored with individuals other than myself, I am the first author in each case and have conducted the research and manuscript preparation. The objectives of each chapter and research approach are based on my initiative in consultation with my thesis supervisor Prof. K. U. Mayer. The actual research work including model development, implementation and numerical simulations has been entirely done by myself. My supervisor provided guidance towards content and structure of the manuscripts. Only limited support was provided by co-authors other than my thesis supervisor, as outlined below. A version of chapter 2 has been published as: P. Rasouli, C. I. Steefel, K. U. Mayer, and M. Rolle, “Benchmarks for multicomponent diffusion and electrochemical migration,” Comput. Geosci., vol. 19, no. 3, pp. 523–533, 2015, doi:10.1007/s10596-015-9481-z. I conducted all verification simulations and wrote most of the manuscript. The “Benchmark 2: tracer isotope diffusion” problem was proposed by Steefel, C. I. and the section on “Benchmark 3: transverse dispersion” was originally drafted by Rolle, M. A version of chapter 3 will be submitted as: P. Rasouli and K. U. Mayer: Numerical simulation of gas generation and reactive solute transport during electrokinetic remediation. I carried out all model development, implementation and numerical simulations and analyzed all data. I wrote the manuscript guided by my thesis supervisor. A version of chapter 4 will be submitted as: P. Rasouli and K. U. Mayer: A process-based reactive transport model for understanding self-potential signals associated with buried ore v bodies. I carried out all model development, implementation and numerical simulations and analyzed all data. I wrote the manuscript guided by my thesis supervisor. The multicomponent diffusion model developed in chapter 2 was used to contribute to another benchmark study investigating the role of electromigration and its feedback of geochemical reactions. This manuscripts was published as: P. Alt-Epping, C. Tournassat, P. Rasouli, C. I. Steefel, K. U. Mayer, A. Jenni, U. Mäder, S. S. Sengor, and R. Fernández, “Benchmark reactive transport simulations of a column experiment in compacted bentonite with multispecies diffusion and explicit treatment of electrostatic effects,” Comput. Geosci., vol. 19, no. 3, pp. 535–550, 2015, doi: 10.1007/s10596-014-9451-x. vi Table of Contents Abstract .................................................................................................................................... ii Preface ...................................................................................................................................... iv Table of Contents ..................................................................................................................... vi List of Tables ............................................................................................................................ ix List of Figures........................................................................................................................... xi Acknowledgements ................................................................................................................ xvii Dedication .............................................................................................................................xviii Chapter 1: Introduction ............................................................................................................ 1 1.1 Background and Motivation .........................................................................................1 1.2 Previous Studies ..........................................................................................................7 1.3 Research Objectives .....................................................................................................9 1.4 Organization .............................................................................................................. 10 Chapter 2: Benchmarks for Multicomponent Diffusion and Electrochemical Migration ... 12 2.1 Introduction ............................................................................................................... 12 2.2 Governing Equations ................................................................................................. 15 2.2.1 Mass Transfer in Electrolytic Systems ................................................................... 15 2.2.2 Nernst-Planck Equation for Multicomponent Systems ............................................ 15 2.3 Participating Codes .................................................................................................... 16 2.4 Benchmark Descriptions ............................................................................................ 17 2.4.1 Benchmark 1: Transient Electromigration ............................................................. 18 2.4.2 Benchmark 2: Tracer Isotope Diffusion ................................................................. 19 2.4.3 Benchmark 3: Transverse Dispersion ..................................................................... 20 vii 2.5 Results and Discussion .............................................................................................. 24 2.5.1 Benchmark 1 .......................................................................................................... 24 2.5.2 Benchmark 2 .......................................................................................................... 25 2.5.3 Benchmark 3: 1D Transverse Dispersion .............................................................. 28 2.5.4 Benchmark 3: 2D Flow and Transverse Dispersion ............................................... 30 2.6 Concluding Remarks.................................................................................................. 33 Chapter 3: Numerical Simulation of Gas Generation and Reactive Solute Transport during Electrokinetic Remediation .................................................................................................... 34 3.1 Introduction ............................................................................................................... 34 3.2 Governing Equations ................................................................................................. 37 3.2.1 Generalized Mass Continuity in Electrolyte Solutions ............................................ 37 3.2.2 Nernst-Planck Equation for Multicomponent Systems ............................................ 39 3.2.3 Generalized Charge Continuity in Electrolyte ........................................................ 40 3.2.4 Geochemical Relationships .................................................................................... 43 3.2.5 Gas Exsolution and Degassing ............................................................................... 43 3.3 Solution Procedure .................................................................................................... 45 3.4 Model Verification and Validation ............................................................................. 45 3.4.1 Free Diffusion ........................................................................................................ 46 3.4.2 Electrokinetic Permanganate Transport .................................................................. 51 3.4.3 Electrokinetic Desalinization ................................................................................. 54 3.4.3.1 Electrode Processes and Boundary Conditions ............................................... 56 3.4.3.2 Simulation Results ......................................................................................... 57 3.5 Evaluation of the Potential for Gas Generation, Exsolution and Consumption ............ 66 viii 3.5.1 Simulation Results ................................................................................................. 69 3.5.2 Sensitivity and Scenario Analysis .......................................................................... 73 3.6 Summary and Conclusions ......................................................................................... 76 Chapter 4: A Process-based Reactive Transport Model for Understanding Self-potential Signals Associated with Buried Ore Bodies ........................................................................... 79 4.1 Introduction ............................................................................................................... 79 4.2 Conceptual Model...................................................................................................... 82 4.3 Governing Equations ................................................................................................. 84 4.4 Model Implementation and Verification ..................................................................... 87 4.5 Experimental Case Study ........................................................................................... 88 4.6 Illustrative Ore Deposit Model ................................................................................. 101 4.7 Summary and Conclusions ....................................................................................... 115 Chapter 5: Summary and Conclusions ................................................................................ 118 5.1 Model Development Summary ................................................................................ 118 5.2 Limitations, Future Works and Contributions........................................................... 122 Bibliography .......................................................................................................................... 125 Appendix A ............................................................................................................................ 146 ix List of Tables Table 2.1. Summary of multicomponent diffusion benchmarks. ................................................ 18 Table 2.2. Boundary conditions, initial conditions, and species dependent diffusion coefficients for Benchmark 1 (transient electromigration problem)............................................................... 19 Table 2.3. Boundary conditions and diffusion coefficients for Benchmark 2 (isotope tracer problem). .................................................................................................................................. 20 Table 2.4. Chemical conditions and transverse dispersion coefficients for Benchmark 3 (transverse dispersion problem). ................................................................................................ 22 Table 3.1. Chemical parameters: boundary conditions and diffusion coefficients for free diffusion problem. ..................................................................................................................... 47 Table 3.2. Physical parameters for free diffusion problem. ........................................................ 48 Table 3.3. Boundary condition, initial condition, and model parameters used for permanganate removal problem. ...................................................................................................................... 52 Table 3.4. Boundary conditions, initial conditions, and species-specific diffusion coefficients for desalination problem. ................................................................................................................ 56 Table 3.5. Physico-chemical model parameters used in desalination problem. ........................... 68 Table 3.6. Reactions network and corresponding coefficients. ................................................... 71 Table 3.7. Dimensionless and composite scaled sensitivities for Scenario 6 calculated by UCODE [141]. .......................................................................................................................... 74 Table 4.1. Reactions network and corresponding coefficients [166]. .......................................... 90 Table 4.2. Boundary conditions, initial conditions, and species-specific diffusion coefficients... 93 Table 4.3. Physical parameters. ................................................................................................. 93 Table 4.4. Initial conditions for three different zones. .............................................................. 103 x Table 4.5. Physical parameters. ............................................................................................... 103 Table 4.6. Reactions network and corresponding coefficients [166]. ........................................ 104 xi List of Figures Figure 2.1. Schematic of the 2D flow and transverse dispersion experiment conducted by Rolle et al. [70]. ..................................................................................................................................... 21 Figure 2.2: Species concentrations after 1 hour simulation time for HNO3 diffusion (Benchmark 1). The left boundary is a fixed concentration (Dirichlet) boundary, while the right boundary is no-flux. ..................................................................................................................................... 25 Figure 2.3: Na+, Cl-, H+, OH- and 22Na+ concentrations after 1500 days for system summarized in Table 2 (Benchmark 2). The left boundary is a fixed concentration (Dirichlet) boundary at 0.5 mM, while the right boundary is a fixed concentration boundary at 0.1 mM for Na+, Cl-. The fixed gradient in NaCl results in a flux of H+, OH- and 22Na+, despite the fact that their concentrations are the same at either end of the column. ............................................................ 26 Figure 2.4: Na+, Cl-, H+, OH- and 22Na+ diffusive, electromigration and total fluxes after 1500 days (representing conditions close to steady state) for the system summarized in Table 2 (Benchmark 2). H+, OH- and 22Na+ diffusive fluxes change direction at 0.65 cm. ....................... 27 Figure 2.5. 1D simulation results of transverse profiles for Cl-, K+ and Mg+2 at the outlet (corresponding to a residence time of 16 hours) demonstrate the effect of species-dependent dispersion and electromigration on the transverse displacement of charged species (1D Benchmark 3 solved with MIN3P). ........................................................................................... 29 Figure 2.6. 1D simulation of transverse multicomponent diffusion for the case of transport of mixed electrolytes (KCl and MgCl2 solution) in pure water described by Rolle et al. [70], comparing CrunchFlow, MIN3P and PHREEQC results. .......................................................... 30 xii Figure 2.7. Simulation results for Benchmark 3 considering flow (uni-directional) and multicomponent transverse dispersion for steady-state conditions, from top to bottom are shown: K+, Mg+2 and Cl- for CrunchFlow and K+, Mg+2 and Cl- for MIN3P. ......................................... 31 Figure 2.8. Comparison of 1D PHREEQC results (no explicit consideration of flow, only following the plume as it moves down the flow path) and transverse profiles derived from 2D CrunchFlow and MIN3P runs for the transverse dispersion problem. The CrunchFlow 2D runs are based on GIMRT and use a first order upwind formulation, along with a backwards Euler time stepping approach, the same numerical methods are used in the MIN3P simulations. ........ 32 Figure 3.1. MIN3P Gauss-Ampère formulation (symbols) and Paz-García Poisson-Nernst-Planck formulation (solid line) [20]: (a) Na+ concentration profile; (b) pH profile and (c) diffusion potential profile (initial conditions for potential is zero for entire domain). ................................ 49 Figure 3.2. Schematic of the 1D pH-isolated electrokinetic experiment conducted by Hodges et al. [129]. Electrode reservoirs are separated by clay cores to prevent pH disturbance of the electrokinetic cell. ..................................................................................................................... 51 Figure 3.3. MIN3P (solid line), measured (symbols) and PHT3D (dashed line) [39]: MnO4_concentration in target reservoir over time. ..................................................................... 53 Figure 3.4. Schematic of the 1D representation of electrokinetic desalination experiment conducted by Kamran et al. [61]. ............................................................................................... 55 Figure 3.5. (a) Na+ concentration profiles for selected times (0, 2, 4, 6, 8, 10, 12, 14 and 16 hours): model (solid sharp lines) vs. experimental measurements by NMR (shadow background lines); the arrows show concentration profile evolution over time; (b) potential profile for 0, 10, 12 and 16 hours: model (lines) vs. experimental measurements (symbols); (c) Na+ and Cl− concentration profile at t = 16 hours (end of experiment): model (solid and dashed sharp lines), xiii experimental measurements of Na+ concentration by NMR (shadow background lines) and experimentally measured Cl− concentrations by ion chromatography (symbols); (d) acid and alkaline fronts progress (pH = 1 and pH = 13 respectively): model (lines) vs. measurements (symbols). Simulation results are for Scenario 1 and 2 and experimental data are taken from Kamran et al. [61]. .................................................................................................................... 59 Figure 3.6. Gas evolution over time at anode (solid line) and cathode (dashed line): (a and c) base case (Scenario 1) including gas generation, but without degassing; (b and d) including gas generation and degassing (Scenario 2). Cl2(g) is plotted separately (c and d) as its contribution to total pressure was minimal with respect to other gas components. Degassing starts from the beginning and maintains total gas pressure by not allowing the sum of all gas pressures to exceed the confining pressure. Degassing strips out N2(g) from the mixture to maintain pressure (cf. a and b). ....................................................................................................................................... 63 Figure 3.7. Simulated total and partial pressures and dissolved gas profiles during electrokinetic desalination for selected times (0, 2, 4, 6, 8, 10, 12, 14 and 16 hours): (a and c) base case (Scenario 1) including gas generation, but without degassing; (b and d) including gas generation and degassing (Scenario 2). Degassing starts after 13 hours and 20 minutes, maintaining total pressure by not allowing the sum of all gas pressures to exceed the confining pressure. Degassing strips out N2(g) from the mixture to maintain total gas pressure at the confining pressure. ......... 65 Figure 3.8. Gas evolution over time at anode (solid lines) and cathode (dashed lines): (a) degassing and gas consumption not active (Scenario 3); (b) degassing active, but no gas consumption (Scenario 4); (c) degassing not active, but including gas consumption (Scenario 5); (d) degassing and gas consumption active (Scenario 6). Degassing maintains total pressure by not allowing the sum of all gas pressures to exceed the confining pressure (b and d). Gas xiv consumption reaction first depletes oxygen at the cathode then lets the pressure increase (c and d). ............................................................................................................................................. 72 Figure 3.9. Total gas pressure at anode (solid lines) and cathode (dashed lines) at the end of simulation and as a function of: (a) electric potential; (b) rate of hydrogen oxidation and (c) rate of degassing. (Note: confining pressure is 1.0097 [atm]). .......................................................... 76 Figure 4.1. Geobattery concept. An electric field is established when a gradient of the redox potential is connected by an electrical conductive body. Oxidation-reduction reactions (oxygen reduction near the surface and metal oxidation e.g. metal or sulfide oxidation at depth) are responsible for the generation of an electric current by electron flux through the ore body. The generated electric field initiates electromigration of ions in the opposite direction in the surrounding aquifer to maintain electroneutrality. Geochemical reactions such as mineral dissolution-precipitation, gas generation and exsolution may occur in the surrounding environment and within the ore body. ........................................................................................ 83 Figure 4.2. (a) Domain and boundary conditions and (b) MIN3P model grid applied for solving iron bar case study..................................................................................................................... 91 Figure 4.3. Distribution of the redox potential: simulated (top) vs. measured (bottom) [12] 42 days after inserting the iron bar on the left hand side of the tank. Measured redox potential from [12]. Copyright © 2008 European Association of Geoscientists & Engineers. Reprinted with permission from John Wiley & Sons. ........................................................................................ 95 Figure 4.4. Self-potential distribution on the surface; line is the model and symbols are measured at the end of the experiment after 42 days. ................................................................................. 96 xv Figure 4.5. Self-potential distribution after 42 days: simulated (top) vs. measured (bottom) [12]. Measured self-potential from [12]. Copyright © 2008 European Association of Geoscientists & Engineers. Reprinted with permission from John Wiley & Sons. ............................................... 97 Figure 4.6. Oxygen partial pressure (atm) after 42 days: depletion of oxygen due to oxygen reduction reactions. ................................................................................................................... 98 Figure 4.7. pH distribution: simulated (top) vs. measured (bottom) [12] after 42 days after inserting the iron bar on the left hand side of the tank. Measured pH from [12]. Copyright © 2008 European Association of Geoscientists & Engineers. Reprinted with permission from John Wiley & Sons. ........................................................................................................................... 99 Figure 4.8. Secondary mineral formation in the direct vicinity of the iron bar after 42 days, shown as volume fractions. ..................................................................................................... 100 Figure 4.9. Initial mineral distribution a) pyrite, b) calcite, and c) K-mica. .............................. 105 Figure 4.10. Redox potential at quasi steady state. Electron drift through the ore body decreases the redox potential gradually. .................................................................................................. 107 Figure 4.11. Oxygen distribution at quasi steady state. ............................................................ 108 Figure 4.12. pH distribution at quasi steady state. pH decreases at depth due to pyrite oxidation and increases near the surface due to oxygen reduction. .......................................................... 109 Figure 4.13. (a) Pyrite dissolution rate, (b) K-mica dissolution rate, (c) pyrite saturation index and (d) calcite dissolution rate at the end of simulation (quasi steady state). ............................ 109 Figure 4.14. Self-potential generated by pyrite oxidation at quasi steady state. ........................ 110 Figure 4.15. Self-potential response for (a) base case scenario, (b) reduced oxygen influx by a factor of 3.5 by reducing upper boundary redox potential and (c) surface signal. ..................... 112 xvi Figure 4.16. Generation and transport of ions: (a) ferrous iron, (b) sulfate, (c) potassium, (d) carbonate, (e) aluminum, (f) silicic acid, (g) redox potential and (h) chloride. ......................... 114 Figure 4.17. Precipitation of secondary minerals: (a) siderite precipitation rate and (b) ferrihydrite precipitation rate. .................................................................................................. 114 xvii Acknowledgements First and foremost I thank Prof. Uli Mayer for his continuous support, leadership, constructive guidance and providing coherent answers to my endless questions. I also would like to express my sincere gratitude for his support during some very difficult times that I went through and his enormous patience in fixing my misuse of definite and indefinite articles in my writings. I would like to extend my sincere appreciation and gratitude to my committee members Prof. Roger Beckie and Prof. Mark Jellinek for their advice during my time at UBC. I would like to thank my wife, Maryam, for her everlasting and unconditional love, support and all sacrifices that she made during all these years. I know that a lot of this work would not have been possible without her love and patience. Funding for this research was provided by the Natural Sciences and Engineering Research Council of Canada (NSERC) in form of a Discovery Grant and a Discovery Accelerator Supplement awarded to K. Ulrich Mayer. xviii Dedication To Maryam 1Chapter 1: Introduction 1.1 Background and Motivation Reactive transport codes can be used to better understand flow, transport and reactions in aquifer systems. These codes help to verify conceptual models, design experiments and increase the knowledge about important parameters and processes in the subsurface. A multicomponent reactive transport model in saturated porous media can be used to address many problems in the field of geology, hydrogeology, engineering and environmental research. The migration of aqueous species is affected by a variety of transport mechanisms. Conventionally, the three common transport mechanisms considered in hydrogeology are advection, dispersion and diffusion. Advection is solute transport by groundwater flow; its direction and magnitude is related to the groundwater velocity. Dispersion is due to mixing and spreading of the solute created by variation in the magnitude and direction of the pore velocity. Molecular diffusion or diffusion is a pore-scale mixing by random molecular movements of solutes driven by Brownian motion. In many transport models, diffusion coefficients appear as proportionality constants in Fick’s law and are treated as adjustable parameters [1]. In order to calculate a specific value of flux and concentration distribution, experimental measurements of these coefficients are required. The common basis for estimating diffusion coefficients in liquids is the Stokes-Einstein equation. The Stokes-Einstein model suggests direct dependency of diffusion coefficients to temperature variations and inverse dependency to solvent viscosity and solute ions radii. The effects of solvent viscosity and solute radius are much larger than the temperature effect [2]. 2 Although these three transport processes suffice in many cases to adequately describe migration behavior, there are systems in which significant interactions between diffusing molecules strongly affect the apparent diffusion coefficients. In some cases, these interactions produce unexpected diffusion behavior such as uphill diffusion [3], [4]. For specific solutes sometimes diffusion coefficients show a strong dependency on concentration [2]. As the measurement of diffusion coefficients is a very difficult task [1] and it is almost impossible to generate apparent diffusion coefficients for different species combinations for varying conditions, it is useful to consider electrochemical interactions affecting diffusion explicitly rather than hiding them as part of the empirically measured diffusion coefficient. Species-specific diffusion is often necessary to describe this behavior, if diffusive transport is the dominant mass transport process. Ions have different charges and mobilities, by considering only concentration gradient a charge imbalance will develop in the electrolyte, which is not physically possible. Therefore, additional transport phenomena are required to counteract charge imbalance. The electric interactions of the ions in an electrolyte affect ion mobility and contribute to a transport phenomenon called “electromigration”. In a general sense, this transport mechanism ties the migration of big cations and small anions together electrostatically [1], [5] to enforce electrical neutrality. Fick’s law does not consider electrostatic bonding (columbic interactions) of charged species. To account for these electric interactions, an explicit consideration of charge and electric potential field is required. This coupling can be represented by the Nernst-Planck equation. This formulation ensures that the electrical field that drives electromigration 3counteracts the charge imbalance. As a result, multicomponent diffusion is controlled by gradients of chemical and electrical potentials. To adequately describe multicomponent transport in diffusion-dominated systems, a formulation is required that includes species-dependent diffusion coefficients and an electrochemical migration term. According to Ohm’s law, movement or migration of charges induce electric current. Based on the way electric current is conducted and the value of the electric conductivity, media and materials are classified as conductors, semiconductors or insulators (or dielectrics). Water is not an electric conductor in itself and serves as dielectric but the presence of charged species (ions) in the solution makes it electrically conductive [6]. Electric charge is defined physically by carriers such as electrons and ions. The conductors are characterized by the nature and concentration of the free charges. Depending on the involved charges an electric current can be electric or ionic. Electrical conduction establishes where there is an electric potential gradient and an electrical conductor is present. In the geochemical systems, this gradient can be the redox potential gradient. Redox potential is measured with respect to standard hydrogen electrode as a reference point and it indicates the tendency for reduced or oxidized species to dominate [7]. On the other hand, self-potential is the electric potential generated by movements of charges (ions and electrons) [8] and measured by non-polarizing electrodes, often placed at the ground surface [9]. Ore bodies have been shown to act as electrical conductors [10]–[12]. Metals and certain other substances like carbon materials (graphite and carbon black), some oxides and inorganic compounds (e.g. sulfide ore deposits and tungsten carbide) and a number of organic substances have also been identified to act as electrical conductors [6]. Ionic conductors (the second type 4conductors) are known as electrolytes. Aqueous solutions of acids, bases and salts are ionic conductors. Certain hydrogeochemical settings can favor conditions in which a mixed electric current establishes consisting of electric and ionic components [6]. Most reactive transport models either neglect explicit electrostatic coupling by using Fick’s law [13], use a simplified form of the Nernst-Planck equation by assuming no externally induced current (null current assumption) [14], [15], or only consider limited and simplified geochemical suites for the reactions [16]. Neglecting the electric coupling in a reactive transport model also limits the ability of the model to describe transport and reactions of chemical components exposed to an external electric current as applied in electrokinetic remediation techniques. Surface measurable self-potential associated with buried metallic ore bodies is another phenomenon that is associated with electromigration and electron transfer in the subsurface. Self-potential measurements have been used for decades to locate and characterize ore deposits without quantitative analysis of geochemical reactions that contribute to generation of electric potential. Therefore, a new development of a multicomponent electrochemical reactive transport model is necessary to quantitatively analyze these systems. Extending a multicomponent reactive transport model to an electrochemical multicomponent reactive transport model requires departure from Fick’s law by including the Nernst-Planck equation. Fick’s law is based on the chemical potential of aqueous species while the Nernst-Planck approach includes also the effect of electrochemical potential by considering electrostatic interactions of species with each other [5], [6], [17], [18]. The movement of ions generates an ionic electric field and consequently electric current in the system. An external electric field can also contribute to the inner electric 5field; as for the charged ions there is no way of understanding the origin of the electric field without considering coupled electromigration processes [5], [17]. As a result, the same formulation can be used for an applied electric field condition applicable for electrokinetic remediation technology. Application of external electric potential can also lead to hydrolysis of water on the electrodes. Hydrolysis of water on electrodes forms O2(g) and H2(g) gas bubble on the electrode surfaces, decreases their electrical conductivity and consequently reduces remediation efficiency [19]–[23]. Therefore, a thorough understanding of geochemical processes including the fate of dissolved species, as well as gases is required to assess the performance of electrokinetic remediation on a process-based level. The gas generation and exsolution processes at the electrodes have not been included in electrokinetic simulations to date. Although the full Nernst-Planck equation is the most accurate approach for describing electrochemical migration processes, a variety of simplified approaches have been reported in the literature, e.g. for multicomponent diffusion in sediment diagenesis [2], [24]–[26]. The full Nernst-Planck equation is based on the gradient of the electrochemical potential that incorporates activity coefficient gradients and an electromigration term [26]–[32]. A common simplification valid for low concentrations and under isothermal conditions is neglecting the gradient of activity coefficients, which is negligible under these conditions [28], [31]. Electroneutrality is an essential condition in multicomponent systems. Many models assume that coupling of the Nernst-Planck equation with electroneutrality condition or constant electric field gives an adequate description of the ionic diffusion. This is a valid assumption in many practical cases [33], [34]; however, coupling the 6Nernst-Planck equation with an electric flux continuity equation is a more rigorous approach to the problem [17], [33], [35], [36]. Coupling mass flux continuity i.e. the Nernst-Planck equation and electric flux continuity i.e. Gauss-Ampère’s law provides a formulation applicable to a wide range of applications. In recent years efforts have been focused on the inclusion of electrokinetic simulation capabilities into reactive transport codes [37]–[40]. Integrating electrokinetic processes into the reactive transport framework utilizes the flexibility provided by the multicomponent reactive transport approach [28], [41] for electrokinetic remediation technologies. One of the motivations for the present contribution is to examine this potential with a specific emphasis on gas evolution in remediation systems. Another motivation for the present research is to study the geobattery concept from a reactive transport point of view. According to present geobattery models [10], [11], [42], [43], infiltrating oxygenated rain water creates a redox gradient between shallow and deep parts of the aquifer; when there is an ore or metallic body that connects these two regions, electron transport though the metallic body may generate an electric current. To maintain electroneutrality in the early stages, an ionic current in opposite direction must be generated. Even if the magnitude of the ionic current declines over time, it will leave behind a geochemical signature when steady state conditions are approached. The movement of ions and electrons in the system contributes to the total electric current flux that generates the electric field which manifests itself as self-potential signals measurable at the ground surface. Although the association of self-potential signals with mineral deposits [44], [45] contaminant plumes [46], [47], redox fronts [48], [49] and corrosion of buried metallic bodies [12], [50] has been investigated, 7so far there is no forward model that includes the wide range of (bio)geochemical (kinetic and equilibrium) reactions, taking into consideration porous media characteristics and geochemical composition along with the various transport phenomena that ultimately control and determine electric currents in geochemical systems and associated potential signals. 1.2 Previous Studies Multicomponent diffusion has been simulated by simplified forms of the Nernst-Planck equation [15], [30], [51] neglecting explicit electric coupling. There are also numerical solutions for coupled mass flux and electric flux system of equations. James et al. [52] developed a numerical solution for steady-state conditions using the finite element method and Galerkin residual weighting to solve the Poisson-Nernst-Planck system of equations where the Poisson equation was used to describe spatial distribution of charges. The model simulates the flow of liquids containing charged particles in a cylinder. Kato [53] also developed a similar solution as James et al. with the difference that its starting point is a potential profile which is used subsequently to calculate concentration distributions. The James et al. [52] and Kato [53] models are lacking geochemical reaction calculation capabilities. Samson et al. [54] developed a 1D model for a variable number of species for transient conditions. The model was further developed by incorporating heterogeneous reactions in variably saturated porous media. Their model was used to simulate the degradation of cement-based materials in aggressive environments [16], [55]. Paz-García et al. [19] used the finite element method to solve the Poisson-Nernst-Planck system of equations with a chemical equilibrium model to simulate clay brick desalination. To simulate electrokinetic in-situ oxidation 8[40], Wu et al. [39] incorporated an electric potential distribution module into an operator-splitting reactive transport code. However, the use of numerical models to investigate geochemical reactions associated with electrokinetic remediation remains limited to problem-specific cases for defined geochemical reaction networks [38], [56]–[58] and the process of gas generation and exsolution at the electrodes has not been studied. Sato and Mooney [42] proposed an electrochemical model describing the geobattery concept for the first time. Its mathematical formulation was later presented by Stoll et al. [10] and Bigalke and Grabner [11]; these authors were able to reproduce surface measurable self-potential signals at a deep continental drilling site in Germany [10], [11]. Stoll et al’s. [10] model is based on electrochemical potential equilibrium between the ore body and an electrolyte solution and relies on a depth-dependent exponential approximation of redox potential in the subsurface without considering the geochemical properties of the environment. Bigalke and Grabner [11] conducted laboratory experiments with synthetic materials to calculate the redox potential field and electrochemical reactions. None of these models incorporates geochemical reactions into a rigorous description of the feedback between geochemical reactions and self-potential signals. There is a lack of studies available in the literature for forward modeling of self-potential signals associated with reactive minerals. The motivation of this research is to provide a multicomponent electrochemical model to evaluate quantitatively the role of electric coupling in multicomponent diffusion and electromigration in a reactive transport framework. The contribution is to provide a 9tool that will help to better understand multicomponent diffusion and electromigration processes, electrokinetic remediation and the geobattery model. 1.3 Research Objectives The objectives of this research project are: 1) to develop an integrated multicomponent electrochemical reactive transport model for porous media that includes electromigration transport based on the Nernst-Planck equation and null current assumption in the aqueous phase with a suite of geochemical reactions; 2) to introduce benchmark problems to evaluate the effect of electric coupling during multicomponent diffusion and electrochemical migration and to facilitate an inter-comparison of existing reactive transport codes; 3) to develop and benchmark the option to couple the Nernst-Planck equation with Gauss-Ampère’s law in the presence or absence of external electric current; 4) to evaluate the potential for gas generation, exsolution and consumption during electrokinetic remediation; and 5) to develop a process-based reactive transport model to describe the geobattery concept with a focus on applications in mineral exploration. This thesis documents the development of the first numerical model for simulation of coupled reactive transport and electrical coupling for a wide range of applications. The identification of the degree and nature of coupling between the various transport and reaction mechanisms is important for the investigation of systems dominated by diffusive transport, the assessment of electrochemical remediation technologies and the quantitative evaluation of the geobattery concept. 101.4 Organization The main body of the dissertation consists of three chapters (chapters 2-4). The model development is based on the existing reactive transport code MIN3P [59] and focuses on the derivation, implementation and presentation of the equations accounting for the contributions of multicomponent diffusion and electromigration, full coupling of electric flux and mass flux, application of external electric current and electrode boundary conditions and electron transfer in electrical conductors. Chapter 2 contains the theoretical background of multicomponent diffusion and electromigration based on the Nernst-Planck equation and null current assumption. Three benchmark problems are introduced to evaluate the effect of electromigration on ion transport and also to facilitate the inter-comparison of existing reactive transport codes. Benchmark problems are essential tools for future developments of electromigration and multicomponent diffusion models. Chapter 3 presents theoretical background, implementation and verification of the Gauss-Ampère and Nernst-Planck coupling. This chapter explores circumstances in which explicit electric coupling is necessary. We investigate the application of the extended model for electrokinetic permanganate transport [60] and electrokinetic desalination [61] experiments. The model application is further extended by including, gas exsolution, degassing, and gas consumption kinetics to assess the effect of an evolving gas phase on the system. Chapter 4 focuses on incorporation of electron transfer in an electrical conductor and its application for simulating geobattery systems. A validation example and an illustrative application example for the geobattery concept are described in order to 11demonstrate the problems that can be simulated with this tool. For this purpose, experimental data from a sandbox experiment by Castermant et al. [12] is used as case study and for model validation. An ore deposit illustrative model is presented to demonstrate the applicability of this tool to simulate the geobattery concept associated with mineral deposits and to study feedback processes between ionic transport, electron transport and geochemical reactions on generation of self-potential signals. Appendix A contains derivation steps of formulations presented in the main chapters with greater detail. 12Chapter 2: Benchmarks for Multicomponent Diffusion and Electrochemical Migration1 2.1 Introduction It is well known that diffusive transport in multicomponent electrolyte systems cannot be fully described by Fickian diffusion alone, but is affected by a variety of processes including the electrostatic interactions between individual ions [1], [5], [62], [63]. Each dissolved species is subject to its own species-dependent diffusion coefficient, affected by parameters such as charge and size of the ion [1] and ionic conductivity [2]. As a result, dissolved species will tend to diffuse at different rates, promoting the development of a charge imbalance in solution. However, positively and negatively charged species are also affected by electric coupling, which ensures that charge balance in solution is maintained. Generally speaking, “large” cations and “small” anions are tied together electrostatically [1], [5] to enforce electroneutrality at the macroscale - an essential condition in electrolyte solutions [32], [64]. This electric coupling leads to an additional mass transport process called “electrochemical migration” or “electromigration” [5], [25]. Fick’s law neglects these interactions, describes ion migration solely based on concentration gradients, and consequently does not consider the electric field generated by electrostatic bonding (Coulombic interactions) of charged species [2], [5], [32], [64], [65]. In a multicomponent system that includes charged species, diffusive ion migration is therefore better described by the Nernst-Planck equation, a formulation that explicitly 1 A version of this chapter has been published as: P. Rasouli, C. I. Steefel, K. U. Mayer, and M. Rolle, “Benchmarks for multicomponent diffusion and electrochemical migration,” Comput. Geosci., vol. 19, no. 3, pp. 523–533, 2015, doi:10.1007/s10596-015-9481-z. 13considers the electric coupling between species and ensures the conservation of charge [2], [5], [14], [30], [32], [64]. In some cases, electrostatic interactions between diffusing species can have a strong effect on ion mobility and can produce unexpected behavior such as uphill diffusion [3], [4]. In addition, apparent diffusion coefficients (i.e. diffusion coefficients derived from Fick’s law) may show a strong dependency on concentrations. Considering that the quantification of diffusion coefficients is labor-intensive [1], [66], it is impractical to determine apparent diffusion coefficients as a function of solution composition for a range of conditions. Instead, it is advantageous to consider electrochemical interactions affecting diffusion explicitly rather than lumping this effect into empirically measured apparent diffusion coefficients. Species-dependent diffusion and electromigration can play a significant role for a variety of problems in the Earth sciences, in particular in the absence of advective transport and when the diffusing species are charged. For example, electromigration can be relevant during early diagenesis [2], [24], [64], in the near and far-field of radioactive waste repositories [31], [67], [68] and in the field of CO2 sequestration [65], [69]. Species-dependent diffusion and electromigration also play a role on a range of scales extending from the pore scale to the field scale [70]–[73]. Reactive transport models are commonly used for the quantitative investigation of flow, transport and reaction processes in porous media. These models aid with the verification of conceptual models, are used to design and evaluate experiments, and assist with the interpretation of field data in the fields of geology, engineering and environmental research [26], [30], [74]–[76]. Traditionally, diffusion has been 14implemented into reactive transport models based on Fick’s law and diffusion coefficients are often treated as adjustable parameters [1]. However, the number of reactive transport models that include electromigration and consider the chemical potential gradient as the driving force of diffusion has been growing in recent years [15], [20], [28], [56], [77], [78]. Although some aspects of electromigration on solute transport have been investigated [4], [28], [79], a direct model inter-comparison that specifically focuses on the role of electromigration and electrostatic effects on ion transport has not been published to date. This contribution was motivated by the need for benchmark problems suited to evaluate the effect of electric coupling during multicomponent diffusion and electrochemical migration and to facilitate an inter-comparison of existing reactive transport codes. The benchmark focuses on the assessment of electromigration and species-dependent diffusion in the context of the Nernst-Planck equation. More complex processes, such as anion exclusion and enhanced cation diffusion in the diffusive double layer, relevant in clays and other low permeability media (e.g. [3], [31], [68], [80]), are not included. These aspects are addressed in a companion contribution by Alt-Epping et al. [67]. The following benchmark problems were specifically designed to highlight effects of electromigration and are characterized by strong electric coupling effects, substantially more pronounced than in the contribution by Alt-Epping et al. [67]. The first two benchmarks are one-dimensional and the third benchmark includes two parts, involving one- and two-dimensional scenarios. Three reactive transport codes were used independently for the inter-comparison, namely CrunchFlow [14], MIN3P [59] and PHREEQC [15]. 152.2 Governing Equations 2.2.1 Mass Transfer in Electrolytic Systems Species-specific diffusion is necessary to describe the behavior of electrolyte systems [79] where diffusive transport is the dominant mass transport process. The most important feature that distinguishes the electrolyte systems from non-electrolyte systems is the electric coupling of the ionic fluxes [5], [17]. In the electrolyte systems, electric interaction of ion-ion, ion-solvent and ion-interface induces an electric field. The treatment of electrolytic diffusion follows naturally from the generalized treatment of diffusion [81]. 2.2.2 Nernst-Planck Equation for Multicomponent Systems The migration of interacting species is described by the Nernst-Planck equation, which can be derived from expressions for the diffusive flux written in terms of the chemical potential [14]. Written in terms of the aqueous flux of an arbitrary species i, the Nernst-Planck equation is given by: ۸୧ୟ = ϕτSୟ۲୧ −∇C୧ୟ − C୧ୟ∇lnγ୧ୟ −FRTC୧ୟz୧∇ψ൨ (2.1) where ϕ is porosity (m3 void m-3 bulk), τ is tortuosity (m void m-1 bulk), Sa is the saturation of the aqueous phase (m3 fluid m-3 void), Di is the species-dependent diffusion coefficient (m2 s-1), Cia is the concentration (mol L-1 H2O), γia is the activity coefficient (-), F is the Faraday constant (96485 C mol-1), R is the gas constant (8.341 J K-1 mol-1), T is the absolute temperature (K), zi is the charge number (-) and ψ is the electric potential (V or J C-1). In the presence of advection with a Darcy’s velocity qa (m s-1), the modified flux term is: 16۸୧ୟ = ϕτ ܙୟC୧ୟ−Sୟ۲୧∇C୧ୟ − Sୟ۲୧C୧ୟ∇lnγ୧ୟ −FRTSୟ۲୧C୧ୟz୧∇ψ൨ (2.2) This expression is known as the extended Nernst-Planck equation and holds, in ideal systems, for all mobile species. It describes the movement of ions in a solution with or without external electric field [6], [17], [62]. In a multicomponent system, the set of Nernst-Planck equations, one for each component, must be solved simultaneously. By assuming small gradients in ionic strength, a dilute solution with low ionic strength and isothermal conditions, the contribution of the flux from the gradients in the logarithms of the activity coefficients can be neglected [28], [79]. With this approximation, the flux of an individual species becomes: ۸୧ୟ = ϕτ ܙୟC୧ୟ−Sୟ۲୧∇C୧ୟ −FRTSୟ۲୧C୧ୟz୧∇ψ൨ (2.3) This equation represents the contributions of advection, diffusion and electromigration to the total mass transfer. Assuming there is no externally induced current (null current assumption), a simplified version of the mass flux can be derived [28], [56]: ۸୧ୟ = ϕτ ܙୟC୧ୟ − Sୟ۲୧∇C୧ୟ − Sୟ۲୧z୧C୧ୟ ܙ܉ ∑ z୩C୩ୟ౪౪୩ − Sୟ ∑ ۲୪z୪∇C୪ୟ౪౪୪Sୟ ∑ ۲୫z୫ଶ C୫ୟ౪౪୫൩ (2.4) where Ntot is the total number of dissolved species in the system. This formulation has the advantage that the electric potential does not appear as a primary unknown and is therefore well suited for implementation in standard reactive transport codes. 2.3 Participating Codes The three reactive transport codes participating in this benchmarking exercise are CrunchFlow [14], MIN3P [59] and PHREEQC [15]. CrunchFlow and MIN3P are 3D block-centered finite difference (finite volume) models using the global implicit method 17(GIM) to solve the fully coupled transport and reaction equations. PHREEQC solves the multicomponent diffusion problem with a 1D finite difference method using the sequential non-iterative approach (SNIA). A detailed description of the formulation and capabilities of the codes is discussed elsewhere [14] and references therein. PHREEQC considers the gradients of the activity coefficients in its implementation [68][4] whereas CrunchFlow and MIN3P neglect this contribution following the formulation of Giambalvo et al. [28] (see Equation A13). In PHREEQC, electromigration is implemented in the TRANSPORT framework and since it is a SNIA approach, transport calculations are separated from speciation calculations (e.g. the determination of OH- concentrations from water equilibrium reactions). CrunchFlow and MIN3P on the other hand use the GIM to implement multicomponent diffusion and electrochemical migration following the formulation of Giambalvo et al. [28]. In this approach, the total concentration terms are weighted by species-specific diffusion coefficients for both primary and secondary species when calculating the transport term. Both formulations account for contributions of both primary and secondary species towards the electrochemical migration term. 2.4 Benchmark Descriptions The three benchmark problems are summarized in Table 2.1 The first benchmark (1) focuses on the role of electromigration in driving the flux of the various charged species to maintain local charge balance and was first presented by Lichtner [82]. This problem considers diffusion of HNO3 from a low pH solution (pH = 4.007) into a circum-neutral reservoir (pH = 6.001) with low HNO3 concentrations, both with the same elevated NaCl background concentrations. 18Table 2.1. Summary of multicomponent diffusion benchmarks. Benchmark Description Processes Remarks 1 HNO3 (pH 4.007) diffusion into a circum-neutral pH reservoir Diffusion/Electromigration 1D 2 Sodium isotope fractionation induced by sodium chloride diffusion in neutral pH water Diffusion/Electromigration 1D 3 Transverse dispersion affected by electromigration Advection/Diffusion/ Electromigration 1D/2D The second benchmark (2) shows the electromigration and subsequent fractionation of the sodium isotope 22Na+ due to diffusion of NaCl under neutral pH conditions. The problem is loosely based on Glaus et al. [3]; however, the benchmark case presented here is set up for a uniform relatively coarse-grained uncharged porous medium and does not include diffusion through charged micropores as would be the case if porous clay were considered. In the first two benchmarks, diffusion and electromigration are the only transport processes and models are set up in one dimension. The third benchmark (3) investigates the effect of electromigration on transverse dispersion and is based on experiments and modeling carried out by Rolle et al. [70]. This benchmark also includes advection and is characterized by a higher level of complexity; it is simulated in one and two dimensions. 2.4.1 Benchmark 1: Transient Electromigration This problem was initially presented by Lichtner [82] and has previously been used as an example to illustrate the multicomponent capabilities of PHREEQC [15]. It is a 1D transient problem with a fixed concentration (Dirichlet) boundary condition on the left (at x = 0), representing the reservoir, and a no-flux (Neumann) boundary condition on the right (x = 0.01 m). The chemical system is composed of four primary (component) 19species (H+, NO3-, Na+ and Cl-) and one secondary species (OH-). The porosity is set to 1.0 and the domain is discretized into 100 equally spaced cells of 100 microns each. The temperature is 25ºC and there is no flow, the only transport process is multicomponent diffusion according to the Nernst-Planck equation (equation 2.4). Activity coefficients are calculated with the extended Debye-Hückel equation. Species-dependent diffusion coefficients, as well as the initial and boundary conditions defining the chemical system are given in Table 2.2. Table 2.2. Boundary conditions, initial conditions, and species dependent diffusion coefficients for Benchmark 1 (transient electromigration problem). Species Boundary Condition at x = 0 (mM) Initial Condition (mM) Diffusion Coefficient (m2 s-1) pH 6.001 4.007 9.31 × 10-9 Na+ 0.1 0.1 1.33 × 10-9 Cl- 0.1 0.1 2.03 × 10-9 NO3- 0.001 0.1 1.90 × 10-9 OH- a1.03 × 10-5 a1.06 × 10-7 5.27 × 10-9 aOH- concentrations are only provided for completeness, calculated from H+ and H2O (Kw = 10-14) The problem is run for 1 hour using a constant time step of 0.001 hour (corresponding to 1,000 time steps). Results are compared along the spatial profile after T = 1 hour for H+, Na+, NO3- and Cl-. 2.4.2 Benchmark 2: Tracer Isotope Diffusion This 1D problem involves three primary (component) species, Na+, Cl-, and H+, along with an isotope of Na that is also treated as a distinct component, 22Na+. In addition, a single secondary species, OH-, is considered. In this case, fixed concentration (Dirichlet) boundary conditions are considered at either end of the domain. The initial condition in the domain is divided into two regions; concentrations in half of the domain are 20equivalent to those at the left boundary, while concentrations in the other domain half are defined by the right boundary condition. However, the initial conditions are not significant since the simulation is run until conditions close to steady state are achieved. The gradients in activity coefficients are not considered in CrunchFlow and MIN3P, despite the significant gradient in ionic strength. The porosity is set to a constant and uniform value of 0.5 and the domain is discretized into 100 equally spaced cells of 100 microns each. The diffusion coefficients of Na+ and 22Na+ are assumed to be identical. A constant time step of 1 hour is used and the simulation is run to 1,500 days representative of conditions close to steady-state. Concentrations at the boundaries and species-dependent diffusion coefficients are described in Table 2.3. The simulation also assumes no flow. Table 2.3. Boundary conditions and diffusion coefficients for Benchmark 2 (isotope tracer problem). Species Left Boundary Condition (mM) Right Boundary Condition (mM) Diffusion Coefficient (m2 s-1) pH 7.0 7.0 9.31 × 10-9 Na+ 0.5 0.1 1.33 × 10-9 22Na+ 10-6 10-6 1.33 × 10-9 Cl- 0.5 0.1 2.03 × 10-9 OH- a1.05 × 10-4 a1.03 × 10-4 5.27 × 10-9 aOH- concentrations are only provided for completeness, calculated from H+ and H2O (Kw = 10-14) 2.4.3 Benchmark 3: Transverse Dispersion Rolle et al. [70] investigate the effect of electromigration on transverse dispersion under steady state flow conditions. In the full 2D case, the problem involves unidirectional flow and transport of a multicomponent tracer plume down the length of a 2D flow-through chamber. 21 Figure 2.1. Schematic of the 2D flow and transverse dispersion experiment conducted by Rolle et al. [70]. Using PHREEQC, Rolle et al. [70] solved the problem numerically by simulating transverse dispersion and electromigration perpendicular to the flow path as a 1D problem. This approach simplifies a 2D steady-state problem into a 1D transient problem by making use of the transformation t = x/v, where x is the distance from the source for the 2D problem, v is the uniform average linear groundwater velocity, and t defines the travel time to reach the location x. At the same time, t defines the simulation time for the 1D transient transverse dispersion problem [70]. Coinciding with experimental conditions, a 1 cm source in the middle of the 12 cm wide cross section at x = 0 m describes the continuous release of the electrolyte solution. The simulation was run for the case of an average linear velocity of 1.5 m day-1. The results of the 1D transient simulations are compared among the three participating codes, whereas fully 2D simulations with explicit treatment of flow were performed with CrunchFlow and MIN3P. 22Table 2.4. Chemical conditions and transverse dispersion coefficients for Benchmark 3 (transverse dispersion problem). Species Tracer Injection Ports (mM) Initial Condition (1D) and Remaining Injection Ports (2D) (mM) Diffusion Coefficient (m2 s-1) Transverse Dispersion Coefficient (m2 s-1) K+ 0.29 10-6 1.77 × 10-9 2. 405 × 10-9 Mg+2 0.29 10-6 6.26 × 10-10 1. 745 × 10-9 Cl- 0.87 3 × 10-6 1.81 × 10-9 2. 425 × 10-9 The dispersion coefficients used in these simulations require some discussion. In fact, the parameterization of the hydrodynamic transverse dispersion coefficient used in Rolle et al. [70] differs from the classical linear model commonly adopted in subsurface applications of solute transport and reads as: D୧ = D୧ + D୧ୟ୯ ቆPeଶPe + 2 + 4δଶቇஒ (2.5) where DiP is the pore diffusion coefficient approximated as the product of the aqueous diffusion coefficient of a species i and the porosity of the medium (0.41). Pe = vd/Diaq is the grain Péclet number where d is the average grain size (1.25 mm). δ = 6.2, and β = 0.47 are empirical parameters determined in previous multitracer experiments and pore-scale simulations [83]. Equation 2.5 explicitly retains a direct dependence of the mechanical dispersion term on the aqueous diffusivity of the transported species; the non-linear dependence on the average flow velocity arises from the incomplete mixing in the pore channels [84], [85]. For this benchmark analysis we considered the mixed electrolyte case described in Rolle et al. [70], where a dilute solution of KCl and MgCl2 was continuously injected in ambient deionized water. The free aqueous diffusion coefficients of the ions at T=20 ºC 23are DK+ = 1.77 × 10-9 m2 s-1, DMg+2 = 6.26 × 10-10 m2 s-1, and DCl- = 1.81 × 10-9 m2 s-1. Using these values in equation (5.2) yield the transverse dispersion coefficients given in the last column of Table 2.4. 1D Benchmark: The 1D benchmark consists of a pure transverse diffusion problem discretized into 48 grid cells of 2.5 mm. In the 1D system, the injection ports constitute initial conditions used at grid cells 23-26, corresponding to a 10 mm wide region in the center of the symmetrical system. The initial condition is used everywhere else in the domain and is intended to represent deionized water. The transverse dispersion coefficients given in Table 2.4 are used. The boundaries at either end of the system are treated as no-flux, but they do not influence the system behavior for the 16 hour simulation time (corresponding to x = 1 m, i.e. the outflow boundary of the domain). The simulation was run with a constant 0.001 hour time step. 2D Benchmark: For the full 2D problem solved with CrunchFlow and MIN3P, the transverse discretization is 50 grid cells with a spacing of 2.4 mm (corresponding to a total width of 0.12 m). At the inlet boundary, grid cells 24, 25, 26 and 27 in the transverse direction are set at the tracer injection port concentrations of 0.29 mM K+, 0.29 mM Mg+2, and 0.87 mM Cl- (see Table 2.4) while the remaining injection ports carry deionized water. The longitudinal discretization is 500 grid cells with a spacing of 2.4 mm thus a total length of 1.2 m; the concentrations are reported at x = 1.0 m, corresponding to the outflow boundary of the experimental setup. The additional length of 0.2 m is considered in the models to avoid any possible boundary effects. In this case, lateral flow can be calculated, or simply prescribed at 1.5 m day-1. A maximum time step of 1 hour is used with an initial minimum time step of 10-6 hours. The simulation time is 2432 hours to ensure that the final results correspond to steady state conditions representative of the experiment. 2.5 Results and Discussion 2.5.1 Benchmark 1 Simulation results for Benchmark 1 [82] depict the diffusion of HNO3 (pH = 4.007) from the solution domain towards the boundary where NO3- concentrations are 100 times lower and pH = 6.001. Results for NO3- and H+ reveal that both ions continue to diffuse towards the left boundary after 1 hour simulation time (Figure 2.2). Because the diffusion coefficient for H+ is much larger than the corresponding value for NO3-, H+ has become substantially more depleted in the domain than NO3-. The discrepancy in diffusion rates of H+ and NO3- triggered electromigration of Na+ and Cl- to maintain local charge balance; Na+ is entering the domain to offset the preferential loss of H+, while Cl- is leaving the system to counterbalance NO3-, which is preferentially retained. Migration of Na+ and Cl- occurs despite the fact that there was no initial concentration gradient for both species (Table 2.2) and takes place even against the developing concentration gradients of Na+ and Cl-. If Fick’s Law were used to describe this multispecies diffusion problem, there would be no change in Na+ and Cl- concentration and consequently electroneutrality would be violated. There is very good agreement between the simulation results of all three codes and demonstrated by near identical outputs. 25 Figure 2.2: Species concentrations after 1 hour simulation time for HNO3 diffusion (Benchmark 1). The left boundary is a fixed concentration (Dirichlet) boundary, while the right boundary is no-flux. 2.5.2 Benchmark 2 The results of the Benchmark 2 simulation visually show conditions close to steady state diffusion with identical concentration profiles and fluxes for Na+ and Cl- from left to right (Figure 2.3 and Figure 2.4). However, it has to be kept in mind that the diffusion coefficient for Cl- is considerably larger than the one for Na+. In fact, considering that the equations are based on the null current assumption, this holds back Cl- migration and accelerates Na+ migration. Although there are no initial concentration gradients for 22Na+, H+ and OH-, these species, present at much lower concentrations, also become affected by the electrostatic coupling. 26 Figure 2.3: Na+, Cl-, H+, OH- and 22Na+ concentrations after 1500 days for system summarized in Table 2 (Benchmark 2). The left boundary is a fixed concentration (Dirichlet) boundary at 0.5 mM, while the right boundary is a fixed concentration boundary at 0.1 mM for Na+, Cl-. The fixed gradient in NaCl results in a flux of H+, OH- and 22Na+, despite the fact that their concentrations are the same at either end of the column. A closer look at the results reveals that H+ migrates from the left to the right to enhance the positive charge flux, while OH- migrates from the right to the left to counteract the negative charge flux from the left to the right dominated by Cl- (Figure 2.4). The sodium isotope is present at very low concentrations and is also subject to a net migration from the left to the right. 27 Figure 2.4: Na+, Cl-, H+, OH- and 22Na+ diffusive, electromigration and total fluxes after 1500 days (representing conditions close to steady state) for the system summarized in Table 2 (Benchmark 2). H+, OH- and 22Na+ diffusive fluxes change direction at 0.65 cm. 28More importantly, the coupling between diffusive and electromigration fluxes induces a substantial isotope fractionation for 22Na+ (Figure 2.3). Diffusive fluxes of H+, OH- and 22Na+ change direction at 0.65 cm associated with a strong non-linear increase of the electromigration flux towards the right boundary (Figure 2.4). Solving this problem with Fick’s law would not predict 22Na+ isotope fractionation, H+ and OH- migration, and would result in a net negative charge flux across the domain. These results suggest that multicomponent diffusion can introduce isotope fractionation, even in the absence of fractionating reactions. Overall, there is very good agreement between the three codes with better agreement between CrunchFlow and MIN3P. Slight differences are observed for the PHREEQC results. It is difficult to decisively determine the reasons for these differences, but it is likely that the discrepancies are due to slight variations in model formulation (i.e. consideration of activity gradients in the PHREEQC formulation, absent in the other two codes). However, all codes show identical trends and concentration differences are small, implying that the residual discrepancies will not affect the interpretation of the results. 2.5.3 Benchmark 3: 1D Transverse Dispersion The transverse concentration profiles for Cl-, K+ and Mg+2 are plotted at the outlet (x = 1.0 m) corresponding to a residence time of 16 hours in the 2D domain. The separation of the three tracer profiles (Figure 2.5) demonstrates the effect of species-dependent dispersion coefficients and electrochemical migration on transverse displacement. The Cl- concentration profile is located between K+ and Mg+2 despite having the largest diffusion coefficient. In fact, DCl- in liberated state is considerably larger than DMg+2 and also slightly larger than DK+ (Table 2.4). These results show that Cl- migration is retarded due 29to electrostatic coupling with the cations and in particular with Mg+2, which diffuses more slowly. The outcomes reported in Figure 2.5 demonstrate the positive contribution of electromigration to transverse displacement of the two cations and the negative contribution of electromigration to transverse displacement of chloride [35]. Figure 2.5. 1D simulation results of transverse profiles for Cl-, K+ and Mg+2 at the outlet (corresponding to a residence time of 16 hours) demonstrate the effect of species-dependent dispersion and electromigration on the transverse displacement of charged species (1D Benchmark 3 solved with MIN3P). There is a very good agreement between the three codes and an excellent match between CrunchFlow and MIN3P (Figure 2.6). Similar to the two previous benchmarks, there are slight differences between the results of CrunchFlow and MIN3P on the one hand and PHREEQC on the other hand. Peak chloride concentrations predicted by PHREEQC are slightly higher than those calculated by CrunchFlow and MIN3P (~ 0.6%). Magnesium and potassium concentration profiles are in very good agreement for all codes. 30 Figure 2.6. 1D simulation of transverse multicomponent diffusion for the case of transport of mixed electrolytes (KCl and MgCl2 solution) in pure water described by Rolle et al. [70], comparing CrunchFlow, MIN3P and PHREEQC results. 2.5.4 Benchmark 3: 2D Flow and Transverse Dispersion Using CrunchFlow and MIN3P it was possible to carry out a full two-dimensional flow and multicomponent transport simulation of the flow-through system. The simulation was run for two pore volumes (32 hours) to ensure that steady state conditions at the outflow were reached. To illustrate the 2D concentration distributions and to provide a means for 31visual comparison of the CrunchFlow and MIN3P results, 2D contour plots are provided for K+, Mg+2 and Cl- (Figure 2.7). Figure 2.7. Simulation results for Benchmark 3 considering flow (uni-directional) and multicomponent transverse dispersion for steady-state conditions, from top to bottom are shown: K+, Mg+2 and Cl- for CrunchFlow and K+, Mg+2 and Cl- for MIN3P. Cross-sections extracted from two-dimensional steady state CrunchFlow and MIN3P results are compared at the outflow to one-dimensional transient PHREEQC 32results, corresponding to a residence time of 16 hours. Overall, there is an excellent agreement between MIN3P and CrunchFlow results (Figure 2.8) and results are also very close to the concentrations computed with the 1D PHREEQC approach. PHREEQC concentration profiles are slightly higher than CrunchFlow and MIN3P (the differences of the peak concentrations are ~ 0.6% for Cl-, ~ 0.7% for Mg2+ and ~ 0.4% for K+). Figure 2.8. Comparison of 1D PHREEQC results (no explicit consideration of flow, only following the plume as it moves down the flow path) and transverse profiles derived from 2D CrunchFlow and MIN3P runs for the transverse dispersion problem. The CrunchFlow 2D runs are based on GIMRT and use a first order upwind formulation, along with a backwards Euler time stepping approach, the same numerical methods are used in the MIN3P simulations. 332.6 Concluding Remarks Three benchmark problems were presented, each with significant effects of multicomponent diffusion and electromigration on transport of solutes in saturated porous media. The benchmarks were specifically designed to be sensitive to the effect of electromigration on diffusion and lateral concentration displacement. Benchmarks 1 and 2 are hypothetical problems that provide opportunities to verify the implementation of multicomponent diffusion and electromigration formulations in reactive transport codes. Benchmark 3 is based on the outcomes of laboratory experiments [70] and provides the opportunity to verify and validate multicomponent diffusion and species-dependent transverse dispersion formulations under flow-through conditions. Three reactive transport codes with the capability of simulating multicomponent diffusion and electrochemical migration participated in this study (CrunchFlow, MIN3P and PHREEQC). For all benchmark problems considered in this work an overall very good agreement between the simulation results obtained with the different codes. Despite some residual discrepancies between the simulation results, all three codes were able to consistently reproduce the same trends and evolution in concentration patterns induced by multicomponent diffusion and by the electrostatic interactions between the charged species. Small discrepancies between the results indicate that different approaches in implementing the governing equations are not a significant source of uncertainties for model applications; uncertainties will rather be dominated by the underlying conceptual model. 34Chapter 3: Numerical Simulation of Gas Generation and Reactive Solute Transport during Electrokinetic Remediation 3.1 Introduction Application of an electric potential gradient that drives an electric current between electrodes to transform and transport chemical species in soils has become popular as a cost-effective remediation technology [86]. This technology, known as electrokinetic remediation, can be used for desalination of sea or brackish water [87], [88], heavy metal remediation of soils [89]–[92], desalination of concrete [61] and to promote the decomposition of organic contaminants to less harmful and less toxic compounds [93]–[96]. Ionic solute transport in groundwater is commonly dominated by advection and diffusion, driven by pressure and concentration gradients, respectively. However, during electrokinetic remediation, other transport processes such as electromigration and electroosmosis, driven by electric potential gradients, are also of importance. The applied electric potential gradient between electrodes leads to oxidative conditions at the anode and reducing conditions at the cathode [97]. This modified redox environment provides thermodynamic conditions favorable for the degradation of contaminants to non-toxic and/or less-toxic compounds [93], [95]. However, the modification of redox conditions also leads to water electrolysis, producing H+ and O2(g) at the anode and OH- and H2(g) at the cathode, therefore creating acidic conditions at the anode and alkaline conditions at the cathode [86]. The generation of product species at the electrodes may have detrimental side effects and may inhibit electrokinetic remediation, potentially reducing its efficiency. For example, electrokinetic desalination 35of reinforced concrete can enhance corrosion processes [61], [98], [99]. In addition, it was shown that the precipitation of metal hydroxides in the alkaline region at the cathode inhibited transport of metals and reduced remediation efficiency [61], [100], [101]. Activation polarisation due to O2(g) and H2(g) gas bubble formation on the surface of electrodes tends to decrease the electric conductivity and consequently reduces remediation efficiency [19]–[22], [102]. These previous studies show that a thorough understanding of geochemical processes including the fate of dissolved species, as well as gases is required to assess the performance of electrokinetic remediation on a process-based level. Reactive transport codes have been used successfully for quantitative investigation of flow, transport and reaction processes in porous media [79]. Specifically, these models were applied to quantitatively assess conceptual models, to design and evaluate experiments, and to assist with the interpretation of field data in the fields of geology, engineering and environmental sciences [26], [30], [74], [103]. Reactive transport codes are capable of modelling a wide range of biogeochemical reaction processes under dynamic flow conditions [14]; thus providing a versatile framework for environmentally related applications. However, most reactive transport codes assume that ion migration is dominated by advective and diffusive solute transport. Fick’s law is commonly used to simulate the process of diffusion and diffusion coefficients are often treated as adjustable parameters [1]. Although the number of reactive transport codes that include electromigration and consider the chemical potential gradient as the driving force of diffusion has been growing in recent years [14], [19], [20], [28], [56], [77], [104], the implementations are 36commonly based on the null current assumption. The null current assumption implies that there is no electric current in the solution, which prevents this formulation from being used for simulating electrokinetic remediation. Integrating electrokinetic processes into the reactive transport framework is necessary in order to fully benefit from the flexibility provided by the multicomponent reactive transport approach [14], [74] and recent efforts have focused on the inclusion of electrokinetic simulation capabilities [19], [39]. For example, Paz-García et al. [19], using the finite element method, coupled a generalized Poisson-Nernst-Planck transport model [20] with a chemical equilibrium model to simulate the desalination of brick. Wu et al. [39] incorporated an electric potential distribution module into PHT3D, an operator-splitting reactive transport modeling code, to simulate electrokinetic in-situ oxidation [40]. However, the use of numerical models to investigate geochemical reactions associated with electrokinetic remediation remains limited to problem-specific cases and relatively simple geochemical reaction networks [38], [56]–[58]. The process of gas generation and exsolution at the electrodes has not been included to date. The present contribution was motivated by the need for incorporating electrokinetic phenomena into reactive transport codes in order to demonstrate their capability for evaluating electrokinetic remediation including an assessment for the potential of gas generation at the electrodes. To this end, the fully coupled 3D reactive transport code MIN3P [105] was extended by implementing the Gauss-Ampère formulation. In this formulation, the electric potential is expressed explicitly and is fully coupled with transport and geochemical reactions with appropriate boundary conditions, enabling the simulations of electrokinetic problems. Using this approach, charge 37continuity and mass continuity for all components are solved simultaneously, maintaining charge balance locally and globally at the same time. In previous works, depending on application, select chemical components were used to maintain the charge [38], [106], [107] or local charge balance was achieved by adjusting the concentration of inert ions [39]. In addition, the new model formulation covers a wide range of geochemical reactions that can potentially affect the electrokinetic remediation process, including the spontaneous exsolution of gases in response to gas generation at the electrodes. The enhanced code is verified against existing experimental observations and models from the literature [20], [39], [61]. The verification cases are focused on coupling mass and charge continuity in the context of the Nernst-Planck equation and in the absence of advective fluxes. The code’s capability to reproduce the evolution of electrochemical systems is demonstrated and further advanced by simulating gas generation and exsolution at the electrodes. Scenario and sensitivity analysis are performed to assess the importance and contribution of selected input parameters on the onset and vigor of gas generation, affecting total dissolved gas pressure and solution chemistry. 3.2 Governing Equations 3.2.1 Generalized Mass Continuity in Electrolyte Solutions In the absence of an externally applied electric potential gradient, the reactive transport problem involves the simultaneous solution of advective-diffusive-dispersive transport of dissolved species, diffusive gas transport, and the contribution from geochemical reactions involving aqueous, gaseous and mineral species [59]. The governing equation describing reactive transport in variably saturated media including electromigration is 38based on the modification of Fick’s law to include the Nernst-Planck equation [28], [104]. The mass continuity equation for component Aj, written in terms of total component concentration Tj, is: ∂∂tൣϕτSୟT୨ୟ൧ +∂∂tቂϕτST୨ቃ + ∇ ∙ ൣϕτܙୟT୨ୟ൧ − ∇ ∙ ൣϕτSୟ۲ୟ∇T୨ୟ൧ −∇ ∙ ϕτFRTSୟ۲ୟ∗ T୨ୟz୨∇ψ൨ − ∇ ∙ ቂϕτS۲∇T୨ቃ −Q୨ୟ,ୟ − Q୨ୟ,ୱ − Q୨ୟ,ୣ୶୲ − Q୨,ୣ୶୲ = 0, j = 1, Nୡ (3.1) where ϕ is porosity (m3 void m-3 bulk), τ is tortuosity (m void m-1 bulk), Sa is the saturation of the aqueous phase (m3 fluid m-3 void), Sg is the saturation of the gaseous phase (m3 gas m-3 void), Tja (mol L-1 H2O) defines the total aqueous component concentration for component Ajc, while Tjg (mol L-1 gas) is the total gaseous concentration for component Ajc. The total aqueous component concentrations Tja provide a relationship between primary (component) and secondary (complexed) aqueous species defined by T୨ୟ = C୨ୡ + ∑ ν୧୨୶ C୶୧౮୧ , where Cjc are the concentrations of the components as species in solution, Cix are the concentrations of secondary (complexed) species, ν୧୨୶ are the stoichiometric coefficients of the components in the secondary species, and Nx define the number of secondary species in the system. The gas concentrations are related to the components by defining the total gaseous component concentrations Tjg based on the law of mass action [108]. qa is the Darcy velocity vector, while Da and Dg (m2 s-1) are the dispersion-diffusion tensors for the aqueous and gaseous phase, respectively. Diffusion coefficients in Da are species-dependent. aD is the electroosmosis-electromigration tensor for the aqueous phase. zj is the generalized charge for component Ajc, and ψ is the electric potential in units of (V). F is the Faraday constant 39(96485 C mol-1), R is the gas constant (8.341 J K-1 mol-1), and T is the absolute temperature (K). Qja,a (mol dm-3 porous medium) and Qja,s (mol dm-3 porous medium) are internal source and sink terms toward the total aqueous component concentrations due to intra-aqueous kinetic reactions and kinetically controlled dissolution-precipitation reactions. Qja,ext (mol dm-3 porous medium) and Qjg,ext (mol dm-3 porous medium) are external source and sink terms for the aqueous and gaseous phase, respectively. A detailed description of the general reactive transport formulation excluding electromigration terms has been previously presented by [59] and is not described further. Here, we focus on the inclusion of electro-chemical migration in electrolyte solutions capable of transferring an electric current. The most important feature that distinguishes electrolyte systems from non-electrolyte systems is the electric coupling of the ionic fluxes [17]. In electrolyte systems, electric ion-ion, ion-solvent and ion-interface interactions induce an electric field. The treatment of electrolytic diffusion follows naturally from the generalized treatment of diffusion [81]. In this study we only consider ion-ion electric interactions; ion-interface electric interactions on the electrodes are considered as boundary conditions at the electrodes. 3.2.2 Nernst-Planck Equation for Multicomponent Systems For Nc components, equation (3.1) has Nc+1 unknowns; in order to solve it either one unknown must be eliminated or another equation must be coupled to it. In the first approach, called the null current assumption, it is assumed that there is no electric current in solution and the Nernst-Planck term (the fifth term in equation 3.1) is rearranged to meet the null-current constraint. In the second approach, Gauss-Ampère’s law (charge 40continuity) is coupled to the mass continuity equations and solved simultaneously. Both approaches are described in the following. The extended Nernst-Planck equation describing the migration of interacting individual species with low ionic strength in dilute solutions under isothermal conditions [28], [104] is defined by: ۸୧ୟ = ϕτ ܙୟC୧ୟ − Sୟ۲୧∇C୧ୟ −FRTSୟ۲୧∗C୧ୟz୧∇ψ൨ (3.2) where Cia (mol L-1 H2O) is the ith aqueous species concentration, Di is species-dependent diffusion coefficient and Jia is the ith species flux (mol m-2 s-1). Equation (3.2) describes the movement of ions in solution in the presence and/or absence of an external electric field [6], [17]. 3.2.3 Generalized Charge Continuity in Electrolyte Null Current: Equation (3.2) represents the contributions of diffusion, electromigration and advection to the total mass transfer. In a linear mass flow system, the total electric current density can be expressed as: ۷ୟ = F z୧۸୧ୟ౪౪୧ (3.3) where ۷ୟ is the total electric current density due to the migration of aqueous species (C s-1 m-2 or A m-2) and Ntot is the total number of aqueous species. Assuming there is no externally induced electric current, equation (3.3) becomes: ۷ୟ = F z୧۸୧ୟ = 0౪౪୧ (3.4) Substituting (3.2) into (3.4) gives an expression for the electric potential gradient as a function of individual concentrations of the species: 41∇ψ = RTFܙ܉ ∑ z୩C୩ୟ౪౪୩ − Sୟ ∑ ۲୪z୪∇C୪ୟ౪౪୪Sୟ ∑ ۲୫zଶ C୫ୟ౪౪୫൩ (3.5) Thus the mass flux of the ith species is: ۸୧ୟ = ϕτ ܙୟC୧ୟ − Sୟ۲୧∇C୧ୟ− Sୟ۲୧∗C୧ୟz୧ܙ܉ ∑ z୩C୩ୟ౪౪୩ − Sୟ ∑ ۲୪z୪∇C୪ୟ౪౪୪Sୟ ∑ ۲୫zଶ C୫ୟ౪౪୫൩ (3.6) In reactive transport codes with electrodiffusion capability equation (3.5) is commonly substituted into equation (3.1) [14]. By eliminating the electric potential term from the mass continuity equations, the electric potential is eliminated as a primary unknown, effectively implying a null current condition. Electroneutrality is maintained due to the fact that the null current assumption and electroneutrality are equivalent and differ only in the solution strategy [30]. Gauss-Ampère’s Law: To maintain charge continuity in the system, two constituent forms of Maxwell’s equations defined on the REV scale are needed [36]. The charge conservation results from combining the constituent form of Gauss and Ampère’s laws [109]. The macroscopic form of Ampère’s law for quasi-steady state conditions for the ith species in aqueous phase a is (following [36]): ϕSୟ۷ୟ = −∂(ϕSୟ۲ୟ )∂t+ સ × (ϕSୟ۶ୟ ) − સ × (ϕSୟ۾ୟ × ܞୟ) + ϕSୟܐመ ୠୟୟஷୠ (3.7) where ۲ୟ is electric displacement (C m-2), ۶ୟ is magnetic field intensity (A m), ۾ୟ is average polarization density (C m-2), ܞ܉ is the mass-average velocity vector (m s-1), ܐመ ୠୟ is external exchange current density between solution (phase a) and electrodes (phase b) (A m-2). By assuming no water polarization and taking the divergence of equation (3.7) to 42maintain charge continuity, the second and third terms on the right hand side of the equation can be eliminated: ∂ ቀ∇ ∙ (ϕSୟ۲ୟ )ቁ∂t+ ∇ ∙ (ϕSୟ۷ୟ ) − ∇ ∙ ൫ϕSୟܐመ ୠୟ ൯ୟஷୠ൩ = 0 (3.8) According to Gauss’ law: ∇ ∙ (ϕSୟ۲ୟ ) = ϕSୟF C୧ୟz୧౪౪୧ (3.9) leading to: ∂൫ϕSୟF ∑ C୧ୟz୧౪౪୧ ൯∂t+ ∇ ∙ (ϕSୟ۷ୟ ) − ∇ ∙ ൫ϕSୟܐመ ୠୟ ൯ୟஷୠ൩ = 0 (3.10) By substituting (3.2) into (3.3) for mass flux and substituting the result into (3.10) an expression for charge continuity in the aqueous phase can be derived (note the volume averaging for the aqueous phase): ∂൫ϕSୟF ∑ C୧ୟz୧౪౪୧ ൯∂t+ ∇ ∙ Fϕτ z୧ ൜ܞୟC୧ୟ − Sୟ۲୧∇C୧ୟ −FRTSୟ۲୧∗C୧ୟz୧∇ψൠ౪౪୧− ∇ ∙ ൫ϕSୟܐመ ୠୟ ൯ୟஷୠ൩ = 0 (3.11) When there is no externally induced electric current, the last term on the left hand side of (3.11) becomes zero. For a system with Nc components, Nc equations (3.1) and one equation (3.11) must be solved simultaneously. Equations (3.1) and (3.11) are functions of species concentrations (aqueous and gaseous phases) and electric potential. 433.2.4 Geochemical Relationships As a database-driven code [110], MIN3P is capable of simulating a wide range of geochemical processes [59], [110]–[114]. The present study required the consideration of aqueous complexation by hydrolysis reactions and ion pairing, gas dissolution-exsolution, and the consumption and production of aqueous components due to kinetically controlled intra-aqueous reactions. The formulation of the relevant geochemical reactions is provided in [59]. 3.2.5 Gas Exsolution and Degassing Production and exsolution of gas bubbles from the surface of the electrodes has been observed in electrokinetic experiments [115]. Degassing of O2(g), H2(g) ), N2(g) and Cl2(g) due to dissociation of the electrolyte [23] has been implemented in the code. Without gas exsolution processes, the total pressure increases continuously which is not physically realistic in an open system. Equilibrium gas dissolution-exsolution reactions is described by the mass action law: C୨ = RT ቀK୨ቁିଵෑ(γ୧ୡC୧ୡ)ౠౙ୧ୀଵ j = 1, N, (3.12) where Kjg is the equilibrium constant. According to the ideal gas law, partial pressures relate to the molar concentrations of gas species (Cjg [mol L-1 gas]): p୨ = RTC୨ (3.13) 44For common ranges of temperatures and pressures in shallow subsurface environments, we can assume that the use of the ideal gas law for real gases does not introduce significant errors [110], [116]. The model assumes that gas exsolution occurs when the sum of all gas pressures exceeds the confining pressure. The degassing rate is calculated using [111]: Rୢ = kୢmax ቌ∑ p୨୨pୟ + 0.9869 ρ୵g൫h୨ − z୨൯− 1ቍ , 0 (3.14) where kd is the rate constant for the degassing reaction (mol L-1 H2O s-1), pa is the atmospheric pressure (atm), ρw is the density of water (kg m-3) and g is the gravitational acceleration (m s-2). We assume that exsolved gas bubbles are instantaneously released into the atmosphere, gas bubble transport via ebullition is therefore not considered [111]. The hydrogen oxidation reaction (HOR) and oxygen reduction reaction (ORR) can occur in the presence of hydrogen-oxidizing bacteria [117]–[120] or catalysts [121], [122]. Hydrogen and oxygen can be consumed to form water: 12Oଶ(aq) + Hଶ(aq) → HଶO (3.15) In this context, it is assumed that the progress of reaction depends on hydrogen and oxygen gas availability described by the following Monod-type rate expression: Rమିୌమ = −kమିୌమ ቆ[Oଶ]Kమ + [Oଶ]ቇ ቆ[Hଶ]Kୌమ + [Hଶ]ቇ (3.16) where kమିୌమ is rate constant (mol L-1 s-1), and Kమ, Kୌమ are half-saturation constants (mol L-1) for oxygen and hydrogen, respectively. 453.3 Solution Procedure Equation (3.1), (3.6) and (3.11) are implemented into MIN3P [59]. With the enhancements, the code is capable of simulating multicomponent diffusion subject to the null current assumption or the Gauss-Ampère model which enables the user to define electric potential or current constraints for simulating electrokinetic remediation including geochemical reactions. MIN3P employs the modified Newton-Raphson method to solve the nonlinear system of equations in terms of concentration increments on a logarithmic scale and electric potential on a linear scale. The spatial discretization of the domain is based on the finite volume method which is locally mass and charge conservative, a requirement for solving the reactive transport equations using the global implicit approach. The temporal discretization of the governing equations is fully implicit with adaptive time stepping scheme. The total aqueous component concentrations used for the advective flux calculations across the interface between neighboring control volumes are determined by upstream or centered weighting, or a by a flux limiter scheme. Interfacial values for aqueous and gaseous saturations, porosities and diffusion coefficients are estimated using harmonic averaging based on parameter values in the adjacent control volumes. The geochemical reactions are accessible through a thermodynamic database and a user-defined input file. The set of algebraic equations yields a large Jacobian matrix which is solved using WATSOLV, a sparse iterative solver package [59], [123]. 3.4 Model Verification and Validation To verify the enhanced code, previously published experimental and modeling studies were used for evaluation of the different components and functionalities of the code. 46Three verification problems are presented. The first test case compares the null current method with the Gauss-Ampère approach for a free diffusion problem. The Gauss-Ampère model is verified against a Poisson-Nernst-Planck example by Paz-García [20]. The second test case compares the results of the MIN3P electrokinetic formulation with PHT3D and experimental measurements. The third problem investigates controlling factors on electrokinetic desalination, electrode kinetics and gas generation driven by electrokinetic processes and is based on experiments and modeling carried out by Kamran et al. [61]. The third problem is further extended by including, gas exsolution, degassing, and gas consumption kinetics to assess the likeliness for the occurrence of gas exsolution and degassing in the system. 3.4.1 Free Diffusion In order to evaluate the applicability of equation (3.6) and equations (3.1, 3.11) for simulating diffusion and electrochemical migration under null current conditions, a free diffusion example previously presented by Paz-García [20] has been adopted. Paz-García [20] employs a Poisson-Nernst-Planck formulation solved numerically with the finite element method. This test case is a fully saturated 1D transient problem with free exit boundaries on both ends of the domain (x = 0.0 and x = 0.1 m). The temperature is set to 25ºC and there is no flow; the only transport process is diffusion and electromigration. Two sets of simulations are performed: a) free diffusion described by the Nernst-Planck equation subject to the null current condition in the domain (equation 3.6) and b) free diffusion described by the Nernst-Planck equation coupled with Gauss-Ampère’s model (fully coupled equations 3.1 and 3.11). When using the Gauss-Ampère approach, it is necessary to set initial and boundary conditions for electric potential. The initial 47condition for electric potential is set to 0 mV for the entire domain; for the boundary conditions, the electric potential is fixed at 0 mV at x = 0.1 m as a reference point while no boundary condition is imposed at the other side (x = 0) of the domain. Initial and boundary conditions for the chemical composition of water are the same for both approaches. The initial electroneutral solution is composed of five primary (component) species (H2O, Na+, Cl-, K+ and H+) and one secondary species (OH-). The porosity is set to 0.3 and the domain is discretized into 100 equally spaced cells of 1 millimeter each. No reaction with the host porous medium is assumed. The Debye-Hückel equation is used as the activity coefficient model. Species-dependent diffusion coefficients [124], as well as the initial condition defining the chemical system are given in Table 3.1. Table 3.1. Chemical parameters: boundary conditions and diffusion coefficients for free diffusion problem. Species Left Initial Condition (M) Right Initial Condition (M) Diffusion Coefficient (m2 s-1)† pH 7. 0 7.0 9.31 × 10-9 Na+ 25.0 1.0 × 10-7 1.33 × 10-9 Cl- 25.0 1.0 × 10-7 2.03 × 10-9 K+ 1.0 × 10-7 1.0 × 10-7 1.90 × 10-9 OH- a1.0 × 10-7 a1.0 × 10-7 5.27 × 10-9 aOH- concentrations are only provided for completeness, calculated from H+ and H2O (Kw = 10-14) †From “Diffusion Coefficients in Liquids at Infinite Dilution” in [124] Based on pore geometry and pore particle shapes, Mohajeri et al. [125] suggests different tortuosity factors for cations and anions. The tortuosity factor are set to τc = 0.8 for cations and τa = 0.9 anions [20]. The physical parameters are summarized in Table 3.2. 48Table 3.2. Physical parameters for free diffusion problem. Parameter Symbol Value Unit Length L 0.1 m Cross-section area A 1.0 m2 Spatial discretization Δx 1.0 × 10-3 m Simulation time t 30 day Maximum time step Δtmax 1 day Porosity n 0.3 [-] Anion tortuosity† τa 0.9 [-] Cation tortuosity† τc 0.8 [-] †From [125] The problem is run for 30 days and with a maximum time step of 24 hours. The profiles are shown for 6, 12, 18, 24 and 30 days. For high concentration components (i.e. Na+ and Cl-), only the Na+ profile is shown; and for low concentration components (i.e. K+ and H+) only the pH profile is presented. Figure 3.1a depicts total mass flux (diffusion and electromigration) of Na+ taking place from left to right; electromigration flux holds back Cl- to maintain electroneutrality. It is seen in Figure 3.1b that at the same time, protons diffuse from right to left while electromigration causes proton migration from left to right. The total proton flux during early stages are from left to right while at the later stages fluxes take place from right to left; however, it is at early times that significant proton mass transfer to the right occurs; the amount of mass transfer from right to left is significantly less at later stages. The comparison of the results from null current and Gauss-Ampère approaches are nearly identical and cannot be distinguished visually; therefore the null current results are not shown in the plots. The results of Gauss-Ampère (this work) and Poisson-Nernst-Planck [66] approaches are very similar. 49 Figure 3.1. MIN3P Gauss-Ampère formulation (symbols) and Paz-García Poisson-Nernst-Planck formulation (solid line) [20]: (a) Na+ concentration profile; (b) pH profile and (c) diffusion potential profile (initial conditions for potential is zero for entire domain). 50Charge separation generates an electric potential known as “diffusion potential”, “liquid junction potential” or “junction potential” [126]–[128]. Figure 3.1c compares diffusion potential profiles using Gauss-Ampère and Poisson-Nernst-Planck methods. Decreasing concentration gradients in time are reflected in a decreasing electric potential (Figure 3.1). Initial steep concentration gradients at the junction yields an initially negative charge in the concentration region and a positive charge to the concentrated region because Cl- diffuses faster than Na+, creating a boundary potential difference. However, the generated electric field slows down Cl- migration and speeds up Na+ until the effective diffusion rates for all ions become equal. This process creates a detectable steady-state potential. It is not an equilibrium process [62] until the concentration gradient disappears. For the null current approach (Section 3.2.3), the potential is not expressed explicitly in the model and electroneutrality is imposed intrinsically. Therefore, electric potential and electric current are zero for the null current model. However, the Gauss-Ampère and Poisson-Nernst-Planck models predict the generation of an electric potential in the domain (Figure 3.1c). The electric current and charge density remain the same during the simulation within a range of +/- 10-13 A and +/- 10-8 C m-3 H2O respectively. The very small magnitude of the electric current in this process explains why the null current assumption still provides a comparable concentration gradient; however, it fails to provide information about electric potential. The Gauss-Ampère method has the advantage of integrating electrostatic physics into the problem based on process-based theoretical considerations. While using the null current approach provides similar results, our simulations demonstrate that in reality charge separation occurs, generating local electric currents, albeit small in magnitude. These results suggest that 51care must be taken when applying the null-current method to low concentration, smaller scale problems (Debye length order) [109]. 3.4.2 Electrokinetic Permanganate Transport In order to verify and validate the implementation of the electrokinetic model, reactive permanganate transport under an applied electric field through a 1D domain was simulated. The details of setup and methodology of the pH-isolated electrokinetic experiment (Figure 3.2) can be found in Hodges et al. [129]. Figure 3.2. Schematic of the 1D pH-isolated electrokinetic experiment conducted by Hodges et al. [129]. Electrode reservoirs are separated by clay cores to prevent pH disturbance of the electrokinetic cell. The three clay cores (1, 2, and 3) separated source and target reservoirs and the electrode reservoir at both ends of the column [129]. A potassium permanganate (KMnO4) solution was introduced to the source reservoir and permanganate (MnO4-) concentrations were monitored in the target reservoir. pH was maintained near-neutral in both electrode chambers by circulating electrode reservoirs fluids; The column was maintained under saturated conditions with no hydraulic gradient across the domain. [129]. Under these conditions, permanganate transport was driven by diffusion, electromigration and electroosmosis according to equations (3.1) and (3.10). The experiment can be conceptualized as a 1D transient problem with a fixed concentration 52boundary condition in the source reservoir (x = 0), and free exit boundary condition in the target reservoir (x = 0.1 m). The chemical system is composed of K+ and MnO4-. The porosity is 0.3 and the domain is discretized into 50 equally spaced cells of 2.0 × 10-3 m each. A voltage gradient of 2.1 V m-1 was applied across core 2, which defines the solution domain (see Figure 3.2). Electroosmotic conductivity is set to 2.14×10-9 m2 V-1s-1. Flow induced by electroosmosis was neglected, considering that it has no major effect on the final outcome. The assumed tortuosity value is 0.2 [39]. The temperature is 25ºC and activity coefficients are calculated with the extended Debye-Hückel equation. Physical and chemical parameters as well as boundary and initial conditions of the model are summarized in Table 3.3. Table 3.3. Boundary condition, initial condition, and model parameters used for permanganate removal problem. Parameter Symbol Value Unit KMnO4- (initial concentration) C0 0.0 M KMnO4- (boundary concentration) Cb 0.19 M Free diffusion coefficient of K+† D0,K+ 1.96×10-9 m2 s-1 Free diffusion coefficient of MnO4-† D0,MnO4- 1.63×10-9 m2 s-1 Length L 0.1 m Cross-sectional area A 8.49×10-3 m2 Spatial discretization Δx 2.0×10-3 m Simulation time t 320 h Maximum time step Δtmax 144 s Potential gradient Δφ 2.1 V m-1 Electro-osmotic conductivity* keo 2.14×10-9 m2 V-1 s-1 Porosity n 0.3 [-] Tortuosity* τ 0.2 [-] †From “Diffusion Coefficients in Liquids at Infinite Dilution” in [124] *From [39] 53The problem is run for 320 hours using a maximum time step of 144 seconds. Permanganate concentrations in the target reservoir are compared with the corresponding PHT3D results and measurements [39] over time. Although there are small differences between the MIN3P and PHT3D results, there is satisfactory agreement between the results of the two codes and the measurements (Figure 3.3). The discrepancies between the model results can be attributed to differences in the modeling approach. PHT3D considers only MnO4- and uses a different free diffusion coefficient (2.0×10-9 m2 s-1) for permanganate, while MIN3P considers the coupled diffusion of K+ and MnO4- using literature values [124] for species-dependent diffusion coefficients. On the other hand, PHT3D accounts for electroosmotic driven flow, which is neglected by MIN3P. Moreover, PHT3D uses a sequential non-iterative operator splitting approach with calculating potential in a separate step after transport while MIN3P is fully coupled with transport and potential equations in the same time step using the global implicit approach. Figure 3.3. MIN3P (solid line), measured (symbols) and PHT3D (dashed line) [39]: MnO4_concentration in target reservoir over time. 543.4.3 Electrokinetic Desalinization Electrokinetic desalinization is a technology that is based on the controlled application of direct current through the soil between appropriately distributed electrodes to mobilize aqueous species in order to remediate contaminated soils. It is well known that during electrochemical remediation, electrolysis of water occurs at the electrodes. At the anode water oxidizes: HଶO − 2eି ↔ 12Oଶ(aq) + 2Hା; E° = −1.229 V (3.17) while at the cathode water reduces: 2HଶO + 2eି ↔ Hଶ(aq) + 2OHି; E° = −0.828 V (3.18) where Eº are the standard reduction electrochemical potentials of the half reactions, which is a measure of the tendency of the reactants to proceed to products in their standard states. The type of electrolysis reactions that take place at the electrodes depends on the availability of the chemical species and the electrochemical potential of these reactions. For example, chloride oxidation is a competitive reaction at the anode that may produce chlorine gas [19], [98]: 2Clି − 2eି ↔ Clଶ(aq); E° = −1.360V (3.19) Reaction (3.17) causes acidification at the anode and reaction (3.18) causes the development of alkaline conditions in the vicinity of the cathode. Aqueous species are transported towards the electrodes by diffusion, electromigration and electroosmosis [86], [130]. Because sodium reduction is much less favorable than water (standard reduction electrochemical potential for sodium is -2.710 V), we neglected this process and assumed only hydrogen gas production at the cathode. 55An experimental and modeling study of water desalination has been carried out by Kamran et al. [61]. Their experimental data has been used as an additional test case for the verification of the enhanced MIN3P code. Specifically, the base case scenario (1) and scenario 2 in this section constitute direct comparisons with the experimental data and modeling results by Kamran et al. [8]. A fired-clay brick specimen with porosity of 0.3 was vacuum saturated with a 4 M sodium chloride solution at 25ºC (room temperature) and was placed between two meshed platinum electrodes [61]. The electrodes on both sides of the sample were in indirect contact with large water reservoirs separated by water soaked sponges. The permeable electrodes allowed the exchange of ions between the reservoirs and the sample. The reservoirs were also in contact with the atmosphere (Figure 3.4). Figure 3.4. Schematic of the 1D representation of electrokinetic desalination experiment conducted by Kamran et al. [61]. The model includes seven aqueous components (Na+, Cl-, H+, O2(aq), H2(aq), N2(aq),Cl2(aq)), two secondary aqueous species (OH-, NaCl(aq)) and four gaseous species (O2(g), H2(g), N2(g) , Cl2(g)). A 1D model with 0.1 m length was set up. Initial tortuosity was set to 0.25 for all components but adjusted later. Electroosmotic 56conductivity was set to 2.14×10-9 m2 V-1s-1. The applied electric potential was 5.5 V. The model was discretized spatially into 50 equal size 2.0×10-3 m cells and run for 16 hours for the base case scenario and scenario 2. There was no hydraulic gradient across the sample. Boundary and initial conditions, as well as species-specific diffusion coefficients are summarized in Table 3.4, Physico-chemical model parameters are presented in Table 3.5 and Table 3.6. Table 3.4. Boundary conditions, initial conditions, and species-specific diffusion coefficients for desalination problem. Species *Initial Condition *Boundary Condition Diffusion Coefficient (m2 s-1) Aqueous species pH 7.0 7.0 †9.31 × 10-9 Na+ 4.0 0.0 †1.33 × 10-9 Cl- 4.0 0.0 †2.03 × 10-9 OH- a1.0 × 10-6 a1.0 × 10-6 †5.26 × 10-9 NaCl a1.26 0.0 2.0 × 10-9 O2(aq) a2.65 × 10-4 a2.65 × 10-4 2.0 × 10-9 N2(aq) a5.15 × 10-4 a5.15 × 10-4 2.0 × 10-9 *Aqueous concentrations are in mol L-1 and gaseous concentrations are pressure in atm. aConcentrations are only provided for completeness; calculated by speciation. †From “Diffusion Coefficients in Liquids at Infinite Dilution” in [124]. 3.4.3.1 Electrode Processes and Boundary Conditions Water electrolysis at the electrodes generates H+ and OH- at a rate proportional to electric current according to Faraday’s Law [115]. In the presence of other species participating in electrolysis reactions, the electric current is apportioned between active species depending on their transference number: ۷୨ୟ =t୨z୨F ϕSୟܐመ ୠୟୟஷୠ (3.20) 57where tj is the transference number of the participating species (e.g. H+, OH-, Cl-) carrying different proportions of the current [131]. In the present model, only water reduces at the cathode producing OH- and H2(aq), therefore tୌష = 1. But at the anode water is being oxidized to H+ and O2(aq), and chloride is oxidized to Cl2(aq). As a result, the applied electric current, ܐመ ୠୟ , must be distributed between H2O and Cl- depending on H+ generation rate and Cl- consumption rate. The transference number for each species depends on the charge transfer rate at the electrodes following the “Frumkin-Butler-Volmer” expression [132]: ۸୨ୟ = K,୨Cୖ,୨ୟ exp ൬−12FRTΔψୱ൰ − Kୖ,୨C,୨ୟ exp ൬12FRTΔψୱ൰ (3.21) where the O and R subscripts denote oxidation and reduction reactions, respectively, K are the oxidation/reduction rate constants, C are the concentrations of the reacting species in the reduced and oxidized state at electrode surface, and Δψୱ is the electric potential drop across electrode surface and electrolyte. Potential drop is on the order of volts per centimeter [133] and it is determined by measuring actual potential across the domain compared to applied potential. Therefore in the model, the initial measured potential is considered as effective potential. The KO,Cl-/ KO,H2O ratio is determined following Kamran et al. [61]. 3.4.3.2 Simulation Results In the base case scenario (1) degassing was not considered. No gas formation was reported during the experiment [61]. In the second scenario (2), the degassing model (equation 3.14) was used to determine if gas formation and degassing has an effect on the results. Electroosmosis was also considered for the base scenario (1) but it did not have any significant effects on the experimental results as the observed water velocity (∼5 × 5810−8 m s-1) was more than one order of magnitude smaller than ion migration velocity (∼10−6 m s-1) [61]. To reduce execution time, electroosmosis was neglected for all other scenarios simulated. Experimental measurements include Na+ concentrations in 2 hours intervals by Nuclear Magnetic Resonance (NMR) Imaging and concentrations measurements of Na+ and Cl- after completion of the experiment by ion chromatography [61]. Base case scenario (1): The comparison of simulation results of MIN3P with experimental measurements is presented in Figure 3.5. To calculate electric potential, Kamran et al. [8] integrate electric current across the domain decoupled from the mass transport equation; their model only includes water self-dissociation reaction [61] and their simulation results are not compared and presented here. The MIN3P results confirm the experimental findings and simulation results by Kamran et al. [8]. During the desalinization process, Na+ is becoming depleted from both sides of the sample (Figure 3.5a and c). Na+ depletion at the anode side occurs more rapidly than on the cathode side. A stagnation point develops near x = 33 mm from the anode (Figure 3.5a). Modeling results confirm the findings by Kamran et al. [61], indicating that the position of the stagnation point depends on hydrogen mobility. The stagnation position cannot be matched using literature values for the H+ diffusion coefficient. It was therefore reduced by a factor of four, in-line with the approach of Kamran et al. [8], who reduced the H+ diffusion coefficient by a factor of three to match the location of the stagnation point [61]. It has been observed previously that proton mobility reduces in concentrated salt solutions [134], [135]. For example, Soudijn [134] reported a fivefold reduction in the 59 Figure 3.5. (a) Na+ concentration profiles for selected times (0, 2, 4, 6, 8, 10, 12, 14 and 16 hours): model (solid sharp lines) vs. experimental measurements by NMR (shadow background lines); the arrows show concentration profile evolution over time; (b) potential profile for 0, 10, 12 and 16 hours: model (lines) vs. experimental measurements (symbols); (c) Na+ and Cl− concentration profile at t = 16 hours (end of experiment): model (solid and dashed sharp lines), experimental measurements of Na+ concentration by NMR (shadow background lines) and experimentally measured Cl− concentrations by ion chromatography (symbols); (d) acid and alkaline fronts progress (pH = 1 and pH = 13 respectively): model (lines) vs. measurements (symbols). Simulation results are for Scenario 1 and 2 and experimental data are taken from Kamran et al. [61]. 60proton effective diffusion coefficient in a 3 M salt solution compared to the average diffusion coefficient in a solution containing no salt [134]. This reduction of mobility is related to suppression of water molecule reorientation by salt ions [136], [137] or proton interaction with the pore surface [61], [138]. Redistribution of electric resistivity changes the electric potential distribution over time (Figure 3.5b). The electric potential drop location coincides with the low concentration region. As resistivity increases the role of electromigration decreases and diffusion plays a more pronounced role, resulting in a competition between diffusion and electromigration in this region. Despite the fact that a voltage of 9 V is applied, the potential drop between electrode and bulk solution reduces the actual measured potential [61]. Based on the measurements (Figure 3.5b) we used a 5.5 V potential difference in the model. The only chloride measurements are by ion chromatography after completion of the experiment ([61], Figure 3.5c). Kamran et al. [61] report that during sample preparation some salt might have been lost, possibly providing an explanation for the differences between observations and simulated results. The competition between the oxidation of chloride and water at the anode is used as an adjustable parameter in the model. The best fit is obtained with KO,Cl-/ KO,H2O = 20. This ratio does not have a significant effect on the location of stagnation point, which is mainly controlled by H+ mobility. By increasing the KO,Cl/KO,H2O ratio, minor increases in pH at the acidic front was observed, as generation of Cl2(g) at the anode inhibits production of H+. The acidic and alkaline fronts originating from the electrodes encounter each other at 1/3 of the distance from the anode corresponding to the stagnation point (x = 33 mm). The encounter happens after 10 hours of treatment (Figure 3.5d). Towards the end of the simulation, a sharp pH change from pH = 1 to pH = 13 can be observed at this 61location. This stagnation point defines the boundary between removal and accumulation regions. On the left side and closer to the anode, Na+ is removed and Cl- accumulates and on the right side and closer to the cathode, the opposite situation occurs. When the alkaline and acidic fronts meet, the fronts stop moving and electrokinetic ion transport becomes inhibited [61]. The excess of positive charge near the anode and negative charge at the cathode attract opposite charges and repel similar charges, leading to the development of low ion concentrations at the stagnation point and reduction of conductivity. At this point, the desalinization process becomes ineffective and remediation stops. This is observed in the experiment starting at around 10 hours [61]. The results presented here demonstrate the capability of the model to reproduce the main processes in the system. Scenario 1, 2 (Gas generation and degassing): These two scenarios are extensions of the interpretations by Kamran et al. [8] to investigate the effect of gas generation by including the degassing model (equation 3.14). Since there was no gas measurement carried out in the experiment [61], the simulations are used as a tool to assess the system dynamics and highlight the capabilities of the code. Scenario 1 is a direct extension of the base scenario by including gas generation, but without degassing, while Scenario 2 considers rapid degassing of generated gases at the electrodes. As there is no indication of pressure build up in the electrode chambers in the Kamran et al. [61] experiment, we assumed open electrode reservoirs, initially in equilibrium with atmospheric pressure. The pressure head in the vicinity of the electrodes has a significant effect, because degassing starts only when total gas pressure exceeds confining pressure (equation 3.14). The maximum degassing rate is set to 1.0 × 10-5 mol L-1 s-1 to represent 62rapid gas exsolution and avoid total gas pressure build-up beyond the confining pressure. The simulated reaction network is presented in Table 3.6. The evolution of generated gases at the electrodes is presented in Figure 3.6 and Figure 3.7. Competition between chloride and water at the anode determines the production rate of oxygen and chlorine gas. Therefore, the consumption of Cl- at the anode by production of chlorine is controlled by the KO,Cl/KO,H2O ratio. Keeping the same ratio (KO,Cl/KO,H2O = 20) as in the base case for Scenarios 1 and 2 yields the same aqueous concentrations and electric potential (Figure 3.5) with no noticeable differences relative to the base case. Decreasing the KO,Cl/KO,H2O ratio results in an increase of oxygen production and elevated total pressures at the anode. However, this modification does not lead to a better fit of simulated aqueous concentration profiles, because chloride concentrations in the aqueous phase are electrostatically linked to other species. This result indicates that competition between H2O and Cl- at the anode is controlling chloride consumption and chlorine generation, and degassing does not play a significant role for this experimental setup. The generation of H2(g) at the cathode and Cl2(g) and O2(g) at the anode increase total gas pressure. 63 Figure 3.6. Gas evolution over time at anode (solid line) and cathode (dashed line): (a and c) base case (Scenario 1) including gas generation, but without degassing; (b and d) including gas generation and degassing (Scenario 2). Cl2(g) is plotted separately (c and d) as its contribution to total pressure was minimal with respect to other gas components. Degassing starts from the beginning and maintains total gas pressure by not allowing the sum of all gas pressures to exceed the confining pressure. Degassing strips out N2(g) from the mixture to maintain pressure (cf. a and b). 64The pressure increase at the cathode is greater than at the anode, because the production rate of oxygen (equation 3.17) is half that of hydrogen (equation 3.18). Moreover, oxygen solubility (1.2×10-3 mol L-1 atm-1) is approximately one order of magnitude higher than that of hydrogen (1.0×10-4 mol L-1 atm-1), thus keeping more O2 in the aqueous phase. Chloride oxidation also produces chlorine gas at the anode (equation 3.19) increasing chlorine gas pressure (Figure 3.6c). The contribution of chlorine gas generation to total pressure is not very significant which can be attributed to its relatively small contribution in relation to water oxidation and its higher solubility in water (9.1×10-2 mol L-1 atm-1) (Figure 3.7c and Figure 3.7d). When total gas pressure exceeds the confining pressure (atmospheric pressure) nitrogen gas starts to exsolve from the solution at the electrodes which is more pronounced in the cathode than at the anode. Nitrogen is stripped out of solution and therefore its pressure declines (Figure 3.6b, Figure 3.7b and Figure 3.7d), maintaining the total pressure at the confining pressure, despite continuous gas generation. Hydrogen also exsolves from solution at the cathode. Hydrogen pressure at the end of process (t = 16 h) in Scenario 2 is 0.176 atm compared to 0.202 atm in Scenario 1. 65 Figure 3.7. Simulated total and partial pressures and dissolved gas profiles during electrokinetic desalination for selected times (0, 2, 4, 6, 8, 10, 12, 14 and 16 hours): (a and c) base case (Scenario 1) including gas generation, but without degassing; (b and d) including gas generation and degassing (Scenario 2). Degassing starts after 13 hours and 20 minutes, maintaining total pressure by not allowing the sum of all gas pressures to exceed the confining pressure. Degassing strips out N2(g) from the mixture to maintain total gas pressure at the confining pressure. 663.5 Evaluation of the Potential for Gas Generation, Exsolution and Consumption The experimental conditions by Kamran et al. [61] involved a scenario with high salt concentrations. As was shown above, ion separation and increasing electric resistivity inhibited electrokinetic desalination and prevented further gas generation. In the absence of the salt or lower salt concentrations, remediation processes continue for longer time periods. On the downside, gas generation at the electrode will continue to progress with more substantial effects than shown in the experiment by Kamran et al. [61]. Such scenarios are relevant for numerous groundwater remediation problems involving trace metals or trace organic contaminants. To investigate this scenario in more detail, Na+ and Cl- were removed and the model was run for a longer period of time up to 64 hours. This scenario provides a more suitable platform to further evaluate the effects of gas generation at the electrodes. Four hypothetical scenarios have been run to evaluate the effects of gas generation, exsolution and also gas consumption in the region where O2 and H2 are present. In the first case, it was assumed that degassing and gas consumption do not occur (Scenario 3). In this scenario generation of hydrogen and oxygen at the cathode and anode, respectively, continue and increase total pressure. The second case assumed instantaneous degassing with a degassing rate of 1.0 × 10-5 mol L-1 s-1, but also without gas consumption (Scenario 4). The rapid degassing is representing gas bubble formation at the electrodes in the early stages of the process. In the third case, a gas consumption reaction (equation 3.15) is included (Scenario 5), in the absence of degassing. The gas consumption model (equation 3.16) is a kinetic Monod-type reaction which accounts for the oxidation of H2(aq) by O2(aq). The effective rate constant of 1.0 × 10-7 mol L-1 s-1 is introduced. In the absence of experimental data, the reaction rate is 67chosen relatively fast deliberately to account for the high energy yield of this reaction and in order to capture the interaction between degassing and reaction during the simulation time. In a follow-up sensitivity analysis, multiple realizations with reaction rates varying by 5 orders of magnitude have been run. The last case includes degassing and gas consumption (Scenario 6). Table 3.5 summarizes model parameters for the different Scenarios and Table 3.6 presents the geochemical reactions considered for each scenario. In all scenarios, molecular diffusion and electromigration are the main transport mechanisms in the aqueous phase. The current model does not explicitly consider the process of ebullition. 68 Table 3.5. Physico-chemical model parameters used in desalination problem. Parameter Symbol Unit Scenarios †1 †2 3 4 5 6 Simulation time t h 16 16 64 64 64 64 Maximum time step Δtmax s 36 36 288 288 288 288 Length L m 0.1 Spatial discretization Δx m 2.0 × 10-3 Initial pressure head ψ m 0.01 Electroosmosis conductivity* keo m2 V-1 s-1 2.14 × 10-9 Electrode permeability‡ β m2 s-1 5.0 × 10-6 Porosity‡ n [-] 0.3 Tortuosity‡ τ [-] 0.25 Potential‡ φ V 5.5 5.5 5.5 5.5 5.5 5.5 Oxidation rate constants ratio†† KO,Cl-/ KO,H+ [-] 20 20 - - - - Degassing - - No Yes No Yes No Yes Gas consumption - - No No No No Yes Yes Degassing rate constant kd mol L-1 s-1 - 1.0 × 10-5 - 1.0 × 10-5 - 1.0 × 10-5 Gas consumption rate constant kO2-H2 mol L-1 s-1 - - - - 1.0 × 10-7 1.0 × 10-7 †Na+, Cl-, NaCl, Cl2(aq) and Cl2(g) are present in base scenario 1 and 2. *From [39]. ‡From [61]. ††Calibrated. 693.5.1 Simulation Results Scenario 3: Generation of hydrogen and oxygen at the cathode and anode increase total gas pressure. At the cathode pressure increases to 2.18 atm and at the anode an increase to 1.36 atm can be observed at the end of simulation (64 h). Figure 3.8a illustrates gas evolution at the electrodes. Scenario 4: In this scenario instantaneous degassing is assumed. Degassing strips nitrogen, hydrogen and oxygen out of the solution to maintain total gas pressure below the confining pressure. At the cathode nitrogen and hydrogen exsolve significantly more out of the solution than oxygen, while at the anode mainly nitrogen is exsolved. Oxygen partial pressure at the end of the simulation at the anode is 0.48 atm compared to 0.57 atm in scenario 3 without degassing (Figure 3.8b). During the simulation, the total gas pressure at the electrodes remains almost constant at near atmospheric pressure, representative of instantaneous degassing conditions. Scenario 5: Degassing is deactivated in this scenario and a Monod-type kinetic (equation 3.16) gas consumption reaction (equation 3.15 or reaction 10 in Table 3.6) with an effective rate constant of 1.0 × 10-7 mol L-1 s-1 is introduced. The final total gas pressures at the electrodes decrease as result of hydrogen and oxygen consumption in relation to Scenario 3. At the cathode, generated hydrogen immediately reacts with dissolved oxygen initially present in the water, until oxygen is depleted at around 57 hours, temporarily dropping total pressure to 0.85 atm. In the following, pressure starts to increase to reach 1.28 atm at the end of the simulation as the H2 consumption reaction is inhibited by a lack of oxygen. At the anode, due to the absence of hydrogen in the water, no reaction occurs and generated oxygen contributes directly to total gas pressure. The 70maximum total gas pressure at the anode is similar to Scenario 3 and reaches 1.35 atm at the end of simulation (Figure 3.8c). The oxidation of hydrogen by oxygen inhibits gas exsolution, as long as total pressure is kept below the confining pressure (demonstrated in the next scenario); and if the groundwater condition are initially aerobic. Scenario 6: This scenario includes degassing and gas consumption reactions at the same time. Gas consumption at the cathode prevents total gas pressure to exceed confining pressure; therefore degassing does not occur at the cathode for most of the simulation period. Increasing pressure at the anode by oxygen generation initiates degassing by stripping nitrogen out of the solution and maintaining confining pressure (Figure 3.8d). At the cathode, hydrogen starts to react with ambient oxygen from the beginning and stops at about 57 hours when oxygen is completely depleted. From this point hydrogen pressure starts to rise until total pressure reaches the confining pressure and gas exsolution eventually initiates at around 61 hours. These scenarios demonstrate the dynamics between gas exsolution and geochemical reactions that consume gases generated at the electrodes. There is a competition between the processes of degassing and the oxidation of hydrogen by oxygen to maintain the total gas pressure. Rapid degassing in the vicinity of the electrodes will inhibit the oxidation of hydrogen by oxygen, in particular considering the fact that degassing occurs by desorption of gas bubbles from the electrode surfaces. Therefore, it becomes more difficult for H2 and O2 generated to migrate away from electrodes. 71Table 3.6. Reactions network and corresponding coefficients. No. Name †Reaction †Reaction coefficients Units Remarks Scenarios 1 H2O dissociation H2O ↔ H+ + OH- log K25 = -13.998 mol L-1 equilibrium 1,2,3,4,5,6 2 *NaCl dissociation NaCl ↔ Na+ + Cl- log K25 = -0.7770 mol L-1 equilibrium 1 and 2 3 H2O oxidation H2O ↔ 2H+ + 2e- + 0.5O2(g) E0 = -1.229 V at anode 1,2,3,4,5,6 4 *Cl- oxidation 2Cl- ↔ Cl2(g) + 2e- E0 = -1.360 V at anode 1 and 2 5 H2O reduction H2O + 2e- ↔ 2OH- + H2(g) E0 = -0.828 V at cathode 1,2,3,4,5,6 6 *Cl2(g) dissolution Cl2(aq) ↔ Cl2(g) kH = 9.1 × 10-2 mol L-1 atm-1 Henry’s law 1 and 2 7 H2(g) dissolution H2(aq) ↔ H2(g) kH = 1.0 × 10-4 mol L-1 atm-1 Henry’s law 1,2,3,4,5,6 8 O2(g) dissolution O2(aq) ↔ O2(g) kH = 1.2 × 10-3 mol L-1 atm-1 Henry’s law 1,2,3,4,5,6 9 N2(g) dissolution N2(aq) ↔ N2(g) kH = 6.5 × 10-4 mol L-1 atm-1 Henry’s law 1,2,3,4,5,6 10 O2(g) / H2(g) consumption H2O ↔ H2(aq) + 0.5O2(aq) kO2-H2 = 1.0 × 10-7 mol L-1 s-1 Monod kinetic 5 and 6 †Henry’s constant taken from Sander 2015 [139]. *Na+, Cl-, Cl2(g) and Cl2(aq) are only considered in scenario 1 and 2. 72 Figure 3.8. Gas evolution over time at anode (solid lines) and cathode (dashed lines): (a) degassing and gas consumption not active (Scenario 3); (b) degassing active, but no gas consumption (Scenario 4); (c) degassing not active, but including gas consumption (Scenario 5); (d) degassing and gas consumption active (Scenario 6). Degassing maintains total pressure by not allowing the sum of all gas pressures to exceed the confining pressure (b and d). Gas consumption reaction first depletes oxygen at the cathode then lets the pressure increase (c and d). 733.5.2 Sensitivity and Scenario Analysis The influence of the different processes on total gas pressures was discussed in the previous section. In order to further investigate the dynamics between the different processes and their implications on electrokinetic remediation design, sensitivity and scenario analyses were used. Sensitivity analysis is a tool to identify the contribution of input parameters on model outputs, specifically total gas pressure and composition at the electrode in this case. For the sensitivity analysis, Scenario 6 was selected, which represents the most comprehensive test case by considering degassing and the oxidation of hydrogen by oxygen. Sensitivities were calculated with respect to nominal values of potential (5.5 V), hydrogen oxidation rate (1.0 × 10-7 mol L-1 s-1) and the rate of degassing (1.0 × 10-5 mol L-1 s-1). The electric potential was varied between 3.5 and 11 V, the hydrogen oxidation rate ranged from 1.0 × 10-8 mol L-1 s-1 to 1.0 × 10-4 mol L-1 s-1 and the degassing rate was varied between 1.0 × 10-10 mol L-1 s-1 and 1.0 × 10-5 mol L-1 s-1. The wide range of degassing and reaction rates was chosen to account for the uncertainty and to cover the full range from slow to fast progress. Dimensionless (dss) and composite scaled sensitivities (css) [140] were calculated using UCODE [141]. It was found that the electric potential and rate of degassing provide the main controls on the total gas pressure and their effects are greater at the cathode, because gas production is greater than at the anode (Table 3.7). The composite scaled sensitivities are computed using model observations of total pressure in vicinity of the electrodes. Electric potential is the most sensitive parameter with the largest effect on total gas pressure at the electrodes, even though there is less uncertainty in determining electrode potential in comparison to other parameters. The rate of degassing is the second most sensitive parameter with a negative sensitivity, implying that increasing the rate of degassing decreases total pressure. The model shows the least sensitivity to the rate of hydrogen 74oxidation with almost no sensitivity at the anode, because there is no hydrogen present to trigger the reaction. Table 3.7. Dimensionless and composite scaled sensitivities for Scenario 6 calculated by UCODE [141]. Parameter Dimensionless Sensitivity (dss) Composite Scaled Sensitivity (css) Anode Cathode Potential 81.8 220.8 578.7 Degassing rate -19.5 -44.8 19.7 Reaction rate 0.0 -0.1 8.5 Scenario analysis provides a tool to evaluate the effect of perturbing sensitive parameters on model output with respect to the reference model with nominal values. For scenario analysis, 380 realizations of scenario 6 were produced using PEST [142] by assuming uniform distributions of electric potential, degassing and reaction rates between upper and lower bounds. The effect of the input parameters on total gas pressure at the end of simulations (t = 64 hours) is depicted in Figure 3.9. Total gas pressure increases at the anode with potential increase (Figure 3.9a). In the low electric potential range between 3.5 V to 4.5 V, pressure at the cathode drops from initial pressure considerably due to the hydrogen oxidation reaction, which consumes both ambient oxygen and produced hydrogen. By increasing potential and consequently the gas production rate, total pressure quickly increases because the rate of hydrogen production at the electrode outpaces the rate of hydrogen oxidation. For potentials greater than 5.5 V, the generation rate exceeds hydrogen oxidation and degassing rates and the total gas pressure exceeds confining pressure. Total gas pressure increase at the anode is directly controlled by the applied potential and the rate of degassing. When the electric potential becomes greater than 6 V, the rate of 75oxygen generation exceeds the rate of degassing and total gas pressure increases above the confining pressure. Hydrogen oxidation does not affect total pressure at the anode significantly (Figure 3.9b) due to a lack of H2(aq) in this region. At the cathode, three different steps are noticeable in total gas pressure curve as a function of the hydrogen oxidation rate. For rate coefficients up to approximately 7.0 × 10-10 mol L-1 s-1, the consumption rate is too slow to affect total gas pressure. The range from 7.0 × 10-10 to 3.0 × 10-8 mol L-1 s-1 is a transition region, showing increasing effects on the total gas pressure. For rate coefficients greater than 3.0 × 10-8 mol L-1s-1 total gas pressure is maintained at the confining pressure (Figure 3.9b). The effect of the rate of degassing on total gas pressures at the electrodes is demonstrated in Figure 3.9c. The rate of degassing controls total gas pressure at both electrodes in a similar fashion. From 1.0 × 10-8 to 1.0 × 10-5 mol L-1 s-1, increasing the rate of degassing decreases total dissolved gas pressure, but still remains above confining pressure; at 1.0 × 10-5 mol L-1 s-1, total gas pressure is maintained at the confining pressure, and a further rate increase does not affect total gas pressure at the electrodes. During electrokinetic remediation, potential is generally controlled and regulated either by applying a constant potential or a constant electric current, which reduces uncertainty with respect to the rate of gas production at the electrodes. However, gas accumulation and adsorption on the electrodes can reduce electric conductivity of the electrodes [143], [144] and create a decline of the effective potential between the electrodes [22], consequently reducing the efficiency of the method. The results of sensitivity and scenario analyses suggest that for high voltage, when the rates of degassing and hydrogen oxidation are outpaced by gas generation, gas accumulation and adsorption can be diminishing method efficiency. 76 Figure 3.9. Total gas pressure at anode (solid lines) and cathode (dashed lines) at the end of simulation and as a function of: (a) electric potential; (b) rate of hydrogen oxidation and (c) rate of degassing. (Note: confining pressure is 1.0097 [atm]). 3.6 Summary and Conclusions The development of a general purpose reactive transport code incorporating electrochemical processes based on the Nernst-Planck system of equations has been presented. Electrode boundary conditions were implemented for application to electrokinetic remediation problems. The resulting model facilitates the study of interactions between electrode reactions, gas generation and consumption, as well as solute transport in electrokinetic systems. 77The necessity of coupling mass transport and electrochemical processes to simulate electrokinetic remediation has been demonstrated. Two physically different approaches, one based on the null current assumption and the other based on the direct integration of an electric field (Gauss-Ampère model) were implemented and compared. For a free diffusion problem without an applied external potential gradient, near-identical results were obtained from both methods in terms of simulated concentration distributions, despite different physical and mathematical basis. However, using the Gauss-Ampère method provides additional information on the generation of electric potential and electric current within the domain, triggered by charge separation. This example also demonstrates that a weak electric current can generate detectable electric potential which is not predictable using the null-current assumption. Using the null-current assumption eliminates the electric potential term from the Nernst-Planck equation and the electric potential field cannot be calculated explicitly. In the presence of an external electric field, it is necessary to use the Gauss-Ampère method. The implementation of the Gauss-Ampère method was verified with two test examples. The first test case included pH-isolated electrokinetic transport of permanganate preventing alkaline and acidic fronts from advancing into the domain. The second verification example focused on the electrokinetic desalination of brick. For the permanganate problem, the fully coupled model was able to reproduce experimental observations and results from PHT3D, a two-step electrokinetic reactive transport code and experimental observations. Good agreement with experimental results was also obtained for the desalination problem. To demonstrate the model capabilities and in order to evaluate the potential for gas generation, exsolution and consumption, the desalination problem was extended to include water electrolysis and four different scenarios were investigated. Water electrolysis generates hydrogen 78and oxygen at the electrodes increasing total gas pressure and possibly driving gas exsolution. The build-up of total gas pressure was controlled by including degassing and combined hydrogen oxidation and oxygen reduction. Scenario and sensitivity analyses provide a quantitative assessment on the roles of electric potential, degassing and hydrogen oxidation on the build-up of total gas pressure. An understanding of these interactions is particularly important for designing electrokinetic remediation technologies for heavy metal removal, where a decrease in current density and removal efficiency [38], [86], [145] has been observed. Utilizing reactive transport codes that include a complex geochemical network can be used to further assess the performance of electrokinetic remediation. 79Chapter 4: A Process-based Reactive Transport Model for Understanding Self-potential Signals Associated with Buried Ore Bodies 4.1 Introduction Self-potential signals are associated with a variety of processes occurring in subsurface environments such as the natural attenuation of contaminant plumes [46], [146], the development of redox fronts [48], [49], as well as the corrosion of buried metallic bodies [12] and the weathering of mineral deposits [147]–[149]. Self-potential signals can typically be characterized by an ohmic potential drop that is measurable at the ground surface. Sulfide minerals, most notably pyrite and pyrrhotite, are largely responsible for self-potential signals related to buried ore bodies. Spatially separated redox reactions associated with the ore body create an electric current through the ore body, generating the self-potential signal. The oxidation-reduction reactions are generally metal or sulfide mineral oxidation at depth leading to the transfer of electrons towards the ground surface, facilitated by the conductivity of the ore body, and oxygen reduction near the ground surface, receiving the electrons. These reactions initiate an array of subsequent geochemical processes that tend to alter the ambient groundwater geochemistry considerably. The generated electric field and induced electrochemical transport processes contribute to further geochemical alterations of groundwater within and in the vicinity of the ore body. In this context, the conductive ore body and the surrounding aquifer effectively act as a “geobattery”. The first electrochemical model of the geobattery concept describing the principles behind self-potential signals of massive mineral deposits was proposed by Sato and Mooney [42]. Later, Stoll et al. [10] and Bigalke and Grabner [11] presented a mathematical formulation of the geobattery model and were able to reproduce self-potential signals measured on the 80surface above a deep continental drilling site in Germany [10], [11] (Stoll et al., 1995; Bigalke and Grabner, 1997). The Stoll et al. and Bigalke and Grabner's mathematical “inert electrode model” is based on electrochemical potential equilibrium between an ore body and an ambient electrolyte, and relies on the approximation of the vertical redox potential gradient [10]. In the “inert electrode model” the chemical composition and properties are constant and only electron transfer is taken into account [10], [11]. Bigalke and Grabner [10] conducted laboratory experiments with synthetic materials to calculate the redox potential field and electrochemical reactions at the mineralization surface. The assumed reactions were oxygen reduction at the cathode and iron and methane oxidation at the anode. The oxygen diffusion coefficient and standard current density were treated as fitting parameters in their study. The Bigalke and Grabner simulations were able to reproduce measured self-potential signals, but did not take into account geochemical changes in the subsurface [11]. To date, many aspects of self-potential signals cannot be explained by the existing models [148]. For example, the inert electrode model is not able to explain the effect of reactive minerals on the electric field [44]. There is no study available in the literature for forward modeling of self-potential signals associated with reactive minerals. Although the association of self-potential signals with contaminant plumes [46]–[48], [146], redox fronts [48], [49], corrosion of buried metallic bodies [12] and ore deposits [44] has been investigated using inverse modeling, so far there is no forward model that includes the wide range of (bio)geochemical (kinetic and equilibrium) reactions accounting for porous media properties, aqueous composition and the various transport phenomena that ultimately control and determine electric currents in the geochemical systems and the associated potential signals. 81Recent investigations of relationships between self-potential and redox potential have focused on the mapping and assessment of contaminant plumes. Self-potentials have been associated with redox potential gradients resulting from the biodegradation of dense nonaqueous phase liquids [150], [151], bacterially-mediated organic carbon oxidation in landfill leachates [43], [47], [48], [146], [150]–[155], and corrosion of buried iron bars [12], [50]. In this study, we present the first process-based reactive transport model to investigate the relationship between self-potential surface signatures and multicomponent reactive transport processes in the subsurface. Processes considered by the model include geochemical reactions such as aqueous complexation, mineral dissolution-precipitation and gas dissolution-exsolution along with mass and charge transport by diffusion, electromigration and electron drift through the ore body that acts as an electrical conductor. There are a limited number of previous studies on this subject and their two major limitations are 1) performing only inverse modeling for self-potential source currents or dipole localization [44], [47]–[50], [146], [156]–[158] and 2) not considering the geochemical reactions generating the redox gradient in the forward model [50], [158]. The present work focuses on a quantitative and process-based description of the geobattery concept related to self-potential signals and geochemical reactions within and in the vicinity of an ore body. Specific objectives of this contribution include: 1) the development of the governing equations describing the geobattery concept in a multicomponent reactive transport framework and implementation of the equation into the reactive transport code MIN3P [59], [104], 2) evaluation of the model based on a published dataset from a controlled sandbox experiment focusing on the corrosion on a buried iron bar and 3) to demonstrate, based on an illustrative example, that a process-based reactive transport model can be used to simulate the 82geochemical and transport processes that generate redox anomalies and self-potential signatures associated with buried ore deposits. In addition, we also investigate the effect of geochemical and physical properties on self-potential response. To serve these objectives, we extended the fully coupled 3D reactive transport code MIN3P [59], [104] by implementing the Gauss-Ampère formulation to couple electric potential implicitly with relevant transport processes and geochemical reactions. 4.2 Conceptual Model In the absence of a conductive ore body, redox potential in the subsurface is affected by the ingress of atmospheric oxygen, which enters with infiltrating rainwater and by diffusion, generating oxidizing conditions in the shallow subsurface. Mineral phases containing metals in reduced state (e.g.: Fe(II), Mn(II)), disseminated sulfides, or organic carbon fragments act as effective redox buffers, consuming oxygen and maintaining reduced conditions in the deeper subsurface. However, the presence of a conductive ore body may allow oxidation reactions to occur at greater depth, effectively leading to the formation of two spatially separated half-cell oxidation-reduction reactions (Figure 4.1). The mobility of electrons through the conducting ore body enables electron transfer between the spatially separated reaction sites [159]. Therefore, when an electrical conductive metallic body is present, an upward electron flow in the ore body is established facilitated by deep anodic reactions releasing electrons as a result of oxidation reactions (e.g. corrosion reactions), and near surface cathodic reactions typically by the reduction of oxygen radicals [43], [160]. Simultaneously, the ionic current in the surrounding electrolyte establishes upward migration of cations from the anodic region at depth to the cathodic region near the ground surface and downward migration of anions from cathode to anode, both processes contributing to the generation of a potential field [10], [11], [160]. Electrokinetic 83processes such as diffusion, electromigration, electroosmosis and streaming potential contribute to the movement of charged species and generate electric potential. The electric current through the conducting ore body, known as “source current” [12], [43] or “forcing current” [50], contributes to the total electric current flux [10], [11], [42], [43] and affects the self-potential measurable at the ground surface. Figure 4.1. Geobattery concept. An electric field is established when a gradient of the redox potential is connected by an electrical conductive body. Oxidation-reduction reactions (oxygen reduction near the surface and metal oxidation e.g. metal or sulfide oxidation at depth) are responsible for the generation of an electric current by electron flux through the ore body. The generated electric field initiates electromigration of ions in the opposite direction in the surrounding aquifer to maintain electroneutrality. Geochemical reactions such as mineral dissolution-precipitation, gas generation and exsolution may occur in the surrounding environment and within the ore body. 84This system acts as a geobattery in which the electric current through the electrical conductive minerals or other forms of conductors, such as bacterial conductive pili, and/or nanowires [43], is balanced by ionic electrochemical migration. In this study we establish the link between geochemical reactions occurring in subsurface environments and the generation of surface-measureable electric potential signals by reactive transport modeling. 4.3 Governing Equations As described in the previous chapter, the migration of interacting individual aqueous species with low ionic strength in dilute solutions under isothermal conditions [28], [104], and in the absence of advection, is defined by the Nernst-Planck equation: ۸୧ୟ = −ϕτSୟ۲୧ ∇C୧ୟ −FRTC୧ୟz୧∇ψ൨ (4.1) where ϕ is porosity (m3 void m-3 bulk), τ is tortuosity (m void m-1 bulk), Sa is the saturation of the aqueous phase (m3 fluid m-3 void), Cia (mol L-1 H2O) is the ith aqueous species concentration, Di is species-dependent diffusion coefficient (m2 s-1), zi is the charge of a species, and ψ is the electric potential (V). F is the Faraday constant (96485 C mol-1), R is the ideal gas constant (8.341 J K-1 mol-1), and T is the absolute temperature (K). Jia is the ith aqueous species flux (mol m-2 s-1). This equation describes the total mass flux of aqueous species in the sediments as a function of total concentrations and electric potential. To complete the geobattery model, the Nernst-Planck equation requires an additional term to account for movement of electrons in the conductor called electron drift. Electron drift in the electrical conductor can be expressed as [161]: ۸܍ = −enμୣF∇ψ (4.2) 85where e is the elementary electron charge (1.6 × 10-19 C), n is number of free electrons per volume (m-3) and μe is electron mobility (m2 V-1 s-1). In the absence of an external electric current, the total electric current associated with diffusion, migration and electron drift is: ۷۸ = F ൭ z୧Jia + Je୧൱ (4.3) where ۷۸ୟ is the total electric current density due to the movements of aqueous species, ۸ܑ܉,(C s-1 m-2 or A m-2) and includes the electron flux in the electrical conductor, ۸܍. N is the total number of aqueous species. In the aquifer, ۸܍, is equal to zero because there is no conductor and electron drift. In the absence of external electric current and electrical conductive media, the total electric current and resulting electric potential are small and commonly termed as "shale potential", "liquid junction potential" or "diffusion potential" [8]. An example of diffusion potential was presented in the previous chapter in a free diffusion verification problem where mixing of two aqueous solutions with different ionic strength generated electric potential across the diffusion cell. According to the geobattery model presented by Revil and co-workers [12], [43], [46], [47], [49], [146], [162], an electric current, termed as source (or driving) electric current density, ۷܁ୟ (A m-2), is established in the presence of an electrical conductor and a redox potential gradient, facilitating electron transfer from the reducing regions to the oxidizing regions. The source electric current density can be defined as: ۷܁ୟ = −σ୫∇ൣ൫E୫_ୟ୬ − Eୌ_ୟ୬൯ − ൫E୫_ୡୟ − Eୌ_ୡୟ൯൧ (4.4) where σm (S m-1) is the electric conductivity of the electrical conductive medium (ore or metallic body) and EH_an and EH_ca (V) are the equilibrium redox potentials of anodic and cathodic 86reactions, respectively. Em_an and Em_ca are the redox potentials of the electrical conductive medium at the anode and cathode, corresponding to the top and bottom of the ore body. The redox potential of the electrical conductive medium eventually becomes constant (Em_an = Em_ca) as the chemical potential of electrons is constant inside the conductor [11], [12]. Here we adopted the “active electrode model” after [12], [46], [163], which assumes participation of the ore body in geochemical reactions and for simplification considers a linear relationship between the redox potential and source current density. The total electric current density, ۷ܜܗܜୟ , is the sum of total electric current density due to the movements of aqueous species in the aquifer and electron drift in the ore body, ۷۸ୟ , and source current density ۷܁ୟ [12], [43], [46], [157], [164]: ۷ܜܗܜୟ = ۷۸ୟ + ۷܁ୟ (4.5) By substituting equations (4.3) and (4.4) into equation (4.5) and applying Gauss-Ampère’s law (i.e. continuity of charge written in terms of electric flux), the temporal and spatial distribution of concentrations and the electric potential will be obtained: સ ∙ ۷ܜܗܜୟ = સ ∙ ൣ۷۸ୟ + ۷܁ୟ൧ (4.6) and: ∂൫ϕSୟF ∑ C୧ୟz୧ା୧ ൯∂t+ સ ∙ ۷ܜܗܜୟ = 0 (4.7) substituting equations (4.3) and (4.4) into equations (4.6) and then (4.7): ∂൫ϕSୟF ∑ C୧ୟz୧ା୧ ൯∂t− સ ∙ ቌFϕτSୟ ۲୧z୧౪౪୧∇C୧ୟ +FRTC୧ୟz୧∇ψ൨ − (1 − ϕ) enμୣ∇ψቍ − સ ∙ ൫ߪ∇ൣ൫E୫ − Eୌ൯ − ൫E୫ౙ − Eୌౙ൯൧൯ = 0 (4.8) 87Since electric conductivity of the surrounding porous media (σm) is very low in comparison to the electrical conductive body, the source term will be negligible for the surrounding sediments. The electron movement is restricted to the conductor according to Gauss's law; in this model the electrical conductor zone transmits electrons from a high potential region to a low potential region (equation 4.2). A more detailed step-by-step derivation of the equations is presented in Appendix A. 4.4 Model Implementation and Verification Equation (4.6) without the addition of source electric current was implemented into previously verified and validated developments of MIN3P described in the preceding chapter [59]. In the previous model for simulation of electrokinetic remediation, the electrode boundary conditions were applied at the electrode-aquifer interface (equation 3.21 in chapter 3). For the present model, equation (4.8) was implemented in MIN3P and electric flux boundary conditions are specified at the ore-aquifer interfaces (i.e. anode-aquifer and cathode-aquifer) where the sum of the electron flux through the conductor is equal to the sum of the electric flux in the aquifer induced by ion migration. The electrical conductive zone is defined as a co-existing medium where electrons are drifting from high potential to low potential regions. A vertical electric flux in the conductor was assumed, which only can enter and exit at the anodic and cathodic ends of the ore body. Outside the electrical conductive zone, electrons are immobile. Since there is no other forward reactive transport numerical model coupling geochemical transport and reactions with self-potential signals, the new model was validated against published results of a controlled sandbox iron corrosion experiment conducted by Castermant et al. [12]. 884.5 Experimental Case Study Castermant et al. [12] studied the relationship between redox potential and self-potential anomalies by inserting a 50 × 2 × 2 cm iron bar in a 200 × 50 × 7 cm sandbox filled with sand composed of 95% silica, 4% orthoclase and less than 1% albite. The sandbox was carefully filled by an electrolyte of known composition (1×10−2 M KCl and 2×10−3 M NaOH) with an initial pH of 7. The use of NaOH was intended to accelerate the process of iron corrosion. The upper boundary of the tank was in contact with air. Surface self-potential signals were measured by non-polarizing electrodes. The redox potential, self-potential and pH were measured in a regular grid 42 days after the iron bar insertion into the sandbox [5]. The iron bar was inserted vertically at the left-hand side of the tank. The results were interpreted in terms of a geobattery model combined with iron oxidation reactions. Here iron acts as the electron donor and electrical conductor to transport electrons released by iron oxidation towards the top of the sandbox, where oxygen reduction takes place. The iron oxidation reactions lead to the corrosion of iron generating ferrous or ferric iron and dissolved gases. The pH increased significantly during the experiment and the redox potential of the ambient water decreased [12]. The pH increase initiated precipitation of secondary minerals on the surface and in the vicinity of the iron bar which partially buffered further pH increases. Visual observations of blue-greenish sand in the vicinity of the iron bar suggested the formation of green rust represented by fougerite (Fe(OH)2). Upon exposure to air, unstable Fe(OH)2 oxidized to an ochre-coloured mineral which was likely goethite or lepidocrocite (FeOOH). Atmospheric oxygen and CO2 gases are diffusing into the tank across its surface. Oxygen plays the role of the electron acceptor in this system. The charge transfer resulting from iron oxidation and oxygen reduction reactions is essential to establish a measurable surface self-potential response. This charge transfer is controlled by the electric 89conductivity of the conductor. In the absence of a conductive medium, only a very small electric potential field can be generated by diffusion potential. Geochemical Relationships: the geochemical reaction network considered for this study can be represented by the following 8 components: H+, Na+, K+, Fe+2, Fe+3, Cl-, CO3–2 and H2O and 3 gases: O2(g), H2(g), CO2(g) and 4 possible secondary minerals: Fe(OH)2(am), FeOOH, Fe(OH)3(am) and FeCO3(am). Fe+2/Fe+3 was considered as an equilibrium redox couple. Electron migration was restricted to the conductive ore body and the electron was also used as a master variable to calculate redox potential. A total number of 9 aqueous complexation reactions from the MIN3P database were also included in the model. Reactions at the electrodes and mineral dissolution-precipitation: electrode and mineral precipitation reactions are assumed to progress in parallel with no inhibitive or competitive interaction. A general rate expression for the surface-controlled irreversible dissolution was used for oxidation of iron [111], [165]. The rate constants were obtained by calibration to match the observed pH and self-potential signal. The precipitation of secondary mineral phases is represented by the following rate expression: R୫ = −kୣ ൬1 −IAP୫K୫൰ (4.9) where keff is an effective rate constant for mineral precipitation, IAPm is the mineral ion activity product and Km is the mineral equilibrium constant. The reaction stoichiometries of electrode reactions, mineral formation reactions, gas exchange and aqueous complexation are summarized in Table 4.1. 90Table 4.1. Reactions network and corresponding coefficients [166]. Name Reaction log K log keff mol L-1 H2O s-1 Electrode reactions Anode Fe(s) → Fe+2 + 2e− − 14.9 − 8.3† Cathode O2(aq) + 4e− + 2H+ → 2H2O − 86.018 - Secondary minerals reactions Fougerite Fe(OH)2(am) ↔ Fe+2 + 2H2O − 2H+ − 13.905 − 9.6 Ferrihydrite Fe(OH)3(am) ↔ Fe+3 + 3H2O − 3H+ − 4.891 − 9.6 Lepidocrocite FeO(OH) ↔ Fe+3 + 2H2O − 3H+ − 1.371 − 9.6 Siderite FeCO3(am) ↔ Fe+2 + CO3‒2 10.450 − 9.9 Gas dissolution-exsolution reactions O2(g) dissolution O2(aq) ↔ O2(g) 2.898 mol L-1 atm-1 H2(g) dissolution H2(aq) ↔ H2(g) − 43.009 mol L-1 atm-1 CO2(g) dissolution CO2(aq) ↔ CO2(g) 18.160 mol L-1 atm-1 Redox couples Fe2+ /Fe3+ Fe2+ → Fe+3 + e− 13.032 - Aqueous complexation OH- H2O ↔ H+ + OH- − 13.998 - O2(aq) O2(aq) + 4e− + 2H+ ↔ 2H2O − 86.018 - H2(aq) H2(aq) ↔ 2e− + 2H+ − 3.124 - HCO3‒ HCO3‒ ↔ CO3‒2 + H+ − 3.617 - H2CO3(aq) H2CO3(aq) ↔ CO3‒2 + 2H+ − 2.247 - FeHCO3+ FeHCO3+ ↔ Fe2+ + CO3‒2 + H+ 12.330 - FeCO3(aq) FeCO3(aq) ↔ Fe2+ + CO3‒2 + H+ 4.380 - FeCl2+ FeCl2+ ↔ Fe3+ + Cl‒ 1.480 - FeCl+ FeCl+ ↔ Fe2+ + Cl‒ 0140 - †Calibrated. Physical Parameters: a 2D model of the sandbox experiment with overall dimensions of 2 m length and 50 cm depth was constructed and discretized into 1184 control volumes. In the experiment, the iron bar had a thickness of 2.0 cm. The sandbox thickness was 7.0 cm with a porosity of 0.34 for the sand; since the iron bar was not present over the entire thickness of the sandbox, the volume fraction of iron was adjusted accordingly in the 2D model and porosity of the porous media in the iron bar section was reduced to 0.24 to account for the presence of the iron bar. The domain was divided into two zones: “background” from 0.05 to 2.0 m representing 91the aquifer with identical material properties and initial conditions, and the iron bar zone from 0.0 to 0.05 m with an iron bar volume fraction of 0.28 m3 m-3 (Figure 4.2). Figure 4.2. (a) Domain and boundary conditions and (b) MIN3P model grid applied for solving iron bar case study. The initial volume fractions of secondary minerals were set to 0.0. The initial geochemical composition of the pore water in the entire domain is listed in Table 4.2. The diffusion coefficients were taken from [124] for available ions. When diffusion coefficients were not available, they were calculated based on the Stokes-Einstein equation [167]. The physical parameters of the porous medium are given in Table 4.3. Equilibrium constants for geochemical reactions were derived from the MINTEQA2 database [166]. Details can be found in [108] and [111]. The temperature was set to 25ºC and there was no flow in the domain; the transport processes were diffusion and electromigration in porous media and electron drift in the iron bar. The Debye-Hückel equation was used as the activity coefficient model. The top boundary 92condition of the domain allowed continuous O2(g), H2(g) and CO2(g) exchange with atmosphere. The initial electric potential condition was set to 0 V for the entire domain. It has to be noted that here the electric potential is the self-potential which is measured by non-polarizing electrodes while redox potential is expressed with respect to the standard hydrogen electrode as a measure of tendency to accept electrons; therefore a more positive redox potential indicates the tendency to accept electrons. For the electric boundary conditions, the electric potential at the top right corner of the domain was fixed at a reference voltage of 0 V. The top boundary condition was determined by a redox potential representative of atmospheric conditions (i.e. equilibrium with atmosphere: Eh = 803 mV). pH is fixed at the top boundary to 7.0. Electric flux boundary conditions were applied to cathode and anode boundary cells at the top and bottom of the iron bar respectively. This boundary condition is adapted over time to account for changes in redox potential gradient. The results of the simulation were compared to experimental measurements to assess the capabilities of the model to capture the electrochemical evolution in the experiment including geochemical reactions. 93Table 4.2. Boundary conditions, initial conditions, and species-specific diffusion coefficients. Species *Initial Condition *Boundary Condition Unit †Diffusion Coefficient (m2 s-1) Aqueous phase pH 7. 0 7.0 - 9.31 × 10-9 Na+ 2.000 × 10-3 - mol L-1 1.33 × 10-9 Cl‒ 1.059 × 10-2 - mol L-1 2.03 × 10-9 K+ 1.000 × 10-2 - mol L-1 1.90 × 10-9 Fe+3 4.720 × 10-12 - mol L-1 1.33 × 10-9 Fe+2 2.530 × 10-13 - mol L-1 2.03 × 10-9 CO2‒3 1.950 × 10-3 - mol L-1 9.55 × 10-10 Eh 400 803 mV - Mineral phase Fe(s) 0.28 - m3 m-3 - Fougerite 0 - m3 m-3 - Ferrihydrite 0 - m3 m-3 - Lepidocrocite 0 - m3 m-3 - Siderite 0 - m3 m-3 - *Aqueous concentrations are in mol L-1. †From “Diffusion Coefficients in Liquids at Infinite Dilution” in [124] unless otherwise indicated. Table 4.3. Physical parameters. Parameter Symbol Value Unit Length × Depth L × D 2.0×0.5 m Spatial discretization Δx 1.0 × 10-3 m Control volume n 1184 - Simulation time T 42 day Maximum time step Δtmax 0.5 day Porosity - sand φp 0.34 [-] Porosity – iron bar φm 0.24 [-] †Electron mobility in iron μe 3.29 × 10-8 (m2 V-1 s-1) ††Electric conductivity of iron σm 450 (S m-1) †Calculated for iron with electric conductivity of 450 [S m-1]. ††Calibrated. 94Results and Discussion: two-dimensional simulations of 42 days length were carried out and compared to the experimental data of Castermant et al [12]. The model was calibrated against available pH, Eh and self-potential measurements. The model was also used to demonstrate the geochemical changes which likely have occurred in the sand box in response to the corrosion of the iron bar. Eh evolution, electric flux and self-potential signals: the occurrence of iron oxidation created a vertical upward electron flux through the iron bar. The oxidation of iron and the initial vertical natural redox gradient initiated the upward electron drift generating an electric current of 25.2 mA across the iron bar. The observations and simulation results show a decrease in redox potential indicative of a reducing environment around the iron bar (Figure 4.3). Oxygen diffusion from the surface created an oxidizing zone in the upper section of the tank. The redox conditions within the tank were affected by two processes: 1) Gas exchange with the atmosphere which determines the vertical distribution of redox conditions within the sand box at distances greater than 30 cm from the iron bar, exhibiting a transition from an oxidizing environment to a reducing environment; 2) the near field of the iron bar, up to a distance of 30 cm from the bar, was strongly disturbed by iron oxidation reactions (Figure 4.3). The initial redox potential at the upper boundary open to air in the entire tank before introducing the iron bar was 803 mV. After inserting the iron bar, electrons generated by iron oxidation drifted to the upper oxygen rich section and initiated oxygen reduction consuming the dissolved oxygen which created a lower redox potential surrounding the iron bar (Figure 4.3). The far-field measured redox potential ranged between 300 mV and -50 mV and between -250 mV and -120 mV in the vicinity of iron bar. The overall reducing nature of the redox potential in the far-field region of the tank was indicative of other unknown reactions than iron oxidation. The difference between computed and 95measured redox potential can possibly be attributed to organic matter degradation and/or interactions with mineral phases present in the sand tank which were not considered in the model due to the unknown nature of these reactions (Figure 4.3). In addition, it is a well-known fact that redox reactions in groundwater do not easily approach equilibrium conditions due to kinetic limitations, whereas the model simulations were based on the equilibrium assumption. Despite these limitations, the modeling adequately captures the redox zonation within the sandbox. Figure 4.3. Distribution of the redox potential: simulated (top) vs. measured (bottom) [12] 42 days after inserting the iron bar on the left hand side of the tank. Measured redox potential from [12]. Copyright © 2008 European Association of Geoscientists & Engineers. Reprinted with permission from John Wiley & Sons. 96As the source current associated with the electron flux through the iron bar increased over time, the resulting self-potential anomaly also increased significantly and reached its peak at the end of the simulation. The best fit for the self-potential signal at the end of the simulation was obtained for an electric conductivity of 450 S m-1 for the iron bar. The iron bar was coated gradually by a thin layer of rust with a measured electric conductivity of 2.5×10-3 S m-1 [12]. The low electric conductivity of the iron bar can be attributed to rust accumulation on its surface. The computed self-potential signal at the surface was -70 mV, and agreed reasonably well with the measured surface self-potential at the end of the experiment (Figure 4.4). By increasing the electrical conductivity, the self-potential signature magnitude also increased. The simulated self-potential distribution showed a negative potential anomaly located on the top of the iron bar and a positive potential anomaly at the bottom which is consistent with observations and the proposed dipolar nature of the geobattery concept (Figure 4.5). Figure 4.4. Self-potential distribution on the surface; line is the model and symbols are measured at the end of the experiment after 42 days. 97 Figure 4.5. Self-potential distribution after 42 days: simulated (top) vs. measured (bottom) [12]. Measured self-potential from [12]. Copyright © 2008 European Association of Geoscientists & Engineers. Reprinted with permission from John Wiley & Sons. Iron oxidation and oxygen reduction: oxygen was the electron acceptor in the system. There were two iron oxidation mechanisms in this system: one was the kinetic oxidation of zero-valent iron to ferrous iron and the other oxidation of ferrous iron to ferric iron. Electrons released from both reactions were transferred through the iron bar and enhanced oxygen reduction (available from atmospheric diffusion) in the shallower part of the tank (Table 4.1). This result also shows that oxygen availability plays an important role in driving the system. The reducing environment around the iron bar (Figure 4.3) can be attributed to oxygen depletion by iron corrosion reactions (Figure 4.6). The simulation results indicate that iron oxidation by oxygen reduction will be stopped if oxygen as electron acceptor depletes. Under steady state conditions, the transport rate of oxygen eventually should be equal to its consumption rate. In this simulation, a fully saturated 98environment was present, where oxygen is only available within the uppermost few centimeters of the tank. However, in variably saturated environments, if the metallic body intercepts the capillary fringe across the oxic/anoxic interface, then stronger iron oxidation would be expected at the anode in the reduced zone due to greater oxygen availability at the cathode. Figure 4.6. Oxygen partial pressure (atm) after 42 days: depletion of oxygen due to oxygen reduction reactions. Mineral formation and pH evolution: the oxidation of iron increases pH by donating electrons to drive oxygen reduction. Figure 4.7 compares the computed pH values obtained by numerical simulations and from experimental observations. The pH of the sandbox increased from 7.0 to 9.2 in close proximity of the iron bar during the experiment. The computed pH in the vicinity of the iron bar was 9.9. The oxidation reaction rate was calibrated to reproduce the observed pH values. The pH distribution depends on the oxidation reaction rate and the presence of oxygen. 99 Figure 4.7. pH distribution: simulated (top) vs. measured (bottom) [12] after 42 days after inserting the iron bar on the left hand side of the tank. Measured pH from [12]. Copyright © 2008 European Association of Geoscientists & Engineers. Reprinted with permission from John Wiley & Sons. The simulation results also showed that iron oxidation and pH increases lead to the precipitation of secondary minerals (Figure 4.8). For the simulations, minerals are assumed to precipitate as ferrous and ferric hydroxides. The model results indicated that all mineral precipitation took place in the direct vicinity of the iron bar confirming laboratory observations by Castermant et al. [49]. Mineral precipitation of the iron oxides emulates the formation of green rust (here represented as fougerite) observed in the experiment (Figure 4.8). 100 Figure 4.8. Secondary mineral formation in the direct vicinity of the iron bar after 42 days, shown as volume fractions. Ionic transport: In the absence of a hydraulic gradient, diffusion and electromigration were the main transport mechanisms. There was a downward transfer of oxygen and CO2 gas from the atmosphere entering the tank via diffusion. The oxygen diffusion was the driving force of the iron corrosion reaction in the upper part of the tank. Iron and the high pH front produced by corrosion diffused through the pore water. Ferric and ferrous iron transport was limited to a few centimeters surrounding the iron bar due to their small diffusion coefficients and rapid immobilization by precipitation into the solid phase as secondary minerals. However, a high pH front propagated through the domain farther from the iron bar. Simulated pH values are generally lower than measured pH (Figure 4.7). The simulation results showed less pH front movement with respect to iron bar compared to the laboratory measurements, but the overall pH distribution is well captured. The observed pH increase in the far field region cannot be attributed to iron oxidation reactions, because transport by diffusion and migration in 42 days cannot reach the far 101field region despite the high diffusion coefficient of OH-. The elevated pH level in the far field zone can be attributed to possible degradation of organic matter observed inside the sand tank, as discussed by Castermant et al [12]. This process was not simulated here. The local precipitation of secondary minerals over time reduced the reactivity of the iron bar by forming a coating of high electric resistivity around the bar surface. This reduced electric conductivity could decrease electron flux from bottom to the top and consequently might have inhibited iron oxidation. We did not consider the effect of increasing electric resistivity in the model. 4.6 Illustrative Ore Deposit Model In order to evaluate the code’s capabilities to simulate the geobattery mechanism associated with mineral deposits, an illustrative numerical model was constructed. The model and its parametrization include the key features of geobattery processes attributed to buried mineral deposits. Based on the conceptual model presented in Figure 4.1, a two dimensional numerical model was constructed for a 20 × 4 m domain. A pyrite ore body with 1 × 2.5 m dimension was assumed in the middle of the domain; the top and bottom of the ore body were 0.5 m and 1 m from upper and lower boundaries respectively. Although the dimensions of this illustrative model are not representative of field conditions, the model allows capturing the relevant mechanisms and demonstrates their interactions. The upper boundary was in contact with the atmosphere, providing a constant source of oxygen diffusing into the domain and the lower boundary was assumed anoxic with reducing conditions maintained by iron-bearing silicate minerals. An initial mineral distribution is assumed to represent the result of long-term weathering and is meant to represent present day conditions. The ore body was electrical conductive, allowing electron drift though the ore. For simplicity, it is assumed that there was no hydraulic gradient. Diffusion and electromigration 102were the main transport mechanisms for solutes and electron migration occurs in the electrical conductive zone controlled by the electric potential gradient between the two ends of the ore body. The domain was discretized into 1302 control volumes (62 in horizontal and 21 in vertical direction). The mesh around the ore body was refined to better capture the processes occurring at the interface. Geochemical Relationships: The following 12 components were considered in the simulations: H+, Ca+2, Mg+2, Al+3, K+, Fe+2, Fe+3, Cl-, CO3–2, HS-, SO4–2, H4SiO4 and H2O In addition, 3 gases (O2(g), H2(g), CO2(g)) and 7 minerals (pyrite [FeS2], K-mica [KAl3Si3O10(OH)2] and calcite [CaCO3] as primary minerals and gypsum [CaSO4.2H2O], iron hydroxide [Fe(OH)2(am)], ferrihydrite [Fe(OH)3(am)] and siderite [FeCO3(am)] as secondary minerals) were included. Fe+2/Fe+3 and SO4‒2/ HS‒ were considered as redox couples. Electron migration took place through the conductive ore body and the electron was also used as master variable to calculate redox potential. A total number of 10 aqueous complexation reactions from the MIN3P database were included in the model to adequately simulate mineral solubilities. The initial conditions for the three different zones are presented in Table 4.4. The ore body is assumed passivated and non-reactive from 1.5 m to 3.5 m elevation. The diffusion coefficients were taken from [124] for available ions or calculated based on the Stokes-Einstein equation [167] for other ions. The upper boundary was fixed at a redox potential to 808 mV to represent oxidizing conditions by contact with the atmosphere. As O2(aq) is a secondary species, it was not possible to assign oxygen flux directly as boundary condition and as an approximation its concentration was calculated based on redox potential. The initial mineral distribution is illustrated in Figure 4.9. The electric potential was set to 0 V from x = 0 to x = 20 m and z = 4 m (top surface) as reference points for electric potential field calculation. 103Table 4.4. Initial conditions for three different zones. Species Upper Domain Lower Domain Ore Body Units Aqueous phase pH 7. 17 7. 17 7. 17 - Cl‒ 1.000 × 10-2 1.000 × 10-2 1.000 × 10-2 mol L-1 CO2‒3 1.950 × 10-3 1.950 × 10-2 1.950 × 10-3 mol L-1 Ca+2 5.825 × 10-3 7.695 × 10-3 5.869 × 10-3 mol L-1 Fe+3 4.720 × 10-12 4.720 × 10-12 4.720 × 10-12 mol L-1 Fe+2 1.000 × 10-6 1.000 × 10-6 1.000 × 10-6 mol L-1 SO4‒2 1.000 × 10-10 1.000 × 10-10 1.000 × 10-10 mol L-1 HS‒ 1.000 × 10-10 1.000 × 10-10 1.000 × 10-10 mol L-1 Mg+2 1.000 × 10-10 1.000 × 10-8 1.000 × 10-10 mol L-1 K+ 7.000 × 10-9 7.000 × 10-9 7.000 × 10-9 mol L-1 Al+3 2.100 × 10-8 2.100 × 10-8 2.100 × 10-8 mol L-1 H4SiO4 2.100 × 10-8 2.100 × 10-8 2.100 × 10-8 mol L-1 Eh 400 75 400 mV Mineral phase Pyrite 0 0 1.0 m3 m-3 Calcite 1.0 0.5 0 m3 m-3 K-mica 0 0.5 0 m3 m-3 Fougerite 0 0 0 m3 m-3 Ferrihydrite 0 0 0 m3 m-3 Gypsum 0 0 0 m3 m-3 Siderite 0 0 0 m3 m-3 Table 4.5. Physical parameters. Parameter Symbol Value Unit Length × Depth L × D 20.0 × 4.0 m Control volume n 1304 - Maximum time step Δtmax 0.01 year Porosity – aquifer φp 0.38 [-] Porosity – ore body φm 0.38 [-] †Electron mobility in pyrite μe 3.765 × 10-8 (m2 V-1 s-1) †Calculated for pyrite with electric conductivity of 152 [S m-1] from [168]. 104Table 4.6. Reactions network and corresponding coefficients [166]. Name Reaction log K log keff mol L-1 H2O s-1 Primary minerals reactions Pyrite FeS2 + 8H2O ↔ Fe+2 + 2SO4–2 + 16H+ + 14e− 85.798 − 10.0 Calcite CaCO3 ↔ Ca+2 + CO3-2 8.480 − 8.0 K-mica KAl3Si3O10(OH)2 + 9H+ ↔ K+ + 3 Al+3 + 3 H4SiO4+ e− − 12.703 − 10.0 Secondary minerals reactions Fougerite Fe(OH)2(am) ↔ Fe+2 + 2H2O − 2H+ − 13.905 − 9.0 Ferrihydrite Fe(OH)3(am) ↔ Fe+3 + 3H2O − 3H+ − 4.891 − 9.0 Gypsum CaSO4.2H2O ↔ Ca+2 + SO4–2 + 2H2O 4.581 − 8.0 Siderite FeCO3(am) ↔ Fe+2 + CO3‒2 10.450 − 10.0 Gas dissolution-exsolution reactions O2(g) O2(aq) ↔ O2(g) 2.898 mol L-1 atm-1 H2(g) H2(aq) ↔ H2(g) − 43.009 mol L-1 atm-1 CO2(g) CO2(aq) ↔ CO2(g) 18.160 mol L-1 atm-1 Redox couples Fe2+ /Fe3+ Fe+3 + e− ↔ Fe+2 13.032 - SO4–2/HS- HS- + 4H2O ↔ SO4–2 + 9H+ + 8e− − 33.660 - Aqueous complexation OH- H2O ↔ H+ + OH- − 13.998 - O2(aq) O2(aq) + 4e− + 2H+ ↔ 2H2O − 86.018 - H2(aq) H2(aq) ↔ 2e− + 2H+ − 3.124 - HCO3‒ HCO3‒ ↔ CO3‒2 + H+ − 3.617 - H2CO3(aq) H2CO3(aq) ↔ CO3‒2 + 2H+ − 2.247 - HSO4‒ HSO4‒ ↔ SO4‒2 + H+ 3.850 - CaSO4(aq) CaSO4(aq) ↔ SO4‒2 + Ca+2 1.470 - CaOH+ CaOH+ + H+ ↔ Ca+2 + H2O 14.535 - CaHCO3+ CaHCO3+ ↔ Ca+2 + CO3‒2 + H+ 5.410 - CaCO3(aq) CaCO3(aq) ↔ Ca+2 + CO3‒2 4.030 - 105 Figure 4.9. Initial mineral distribution a) pyrite, b) calcite, and c) K-mica. 106Results and discussion: assuming present day mineral distribution, the model was run until aqueous concentrations, redox distribution and self-potential distribution did not show significant change, representing a quasi-steady-state solution. The initial aqueous concentrations were set to ensure global electroneutrality in the domain. The model was used to illustrate the major geochemical and physical changes which are likely to occur in a buried ore deposit. The ore body straddles two regions of reducing condition at depth and oxidizing condition in the shallower region. The redox gradient and pyrite oxidation at depth induced upward electron transfer through the ore body as well as migration of chloride toward the bottom to compensate excessive negative charge on top and maintain electroneutrality in the early stages before approaching steady state. As the process progresses toward steady state, the net movement of chloride diminishes and only a signature of chloride depletion and enrichment remains. At steady state, the Fickian diffusion compensates for electromigration and maintains electroneutrality. The redox potential within the ore is more reducing than in the surrounding aquifer near the top (Figure 4.10). The electrical conductivity of the ore body is the key parameter that enables electron transfer from the anode (oxidation of unweathered sulfides at depth) to the cathode (oxygen reduction near the surface). The passivation of reaction surfaces can lead to persistence of sulphides in regions with elevated redox potential, as was for example observed in the Kambalda nickel deposits [169]. This passivation can be due to deposition of goethite and ferric hydroxide in early stages on reactive surfaces of the shallower portion of the deposits [170]. Therefore, the upper 2m of the ore was assumed passivated due to weathering and coating with insoluble minerals and only the lower 0.5 m of the ore body was reactive and underwent oxidation releasing electrons into the conductive ore according to pyrite oxidative dissolution reaction (FeS2 + 8H2O → Fe+2 + 2SO4–2 + 16H+ + 14e−). As a result, electrons drifted upward 107along the redox potential gradient through the conducting ore body. Simultaneously, strong oxidizing conditions propagate to greater depth outside of the ore body due to the absence of reducing mineral phases in this region (Figure 4.11). The encounter of dissolved oxygen and electrons at top of the ore body leads to oxygen reduction according to the reaction (O2(aq) + 4e− + 2H+ → 2H2O) and resulted in oxygen depletion near the top of the ore (Figure 4.11). In this model, the upper boundary redox potential controls oxygen influx into the domain. Figure 4.10. Redox potential at quasi steady state. Electron drift through the ore body decreases the redox potential gradually. 108 Figure 4.11. Oxygen distribution at quasi steady state. The electrons released by pyrite oxidation are balanced by the electrons accepted by oxygen reduction under steady state conditions. The source current inside the conducting ore body (equation 4.4) is the result of this electron flux. In order to calculate the rate of electron transfer by pyrite, normalized pyrite oxidation (with respect to electrons) was integrated over the ore body and compared with the fraction of electron intake from the ore-aquifer interface, which was consumed by oxygen. The electric conductivity of the conductor was calibrated to balance the two rates for quasi-steady state conditions. Pyrite oxidation reduces pH in the deeper part creating acidic condition and oxygen reduction increases pH in the shallower part (Figure 4.12). However the presence of calcite in the system buffers the pH in the surrounding areas outside the ore body and maintains the pH near neutral pH values. 109 Figure 4.12. pH distribution at quasi steady state. pH decreases at depth due to pyrite oxidation and increases near the surface due to oxygen reduction. At quasi-steady state the rates of pyrite oxidation and oxygen reduction became equal. Pyrite dissolution led to conditions very close to equilibrium in the deeper region (Figure 4.13a and c). The protons generated by pyrite oxidation were consumed by calcite and K-mica dissolution outside the ore body (Figure 4.13b and d). Figure 4.13. (a) Pyrite dissolution rate, (b) K-mica dissolution rate, (c) pyrite saturation index and (d) calcite dissolution rate at the end of simulation (quasi steady state). 110The self-potential signal at quasi steady state is presented in Figure 4.14. According to Sato and Mooney [42], the maximum self-potential anomaly would be one-half of the total potential difference. The maximum self-potential anomaly for pyrite would be on the order of 356 mV; in our simulations a maximum anomaly of 300 mV was reached. Figure 4.14. Self-potential generated by pyrite oxidation at quasi steady state. Pyrite oxidation produced protons and led to dissolution of calcite; near the interface of the ore body, resulting in the depletion of carbonates and, increased CO2(aq) and increased ferrous iron concentration in the pore water. Dissolution of calcite and pyrite caused the formation of siderite 111at depth. Ferrous iron was transported upward, and was oxidized to ferric iron in the shallower part, where it precipitated as ferrihydrite. The maximum oxidation rate of the ore body that can be supported is limited by oxygen diffusion in water to the top of the ore body. In order to investigate the effect of oxygen influx on self-potential response, the oxygen influx rate was reduced by a factor of 3.5 by decreasing the redox potential at the upper boundary. The self-potential results at quasi steady state compared to the base case scenario are presented in Figure 4.15. Decreasing the oxygen influx rate decreases source current and consequently self-potential response. 112 Figure 4.15. Self-potential response for (a) base case scenario, (b) reduced oxygen influx by a factor of 3.5 by reducing upper boundary redox potential and (c) surface signal. 113In the absence of a hydraulic gradient, diffusion, electromigration and electron transfer were the main transport mechanisms. In this simulation pyrite oxidation released ferrous iron (Figure 4.16a) and sulfate (Figure 4.16b) at depth. Ferrous iron and carbonate from calcite dissolution (Figure 4.16d) led to precipitation of siderite (Figure 4.17a). Ferrous iron was transported upward by diffusion and electromigration (Figure 4.16a) and when encountering more oxidizing conditions in the shallower part it oxidized to ferric iron and eventually led to precipitation of ferrihydrite in the upper part (Figure 4.17b). The acidic environment surrounding the reactive ore body contributed to dissolution of K-mica in the vicinity of pyrite (Figure 4.13b) releasing potassium and aluminum (Figure 4.16c and e). The transport of potassium and aluminum around the ore body are dominantly driven by diffusion. Electromigration contributes to cation transport including ferrous iron, potassium and aluminum toward the cathode (Figure 4.16a, c and e) and anions such as sulfate and carbonate toward the anode (Figure 4.16b and d), but is not very significant. In case of silicic acid, since it has no electric charge, its transport is only by diffusion (Figure 4.16e). The effect of electromigration can be seen most clearly in the chloride concentration distribution (Figure 4.16h). Chloride unlike other ions did not participate in any heterogeneous and/or homogenous reaction and its initial concentration was uniformly distributed in the entire domain, therefore any disturbance to chloride concentration and consequently its transport can only be attributed to electromigration to maintain electroneutrality in the system. When electrons drift to the upper part of the ore, a negative net charge in the shallower part and a positive net charge in the deeper part starts to build up. To offset the charge imbalance, chloride migrated toward the anode and its concentration around the cathode decreased as its concentration around the anode increased. Chloride migration is the direct result of electric coupling and the electric field is the 114driving force for its migration. It should be noted that at the steady state only signatures of chloride depletion and enrichment remains and electromigration will be compensated by diffusion. Figure 4.16. Generation and transport of ions: (a) ferrous iron, (b) sulfate, (c) potassium, (d) carbonate, (e) aluminum, (f) silicic acid, (g) redox potential and (h) chloride. Figure 4.17. Precipitation of secondary minerals: (a) siderite precipitation rate and (b) ferrihydrite precipitation rate. 115The simulations are able to reproduce the major mechanisms associated with the geobattery conceptual model shown in Figure 4.1 and results are consistent with the observations from a field study involving the weathering of a pyritic quartz-feldspar porphyry deposit. This field study also reported a correlation between self-potential and redox potential response [171]. 4.7 Summary and Conclusions An electrochemical reactive transport model was developed to incorporate generation of surface-measurable self-potential signals associated with redox reactions (e.g. metal or sulfide oxidation and oxygen reduction) based on the Nernst-Planck and Gauss-Ampère equations and including a full suite of geochemical reactions. The model was constructed according to the geobattery conceptual model developed by Sato and Mooney [42], Sivenas and Beales [160] and Revil et al. [43] where electron transfer through an electrical conductive medium creates a source current in a metallic body and leads to self-potential response at the surface. To account for electron flux through conducting media, the MIN3P code was enhanced by implementing electron drift. The associated electric current by electron transfer (source current) was applied as an electric flux boundary conditions at anode-aquifer and cathode-aquifer interfaces. This assumption required a simplification on pre-defining the location of anode and cathode. Reactive transport modelling was performed to describe geochemical and electrochemical processes occurring during a controlled sand tank experiment involving the oxidation of an iron bar conducted by Castermant et al. [12]. This study served as a validation of the model’s ability to simulate the geobattery concept. The simulation results demonstrated that 116kinetic iron oxidation along with oxygen reduction were controlling the redox, pH and water chemistry. Electron transfer towards the surface where oxygen is present in response to iron oxidation at depth created a source current along the redox potential gradient and produced the self-potential signals. The simulation results indicated that metal oxidation reactions can lead to kinetically controlled formation of secondary minerals with higher electric resistivity in vicinity of an iron bar confirming laboratory observations. The study also demonstrated how the connection between regions of different redox potentials can be established by a metallic body [152], [156], [172]–[175]. The electrochemical signature of a hypothetical pyritic ore deposit and its ambient environment was investigated to illustrate the simulation capabilities of the code in capturing geochemical reactions and self-potential response. To produce surface measurable self-potential signals the electrical conductive media should intercept the regions of strong redox potentials (oxic/anoxic interface). To represent present day conditions of the deposit, the shallower part of the ore was assumed passivated due to prior deposition of inert minerals on reactive surfaces and only the deeper portion of the ore body was considered geochemically reactive. The model was run to reach quasi-steady state conditions. The rate of electron transfer through the conducting ore body is controlled by its electric conductivity; decreasing electric conductivity reduced self-potential response. Overall, it could be shown that the model was able to produce key elements of the geobattery conceptual model in a process-based fashion. The simulation results were also 117consistent with field observations at a sulphide mineral deposit, including self-potential response and geochemical alterations. Generally, despite simplifying assumptions, the simulation results agreed reasonably well with available laboratory and field observations and principles of the geobattery conceptual model. The geochemical reactions could only be investigated hypothetically and qualitatively due to limited field data; more quantitative assessments of geobattery systems using reactive transport modelling, requires detailed observations and measurements of water chemistry and mineralogical studies. Lab scale experiments in controlled environment are very crucial to obtain comprehensive and reliable data for further studies of geobattery systems. Further model development is necessary to account for early stages of system development, when the upper portion of the deposit is not yet passivated. However, because of the strong coupling between geochemical reactions and the self-potential anomaly, the developed model can be used to investigate the geochemical process associated with surface-measurable self-potential signals. The model helps to understand the interactions between geochemical reactions and self-potential signatures. Possible applications for using this novel approach include contaminant plume migration and degradation, and identification of buried ore deposits for mineral exploration. 118Chapter 5: Summary and Conclusions 5.1 Model Development Summary The MIN3P reactive transport code [105] was used as a platform to implement new enhancements to facilitate the simulation of reactive transport problems characterized by electric coupling with mass transfer. The three major enhancements implemented in MIN3P are: 1. Implementation of the Nernst-Planck equation using the null current assumption to account for electric coupling for multicomponent diffusion problems, 2. Full coupling of mass flux and electric flux continuity equations i.e. Nernst-Planck and Gauss-Ampère equations and electrode boundary conditions for electrokinetic remediation problems, 3. Implementation of electron transfer in an electrical conductor and consideration of source current term as the driving force behind the geobattery mechanism. Species-dependent diffusion and electromigration have applications in early diagenesis [2], [24], [64], radioactive waste disposal repositories [31], [67], [68] and CO2 sequestration [65], [69]. Implementation of electromigration based on the null current assumption and the Nernst-Planck equation enhanced MIN3P’s capabilities for these applications. Three benchmark problems specifically designed to be sensitive to electric coupling were considered, each with significant effects of electromigration on diffusion and lateral concentration displacement. The implementation and formulation in MIN3P were tested based on these benchmark problems via an inter-comparison with other reactive transport codes (CrunchFlow and PHREEQC). Two of the benchmark problems were hypothetical problems and one of them was based on laboratory 119flow-through experiments [78]. The benchmark problems contribute to the verification of the implementation of multicomponent diffusion and electromigration in all three reactive transport codes. Overall very good agreement between codes was achieved. Nearly identical concentration distributions were produced by all codes, in the presence of electrostatic interactions between the charged species. Miniscule residual discrepancies between results were indicative of different implementation schemes and were not considered a significant source of uncertainty. The benchmark problems and various applications [2]–[4], [65] have shown the importance of electromigration and electric coupling for reactive transport problems under some circumstances. This thesis has shown that full coupling of mass transport and electric flux is necessary for development of a general purpose reactive transport code incorporating electrochemical processes. For application of reactive transport simulation in electrokinetic remediation technologies, the Nernst-Planck and Gauss-Ampère equations were coupled and electrode boundary conditions were implemented. The completed model was shown to be able to simulate interactions between electrode reactions, gas generation and consumption, in electrokinetic systems. The electric coupling was tested against three examples. The first verification was against a hypothetical free diffusion example in the absence of an external electric field, as introduced by Paz-García et al. [20]. These authors used the finite element method to solve the Poisson-Nernst-Planck system of equations numerically. A reasonable agreement between Gauss-Ampère and Poisson-Nernst-Planck approaches was achieved in the inter-comparison. Another example was a pH-isolated permanganate transport experiment [60] in the presence of an induced external electric field. The pH-isolation prevented alkaline and acidic fronts from advancing into the domain. A comparison between MIN3P, PHT3D [39] and experimental results showed a satisfactory agreement. The discrepancies between two codes can 120be related to the differences in modeling approaches as PHT3D considered only MnO4- but MIN3P considered coupled diffusion of K+ and MnO4-. On the other hand, PHT3D also accounted for electroosmotic driven flow, which was not included in the MIN3P simulations. Moreover, MIN3P and PHT3D use different numerical solution strategies. Finally the third verification example was the electrokinetic desalination of brick. Good agreement with experimental results was obtained for the desalination problem. The desalination model was extended by including gas generation, exsolution and consumption, as well as electrode reactions. The simulations demonstrated that total gas pressure build-up by hydrogen and oxygen generation at the electrodes are controlled by degassing, as well as combined hydrogen oxidation and oxygen reduction. Scenario and sensitivity analysis were used to assess the role of electric potential, degassing and hydrogen oxidation on build-up of total gas pressure. This analysis showed that gas accumulation and adsorption can decrease remediation efficiency when degassing and hydrogen oxidation rates are outpaced by gas generation. The new model can be used to assess and design electrokinetic remediation technologies especially since it benefits from an existing complex geochemical calculation framework already in place and suitable to assess the performance of electrokinetic remediation. Based on geobattery conceptual models, surface measurable self-potential signals are typically generated when an electric source current establishes in an electrical conductor as a result of electron drift between regions of redox potential differences. MIN3P was enhanced by incorporating electron drift through a metallic conductive body to account for electric source current. At anode-aquifer and cathode-aquifer interfaces, electric flux boundary conditions equal to source current was applied. Metal oxidation and oxygen reduction are major redox reactions that drive the geobattery system. For concept validation an iron oxidation experiment in a 121controlled sand box was adopted. Simulation results for transient conditions demonstrated reasonable agreement with experimental observations and reproduced key mechanisms of the geobattery concept. Iron oxidation and oxygen reduction reactions were controlling ambient water chemistry and redox potential distribution. Increasing source current associated with electron flux through the iron bar increased the self-potential signal until reaching its maximum value at the end of the simulation. The dipolar nature of the geobattery was demonstrated by a negative potential anomaly on top and a positive potential anomaly at depth. The model also adequately captured the distribution of redox conditions despite some differences in absolute values, which can be attributed to kinetic limitations of redox reactions and measurements. An illustrative example was developed in order to evaluate the code’s capability in capturing geobattery model behavior of a buried ore body under long-term quasi-steady state conditions. The upper part of the ore was assumed passivated due to early time weathering and more stable mineral precipitation. Pyrite oxidation at depth provided the electrons required for oxygen reduction at shallow depth. The electron transfer rate through the conducting ore body was controlled by its electric conductivity. Redox reactions initiated subsequent mineral dissolution reactions and precipitation of secondary minerals. The generated electric field and electron drift also caused electromigration of chloride from top to bottom to maintain overall electroneutrality, leaving a lasting signature under quasi-steady state conditions. Simulation results of an illustrative example were compared qualitatively with a hydrogeochemical and geophysical study of oxidation of a massive sulfide mineral deposit. Observations indicated association of self-potential signals with pyrite oxidation and oxygen reduction. pH increased in the shallower part of the deposit and increased acidity in the deeper part leading to subsequent carbonate 122mineral dissolution and secondary carbonate deposition. The modeling results were overall consistent with these observations. 5.2 Limitations, Future Works and Contributions Despite implementation of full electric coupling, the model has certain limitations due to simplifying assumptions; there are three major simplifications that could be addressed in the future: 1) advective flux was not considered in the geobattery model; therefore the model is not applicable for scenarios which involve advection and its contribution to self-potential signal known as streaming potential. 2) The location of anode-aquifer and cathode-aquifer interfaces are user defined and it was assumed that electron transfer only occurs at the top-end and bottom-end of the ore body; consequently, internal electric flux boundary conditions were defined at these interfaces. This limits the applicability of the code for early time ore deposit evolution when oxidation starts from the shallow part and advances through the ore body, passivating the upper part by secondary mineral depositions. For the same reason, the simulations assumed that the present day ore deposit contains a passivated upper portion of the ore body, representing weathering over geologic time. 3) To balance the rate of electron transfer from pyrite to the reduction of oxygen, the electric conductivity of the ore body was calibrated. The Nernst-Planck equation with null current assumption used for multicomponent diffusion and full coupling of Nernst-Planck and Gauss-Ampère equations are general propose formulations coupled with geochemical reactions. This feature makes the new model applicable for a wide variety of problems like nuclear waste repositories, CO2 sequestration, soil remediation by electrokinetic technology and self-potential signatures associated with buried ore bodies or contaminant plumes. 123Another potential addition to the multicomponent diffusion model and the Nernst-Planck and Gauss-Ampère coupling includes modelling the diffuse double layer to account for surface complexation, cation surface migration and anion exclusion in low permeability media. The new model containing the fully coupled Nernst-Planck and Gauss-Ampère equations incorporates induced external electric current as applied in electrochemical remediation techniques. Currently there is no other model available to simulate gas phase evolution on electrodes in electrokinetic electrochemical remediation scenarios. The capability of the model to simulate gas generation, exsolution and consumption on electrodes can be utilized to aid in the understanding of electrode processes and improve the design and efficiency of remediation strategies. However, more experimental scale studies are required to obtain reliable data on gas phase dynamics on the electrodes and its effect on technology efficiency. Finally, the geobattery model formulation is capable of simulating the self-potential signals and is the first geochemical/reactive transport-based model for assessing the geobattery concept. The previous models by Stoll [10] and Bigalke and Grabner [11] and models based on Revil et al’s [43] formulation neglects geochemical processes that affect the redox regime, while focusing on the representation of the electric flux in the conductor. The model developed in this thesis includes full geochemical coupling and has the potential to contribute significantly to the understanding of generation of self-potential anomalies related to mineral exploration studies. Detailed quantitative study of geobattery systems requires laboratory scale experiments in controlled environments with a focus on water and mineral chemistry analysis to collect reliable data. The model presented in this project can be used for future developments; MIN3P provides a suitable development platform and it already offers various options for flow (e.g. 124saturated, unsaturated and density-driven) and geochemistry (Pitzer model and diffuse double layer model in batch mode). 125Bibliography [1] E. L. Cussler, Diffusion: Mass Transfer in Fluid Systems, 2nd ed. Cambridge: Cambridge University Press, 2003. [2] A. C. Lasaga, “Treatment of multicomponent diffusion and ion-pairs in diagenetic fluxes,” Am. J. Sci., vol. 279, no. 3, pp. 324–346, 1979. [3] M. A. Glaus, M. Birgersson, O. Karnland, and L. R. Van Loon, “Seeming steady-state uphill diffusion of 22Na+ in compacted montmorillonite,” Environ. Sci. Technol., vol. 47, no. 20, pp. 11522–11527, 2013. [4] E. H. Oelkers, “Physical and chemical properties of rocks and fluids for chemical mass transport calculations,” Rev. Mineral. Geochemistry, vol. 34, no. 1, pp. 131–191, 1996. [5] J. Newman and K. E. Thomas-Alyea, Electrochemical Systems, 3rd ed. New Jersey: John Wiley & Sons, Inc., 2012. [6] V. S. Bagotsky, Fundamentals of Electrochemistry, 2nd ed. New Jersey: John Wiley and Sons, 2005. [7] C. A. J. Appelo and D. Postma, Geochemistry, Groundwater and Pollution, 2nd ed. Leiden: A.A. Balkema Publishers, 2005. [8] W. M. Telford, L. P. Geldart, and R. E. Sheriff, Applied Geophysics, 2nd ed. Cambridge: Cambridge University Press, 1990. [9] K. H. Williams, S. S. Hubbard, and J. F. Banfield, “Galvanic interpretation of self-potential signals associated with microbial sulfate-reduction,” J. Geophys. Res., vol. 112, p. G03109, Sep. 2007. [10] J. Stoll, J. Bigalke, and E. W. Grabner, “Electrochemical modeling of self-potential anomalies,” Surv. Geophys., vol. 16, no. 1, pp. 107–120, 1995. 126[11] J. Bigalke and E. W. Grabner, “The Geobattery model: a contribution to large scale electrochemistry,” Electrochim. Acta, vol. 42, no. 23–24, pp. 3443–3452, 1997. [12] J. Castermant, C. A. Mendonca, A. Revil, F. Trolard, G. Bourrie, and N. Linde, “Redox potential distribution inferred from self-potential measurements associated with the corrosion of a burden metallic body,” Geophys. Prospect., vol. 56, no. 2, pp. 269–282, 2008. [13] T. Xu, E. Sonnenthal, N. Spycher, and L. Zheng, “TOUGHREACT V3.0-OMP Reference Manual: A Parallel Simulation Program for Non-Isothermal Multiphase Geochemical Reactive Transport,” Berkeley, 2014. [14] C. I. Steefel, C. A. J. Appelo, B. Arora, D. Jacques, T. Kalbacher, O. Kolditz, V. Lagneau, P. C. Lichtner, K. U. Mayer, J. C. L. Meeussen, S. Molins, D. Moulton, H. Shao, J. Šimunek, N. Spycher, S. B. Yabusaki, and G. T. Yeh, “Reactive transport codes for subsurface environmental simulation,” Comput. Geosci., vol. 19, no. 3, pp. 445–478, 2015. [15] D. L. Parkhurst and C. A. J. Appelo, “User’s Guide To PHREEQC (version 2) — a Computer Program for Speciation, and Inverse Geochemical Calculations,” Denver, 1999. [16] E. Samson and J. Marchand, “Modeling the transport of ions in unsaturated cement-based materials,” Comput. Struct., vol. 85, no. 23–24, pp. 1740–1756, 2007. [17] F. Helfferich, Ion Exchange, 2nd ed. New York: McGraw-Hill, 1962. [18] A. C. Lasaga, Kinetic Theory in the Earth Sciences. Priceton: Princeton University Press, 2014. [19] J. M. Paz-García, B. Johannesson, L. M. Ottosen, A. N. Alshawabkeh, A. B. Ribeiro, and J. M. Rodríguez-Maroto, “Modeling of electrokinetic desalination of bricks,” Electrochim. 127Acta, vol. 86, pp. 213–222, 2012. [20] J. M. Paz-García, B. Johannesson, L. M. Ottosen, A. B. Ribeiro, and J. M. Rodríguez-Maroto, “Modeling of electrokinetic processes by finite element integration of the Nernst-Planck-Poisson system of equations,” Sep. Purif. Technol., vol. 79, no. 2, pp. 183–192, 2011. [21] A. B. Ribeiro and J. M. Rodríguez-Maroto, “Electroremediation of Heavy Metal Contaminated Soils -Processes and Applications,” in Trace Elements in the Environment. Biogeochemistry, Biotechnology, and Bioremediation, M. N. V. Prasad, K. S. Sajwan, and R. Naidu, Eds. Florida: Taylor & Francis, CRC Press, 2006, pp. 341–368. [22] J. Virkutyte, M. Sillanpää, and P. Latostenmaa, “Electrokinetic soil remediation - Critical overview,” Sci. Total Environ., vol. 289, no. 1–3, pp. 97–121, 2002. [23] J. G. Sah and J. Y. Chen, “Study of the electrokinetic process on Cd and Pb spiked soils,” J. Hazard. Mater., vol. 58, pp. 301–315, 1998. [24] R. E. McDuff and R. E. Ellis, “Determining diffusion coefficients in marine sediments: a laboratory study of the validity of resistivity techniques,” Am. J. Sci., vol. 279, no. 6, pp. 666–675, 1979. [25] S. Ben-Yaakov, “Diffusion of seawater ions - significance and consequences of cross coupling effects,” Am. J. Sci., vol. 281, no. 7, pp. 974–980, 1981. [26] Y. Wang and P. Van Cappellen, “A multicomponent reactive transport model of early diagenesis: application to redox cycling in coastal marine sediments,” Geochim. Cosmochim. Acta, vol. 60, no. 16, pp. 2993–3014, 1996. [27] O. Truc, J. P. Ollivier, and L. O. Nilsson, “Numerical simulation of multi-species transport through saturated concrete during a migration test - MsDiff code,” Cem. Concr. 128Res., vol. 30, no. 10, pp. 1581–1592, 2000. [28] E. R. Giambalvo, C. I. Steefel, A. T. Fisher, N. D. Rosenberg, and C. G. Wheat, “Effect of fluid-sediment reaction on hydrothermal fluxes of major elements, eastern flank of the Juan de Fuca Ridge,” Geochim. Cosmochim. Acta, vol. 66, no. 10, pp. 1739–1757, 2002. [29] S. Lorente, M. Carcasses, and J. P. Ollivier, “Penetration of ionic species into saturated porous media: the case of concrete,” Int. J. Energy Res., vol. 27, no. 10, pp. 907–917, 2003. [30] B. P. Boudreau, F. J. R. Meysman, and J. J. Middelburg, “Multicomponent ionic diffusion in porewaters: Coulombic effects revisited,” Earth Planet. Sci. Lett., vol. 222, no. 2, pp. 653–666, 2004. [31] C. A. J. Appelo and P. Wersin, “Multicomponent diffusion modeling in clay systems with application to the diffusion of tritium, iodide, and sodium in opalinus clay,” Environ. Sci. Technol., vol. 41, no. 14, pp. 5002–5007, 2007. [32] P. C. Lichtner, “Continuum formulation of multicomponent - multiphase reactive transport,” in Reactive Transport in Porous Media, P. C. Lichtner, C. I. Steefel, and E. H. Oelkers, Eds. Washington DC: Mineralogical Society of America, 1996, p. 438. [33] A. D. MacGillivray and D. Hare, “Applicability of Goldman’s constant field assumption to biological systems,” J. Theor. Biol., vol. 25, pp. 113–126., 1969. [34] E. Samson, J. Marchand, J. L. Robert, and J. P. Bournazel, “Modelling ion diffusion mechanisms in porous media,” Int. J. Numer. Methods Eng., vol. 46, no. 12, pp. 2043–2060, 1999. [35] J. M. Galíndez and J. Molinero, “On the relevance of electrochemical diffusion for the modeling of degradation of cementitious materials,” Cem. Concr. Compos., vol. 32, no. 5, 129pp. 351–359, 2010. [36] L. S. Bennethum and J. H. Cushman, “Multicomponent, multiphase thermodynamics of swelling porous media with electroquasistatics: I. macroscale field equations,” Transp. Porous Media, vol. 47, no. 3, pp. 309–336, 2002. [37] E. D. Mattson, R. S. Bowman, and E. R. Lindgren, “Electrokinetic ion transport through unsaturated soil: 1. Theory, model development, and testing,” J. Contam. Hydrol., vol. 54, pp. 99–120, 2002. [38] A. N. Alshawabkeh, “Electrokinetic remediation: II. theoretical model,” J. Geotech. Eng., vol. 122, no. 3, pp. 186–196, 1996. [39] M. Z. Wu, D. A. Reynolds, H. Prommer, A. Fourie, and D. G. Thomas, “Numerical evaluation of voltage gradient constraints on electrokinetic injection of amendments,” Adv. Water Resour., vol. 38, pp. 60–69, 2012. [40] M. Z. Wu, D. A. Reynolds, A. Fourie, H. Prommer, and D. G. Thomas, “Electrokinetic in situ oxidation remediation: Assessment of parameter sensitivities and the influence of aquifer heterogeneity on remediation efficiency,” J. Contam. Hydrol., vol. 136–137, pp. 72–85, 2012. [41] M. T. Harris, D. W. DePaoli, and M. R. Ally, “Modeling the electrokinetic decontamination of concrete,” Sep. Sci. Technol., vol. 32, no. 1–4, pp. 827–848, 1997. [42] M. Sato and H. M. Mooney, “The electrochemical mechanism of sulfide self-potentials,” Geophysics, vol. 25, no. 1, pp. 226–249, 1960. [43] A. Revil, C. A. Mendonca, E. A. Atekwana, B. Kulessa, S. S. Hubbard, and K. J. Bohlen, “Understanding biogeobatteries: Where geophysics meets microbiology,” J. Geophys. Res., vol. 115, p. G00G02, Feb. 2010. 130[44] C. A. Mendonca, “Forward and inverse self-potential modeling in mineral exploration,” Geophysics, vol. 73, no. 1, pp. F33–F43, Jan. 2008. [45] S. M. Hamilton and K. H. Hattori, “Spontaneous potential and redox responses over a forest ring,” Geophysics, vol. 73, no. 3, pp. B67–B75, 2008. [46] T. Arora, N. Linde, A. Revil, and J. Castermant, “Non-intrusive characterization of the redox potential of landfill leachate plumes from self-potential data,” J. Contam. Hydrol., vol. 92, pp. 274–292, 2007. [47] V. Naudet, A. Revil, and J. Y. Bottero, “Relationship between self-potential (SP) signals and redox conditions in contaminated groundwater,” Geophys. Res. Lett., vol. 30, no. 21, 2003. [48] V. Naudet and A. Revil, “A sandbox experiment to investigate bacteria-mediated redox processes on self-potential signals,” Geophys. Res. Lett., vol. 32, no. 11, Jun. 2005. [49] A. Revil, F. Trolard, G. Bourrié, J. Castermant, A. Jardani, and C. a. Mendonça, “Ionic contribution to the self-potential signals associated with a redox front,” J. Contam. Hydrol., vol. 109, pp. 27–39, 2009. [50] J. B. Rittgers, A. Revil, M. Karaoulis, M. A. Mooney, L. D. Slater, and E. A. Atekwana, “Self-potential signals generated by the corrosion of buried metallic objects with application to contaminant plumes,” Geophysics, vol. 78, no. 5, pp. EN65–EN82, Sep. 2013. [51] C. I. Steefel, S. B. Yabusaki, and K. U. Mayer, “Reactive transport benchmarks for subsurface environmental simulation,” Comput. Geosci., vol. 19, no. 3, pp. 439–443, 2015. [52] A. E. James, J. D. Stillman, and D. J. A. Williams, “Finite element solution of the 131equations governing the flow of electrolyte in charged microporous membranes,” Int. J., vol. 20, pp. 1163–1178, 1995. [53] M. Kato, “Numerical analysis of the Nernst-Planck-Poisson system,” J. Theor. Biol., vol. 177, no. 3, pp. 299–304, 1995. [54] E. Samson and J. Marchand, “Numerical solution of the extended Nernst-Planck model,” J. Colloid Interface Sci., vol. 215, no. 1, pp. 1–8, 1999. [55] E. Samson, J. Marchand, K. A. Snyder, and J. J. Beaudoin, “Modeling ion and fluid transport in unsaturated cement systems in isothermal conditions,” Cem. Concr. Res., vol. 35, no. 1, pp. 141–153, 2005. [56] S. Shiba, Y. Hirata, and T. Seno, “Mathematical model for hydraulically aided electrokinetic remediation of aquifer and removal of nonanionic copper,” Eng. Geol., vol. 77, pp. 305–315, 2005. [57] B. Narasimhan and R. Sri Ranjan, “Electrokinetic barrier to prevent subsurface contaminant migration: theoretical model development and validation,” J. Contam. Hydrol., vol. 42, no. 1, pp. 1–17, Mar. 2000. [58] J. Kubo, S. Sawada, C. L. Page, and M. M. Page, “Electrochemical injection of organic corrosion inhibitors into carbonated cementitious materials: Part 2. Mathematical modelling,” Corros. Sci., vol. 49, no. 3, pp. 1205–1227, Mar. 2007. [59] K. U. Mayer, E. O. Frind, and D. W. Blowes, “Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions,” Water Resour. Res., vol. 38, no. 9, p. 1174, 2002. [60] D. Hodges, “Permanganate electromigration in low permeability media,” The University of Western Australia, 2010. 132[61] K. Kamran, M. van Soestbergen, H. P. Huinink, and L. Pel, “Inhibition of electrokinetic ion transport in porous materials due to potential drops induced by electrolysis,” Electrochim. Acta, vol. 78, pp. 229–235, 2012. [62] A. J. Bard and L. R. Faulkner, Electrochemical Methods: Fundamentals and Applications, 2nd ed. John Wiley & Sons, 2001. [63] J. R. Vinograd and J. W. McBain, “Diffusion of electrolytes and of the ions in their mixtures,” J. Am. Chem. Soc., vol. 63, no. 7, pp. 2008–2015, Jul. 1941. [64] P. Van Cappellen and J. F. Gaillard, “Biogeochemical dynamics in aquatic sediments,” in Reactive Transport In Porous Media, vol. 34, P. C. Lichtner, C. I. Steefel, and E. H. Oelkers, Eds. Mineralogical Society of America, 1996, p. 438. [65] S. Molins, D. Trebotich, C. I. Steefel, and C. Shen, “An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation,” Water Resour. Res., vol. 48, p. W03527, 2012. [66] H. J. V. Tyrrell, Diffusion and Heat Flow in Liquids. London: Butterworths, 1961. [67] P. Alt-Epping, C. Tournassat, P. Rasouli, C. I. Steefel, K. U. Mayer, A. Jenni, U. Mäder, S. S. Sengor, and R. Fernández, “Benchmark reactive transport simulations of a column experiment in compacted bentonite with multispecies diffusion and explicit treatment of electrostatic effects,” Comput. Geosci., vol. 19, no. 3, pp. 535–550, 2015. [68] C. A. J. Appelo, L. R. Van Loon, and P. Wersin, “Multicomponent diffusion of a suite of tracers (HTO, Cl, Br, I, Na, Sr, Cs) in a single sample of Opalinus Clay,” Geochim. Cosmochim. Acta, vol. 74, no. 4, pp. 1201–1219, Feb. 2010. [69] S. Ovaysi and M. Piri, “Pore-scale dissolution of CO2+SO2 in deep saline aquifers,” Int. J. Greenh. Gas Control, vol. 15, pp. 119–133, 2013. 133[70] M. Rolle, M. Muniruzzaman, C. M. Haberer, and P. Grathwohl, “Coulombic effects in advection-dominated transport of electrolytes in porous media: Multicomponent ionic dispersion,” Geochim. Cosmochim. Acta, vol. 120, pp. 195–205, 2013. [71] J. Carrera, X. Sánchez-Vila, I. Benet, A. Medina, G. Galarza, and J. Guimerà, “On matrix diffusion: formulations, solution methods and qualitative effects,” Hydrogeol. J., vol. 6, no. 1, pp. 178–190, 1998. [72] G. Chiogna, O. A. Cirpka, P. Grathwohl, and M. Rolle, “Relevance of local compound-specific transverse dispersion for conservative and reactive mixing in heterogeneous porous media,” Water Resour. Res., vol. 47, no. 7, p. W07540, 2011. [73] E. M. LaBolle and G. E. Fogg, “Role of molecular diffusion in contaminant migration and recovery in an alluvial aquifer system,” Transp. Porous Media, vol. 42, no. 1, pp. 155–179, 2001. [74] K. T. B. MacQuarrie and K. U. Mayer, “Reactive transport modeling in fractured rock: A state-of-the-science review,” Earth-Science Rev., vol. 72, no. 3–4, pp. 189–227, 2005. [75] C. I. Steefel, S. Carroll, P. Zhao, and S. Roberts, “Cesium migration in Hanford sediment: a multisite cation exchange model based on laboratory transport experiments,” J. Contam. Hydrol., vol. 67, pp. 219–246, Dec. 2003. [76] Q. Kang, P. C. Lichtner, and D. Zhang, “Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media,” J. Geophys. Res. Solid Earth, vol. 111, p. B05203, May 2006. [77] B. Johannesson, K. Yamada, L.-O. Nilsson, and Y. Hosokawa, “Multi-species ionic diffusion in concrete with account to interaction between ions in the pore solution and the cement hydrates,” Mater. Struct., vol. 40, no. 7, pp. 651–665, 2007. 134[78] M. Muniruzzaman, C. M. Haberer, P. Grathwohl, and M. Rolle, “Multicomponent ionic dispersion during transport of electrolytes in heterogeneous porous media: Experiments and model-based interpretation,” Geochim. Cosmochim. Acta, vol. 141, pp. 656–669, 2014. [79] C. I. Steefel and K. Maher, “Fluid-Rock interaction: A reactive transport approach,” Rev. Mineral. Geochemistry, vol. 70, no. 1, pp. 485–532, 2009. [80] C. A. J. Appelo, “Multicomponent diffusion in clays,” in Water pollution in natural porous media at different scales: assessment of fate, impact and indicators (WAPO2), 2007, pp. 3–13. [81] R. Taylor and R. Krishna, Multicomponent Mass Transfer. John Wiley & Sons, Ltd, 1993. [82] P. C. Lichtner, “Principles and practice of reactive transport modeling,” MRS Proc., vol. 353, pp. 117–130, 1994. [83] M. Rolle, D. Hochstetler, G. Chiogna, P. K. Kitanidis, and P. Grathwohl, “Experimental investigation and pore-scale modeling interpretation of compound-specific transverse dispersion in porous media,” Transp. Porous Media, vol. 93, no. 3, pp. 347–362, 2012. [84] D. L. Hochstetler, M. Rolle, G. Chiogna, C. M. Haberer, P. Grathwohl, and P. K. Kitanidis, “Effects of compound-specific transverse mixing on steady-state reactive plumes: insights from pore-scale simulations and Darcy-scale experiments,” Adv. Water Resour., vol. 54, pp. 1–10, Apr. 2013. [85] M. Rolle and P. K. Kitanidis, “Effects of compound-specific dilution on transient transport and solute breakthrough: A pore-scale analysis,” Adv. Water Resour., vol. 71, pp. 186–199, Sep. 2014. [86] Y. B. Acar and A. N. Alshawabkeh, “Principles of electrokinetic remediation,” Environ. 135Sci. Technol., vol. 27, no. 13, pp. 2638–2647, Dec. 1993. [87] K. Kamran, M. van Soestbergen, and L. Pel, “Electrokinetic salt removal from porous building materials using ion exchange membranes,” Transp. Porous Media, vol. 96, no. 2, pp. 221–235, 2013. [88] K. J. Kim, J. M. Cho, K. Baek, J. S. Yang, and S. H. Ko, “Electrokinetic removal of chloride and sodium from tidelands,” J. Appl. Electrochem., vol. 40, no. 6, pp. 1139–1144, 2010. [89] P. Lu, Q. Feng, Q. Meng, and T. Yuan, “Electrokinetic remediation of chromium- and cadmium-contaminated soil from abandoned industrial site,” Sep. Purif. Technol., vol. 98, pp. 216–220, 2012. [90] K. Baek, D. H. Kim, S. W. Park, B. G. Ryu, T. Bajargal, and J. S. Yang, “Electrolyte conditioning-enhanced electrokinetic remediation of arsenic-contaminated mine tailing,” J. Hazard. Mater., vol. 161, no. 1, pp. 457–462, 2009. [91] S. O. Kim, W. S. Kim, and K. W. Kim, “Evaluation of electrokinetic remediation of arsenic-contaminated soils,” Environ. Geochem. Health, vol. 27, no. 5–6, pp. 443–453, 2005. [92] S. O. Kim, S. H. Moon, and K. W. Kim, “Removal of heavy metals from soils using enhanced electrokinetic soil processing,” Water. Air. Soil Pollut., vol. 125, no. 1–4, pp. 259–272, 2001. [93] D. B. Gent, A. H. Wani, J. L. Davis, and A. Alshawabkeh, “Electrolytic redox and electrochemical generated alkaline hydrolysis of hexahydro-1,3,5-trinitro-1,3,5 triazine (RDX) in sand columns,” Environ. Sci. Technol., vol. 43, no. 16, pp. 6301–6307, 2009. [94] S. Hwang, D. R. Felt, E. J. Bouwer, M. C. Brooks, S. L. Larson, and J. L. Davis, 136“Remediation of RDX-contaminated water using alkaline hydrolysis,” J. Environ. Eng., vol. 132, no. 2, pp. 256–262, 2006. [95] D. M. Gilbert and T. C. Sale, “Sequential electrolytic oxidation and reduction of aqueous phase energetic compounds,” Environ. Sci. Technol., vol. 39, no. 23, pp. 9270–9277, 2005. [96] J. D. Rodgers and N. J. Bunce, “Electrochemical treatment of 2,4,6-trinitrotoluene and related compounds,” Environ. Sci. Technol., vol. 35, no. 2, pp. 406–410, Jan. 2001. [97] K. R. Reddy, C. Cameselle, and P. Ala, “Integrated electrokinetic-soil flushing to remove mixed organic and metal contaminants,” J. Appl. Electrochem., vol. 40, no. 6, pp. 1269–1279, 2010. [98] L. M. Ottosen and I. Rörig-Dalgaard, “Desalination of a brick by application of an electric DC field,” Mater. Struct., vol. 42, no. 7, pp. 961–971, 2009. [99] R. E. Beddoe and H. W. Dorner, “Modelling acid attack on concrete: Part I. the essential mechanisms,” Cem. Concr. Res., vol. 35, no. 12, pp. 2333–2339, 2005. [100] R. A. Jacobs, M. Z. Sengun, R. E. Hicks, and R. F. Probstein, “Model and experiments on soil remediation by electric fields,” J. Environ. Sci. Heal., vol. 29, no. 9, pp. 1933–1955, 1994. [101] R. F. Probstein and R. E. Hicks, “Removal of contaminants from soils by electric fields,” Science (80-. )., vol. 260, no. 5107, pp. 498–503, 1993. [102] J. G. G. Sah and J. Y. Y. Chen, “Study of the electrokinetic process on Cd and Pb spiked soils,” J. Hazard. Mater., vol. 58, no. 1–3, pp. 301–315, Feb. 1998. [103] Q. Kang, P. C. Lichtner, and D. Zhang, “An improved lattice Boltzmann model for multicomponent reactive transport in porous media at the pore scale,” Water Resour. Res., 137vol. 43, no. 12, p. W12S14, 2007. [104] P. Rasouli, C. I. Steefel, K. U. Mayer, and M. Rolle, “Benchmarks for multicomponent diffusion and electrochemical migration,” Comput. Geosci., vol. 19, no. 3, pp. 523–533, 2015. [105] K. U. Mayer, M. Xie, D. Su, and K. T. B. MacQuarrie, “MIN3P-THCm: a three-dimensional numerical model for multicomponent reactive transport in variably saturated porous media user guide,” Vancouver, Jul. 2014. [106] S. O. Kim, J. J. Kim, S. T. Yun, and K. W. Kim, “Numerical and experimental studies on Cadmium(II) transport in kaolinite clay under electrical fields,” Water, Air Soil Pollut., vol. 150, pp. 135–162, 2003. [107] A. T. Yeung and S. Datla, “Fundamental formulation of electrokinetic extraction of contaminants from soil,” Can. Geotech. J., vol. 32, no. 4, pp. 569–583, Aug. 1995. [108] K. U. Mayer, “A numerical model for multicomponent reactive transport in variably saturated porous media,” University of Waterloo, 1999. [109] B. Johannesson, “Comparison between the Gauss’ law method and the zero current method to calculate multi-species ionic diffusion in saturated uncharged porous materials,” Comput. Geotech., vol. 37, no. 5, pp. 667–677, 2010. [110] S. Molins and K. U. Mayer, “Coupling between geochemical reactions and multicomponent gas and solute transport in unsaturated media: A reactive transport modeling study,” Water Resour. Res., vol. 43, no. 5, pp. 1–16, May 2007. [111] K. U. Mayer, D. W. Blowes, and E. O. Frind, “Reactive transport modeling of an in situ reactive barrier for the treatment of hexavalent chromium and trichloroethylene in groundwater,” Water Resour. Res., vol. 37, no. 12, pp. 3091–3103, 2001. 138[112] K. U. Mayer, S. G. Benner, E. O. Frind, S. F. Thornton, and D. N. Lerner, “Reactive transport modeling of processes controlling the distribution and natural attenuation of phenolic compounds in a deep sandstone aquifer,” J. Contam. Hydrol., vol. 53, no. 3–4, pp. 341–368, Dec. 2001. [113] J. Jurjovec, D. W. Blowes, C. J. Ptacek, and K. U. Mayer, “Multicomponent reactive transport modeling of acid neutralization reactions in mine tailings,” Water Resour. Res., vol. 40, no. 11, p. W11202, Nov. 2004. [114] K. U. Mayer, S. G. Benner, and D. W. Blowes, “Process-based reactive transport modeling of a permeable reactive barrier for the treatment of mine drainage.,” J. Contam. Hydrol., vol. 85, no. 3–4, pp. 195–211, May 2006. [115] A. T. Yeung, “Fundamental aspects of prolonged electrokinetic flows in kaolinites,” Geomech. Geoengin., vol. 1, no. 1, pp. 13–25, 2006. [116] R. W. Falta, I. Javandel, K. Pruess, and P. A. Witherspoon, “Density-driven flow of gas in the unsaturated zone due to the evaporation of volatile organic compounds,” Water Resour. Res., vol. 25, no. 10, pp. 2159–2169, Oct. 1989. [117] J. M. Darmon, N. Kumar, E. B. Hulley, C. J. Weiss, S. Raugei, R. M. Bullock, and M. L. Helm, “Increasing the rate of hydrogen oxidation without increasing the overpotential: a bio-inspired iron molecular electrocatalyst with an outer coordination sphere proton relay,” Chem. Sci., vol. 00, no. 5, pp. 1–9, 2015. [118] J. A. Cracknell, A. F. Wait, O. Lenz, B. Friedrich, and F. A. Armstrong, “A kinetic and thermodynamic understanding of O2 tolerance in [NiFe]-hydrogenases.,” Proc. Natl. Acad. Sci., vol. 106, no. 49, pp. 20681–20686, 2009. [119] B. Schink and H. G. Schlegel, “Hydrogen metabolism in aerobic hydrogen-oxidizing 139bacteria,” Biochimie, vol. 60, no. 3, pp. 297–305, 1978. [120] L. Bongers, “Energy generation and utilization in hydrogen bacteria.,” J. Bacteriol., vol. 104, no. 1, pp. 145–151, 1970. [121] D. Strmcnik, M. Uchimura, C. Wang, R. Subbaraman, N. Danilovic, D. van der Vliet, A. P. Paulikas, V. R. Stamenkovic, and N. M. Markovic, “Improving the hydrogen oxidation reaction rate by promotion of hydroxyl adsorption.,” Nat. Chem., vol. 5, no. 4, pp. 300–306, 2013. [122] L. Johnson, A. Ejigu, Licence, and D. A. Walsh, “Hydrogen oxidation and oxygen reduction at platinum in protic ionic liquids,” J. Phys. Chem. C, vol. 116, no. 34, pp. 18048–18056, 2012. [123] K. U. Mayer and K. T. B. MacQuarrie, “Solution of the MoMaS reactive transport benchmark with MIN3P-model formulation and simulation results,” Comput. Geosci., vol. 14, no. 3, pp. 405–419, 2010. [124] D. R. Lide, CRC Handbook of Chemistry and Physics, Internet Version 2005. Boca Raton: CRC Press, 2005. [125] A. Mohajeri, G. A. Narsilio, P. Pivonka, and D. W. Smith, “Numerical estimation of effective diffusion coefficients for charged porous materials based on micro-scale analyses,” Comput. Geotech., vol. 37, no. 3, pp. 280–287, Apr. 2010. [126] E. J. F. Dickinson, L. Freitag, and R. G. Compton, “Dynamic Theory of Liquid Junction Potentials,” J. Phys. Chem. B, vol. 114, no. 1, pp. 187–197, Dec. 2009. [127] A. Revil and D. Jougnot, “Diffusion of ions in unsaturated porous materials,” J. Colloid Interface Sci., vol. 319, no. 1, pp. 226–235, 2008. [128] H. J. Hickman, “The liquid junction potential - the free diffusion junction,” Chem. Eng. 140Sci., vol. 25, no. 3, pp. 381–398, 1970. [129] D. Hodges, A. Fourie, D. A. Reynolds, and D. Thomas, “Development of an apparatus for pH-isolated electrokinetic in situ chemical oxidation,” J. Environ. Eng., vol. 137, no. 9, pp. 809–816, 2011. [130] M. Pazos, E. Rosales, T. Alcántara, J. Gómez, and M. A. Sanromán, “Decontamination of soils containing PAHs by electroremediation: A review,” J. Hazard. Mater., vol. 177, pp. 1–11, 2010. [131] I. Fried, The Chemistry of Electrode Processes, 1st ed. London: Academic Press Inc., 1973. [132] P. M. Biesheuvel, M. van Soestbergen, and M. Z. Bazant, “Imposed currents in galvanic cells,” Electrochim. Acta, vol. 54, no. 21, pp. 4857–4871, Aug. 2009. [133] M. N. V. Prasad, K. S. Sajwan, and R. Naidu, Trace Elements in the Environment: Biogeochemistry, Biotechnology and Bioremediation. Taylor & Francis, CRC Press, 2006. [134] M. L. Soudijn, “Proton transport in aqueous ionic solutions,” University of Amsterdam, 2012. [135] M. N. Nandanwar and S. Kumar, “Modelling of effect of non-uniform current density on the performance of soluble lead redox flow batteries,” J. Electrochem. Soc., vol. 161, no. 10, pp. A1602–A1610, 2014. [136] K. J. Tielrooij, N. Garcia-Araez, M. Bonn, and H. J. Bakker, “Cooperativity in ion hydration,” Science (80-. )., vol. 328, no. 5981, pp. 1006–1009, May 2010. [137] K. J. Tielrooij, R. L. A. Timmer, H. J. Bakker, and M. Bonn, “Structure dynamics of the proton in liquid water probed with terahertz time-domain spectroscopy,” Phys. Rev. Lett., vol. 102, no. 19, p. 198303, May 2009. 141[138] J. Crank, The Mathematics of Diffusion, 2nd ed. Oxford: Oxford University Press, 1975. [139] R. Sander, “Compilation of Henry’s law constants (version 4.0) for water as solvent,” Atmos. Chem. Phys., vol. 15, no. 8, pp. 4399–4981, Apr. 2015. [140] G. Barth and M. C. Hill, “Numerical methods for improving sensitivity analysis and parameter estimation of virus transport simulated using sorptive-reactive processes.,” J. Contam. Hydrol., vol. 76, pp. 251–77, Feb. 2005. [141] E. P. Poeter, M. C. Hill, E. R. Banta, S. Mehl, and S. Christensen, “UCODE_2005 and six other computer codes for universal sensitivity analysis, calibration, and uncertainty evaluation,” Reston, 2005. [142] J. Doherty, “PEST: Model-independent parameter estimation user manual,” 2010. [143] B. E. Bongenaar-Schlenter, “Ohmic resistance and current distribution at gas-evolving electrodes,” Technische Hogeschool Eindhoven, 1984. [144] J. C. Maxwell, A Treatise on Electricity and Magnetism, vol. 1. Oxford: Clarendon Press, 1873. [145] Y. B. Acar, R. J. Gale, A. N. Alshawabkeh, R. E. Marks, S. Puppala, M. Bricka, and R. Parker, “Electrokinetic remediation: basics and technology status,” J. Hazard. Mater., vol. 40, no. 2, pp. 117–137, 1995. [146] V. Naudet, A. Revil, E. Rizzo, J. Y. Bottero, and P. Bégassat, “Groundwater redox conditions and conductivity in a contaminant plume from geoelectrical investigations,” Hydrol. Earth Syst. Sci., vol. 8, no. 1, pp. 8–22, 2004. [147] S. P. Gay, “A 1800 millivolt self-potential anomaly near Hualgayoc, Peru,” Geophys. Prospect., vol. 15, no. 2, pp. 236–245, 1967. [148] C. E. Corry, “Spontaneous polarization associated with porphyry sulfide mineralization,” 142Geophysics, vol. 50, no. 6, pp. 1020–1034, 1985. [149] M. Goldie, “Self-potentials associated with the Yanacocha high-sulfidation gold deposit in Peru,” Geophysics, vol. 67, no. 3, pp. 684–689, May 2002. [150] B. J. Minsley, J. Sogade, and F. D. Morgan, “Three-dimensional source inversion of self-potential data,” J. Geophys. Res. Solid Earth, vol. 112, no. 2, p. B02202, 2007. [151] B. J. Minsley, J. Sogade, and F. D. Morgan, “Three-dimensional self-potential inversion for subsurface DNAPL contaminant detection at the Savannah River Site, South Carolina,” Water Resour. Res., vol. 43, no. 4, p. W04429, 2007. [152] D. Ntarlagiannis, E. A. Atekwana, E. A. Hill, and Y. Gorby, “Microbial nanowires: Is the subsurface ‘hardwired’?,” Geophys. Res. Lett., vol. 34, no. 17, p. L17305, 2007. [153] R. Doherty, B. McPolin, B. Kulessa, A. Frau, A. Kulakova, C. C. R. Allen, and M. J. Larkin, “Microbial ecology and geoelectric responses across a groundwater plume,” Interpretation, vol. 3, no. 4, pp. SAB9–SAB21, 2015. [154] R. Doherty, B. Kulessa, A. S. Ferguson, M. J. Larkin, L. A. Kulakov, and R. M. Kalin, “A microbial fuel cell in contaminated ground delineated by electrical self-potential and normalized induced polarization data,” J. Geophys. Res, vol. 115, p. G00G08, 2010. [155] S. J. S. Fachin, E. L. Abreu, C. A. Mendonça, A. Revil, G. C. Novaes, and S. S. Vasconcelos, “Self-potential signals from an analog biogeobattery model,” Geophysics, vol. 77, no. 4, pp. EN29–EN37, 2012. [156] S. A. Forte and L. R. Bentley, “Mapping degrading hydrocarbon plumes with self potentials: investigation on causative mechanisms using field and modeling data,” J. Environ. Eng. Geophys., vol. 18, no. 1, pp. 27–42, Mar. 2013. [157] A. Revil, L. Ehouarne, and E. Thyreault, “Tomography of self-potential anomalies of 143electrochemical nature,” Geophys. Res. Lett., vol. 28, no. 23, pp. 4363–4366, Dec. 2001. [158] P. Martinez-Pagan, A. Jardani, A. Revil, and A. Haas, “Self-potential monitoring of a salt plume,” Geophysics, vol. 75, no. 4, pp. WA17–WA25, Jul. 2010. [159] J. E. Katz, X. Zhang, K. Attenkofer, K. W. Chapman, C. Frandsen, P. Zarzycki, K. M. Rosso, R. W. Falcone, G. a Waychunas, and B. Gilbert, “Electron small polarons and their mobility in iron (oxyhydr)oxide nanoparticles.,” Science (80-. )., vol. 337, no. 6099, pp. 1200–1203, 2012. [160] P. Sivenas and F. W. Beales, “Natural geobatteries associated with sulphide ore deposits, I. theoretical studies,” J. Geochemical Explor., vol. 17, pp. 123–143, 1982. [161] D. J. Griffiths, Introduction to Electrodynamics, 3rd ed. New Jersey: Prentice Hall, 1997. [162] A. Revil, N. Linde, A. Cerepi, D. Jougnot, S. Matthäi, and S. Finsterle, “Electrokinetic coupling in unsaturated porous media,” J. Colloid Interface Sci., vol. 313, no. 1, pp. 315–327, 2007. [163] N. Linde and A. Revil, “Inverting self-potential data for redox potentials of contaminant plumes,” Geophys. Res. Lett., vol. 34, no. 14, p. L14302, 2007. [164] A. Revil, V. Naudet, J. Nouzaret, and M. Pessel, “Principles of electrography applied to self-potential electrokinetic sources and hydrogeological applications,” Water Resour. Res., vol. 39, no. 5, p. 1114, 2003. [165] E. J. Reardon, “Anaerobic corrosion of granular iron: measurement and interpretation of hydrogen evolution rates.,” Environ. Sci. Technol., vol. 29, no. 12, pp. 2936–2945, 1995. [166] J. D. Allison, D. S. Brown, and K. J. Novo-Gradac, “MINTEQA2/PRODEFA2: a geochemical assessment model for environmental systems: version 3.0 user’s manual,” U. S. Environ. Prot. Agency EPA/600/3-91/021, no. March, p. 115, 1991. 144[167] C. R. Wilke and P. Chang, “Correlation of diffusion coefficients in dilute solutions,” AIChE J., vol. 1, no. 2, pp. 264–270, 1955. [168] R. Schieck, A. Hartmann, S. Fiechter, R. Könenkamp, and H. Wetzel, “Electrical properties of natural and synthetic pyrite (FeS2) crystals,” J. Mater. Res., vol. 5, no. 07, pp. 1567–1572, 1990. [169] M. R. Thornber, “Supergene alteration of sulphides, II. A chemical study of the Kambalda nickel deposits,” Chem. Geol., vol. 15, no. 2, pp. 117–144, 1975. [170] M. R. Thornber, “Supergene alteration of sulphides, I. A chemical model based on massive nickel sulphide deposits at Kambalda, Western Australia,” Chem. Geol., vol. 15, no. 1, pp. 1–14, 1975. [171] M. I. Leybourne, S. M. Hamilton, C. J. Mwenifumbo, W. D. Goodfellow, and D. R. Boyle, “Hydrogeochemical and geophysical evidence for borehole oxidation of massive sulphides below the water table and generation of self-potential anomalies,” Geochemistry Explor. Environ. Anal., vol. 9, no. 1, pp. 39–50, 2009. [172] G. Reguera, K. D. McCarthy, T. Mehta, J. S. Nicoll, M. T. Tuominen, and D. R. Lovley, “Extracellular electron transfer via microbial nanowires.,” Nature, vol. 435, no. 7045, pp. 1098–1101, 2005. [173] G. Reguera, K. P. Nevin, J. S. Nicoll, S. F. Covalla, T. L. Woodard, and D. R. Lovley, “Biofilm and nanowire production leads to increased current in Geobacter sulfurreducens fuel cells,” Appl. Environ. Microbiol., vol. 72, no. 11, pp. 7345–7348, 2006. [174] Y. A. Gorby, S. Yanina, J. S. McLean, K. M. Rosso, D. Moyles, A. Dohnalkova, T. J. Beveridge, I. S. Chang, B. H. Kim, K. S. Kim, D. E. Culley, S. B. Reed, M. F. Romine, D. A. Saffarini, E. A. Hill, L. Shi, D. A. Elias, D. W. Kennedy, G. Pinchuk, K. Watanabe, S. 145Ishii, B. Logan, K. H. Nealson, and J. K. Fredrickson, “Electrically conductive bacterial nanowires produced by Shewanella oneidensis strain MR-1 and other microorganisms,” Proc. Natl. Acad. Sci., vol. 103, no. 30, pp. 11358–11363, Jul. 2006. [175] L. P. Nielsen, N. Risgaard-Petersen, H. Fossing, P. B. Christensen, and M. Sayama, “Electric currents couple spatially separated biogeochemical processes in marine sediment,” Nature, vol. 463, no. 7284, pp. 1071–1074, Feb. 2010. 146Appendix A Extended Nernst-Planck formulation: This equation describes total mass transfer of interacting species by advection, diffusion and electromigration [28], [79]: ۸୧ୟ = ϕτ ܙୟC୧ୟ−Sୟ۲୧∇C୧ୟ −FRTSୟ۲୧C୧ୟz୧∇ψ൨ (A.1) where ϕ is porosity (m3 void m-3 bulk), τ is tortuosity (m void m-1 bulk), Sa is the saturation of the aqueous phase (m3 fluid m-3 void), Di is the species-dependent diffusion coefficient (m2 s-1), qa is Darcy’s velocity (m s-1), Cia is the species concentration (mol L-1 H2O), F is the Faraday constant (96485 C mol-1), R is the gas constant (8.341 J K-1 mol-1), T is the absolute temperature (K), zi is the charge number (-) and ψ is the electric potential (V or J C-1). Assuming zero external electric current (null current assumption), the sum of all electric currents induced by species migration should be equal to zero. Therefore, the total electric current density can be expressed as: ۷ୟ = F z୧۸୧ୟ = 0౪౪୧ (A.2) where ۷ୟ is the total electric current density (C s-1 m-2 or A m-2) and Ntot is the total number of aqueous species. Substituting (A.1) into (A.2) leads to: ۷ୟ = F z୧۸୧ୟ౪౪୧= Fϕτ z୧౪౪୧ܙୟC୧ୟ−Sୟ۲୧∇C୧ୟ −FRTSୟ۲୧C୧ୟz୧∇ψ൨ = 0 (A.3) Here, the sums are taken independently for each term: Fϕτ ܙୟ z୩౪౪୩C୩ୟ−Sୟ z୪౪౪୪۲୪∇C୪ୟ −FRTSୟ∇ψ ۲୫C୫ୟ z୫ଶ౪౪୫ = 0 (A.4) Rearranging equation (A.4) leads to: 147ϕτFଶRTSୟ∇ψ ۲୫C୫ୟ z୫ଶ౪౪୫= Fϕτܙୟ z୩౪౪୩C୩ୟ− FϕτSୟ z୪౪౪୪۲୪∇C୪ୟ (A.5) When solving for the electric potential term: ∇ψ =Fϕτܙୟ ∑ z୩౪౪୩ C୩ୟ− FϕτSୟ ∑ z୪౪౪୪ ۲୪∇C୪ୟϕτFଶRT Sୟ ∑ ۲୫C୫ୟ z୫ଶ౪౪୫ (A.6) similar terms cancel out: ∇ψ =ܙୟ ∑ z୩౪౪୩ C୩ୟ− Sୟ ∑ z୪౪౪୪ ۲୪∇C୪ୟFRT Sୟ ∑ ۲୫C୫ୟ z୫ଶ౪౪୫ (A.7) By rearranging we obtain a new expression for electric potential gradient as a function of concentrations of the species: ∇ψ =RTFܙୟ ∑ z୩౪౪୩ C୩ୟ− Sୟ ∑ z୪౪౪୪ ۲୪∇C୪ୟSୟ ∑ ۲୫C୫ୟ z୫ଶ౪౪୫ (A.8) Substituting equation (A.8) into (A.1) yields: ۸୧ୟ = ϕτ ܙୟC୧ୟ−Sୟ۲୧∇C୧ୟ − Sୟ۲୧C୧ୟz୧ܙୟ ∑ z୩౪౪୩ C୩ୟ− Sୟ ∑ z୪౪౪୪ ۲୪∇C୪ୟSୟ ∑ ۲୫C୫ୟ z୫ଶ౪౪୫൩ (A.9) This formulation is valid when there is no external electric potential. This is the equation implemented in MIN3P and used in chapter 2 for benchmarking multicomponent diffusion problems. Gauss-Ampère’s Law: On the REV scale and for quasi-steady state conditions, Gauss-Ampère’s Law for the ith species in the aqueous phase a can be written as [36]: ϕSୟ۷ୟ = −∂(ϕSୟ۲ୟ )∂t+ સ × (ϕSୟ۶ୟ ) − સ × (ϕSୟ۾ୟ × ܞୟ) + ϕSୟܐመ ୠୟୟஷୠ (A.10) 148where ۲ୟ is electric displacement (C m-2), ۶ୟ is magnetic field intensity (A m), ۾ୟ is average polarization density (C m-2), ܞ܉ is the mass-average velocity vector (m s-1), ܐመ ୠୟ is external exchange current density between solution (phase a) and electrodes (phase b) (A m-2). To maintain continuity of charge, we take the divergence of equation (A.10) and assuming no water polarization; the second and third terms on the right hand side of the equation will be eliminated (Reminder: the divergence of the curl of any vector field is always zero): ∂ ቀ∇ ∙ (ϕSୟ۲ୟ )ቁ∂t+ ∇ ∙ (ϕSୟ۷ୟ ) − ∇ ∙ ൫ϕSୟܐመ ୠୟ ൯ୟஷୠ൩ = 0 (A.11) According to Gauss’ law charge density is equal to change in electric displacement field; so the first term on the left-hand side becomes: ∇ ∙ (ϕSୟ۲ୟ ) = ϕSୟF C୧ୟz୧౪౪୧ (A.12) By substituting (A.12) into (A.11), we obtain: ∂൫ϕSୟF ∑ C୧ୟz୧౪౪୧ ൯∂t+ ∇ ∙ (ϕSୟ۷ୟ ) − ∇ ∙ ൫ϕSୟܐመ ୠୟ ൯ୟஷୠ൩ = 0 (A.13) Substituting (A.3) into (A.13) gives an expression for charge continuity in the aqueous phase: ∂൫ϕSୟF ∑ C୧ୟz୧౪౪୧ ൯∂t+ ∇ ∙ ൦Fϕτ z୧౪౪୧ܙୟC୧ୟ−Sୟ۲୧∇C୧ୟ −FRTSୟ۲୧C୧ୟz୧∇ψ൨൪ − ∇ ∙ ൫ϕSୟܐመ ୠୟ ൯ୟஷୠ൩ = 0 (A.14) When there is no externally induced electric current (ܐመ ୠୟ = ), the last term on the left hand side becomes zero. This is the equation implemented in MIN3P for simulation of electrokinetic remediation presented in chapter 3. 149For the geobattery model presented in chapter 4, advection is neglected. The Nernst-Planck equation without the contribution of advection is: ۸୧ୟ = −ϕτSୟ۲୧ ∇C୧ୟ +FRTC୧ୟz୧∇ψ൨ (A.15) Electron Drift: This equation defines the movement of electrons in the electrical conductor [161]: ۸܍ = −enμୣܨ∇ψ (A.16) where e is the elementary electron charge (1.6 × 10-19 C), n is number of free electrons per volume (m-3) and μe is electron mobility (m2 V-1 s-1). Total electric current flux in porous media is the sum of electric current attribution from aqueous species in the aqueous phase and electron drift in the electrical conductor: ۷۸ୟ = F ቌ z୧Jia + Je౪౪୧ቍ (A.17) The first term on the right-hand side corresponds to equation (A.3) without the advection term. By substituting (A.3) and (A.16) into (A.17), we obtain: ۷۸ୟ = −FϕτSୟ ۲୧z୧౪౪୧∇C୧ୟ +FRTC୧ୟz୧∇ψ൨ − (1 − ϕ) enμୣ∇ψ (A.18) According to the geobattery model [12], [43], [46], [47], [49], [146], [162] a source (or driving) electric current density, ۷܁ୟ (A m-2), is established in the presence of an electrical conductor and a redox potential gradient, facilitating electron transfer from the reducing regions to the oxidizing regions, defined as: 150۷܁ୟ = −ߪ∇ൣ൫E୫_ୟ୬ − Eୌ_ୟ୬൯ − ൫E୫_ୡୟ − Eୌ_ୡୟ൯൧ (A.19) where σm (S m-1) is the electric conductivity of the electrical conductive medium (ore or metallic body) and EH_an and EH_ca (V) are the equilibrium redox potentials for the anodic and cathodic reactions, respectively. Em_an and Em_ca are the redox potentials of the electrical conductive medium at anode and cathode ends. The total electric current density, ۷ܜܗܜୟ , is the sum of electric current density by aqueous species in the aquifer and electron drift in the ore body, ۷۸ୟ , and source current density ۷܁ୟ: ۷ܜܗܜୟ = ۷۸ୟ + ۷܁ୟ (A.20) Applying Gauss-Ampère’s law (i.e. continuity of charge written in terms of electric flux) for (A.20): ∂൫ϕSୟF ∑ C୧ୟz୧ା୧ ൯∂t+ સ ∙ ۷ܜܗܜୟ = 0 (A.21) or: ∂൫ϕSୟF ∑ C୧ୟz୧ା୧ ൯∂t+ સ ∙ ൫۷۸ୟ + ۷܁ୟ൯ = 0 (A.22) Substituting equations (A.18) and (A.19) in to (A.22) yields: ∂൫ϕSୟF ∑ C୧ୟz୧ା୧ ൯∂t− સ ∙ ቌFϕτSୟ ۲୧z୧౪౪୧∇C୧ୟ +FRTC୧ୟz୧∇ψ൨ − (1 − ϕ) enμୣ∇ψቍ − સ ∙ ൫ߪ∇ൣ൫E୫ − Eୌ൯ − ൫E୫ౙ − Eୌౙ൯൧൯ = 0 (A.23) This is the equation implemented in MIN3P for modeling of the geobattery concept, presented in chapter 4.