This is the course outline.
This lesson introduces Reactive Transport Models (RTMs), primarily focusing on a brief history of RTM development, governing equations, and key concepts. It includes a lightboard video about the governing RT equations, a video that introduces the code that we will use in this class, CrunchFlow, and the reading materials here. The idea here is to give you an overview of RTM. Many concepts introduced here will be detailed in later lessons so it is OK if you do not fully grasp them in this lesson.
Please note that this is not a course that teaches how to numerically solve for reactive transport equations, which deserve a separate course by itself. Instead, this is a course that teaches fundamental reactive transport concepts and how to use an existing software CrunchFlow to solve and answer specific questions. So this is a model application course, not a numerical method course.
By the end of this lesson, you should:
There are example files and hw files in each lesson (almost). If you would like extra exercise files, click here [1].
(Optional) Reading |
Note: these are your references. You can skip through quickly in this lesson to get some overview idea. You will use these again and again later. |
---|---|
To Do |
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
Reactive transport models have been applied to understand biogeochemical systems for more than three decades (Beaulieu et al., 2011; Brown and Rolston, 1980; Chapman, 1982; Chapman et al., 1982; Lichtner, 1985; Regnier et al., 2013; Steefel et al., 2005; Steefel and Lasaga, 1994). Multi-component RTMs originated in the 1980s based on the theoretical foundation of the continuum model (Lichtner, 1985; Lichtner, 1988). RTM development advanced substantially in the 1990s with the emergence of several extensively-used RTM codes, including Hydrogeochem (Yeh and Tripathi, 1991; Yeh and Tripathi, 1989), CrunchFlow (Steefel and Lasaga, 1994), Flotran (Lichtner et al., 1996), Geochemist’s Workbench (Bethke, 1996), Phreeqc (Parkhurst and Appelo, 1999), Min3p (Mayer et al., 2002), STOMP (White and Oostrom, 2000), TOUGHREACT (Xu et al., 2000), among others (Ortoleva et al., 1987, Bolton, 1996). Several early diagenetic models were developed at similar times, including STEADYSED (Van Cappellen and Wang, 1996), CANDI (Boudreau, 1996), and OMEXDIA (Soetaert et al., 1996). These codes can be easily found through google their names.
RTMs are distinct from geochemical models that primarily calculate geochemical equilibrium, speciation, and thermodynamic state of a system (Wolery et al., 1990). RTMs also differ from reaction path models (Helgeson, 1968; Helgeson et al., 1969) that represent closed or batch systems without diffusive or advective transport. The major advance of modern RTMs was to couple flow and transport within a full geochemical thermodynamic and kinetic framework (Steefel et al., 2015).
RTMs have been used across an extensive array of environments and applications (as reviewed in MacQuarrie and Mayer, 2005; Steefel et al., 2005). One primary focus has been in the low-temperature (ca < 100˚C) surface and near-surface environment where “rock meets life”, a region often referred to as the Critical Zone (CZ). Within the critical zone, water, atmosphere, rock, soil, and life interact creating the potential for complex chemical, physical, and biological interactions and responses to external forcing. As illustrated in Figure 0.1, RTMs can simulate a wide range of processes in this environment, including fluid flow (single or multiphase), solute transport (advective, dispersive, and diffusive transport), geochemical reactions (e.g., mineral dissolution and precipitation, ion exchange, surface complexation), and biogeochemical processes (e.g., microbe-mediated redox reactions, biomass growth and decay).
RTMs with these capabilities have been applied to understand chemical weathering and soil formation in response to various biological, climatic and physical drivers. RTMs have also been essential to address a wide range of questions at the nexus of energy and the environment, including, for example, environmental bioremediation (Bao et al., 2014), natural attenuation (Mayer et al., 2001), geological carbon sequestration (Apps et al., 2010; Brunet et al., 2013; Navarre-Sitchler et al., 2013), and nuclear waste disposal (Saunders and Toran, 1995; Soler and Mader, 2005). Model frameworks have advanced to incorporate heterogeneous characteristics of natural systems to begin to understand the role of spatial heterogeneities in controlling flow and the interaction between water and reacting components (Scheibe et al., 2006; Yabusaki et al., 2011; Liu et al., 2013). With the expansion of isotopes as tracers of mineral-fluid and biologically mediated reactions, recent advances include development of RTMs that allow an explicit treatment of isotopic partitioning due to both kinetic and equilibrium process.
The RTM approach has also been used to investigate subsurface processes at spatial scales ranging from single pores (Kang et al., 2006; Li et al., 2008; Molins et al., 2012) and single cells (Scheibe et al., 2009; Fang et al., 2011), to pore networks and columns (1 -10s centimeters) (Knutson et al., 2005; Li et al., 2006; Yoon et al., 2012; Druhan et al., 2014), and to field scales (1-10’s of meter) (Li et al., 2011), with a few studies at the watershed or catchment scale (100s of meters) (Atchley et al., 2014). Recent weathering studies have linked regional scale reactive transport models (WITCH) to global climate models to understand the role of climate change in controlling weathering (Godderis et al., 2006; Roelandt et al., 2010). Recent model development also includes full coupling between subsurface biogeochemical processes and surface hydrology, land-surface interactions, meteorological and climatic forcings (Bao et al., 2017; Li et al., 2017a). Such coupling has been argued to be important in understanding the complex interactions between processes of interests in different disciplines (Li et al., 2017b).
Most reaction transport codes solve equations of mass, momentum, and energy conservation (Steefel et al., 2005). For mass conservation, reactive transport models usually partition aqueous species into primary and secondary species (Lichtner, 1985). The primary species are the building blocks of chemical systems of interest, upon which concentrations of secondary species are written through laws of mass action for reactions at thermodynamic equilibrium. The partition between primary and secondary species allows the reduction of computational cost by only solving for mass conservation equations for primary species and then calculating secondary species through thermodynamics. Detailed discussion on primary and secondary species will be in lesson 1 on Aqueous Complexation.
Please watch the following video: Reactive Transport Reactions (7:25)
PRESENTER: What I'm going to do today is really to introduce the general idea of reactive transport equation and the overview of it. This equation that I put here is we call mass conservation equation for chemical species in aqueous phase, for one of the representative species. I don't expect you to know every term, or the details of every term in this equation. Because we will be talking more about each term later. But this is really to give you a general idea and overview before we start about everything.
And the importance of this is that we run these reactive transport codes. And there are a lot of built up architecture behind the code, in terms of what they solve. And it's always a good idea to know some of these, like what they solve, and what are the things behind these codes. So what I will do today is really talk through each of these terms and talk about the [? fake ?] meaning of each term, so that you get a general idea what are they really for.
So the first term we call the mass conservation term. It's called mass accumulation rate. And it should have the units of mol per length cubed, which is volume, of porous media per time, whatever time you pick. But all the terms have to be consistent, have the same length and time units. So this equation really says mass accumulation rate depends on several different processes, right?
So the first term is the rate itself. And the overall rate, the second term, what we call dispersive and diffusive transport. So you think about how a chemical species in water has a change over time, which is this. So some of these rates coming from, for example, the chemical species gold to get different concentrations in different locations. For example, you think about a dye put in a cup of water. And over time, they tend to have the same color everywhere. So this is one of the driving forces, in terms of what we call diffusive or dispersive transport. And they should have the same unit as the first term. Every term should have the same unit. So that's a second term.
And then, the third term is what we call the advective transport. And this is a process where, for example, you think about rivers, right? And the chemical species will flow together with the water. And so essentially, the water brings the chemical species to different places. So this advective transport, in this term, you have the hue, which is we call Darcy velocity. And then, the concentration actually I probably should explain here, the Ci will be the concentration of one chemical species of species i, a representative species, i. And the Ci everywhere is the same. So this is advective transport.
But also in a lot of systems, you have reactions, right? So the last term, four, is for the total reaction rates. And again overall, it has the same units as the first term. But essentially, this could it be a summation of several different-- let's call this equal to summation of i, which is for the chemical species, i. But this i could be involved in, let's say, ik, different numbers of reactions. So essentially, you will be adding all these reactions that this can chemical species, i, is participating in. And this ik would be a total number of reactions.
OK, so this is one representative equation. We write this equation form for the species, i. But I put the i from 1 to n here, meaning you can have any arbitrary number of what we call a primary species. So if you have, let's say, 10 different chemical species, then the primary species you have n equal to 10. And you will be writing 10 of these equations.
And these equations, essentially, if we solve these equations, you get the concentration or different chemical species as a function of time and space. So essentially, the outcome of this is the temporal and spatial distribution of chemical species. So essentially, you can tracing after each chemical species and look at how they change as a function of time and space and how, in different parts, they have different rates, and all that. Species i and i from 1 to n.
OK, so that's what you hope you guys will be exposed to in using this code. You will be solving this equations for particular, specific questions, problem applications. And then, you're usually given a set of initial conditions, like where is the concentration of different species at time 0 at different locations, and then over time, how these countries have different species evolve over time. And you will see you will learn a lot about these, using this code. And we'll be talking more about this in each of these terms, what are they, over time in different lessons that will follow.
The following is a representative mass conservation equation for a primary species I in the aqueous phase:
Here Ci is the total concentration of species i (mol/m3 pore volume), t is the time (s), n is the number of primary species, D is the combined dispersion–diffusion tensor (m2/s), u (m/s) is the Darcy flow velocity vector and can be decomposed into ux and uz in the directions parallel and transverse to the main flow direction. Nr is the total number of kinetic aqueous reactions that involve species i,
Equation (1) implies that the mass change rate of species i depends on physical and chemical processes: the diffusion/dispersion processes that are accounted for by the first term of the right hand side of the equation, the advection process that is taken into account by the second term of the right hand side, and reactions that are represented by the last term of the equation. The last term is the summation of multiple reaction rates, the form of which depend on the number and type of kinetic reactions that species i is involved in. The reaction terms include the rates of kinetically controlled reactions including microbe-mediated bioreduction reactions, mineral dissolution and precipitation reactions, and redox reactions. The reactions also include fast reactions that are considered at thermodynamics equilibrium, including aqueous complexation, ion exchange, and surface complexation. These fast reactions that are at equilibriums however do not show up in the above governing equation (1). Instead they exhibit themselves in the non-linear coupling of primary and secondary species through the expression of equilibrium constants (laws of mass action), as will be detailed later.
The dispersion-diffusion tensor D is defined as the sum of the mechanical dispersion coefficient and the effective diffusion coefficient in porous media D*(m2/s). At any particular location (grid block) with flow velocities in longitudinal (L) and transverse (T) directions, their corresponding diffusion / dispersion coefficients DL (m2/s) and DT (m2/s) are calculated as follows:
Here αL and αT are the longitudinal and transverse dispersivity (m). The dispersion coefficients vary spatially due to the non-uniform distribution of the permeability values.
As will be discussed in lesson 1, in a system with N total number of species and m fast reactions, the total number of primary species is n = N – m. With specific initial and boundary conditions, reactive transport codes solve a suite of n equation (1) with explicit coupling of the physical processes (diffusive/dispersive + advective transport) together with m algebraic equations defined by the laws of mass action of fast reactions. The output is the spatial and temporal distribution of all N species. This type of process-based modeling allows the integration of different processes as a whole while at the same time differentiation of individual process contribution in determining overall system behavior.
The focus of this course is on how to use an existing reactive transport code, not on how to numerically solve reactive transport equations, which deserves a separate course by itself. Readers who are interested in numerical methods of RTM are referred to book chapters and literature for details of discretization of the equations and numerical solution.
The rest of the course is structured as follows. Unit 1 includes lessons that teach principles and set up of geochemical reactions in well-mixed systems (zero-dimension in space, well-mixed systems). In this unit, the equations solved are ordinary differential equations (ODEs) with time as the independent variable however without space dimensions. Lesson 1 focuses on the concepts of primary and secondary species, reaction thermodynamics, and aqueous complexation reactions. Lesson 2 teaches mineral dissolution and precipitation reactions and Transition State Theory (TST) rate law. Lesson 3 is on surface complexation reactions. Lesson 4 teaches ion exchange reactions. Lesson 5 teaches microbe-mediated reactions.
Unit 2 teaches principles and set up of solute transport processes. This is where we introduce the space dimension. Here we solve Advection-Dispersion equation (ADE) without reactions in a one-dimensional system (lesson 6). Lesson 7 teaches principles and set up of heterogeneous, two dimensional domains.
Unit 3 combines biogeochemical reactions in Unit 1 and solute transport processes in Unit 2. Lesson 9 is an example of 1D transport with multiple mineral dissolution and precipitation reactions. If you run the simulation sufficiently long and allow the properties change over time, it becomes a chemical weathering simulation. Lesson 9' intends to have you combine solute transport with microbe-mediated reactions based on lessions 5 and 6. Lesson 10 introduces a 2D system with both physical and geochemical heterogeneities.
CrunchFlow was developed by Carl I. Steefel [2] in the 1990s (Steefel and Lasaga, 1994). It has been continuously evolving with new capabilities and features. Please refer to the CrunchFlow webpage for detailed introduction of code capabilities [3]. Among various existing reactive transport code, we choose to teach CrunchFlow owing to its features that allows each incorporation, computational efficiency with the use of fast solvers from the petsc (Portable, Extensible Toolkit for Scientific Computation) library, and its flexibility in simulating spatailly heterogeneous systems.
CrunchFlow executables: There are several different versions of CrunchFlow executables in Canvas. They are for different operation systems: windows – 32 bit, windows – 64 bit, and Mac. Please download the version that is compatible to the operating system on your PC. For windows executables, you need a dynamic library (.dll file) to run the exe, which is also in the corresponding Canvas folder. The Mac executables do not need a dynamic running library.
CrunchFlow Orientation:
In input file:
Please watch the video (32:46 minutes) Lesson 0 Crunch Flow Orientation that introduces CrunchFlow and its database and input file structure.
Apps, J. A., L. Zheng, Y. Zhang, T. Xu and J. T. Birkholzer (2010). "Evaluation of Potential Changes in Groundwater Quality in Response to CO2 Leakage from Deep Geologic Storage." Transport in Porous Media 82(1): 215-246.
Atchley, A. L., A. K. Navarre-Sitchler and R. M. Maxwell (2014). "The effects of physical and geochemical heterogeneities on hydro-geochemical transport and effective reaction rates." Journal of Contaminant Hydrology 165(0): 53-64.
Bao, C., H. Wu, L. Li, D. Newcomer, P. E. Long and K. H. Williams (2014). "Uranium Bioreduction Rates across Scales: Biogeochemical Hot Moments and Hot Spots during a Biostimulation Experiment at Rifle, Colorado." Environmental Science & Technology 48(17): 10116-10127.
Bao, C., Li, L., Shi, Y. and Duffy, C. (2017) Understanding watershed hydrogeochemistry: 1. Development of RT-Flux-PIHM. Water Resources Research 53, 2328-2345.
Beaulieu, E., Y. Goddéris, D. Labat, C. Roelandt, D. Calmels and J. Gaillardet (2011). "Modeling of water-rock interaction in the Mackenzie basin: Competition between sulfuric and carbonic acids." Chemical Geology 289(1–2): 114-123.
Bethke, C. (1996). Geochemical reaction modeling : concepts and applications. New York, Oxford University Press.
Bolton, E. W., A. C. Lasaga and D. M. Rye (1996). "A model for the kinetic control of quartz dissolution and precipitation in porous media flow with spatially variable permeability: Formulation and examples of thermal convection." Journal of Geophysical Research: Solid Earth 101(B10): 22157-22187.
Boudreau, B. P. (1996). "A method-of-lines code for carbon and nutrient diagenesis in aquatic sediments." Computers & Geosciences 22(5): 479-496.
Brown, B. and D. Rolston (1980). "Transport and transformation of methyl bromide in soils." Soil science 130(2): 68-75.
Brunet, J. L., L. Li, Z. T. Karpyn, B. Kutchko, B. Strazisar and G. Bromhal (2013). "Dynamic evolution of compositional and transport properties under conditions relevant to geological carbon sequestration." Energy & Fuels.
Chapman, B. M. (1982). "Numerical simulation of the transport and speciation of nonconservative chemical reactants in rivers." Water Resources Research 18(1): 155-167.
Chapman, B. M., R. O. James, R. F. Jung and H. G. Washington (1982). "Modelling the transport of reacting chemical contaminants in natural streams." Marine and Freshwater Research 33(4): 617-628.
Druhan, J. L., C. I. Steefel, M. E. Conrad and D. J. DePaolo (2014). "A large column analog experiment of stable isotope variations during reactive transport: I. A comprehensive model of sulfur cycling and delta S-34 fractionation." Geochimica Et Cosmochimica Acta 124: 366-393.
Fang, Y., T. D. Scheibe, R. Mahadevan, S. Garg, P. E. Long and D. R. Lovley (2011). "Direct coupling of a genome-scale microbial in silico model and a groundwater reactive transport model." Journal of Contaminant Hydrology 122(1-4): 96-103.
Godderis, Y., L. M. Francois, A. Probst, J. Schott, D. Moncoulon, D. Labat and D. Viville (2006). "Modelling weathering processes at the catchment scale: The WITCH numerical model." Geochimica et Cosmochimica Acta 70: 1128-1147.
Helgeson, R. C. (1968). "Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions: I. Thermodynamic relations." Geochim. Cosmochim. Acta 32: 853–877.
Helgeson, R. C., R.M. Garrels and F. T. MacKenzie (1969). "Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions: II. Applications,." Geochim. Cosmochim. Acta 33: 455– 482.
Kang, Q. J., P. C. Lichtner and D. X. Zhang (2006). "Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media." Journal of Geophysical Research-Solid Earth 111.
Knutson, C. E., C. J. Werth and A. J. Valocchi (2005). "Pore-scale simulation of biomass growth along the transverse mixing zone of a model two-dimensional porous medium." Water Resources Research 41(7).
Li, L., N. Gawande, M. B. Kowalsky, C. I. Steefel and S. S. Hubbard (2011). "Physicochemical heterogeneity controls on uranium bioreduction rates at the field scale." Environmental Science and Technology 45(23): 9959-9966.
Li, L., C. A. Peters and M. A. Celia (2006). "Upscaling geochemical reaction rates using pore-scale network modeling." Advances in Water Resources 29(9): 1351-1370.
Li, L., C. I. Steefel and L. Yang (2008). "Scale dependence of mineral dissolution rates within single pores and fractures." Geochimica Et Cosmochimica Acta 72(2): 360-377.
Li, L., Bao, C., Sullivan, P.L., Brantley, S., Shi, Y. and Duffy, C. (2017a) Understanding watershed hydrogeochemistry: 2. Synchronized hydrological and geochemical processes drive stream chemostatic behavior. Water Resources Research 53, 2346-2367.
Li, L., Maher, K., Navarre-Sitchler, A., Druhan, J., Meile, C., Lawrence, C., Moore, J., Perdrial, J., Sullivan, P., Thompson, A., Jin, L., Bolton, E.W., Brantley, S.L., Dietrich, W.E., Mayer, K.U., Steefel, C.I., Valocchi, A., Zachara, J., Kocar, B., McIntosh, J., Tutolo, B.M., Kumar, M., Sonnenthal, E., Bao, C. and Beisman, J. (2017b) Expanding the role of reactive transport models in critical zone processes. Earth-Sci. Rev. 165, 280-301.
Lichtner, P. C. (1985). "Continuum Model for Simultaneous Chemical-Reactions and Mass- Transport in Hydrothermal Systems." Geochimica et Cosmochimica Acta 49(3): 779-800.
Lichtner, P. C. (1988). "The Quasi-Stationary State Approximation to Coupled Mass- Transport and Fluid-Rock Interaction in a Porous-Medium." Geochim. Cosmochim. Acta. 52(1): 143-165.
Lichtner, P. C., C. I. Steefel and E. H. Oelkers, Eds. (1996). Reactive transport in porous media. Reviews in Mineralogy, MIneralogical society of America.
Liu, C., J. Shang, H. Shan and J. M. Zachara (2013). "Effect of Subgrid Heterogeneity on Scaling Geochemical and Biogeochemical Reactions: A Case of U(VI) Desorption." Environmental Science & Technology 48(3): 1745-1752.
MacQuarrie, K. T. B. and K. U. Mayer (2005). "Reactive transport modeling in fractured rock: A state-of-the-science review." Earth-Science Reviews 72(3-4): 189-227.
Mayer, K. U., S. G. Benner, E. O. Frind, S. F. Thornton and D. N. Lerner (2001). "Reactive transport modeling of processes controlling the distribution and natural attenuation of phenolic compounds in a deep sandstone aquifer." Journal of Contaminant Hydrology 53(3-4): 341-368.
Mayer, K. U., E. O. Frind and D. W. Blowes (2002). "Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions." Water Resources Research 38(9).
Molins, S., D. Trebotich, C. I. Steefel and C. Shen (2012). "An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation." Water Resources Research 48.
Navarre-Sitchler, A. K., R. M. Maxwell, E. R. Siirila, G. E. Hammond and P. C. Lichtner (2013). "Elucidating geochemical response of shallow heterogeneous aquifers to CO2 leakage using high-performance computing: Implications for monitoring of CO2 sequestration." Advances in Water Resources 53(0): 45-55.
Ortoleva, P., E. Merino, C. Moore and J. Chadam (1987). "Geochemical self-organization, I. Feedbacks, quantitative modeling." Am J Sci 287: 979–1007.
Parkhurst, D. L. and C. A. J. Appelo (1999). User’s guide to PHREEQC (Version 2)—A computer program for speciation, batch reaction, 1 dimensional transport, and inverse geochemical calculation: 310.
Regnier, P., S. Arndt, N. Goossens, C. Volta, G. G. Laruelle, R. Lauerwald and J. Hartmann (2013). "Modelling Estuarine Biogeochemical Dynamics: From the Local to the Global Scale." Aquatic geochemistry 19(5-6): 591-626.
Roelandt, C., Y. Goddéris, M. P. Bonnet and F. Sondag (2010). "Coupled modeling of biospheric and chemical weathering processes at the continental scale." Global Biogeochemical Cycles 24(2): GB2004.
Saunders, J. A. and L. E. Toran (1995). "Modeling of radionuclide and heavy metal sorption around low- and high-{pH} waste disposal sites at {Oak Ridge, Tennessee}." Appl. Geochem. 10(6): 673-684.
Scheibe, T. D., Y. Fang, C. J. Murray, E. E. Roden, J. Chen, Y.-J. Chien, S. C. Brooks and S. S. Hubbard (2006). "Transport and biogeochemical reaction of metals in a physically and chemically heterogeneous aquifer." Geosphere 2(4): doi: 10.1130/GES00029.00021.
Scheibe, T. D., R. Mahadevan, Y. Fang, S. Garg, P. E. Long and D. R. Lovley (2009). "Coupling a genome-scale metabolic model with a reactive transport model to describe in situ uranium bioremediation." Microbial Biotechnology 2(2): 274-286.
Soetaert, K., P. M. J. Herman and J. J. Middelburg (1996). "A model of early diagenetic processes from the shelf to abyssal depths." Geochimica et Cosmochimica Acta 60(6): 1019-1040.
Soler, J. M. and U. K. Mader (2005). "Interaction between hyperalkaline fluids and rocks hosting repositories for radioactive waste: Reactive transport simulations." Nuclear Science and Engineering 151(1): 128-133.
Steefel, C. I., 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. Šimůnek, N. Spycher, S. B. Yabusaki and G. T. Yeh (2015). "Reactive transport codes for subsurface environmental simulation." Computational Geosciences 19(3): 445-478.
Steefel, C. I., D. J. DePaolo and P. C. Lichtner (2005). "Reactive transport modeling: An essential tool and a new research approach for the Earth Sciences." Earth and Planetary Science Letters 240: 539-558.
Steefel, C. I. and A. C. Lasaga (1994). "A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems." Am. J. Sci. 294: 529-592.
Steefel, C. I. and A. C. Lasaga (1994). "A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems." American Journal of Science 294(5): 529-592.
Van Cappellen, P. and Y. F. Wang (1996). "Cycling of iron and manganese in surface sediments: a general theory for the coupled transport and reaction of carbon, oxygen, nitrogen, sulfur, iron, and manganese." Am J Sci 296: 197-243.
White, M. D. and M. Oostrom (2000). STOMP: Subsurface Transport Over Multiple Phases (Version 2.0): Theory Guide. Richland, WA.
Wolery, T. J., K. J. Jackson, W. L. Bourcier, C. J. Bruton, B. E. Viani, K. G. Knauss and J. M. Delany (1990). "CURRENT STATUS OF THE EQ3/6 SOFTWARE PACKAGE FOR GEOCHEMICAL MODELING." Acs Symposium Series 416: 104-116.
Xu, T. F., S. P. White, K. Pruess and G. H. Brimhall (2000). "Modeling of pyrite oxidation in saturated and unsaturated subsurface flow systems." Transport Porous Med 39(1): 25-56.
Yabusaki, S. B., Y. Fang, K. H. Williams, C. J. Murray, A. L. Ward, R. D. Dayvault, S. R. Waichler, D. R. Newcomer, F. A. Spane and P. E. Long (2011). "Variably saturated flow and multicomponent biogeochemical reactive transport modeling of a uranium bioremediation field experiment." Journal of Contaminant Hydrology 126(3-4): 271-290.
Yeh, G.-T. and V. S. Tripathi (1991). "A Model for Simulating Transport of Reactive Multispecies Components: Model Development and Demonstration." Water Resources Research 27(12): 3075-3094.
Yeh, G. T. and V. S. Tripathi (1989). "A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components." Water Resources Research 25(1): 93-108.
Yoon, H., A. J. Valocchi, C. J. Werth and T. Dewers (2012). "Pore-scale simulation of mixing-induced calcium carbonate precipitation and dissolution in a microfluidic pore network." Water Resources Research 48.
In this orientation we have an overview of the history of RTM development and the applications of RTMs in the past decades. We also discussed general reactive transport equations and the physical meaning of different terms. At the end we introduce CrunchFlow, the code we will learn to use in this class.
You have reached the end of Lesson 0! Double-check the to-do list on the Lesson 0 Overview page [15] to make sure you have completed all of the activities listed there. Please download CrunchFlow compatible with your PC, and run the example files following the instructions in the video. Please report any problems if it does not work so we can help you set up.
Chemical species form complexes in the water phase through aqueous complexation reactions. These reactions occur fast and are controlled by reaction thermodynamics. This lesson introduces reaction thermodynamics, equilibrium constants, and how to set up simulations for aqueous speciation / metal complexation reactions in CrunchFlow. We will discuss thermodynamics principles, followed by the concepts of primary and secondary species. A few examples will be illustrated about how to choose primary and secondary species, and setup simulations for aqueous complexation in CrunchFlow.
By the end of this lesson, you should be able to:
To Read (Optional) |
|
|
---|---|---|
To Do |
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
Example 1.1: A closed carbonate system. Imagine we have a closed bottle with clean water except with a head space filled with $\text{C}\text{O}_2$ gas. The system therefore only has carbonate species in the water phase, the pertinent reactions are as follows:
$$\mathrm{H}_2\mathrm{CO}_3^0\Leftrightarrow\mathrm{H}^++\mathrm{HCO}_3^-\quad K_{a1}=\frac{a_{H^+}a_{HCO_3^-}}{a_{H_2CO_3^0}}=10^{-6.35}$$
$$\mathrm{HCO}_3^-\Leftrightarrow H^++\mathrm{CO}_3^{2-}\quad K_{a2}=\frac{a_{H^+}a_{CO_3^{2-}}}{a_{\mathrm{_{ }HCO}_3^-}}=10^{-10.33}$$
$$\mathrm{H}_{2} \mathrm{O} \Leftrightarrow \mathrm{H}^{+}+\mathrm{OH}^{-} \quad K_{w}=a_{H^{+}} a_{O H^{-}}=10^{-14.00}$$
These are the speciation reactions between carbonate species and water species such as $\text{H}^+$ and $\text{OH}^-$. All these reactions are aqueous speciation reactions and occur fast. All Ks are equilibrium constants at standard temperature and pressure conditions. Note that here we are not imposing the condition for charge balance. Reactions in Example 1.1 are examples of aqueous complexation reactions. These reactions occur ubiquitously in water systems and have significant impacts on water chemistry. For example, ligands, dissolved organic carbon (DOC), and other anions can form complexes with cations, including metals, to prevent cations from precipitation, therefore increasing the mobility of cations in water systems. The carbonate complexation and dissociation reactions with $\text{H}^+$ in example 1.1 works as a buffering mechanism to prevent the pH of water from changing rapidly.
Reaction equilibrium constants $K_{eq}$. Here we briefly cover fundamental concepts in reaction thermodynamics. For a more comprehensive coverage, readers are referred to books listed in reading materials [e.g., (Langmuir et al., 1997)]. A chemical reaction transforms one set of chemical species to another. Here is an example,
The reactants A and B combine to form the products C and D. The symbols $\alpha,\beta,\chi,\text{ and }\delta$ are the stoichiometric coefficients that quantify the relative quantity of different chemical species during the reaction in mole basis. That is, $\alpha$ mole of A and $\beta$ mole of B transform into $\chi$ mole of C and $\delta$ mole of D. The Symbol “$\text { É }$” indicates that the reaction is reversible and can occur in both forward and backward directions. According to reaction thermodynamics, the change in the Gibbs free energy $\Delta G_{R}$ during the reaction can be expressed as:
Here, $\Delta G_{R}^{0}$ is the change in reaction free energy under the standard condition (293.15 K and $10^5$Pa); R is the ideal gas constant (1.987 cal/(mol·K)); T is the absolute temperature (K); and $a_{A}, a_{B}, a_{C},a_{D}$ are the activities of the species A, B, C and D, respectively. The derivation of this equation goes to the heart of reaction thermodynamics theory. In any particular aqueous solution, we define ion activity product (IAP) for the reaction (1) as follows:
The reaction reaches equilibrium when the forward and backward reactions occur at the same rate, the point at which the change in the reaction Gibbs free energy $\Delta G_{R}$ reaches zero. At reaction equilibrium, the IAP equals to the equilibrium constant $K_{eq}$:
where $a_{A}, a_{B}, a_{C}, a_{D}$ are the activities of the species A, B, C and D at equilibrium, respectively. In general, larger $K_{eq}$ means larger reaction tendency to the right direction in the equation as written in (1). The $K_{eq}$ is a constant for a particular reaction at any specific temperature and pressure conditions. Although they look very similar, IAP is the ion activity product under any point along the reaction progress, while $K_{eq}$ is the IAP only at one particular point during the reaction, i.e., at equilibrium. We use the saturation index (SI) to compare IAP and $K_{eq}$:
When SI =0, the reaction is at equilibrium; when SI < 0, the reaction proceeds to the right (forward); when SI > 0, the reaction proceeds to the left (backward).
Dependence of $K_{eq}$ on temperature and pressure. Values of $K_{eq}$ are a function of temperature and pressure. According to reaction thermodynamics, the $K_{eq}$ dependence on temperature follows the Van’t Hoff equation under constant pressure:
Here, $\Delta V_{r}^{o}$ is the molar volume change of the reactions under the standard condition (cm3/mol). Typically, the effect of pressure is important when the reaction involves gas (e.g. $\mathrm{CO}_{2(g)}$) and when there is a significant difference in the molar volume of reactants and products.
If the $\Delta V_{r}^{o}$is a constant and does not depend on pressure, the integrated form the above equation is:
Where $K_{eq,1}$ and $K_{eq,2}$ are the equilibrium constants at P1 and P2, respectively. Typically, P1 is at atmospheric condition at 1 bar.
The standard thermodynamic properties (e.g. $H_{r}^{o}$ and $V_{r}^{o}$) of most compounds can be found or calculated from CRC handbooks (e.g., Handbook of Chemistry and Physics (Haynes, 2012), and the NIST Chemistry WebBook [16]. Readers are also referred to the standard geochemical database Eq3/6 for values of equilibrium constants (Wolery et al., 1990).
Biogeochemical reaction systems typically include both slow reactions with kinetic rate laws and fast reactions that are governed by reaction thermodynamics. Kinetic reactions include, for example, mineral dissolution and precipitation and redox reactions. Thermodynamically-controlled reactions are those with rates so fast that the kinetics does not matter for the problem of interest. In geochemical systems, these include, for example, aqueous complexation reactions that reach equilibrium at the time scales of milli-seconds to seconds. For these reactions, the activities of reaction species are algebraically related through their equilibrium constants, or laws of mass action, as shown in equation (4). As such, their concentrations are not independent of each other and should not be numerically solved independently. This necessitates the classification of aqueous species into primary and secondary species. The primary species are essentially the building blocks of chemical systems, whereas the concentrations of secondary species depend on those of primary species through the laws of mass action. As such, with the definition of primary and secondary species, a reactive transport code only needs to numerically solve the number of equations for the species that are independent of each other, which is essentially the number of primary species. The number of primary species is equivalent to the number of components in a system. The concentrations of secondary species can be calculated based on the concentrations of primary species using laws of mass action.
Here we will go through a few examples on how we choose primary and secondary species using the tableau method. Details are referred to in chapter 2, Morel and Hering, 1993.
How do we categorize the primary and secondary species of the system? Here we go through several steps for the choice of primary and secondary species.
We have 3 fast aqueous reactions, which means that we have 3 laws of mass action, as shown in three expressions defining the equilibrium constant for each reaction.
Species | H+ | CO32- | (H2O) |
---|---|---|---|
H+ | 1< | 0 | 0 |
OH- | -1 | 0 | 1 |
H2CO30 | 2 | 1 | >0 |
HCO3- | 1 | 1 | 0 |
CO32- | 0 | 1 | 0 |
We can write all species in terms of H+ and CO32-:
$$\begin{array}{l}
\mathrm{H}^{+}=\mathrm{H}^{+} \\
\mathrm{CO}_{3}^{2-}=\mathrm{CO}_{3}^{2-} \\
\mathrm{OH}^{-}=\left(\mathrm{H}_{2} \mathrm{O}\right)-\mathrm{H}^{+} \\
\mathrm{H}_{2} \mathrm{CO}_{3}=2 \mathrm{H}^{+}+\mathrm{CO}_{3}^{2-} \\
\mathrm{HCO}_{3}^{-}=\mathrm{H}^{+}+\mathrm{CO}_{3}^{2-}
\end{array}$$
Here ${OH}^{-}$, $\mathrm{H}_{2} \mathrm{CO}_{3}{}^{0}$, and ${HCO}_{3}^{-}$ are secondary species.
Species | OH- | HCO3- | (H2O) |
---|---|---|---|
H+ | -1 | 0 | 1 |
OH- | 1 | 0 | 0 |
H2CO3 | -1 | 1 | 1 |
HCO3- | 0 | 1 | 0 |
CO32- | 1 | 1 | -1 |
We can then write all species in terms of primary species:
$$\begin{array}{l} \mathrm{H}^{+}=\left(\mathrm{H}_{2} \mathrm{O}\right)-\mathrm{OH}^{-} \\ \mathrm{OH}^{-}=\mathrm{OH}^{-} \\ \mathrm{H}_{2} \mathrm{CO}_{3}=\mathrm{HCO}_{3}^{-}+\mathrm{H}_{2} \mathrm{O}-\mathrm{OH}^{-} \\ \mathrm{HCO}_{3}^{-}=\mathrm{HCO}_{3}^{-} \\ \mathrm{CO}_{3}{ }^{2-}=\mathrm{HCO}_{3}^{-}-\mathrm{H}_{2} \mathrm{O}+\mathrm{OH}^{-} \end{array}$$Similarly, (H+, HCO3-), (H+, H2CO3), (OH-, CO32-) are also legitimate choices for primary species. However (H2CO3, CO32-), (HCO3-, CO32-), are not. You can practice using these species to write the expression of secondary species.
If we impose the charge balance condition in example 1, we will then an additional algebraic relationship: $$\mathrm{H}^{+}=\mathrm{OH}^{-}+\mathrm{HCO}_{3}^{-}+2 \mathrm{CO}_{3}^{2-}$$ In this case how many primary species will we have? How many secondary species? What are they?
If we add an additional species Ca2+ in the closed carbonate system in example 1, we then have the following reactions in addition to those in Example 1:
$$\begin{array}{l}
\mathrm{CaCO}_{3}^{0} \Leftrightarrow \mathrm{Ca}^{2+}+\mathrm{CO}_{3}^{2-} \\
\mathrm{CaHCO}_{3}^{+} \Leftrightarrow \mathrm{Ca}^{2+}+\mathrm{HCO}_{3}^{-} \\
\mathrm{CaOH}^{+} \Leftrightarrow \mathrm{Ca}^{2+}+\mathrm{OH}^{-}
\end{array}$$
In this subsection we will discuss setting up aqueous complexation calculations in CrunchFlow.
We have a closed carbonate system as an example 1.1 with the total inorganic carbonate concentration (TIC) being 10-3 mol/L. Please answer the following questions:
Assuming this is a dilute system so that the values of activities are the same as concentrations. If we do not impose the charge balance condition, the code solves the following 5 equations for the 5 species questions 0):
$$\begin{array}{l}
K_{a 1}=\frac{C_{H^{+}} C_{H C O_{3}^{-}}}{C_{H_{2} C O_{3}^{0}}}=10^{-6.35} \\
K_{a 2}=\frac{C_{H^{+}} C_{C O_{3}^{2-}}}{C_{H C O_{3}^{-}}}=10^{-10.33} \\
K_{w}=C_{H^{+}} C_{O H^{-}}=10^{-14.00} \\
-\log _{10} C_{H^{+}}=7.0(\mathrm{pH} \text { is } 7) \\
C_{T}=1.0 \times 10^{-3}=C_{H_{2} C O_{3}^{0}}+C_{H C O_{3}^{-}}+C_{C O_{3}^{2}}
\end{array}$$
If charge balance is imposed, with the following equation:
$$C h \text { arg } e \text { balance: } C_{H^{+}}=C_{H C O_{3}^{-}}+2 C_{C O_{3}^{2-}}+C_{O H^{-}}$$
Then either the pH condition or the TIC condition should disappear so we do not over condition they system (6 equations for 5 unknowns).
Please watch the following 37 minute video, Lesson 1, Example 1
We have a closed system with total inorganic carbonate concentration (TIC) being 10-3 mol/L and the total Ca(II) concentration (summation of all Ca-containing species) being 10-4 mol/L. This is the same system as in the example with addition Ca2+, CaHCO3-, CaCO30 species.
Metal complexation in seawater. The seawater composition is as follows:
Component | Conc. (mol/kg) |
---|---|
pH [19] | 8.1 |
Cl− [20] | 0.546 |
Na+ [21] | 0.469< |
Mg2+ [22] | 0.0528 |
SO42− [23] | 0.0282 |
Ca2+ [24] | 0.0103 |
K+ [25] | 0.0102 |
Total Inorganic carbon (TIC) [26] | 0.00206 |
Br− [27] | 0.000844 |
Ba2+ [28] | 0.000416 |
Sr2+ [29] | 0.000091 |
F− [30] | 0.000068 |
Here we assume we have an open carbonate system instead of a closed system. This means H2CO3, or CO2(aq), is in equilibrium with atmospheric CO2 concentration at PCO2 of 10-3.5 atm. Henry’s law constant for CO2 dissolution is KH = 10-1.47 mol/L/atm. Suppose we have a solution with a given pH and H2CO3 is solely from gas dissolution, please answer the following questions:
(Hint: you will need to look into the CrunchFlow manual to know how to set up a solution in equilibrium with a gas phase at given pressure. Look up the table for "Types of Constraints: Aqueous Species" on page 65-66).
Imagine you have two closed bottles. One bottle has pure water at a pH of 7.0 (with only H+, OH-, and water). In the other bottle, you have carbonate water as those in example 1.1 at a pH of 7.0 and a total inorganic carbon concentration of 0.001 mol/L. Both systems are charge balanced. If I add 0.001 mol/L of Ca(OH)2, what is the new pH of the two systems when each system reaches their new equilibrium? How does the presence of carbonate species influence pH changes? Think ahead how the two systems might be different before you do the calculation, and check if your calculation confirms your hypothesis.
HW1 files and solution package [31]
Haynes, W.M. (2012) CRC handbook of chemistry and physics. CRC press.
Langmuir, D., Hall, P. and Drever, J. (1997) Environmental Geochemistry. Prentice Hall, New Jersey.
Wolery, T.J., Jackson, K.J., Bourcier, W.L., Bruton, C.J., Viani, B.E., Knauss, K.G. and Delany, J.M. (1990) CURRENT STATUS OF THE EQ3/6 SOFTWARE PACKAGE FOR GEOCHEMICAL MODELING. Acs Symposium Series 416, 104-116.
Haynes, W.M., 2012. CRC handbook of chemistry and physics. CRC press.
Langmuir, D., Hall, P., Drever, J., 1997. Environmental Geochemistry. Prentice Hall, New Jersey.
Wolery, T.J., Jackson, K.J., Bourcier, W.L., Bruton, C.J., Viani, B.E., Knauss, K.G., Delany, J.M., 1990. Current Status of the EQ3/6 Software Package for Geochemical Modeling. Acs Symposium Series 416, 104-116.
This lesson introduces reaction kinetics of mineral dissolution and precipitation and how to set up simulation in well-mixed batch reactors in CrunchFlow.
By the end of this lesson, you should be able to:
References (Optional) |
|
---|---|
To Do |
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
The natural subsurface is composed of rocks, soils, as well as other forms of porous materials that contain various types of minerals. Minerals dissolve and precipitate when interacting with water. As minerals dissolve, chemicals in the solid phase transform into ions in water, resulting in a decrease in mineral mass and volume. For example, calcite (CaCO3) dissolves into carbonate species (H2CO3, HCO3-, CO32-) and Ca(II)-containing species in water. In contrast, mineral precipitation occurs when aqueous species transform to become solid phases, therefore leading to mass and volume increase in solid phases while decrease in aqueous concentrations.
Mineral dissolution and precipitation are important in both natural and engineered processes. They influence water chemistry, soil formation, contaminant transport and fate, acid stimulation in reservoir engineering, environmental remediation, and global carbon cycling. For example, acid stimulation accelerates mineral dissolution, therefore increasing reservoir porosity and permeability and enhancing oil production. On the other hand, mineral precipitation during hydrocarbon production results in wellbore clogging. The dissolution of carbonate caprocks can potentially result in CO2 leakage from CO2 storage reservoirs. Alteration of flow fields induced by mineral dissolution and precipitation in geothermal systems could affect the long-term production of geothermal energy. Over geological time scale, the CO2 consumption during chemical weathering (through mineral dissolution) helps maintain the clement conditions for life on earth.
Calcite is one of the most common minerals on Earth's surface. Here we take calcite dissolution as an example to illustrate the Transition State Theory (TST) based kinetic rate law. The dissolution of calcite has been proposed to occur via three parallel reactions (Plummer, Wigley et al. 1978, Chou, Garrels et al. 1989):
Each reaction pathway has its own reaction rate dependence on different chemicals. The rate of reaction (1) dominates under acidic conditions and depends on the activity of hydrogen ion. The rate of reaction (2) dominates under CO2- rich conditions, while the rate of reaction (3) prevails under neutral pH conditions. The overall dissolution rate R (mol/s) is the summation of the rates of all three parallel reactions:
Here , A is the reactive surface area A (m2) k1, k2, and k3 are reaction rate constants (mol/m2/s) for the three parallel reactions, respectively; and are the activities of hydrogen ion and carbonic acid: IAP and Keq are the ion activity products and the reaction equilibrium constants for reactions (1), (2), and (3), respectively.
The equilibrium constants determine mineral solubility and depend on temperature, pressure, and salinity in ways similar to those discussed for aqueous complexation reactions. As indicated by Equation (4), the reaction rates depend on several factors, including the intrinsic mineral properties such as the amount of reactive surface area A and the intrinsic rate constants k, as well as external conditions such as the concentration of “catalyzing” species including H+ and H2CO30, and how far away the reactions are from equilibrium (IAP/Keq).
The rate dependence on pH is illustrated in Figure 1. The rates in the figure are normalized by the amount of surface area and therefore are in the units of mol/m2/s. Under low pH conditions with very fast dissolution rates, the dissolution rates can be transport-controlled even in well-mixed reactors, because the speed of mixing may not be as fast as dissolutin. In contrast, under closer to neutral pH conditions, the rates are much slower and the reactions are often kinetics controllled.
A general TST rate law is as follows:
Here nk is the total number of parallel reactions, kj is the rate constant of the parallel reaction , the term describes the rate dependence on pH, and the term describes the rate dependence on other aqueous species that potentially accelerate or limit reactions. In Equation (4) for calcite dissolution rates, also accelerates calcite dissolution. The affinity term quantifies the distance of solution from the equilibrium state of the mineral reaction j. The exponents m1,j and m2,j describe the nonlinear rate dependency on the affinity term and are normally measured in experiments. Note that Equation (4) is a special case of the general reaction rate law in Equation (5).
Various factors affect reaction rates, including, for example, temperature, salinity, pH, organic ligands. Mineral dissolution rates are often accelerated by the presence of H+ or OH-. As a result, many minerals dissolve fast in acidic and alkaline conditions and slow down under neutral conditions. Figure 2 shows silicate dissolution rates as a function of pH.
The intrinsic dissolution rates of minerals vary significantly. Calcite dissolution is in the fast end of dissolution rate spectrum, while quartz dissolves the slowest. Table 1 shows the lifetime of a 1mm crystal of different minerals, indicating more than 5 orders of magnitude difference in dissolution rates of silicates and quartz.
Mineral | Lifetime |
---|---|
Quartz | 34,000,000 |
Muscovite | 2,700,000 |
Forsterite | 600,000 |
K-feldspar | 520,000 |
Albite | 80,000 |
Enstatite | 8,800 |
Diopside | 6,800 |
Nepheline | 211 |
Anorthite | 112 |
The Arrhenius equation is used to quantify the reaction rate dependence on temperature.
Here AF is the pre-exponential factor and is usually considered as a constant that is independent of temperature: the activation energy is always positive; here R is the ideal gas constant , T is the absolute temperature (K). If we take the logarithm of this equation, we obtain:
This means if we draw lnk versus , we get a straight line, as shown in Figure 2. For a particular reaction, if and at the temperature T1 is known, at another temperature T2 can be calculated by
Figure 3 shows the dependence of calcite rate constants on temperature under different CO2 partial pressure.
To model the mineral dissolution/precipitation, we need to inform the code both reaction thermodynamics and kinetics. Here we introduce how to set up simulations for mineral dissolution/precipitation in a well-mixed batch reactor, where the concentrations of all species and temperature are assumed to be spatially uniform. That is, the rate of mixing is faster than the rate of dissolution or precipitation, which is in general true for most minerals except carbonates under very low pH conditions. No transport processes and boundary conditions need to be defined because of the assumption of uniform concentration. The example 6 included in the CrunchFlow package is also a good reference for kinetic reaction setup.
Calcite dissolution proceeds through three parallel pathways as shown in reactions (1)-(3). In addition, a series of fast aqueous complextion / speciation reactions occur simultaneously. These fast reactions are similar to those in the examples in the lesson on Aqueous Complexation.
In a batch reactor of 200 ml, the initial solution is at a pH of 5.0 with close to zero salinity. The initial total carbonate concentration is . The kinetic rate constant for three parallel reactions are , and respectively. The volume fraction of the calcite grains in the batch reactor is 0.5%, with a specific surface area (SSA) of 0.24 m2/g. Please set up the simulation and plot the following quantities as a function of time:
Here is a light board video (12:44) of what equations the code solves for this system.
PRESENTER: So in this video we're going to talk about mineral dissolution and precipitation, the equation involved and interaction involved. This is different from the previous lesson, aqueous [INAUDIBLE], because the previous reaction only involved reaction [INAUDIBLE] dynamics. Meaning we care about the end point of the reaction.
When we have mineral dissolution precipitation, because reaction usually occur at the interface of solid phase and water, reaction occur much slower. So we often need to consider the kinetics of the reaction. Meaning we care about how long it takes for the reaction to occur to reach the end point.
So what I'm going to do today here is using the carbonate dissolution system in a battery reactor. So battery reactor, meaning it's closed while mixed, so that it does not have constitution gradient in different parts of reactor. So essentially you're solving for one constitution for each species for the whole reactor. So that's our systems.
If you think about if you want to draw something relevant to what you do in the lab, but you think about these water, and then you have these calcite grains, and you have these mixers that will keep the system well-mixed. And again, the reaction involving this system is similar to what we talked about last time for the aqueous speciation. The carbon and water-- for example, all these three reactions are the same as what I wrote before. Invest in one. But the key thing is we are adding this calcite dissolution, which is calcite dissolving out to become calcium and carbonate.
Calcite or carbonate in general is a very common type of rock or mineral fissure or surface. So it's a reaction. Compared to the aqueous speciation or [INAUDIBLE], it's much slower. But compared to other type of mineral dissolution, it's actually reactive fast.
So we have these four reactions here. You can think about the calcite dissolution essentially kind of releasing out the elements from the solid phase the water phase. So it dissolve as an add to the water with calcium species and carbonate species. But then once carbonate is released into water, it quickly goes through the speciation reaction to become either bicarbonate or carbonic acid, depending on the pH a system, as we talked about last time.
So if you think about the species in this system, now that the species we have, in addition to the five species, carbonic acid, bicarbonate, carbonate-- this is wrong. So it should be carbonate. And this should be two.
So you have these three carbonate species. You have hydrogen ion, which minus acid five, as before. But on top of that, you are adding calcium.
So essentially, now you have six species. You would have Ca2 plus carbonic acid carbonate, bicarbonate carbonate. And then you also have H plus or H minus. So that's six species.
So what does that mean, is that we essentially have six unknowns to solve for. And we will need six equations to solve for this, to solve for the species. But notice that one things that we need to pay attention to is that we have these three reactions that occur fast.
So you have these three equations already there. We can think about as one, two, three. That is as a three relationship. But we cannot use this one, because this is not a faster reaction. These are fast reaction.
This is a slow reaction, which means we need to care about its kinetics. Like how far it release as these chemical species over time. So when we solve these equations, we actually need to consider these time component.
In this case, we will need to think about the mass balance of system. So these calcite dissolving out and release out. For example, Ca2 plus, and releasing out CO3 bicarbonate. And then this will be exchanged into bicarbonate and exchange into carbonic acid. So these three specie are exchangeable.
So when we think about mass balance of system, you would think about the constitution of a single body calcium first. Volume-- this is a volume of the water in the system. And you have the constitution of calcium.
Ultimately, we are solving for contrition, not activity. So this will give you the mass. So contrition of the unit and mass per volume. And then we will need to solve for order and differential equation now because just time component here. So in this kind of system we don't really have fro transport so the only source of zinc of this calcium specie is through these reactions.
So let's say we have a rate constant of this reaction is Ksp. So that's the rate constant. And we talk about the tst rate, which is transition state rate theory in the online material. So you can go look this up.
So this Ksp and this is rate constant times a surface area of the mineral. And then times for example, how far away the reaction is from equilibrium, which will depends on the activity of Ca2 plus activity of carbonate, and then divided by Ksp, which is how far away this is from equilibrium. So this is Iap.
Now when we have dilute system-- if we have dilute system, then this activity would be equal to c. So you can replace every a is c for simplicity. So essentially, is this equation saying mass of calcium is added to system, to the water, by having this rate lower for calcium. But we also will have adding total carbon to the system. So similarly, you will have V, d, and then we call it total carbon cT, dt.
It will have the same expression here as Ksp A or minus activity ca. Because, essentially, it's through the same reaction. So this rack essentially adding calcium total carbonate-- this CT will be equals to again, constitution of carbonic acid plus concentration of bicarbonate and concentration of carbonate.
But then we need one more equation. One, two, three. And then this is four, this is five. We need a six equation, which would it be-- if you think about it, when this mineral dissolving out to form calcium and carbonate, it's actually changing the pH of the system. So the hydrogen ion-- or somehow we also call it's a measure of acidity of system will also change.
I'm not going through the detailed derivation of the whole process how we come up with. But let's call this dC, which is-- so this a constitutional cT-- concentration of hydrogen ion total dt. And by dissolving out calcite, actually decrease acidity of system. So it will actually have minus Ksp. And then activity one minus a Ca2 plus.
So again, it's the same rate law, but it do actually have the minus signs. That meaning it's decreasing the total acidity of a system. So the expression for the total acidity will be different if we define different to comnination primer specie.
But in general, if we define this carbonate as a primary species, we should have HT equal to hydrogen, concentration hydrogen ion. This should be CHt. Hydrogen ion plus concentration of H2CO3 minus concentration of bicarbonate and minus concentration of OH minus. You can actually derive this equation from the [INAUDIBLE] method. And these expressions is actually-- you can look at [INAUDIBLE] And by combining different terms, you come up with this. So essentially, we have six unknowns, and you have six equation to solve with that. But on top of that, one thing we need to pay attention is that these are OD equations.
So we call it ordinary differential equation. Differential equation, meaning we have one independent variable, which is time in this equation. So typically, the code will be solving these equation first using time stepping.
And then once we get the concentration calcium, concentration of CT, concentration of CHt, you have three kind of variable in particular time step. And then we solve for the algebraic and relationship with these numbers together with these three relationship. And by doing that, we essentially get the concentration of all the six species as a function time in these aqueous species. And what do you see? For example, your homework is essentially generated by solving these equations.
Here is a video demo in CrunchFlow on how to set up the simulation in CunchFlow (43:51 minutes).
Thermodynamics and kinetic parameters that affect mineral dissolution include specific surface area (SSA), kinetic rate constants (k), and equilibrium constants. In addition, geochemical conditions can also have a large impact on mineral dissolution rates, including salinity (e.g. NaCl) and pH. With the provided CrunchFlow template files in example 1 as a starting point, please do the following analysis comparing Ca(II) concentration evolution under different parameters and geochemical conditions. In each question, please only change the parameter that is discussed and keep all other parameters the same as those in Example 1.
1) Calcite specific surface area (SSA). In three different simulations, run the code using SSA being 0.024, 0.24, and 2.4 m2/g. Plot the total Ca(II) as a function of time under these surface area values in one figure. Discuss how SSA affects calcite dissolution kinetics.
2) Kinetic rate constant: increase and decrease the original three rate constant values in Example 1 by an order of magnitude. Compare the Ca(II) concentration evolution under the three k values in one figure. Also compare this figure here with the figure in 1) with the SSA of 0.24 m2/g. Do changing k and A have the same impact on dissolution rates?
3) Equilibrium constant Keq quantifies how much a mineral can dissolve in aqueous phase. Increase and decrease the equilibrium constant by an order of magnitude. Please plot the three curves (total Ca(II) ~ t) with different Keq values in one figure. What do the rate kinetics change with the total Ca(II) evolution figure?
4) Initial pH condition. Compare the base case with two more cases where you have the initial pH being 4.0 and 6.0.
5) The role of speciation. In this question you remove all secondary species such as CaOH+, CaCO3(aq), and CaHCO3+. Draw the same figures as those in Example 1 and solution from this question if the same figure. Compare and describe your observations. Does the reaction system evolve in the same way in these two systems? What do you think are the effects of speciation?
6) Salinity: in the base case scenario there is no NaCl in the solution. Simulate two more cases with NaCl at concentrations of 0.01 and 0.1 mol/kgw, respectively. Please compare total Ca(II) evolution under these three salinity conditions.
7) Go back to your answer to each question 1)-6). Summarize your observations in these questions. Which factors have the most significant control on calcite dissolution?
Feldspars dissolve following the TST rate law and have very different dissolution rates depending on their composition (figure 7 in [Blum and Stillings, 1995]). Pick a feldspar mineral from the paper and/or the cited reference that has a complete TST rate law from acidic to alkaline conditions. Set it up in CrunchFlow to simulate its dissolution with initial pH conditions varying from acidic (for example, pH =4.0), to neutral (pH = 7), and to alkaline condition (pH = 10.0). Compare the reaction production evolution under these three conditions. Please make sure that you use reasonable rate constant and specific surface area values from literature.
In this lesson, we learned about mineral dissolution and precipitation reactions. We discussed two aspects of the reactions: 1) reaction thermodynamics that control how much minerals can dissolve in aqueous solutions; 2) mineral reaction kinetics (TST rate law) and parameters that control the rates of mineral dissolution and precipitation. We also learned how to set up running simulations in CrunchFlow for mineral dissolution in well-mixed reactors.
You have reached the end of Lesson 2! Double-check the to-do list on the Lesson 2 Overview page to make sure you have completed all of the activities listed there before you begin Lesson 3.
Homework assignments for lesson 2 is due on Tuesday before midnight.
Appelo, C. A. J., and D. Pstma (2005), Geochemistry, groundwater, and pollution, 2nd ed., 649 pp., CRC Press, Taylor and Francis group, Boca Raton, FL.
Blum, A. E., and L. L. Stillings (1995), Feldspar dissolution kinetics, in Chemical Weathering Rates of Silicate Minerals, edited by A. F. White and S. L. Brantley, pp. 291-351, Mineralogical Soc America, Washington, D. C.
Brantley, S. L. (2008), Kinetics of mineral dissolution, in Kinetics of water-rock interaction, edited, pp. 151-210, Springer.
Chou, L., R. M. Garrels, and R. Wollast (1989), Comparative study of the kinetics and mechanisms of dissolution of carbonate minerals, Chem Geol, 78(3), 269-282.
Lasaga, A. C. (1984), Chemical kinetics of water-rock interactions, Journal of Geophysical Research, 89(B6), 4009-4025.
Plummer, L., T. Wigley, and D. Parkhurst (1978), The kinetics of calcite dissolution in CO 2-water systems at 5 degrees to 60 degrees C and 0.0 to 1.0 atm CO 2, Am J Sci, 278(2), 179-216.
Pokrovsky, O. S., S. V. Golubev, J. Schott, and A. Castillo (2009), Calcite, dolomite and magnesite dissolution kinetics in aqueous solutions at acid to ciarcumneutral pH, 25 to 150° C and 1 to 55 atm p CO2: New constraints on CO2 sequestration in sedimentary basins, Chem Geol, 265(1), 20-32.
This lesson introduces general principles of surface complexation reactions, as well as how to set up surface complexation models in well-mixed batch reactors in CrunchFlow.
By the end of this lesson, you should be able to:
To Read |
|
---|---|
To Do |
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
Sorption is the adhesion of chemicals to solid surfaces. Adsorption process occurs in many natural and engineered systems. Studies of contaminated systems have shown that sorption–desorption is an important geochemical process that regulates transport and fate of inorganic and organic contaminants in natural subsurface systems. For example, metals (Cd2+, Cr3+, Co2+, Cu2+, Fe3+, Pb2+ or Zn2+) can become immobilized by sorbing on sediments and soils. They can also become mobilized through desorption from the solid surface and re-enter the aqueous phase when geochemical conditions allow. Sorption-desorption is widely used in industrial applications including charcoal activation, air conditioning, water purification, among others.
Sorption can occur either specifically or non-specifically as shown in Figure 1. Chemical sorption (Specific adsorption) is highly selective and occurs only between certain adsorptive and adsorbent species. A chemical bond involves sharing of electrons between the adsorbate and adsorbent and may be regarded as the formation of inner-sphere surface complexes. Chemical adsorption is difficult to reverse because of the strength of the formed bond. Physical adsorption (nonspecific adsorption). A physical attraction resulting from nonspecific, relatively weak Van der Waal's forces. Being only weakly bound, physical adsorption is easily reversed. Multiple layers form through outer-sphere surface complexes during physical adsorption [Goldberg 1991; Webb, 2003].
Sorption via surface complexation has been extensively studied. Surface complexation is the process where species in the aqueous phase form complexes with functional groups on solid surfaces, similar to aqueous complexation in lesson 1. Surface complexation function occurs between aqueous species and functional groups on solid surface, instead of the formation of complexes between aqueous species in aqueous complexation reactions. Surface complexation models use mass action laws that are analogous to aqueous geochemical conditions and solid phase properties.
Surface complexation models describe sorption based on surface reaction equilibrium. Similar to aqueous complexation, surface complexation reactions are considered fast reactions and are controlled by reaction thermodynamics.
There are three commonly used SCMs, the constant capacitance model (CCM), the diffuse layer model (DLM), and the triple layer model (TLM). These models differ in complexity from the simplest CCM that has three adjustable model parameters, to the most complex TLM that has seven adjustable parameters [Hayes et al., 1991]. The double layers exist in practically all heterogeneous fluid-based systems. Here we introduce the principle and thermodynamics of DLM.
The double layer refers to two parallel layers of charge surrounding the solid surface. The first layer, the surface charge (either positive or negative), comprises ions sorbed onto the solid due to chemical interactions. The second layer (“diffuse” layer) is composed of counter ions attracted to the surface charges via the coulomb force, electrically screening the first layer. The schematic of double layer is shown in Figure 2.
In traditional SCM models, all reactions are considered as at equilibrium. As an example, the surface protolysis reactions, where H+ transfers among chemicals, are given by:
Here the $\equiv S O H$ represents a species or functional group on solid surface. In the first reaction, $\equiv S O H$ gains an H+ and becomes positively charged. In the second reaction, $\equiv S O H$ loses an H+ and becomes negatively charged. The apparent equilibrium constants Kapp describe the relationship between activities of different species, written in the same format as we do for aqueous complexations except now we include activities of solid species.
Similarly, for a metal ion M with a positive charge m, the reactions are represented by:
For an anionic ligand L with a negative charge of l, the reactions are represented by:
In these reactions, $\equiv S O H$ represents a surface site for the sorbent functional group, $\equiv S O H_{2}^{+}, \equiv S O^{-}, \equiv S O M^{(\mathrm{m}-1)}, \equiv(S O)_{2} M^{(\mathrm{m}-2)}, \equiv S O H_{2}^{+}-L^{l-} \text { and } \equiv S O H_{2}^{+}-L H^{(l-1)-}$ are surface complexes, [ ] represents the activity of each species or surface complex, $Mm+$ represents a metal ion of charge m+ and Ll- represents an anionic ligand of charge l−.
The apparent equilibrium constant of surface complexation reactions, Kapp, is an important parameter because it determines the ion partition between aqueous and solid phases. Large Kapp values indicate high affinity of the ions to the solid surface. The relationship between total Gibbs free energy $\Delta G_{\text {tot }}$ and Kapp is as follows:
Here the total Gibbs free energy $\Delta G_{\text {tot }}$ can be further expressed as follows:
Here $\Delta G_{c h e m}^{0}$ is the intrinsic free energy of the chemical reactions at the surface; $\Delta G_{\text {coul }}^{0}$ is the electrostatic or Coulombic term that accounts for the electrostatic interactions:
Here Z is the charge of the ion, F is the Faraday constant (96485 C/mol), $\psi_{0}$ is the average potential of the surface plane (V). Therefore,
Where $K^{\mathrm{int} r}=\exp \left(-\frac{\Delta G_{c h e m}^{0}}{R T}\right)$, R is ideal gas constant $(1.987 \mathrm{cal} /(\mathrm{mol} \cdot \mathrm{K}))$; T is the absolute temperature (K). Take equation (1) as an example,
Where $K_{1}^{\mathrm{int} r}=\frac{\left[\equiv \mathrm{SOH}_{2}^{+}\right]}{[\equiv \mathrm{SOH}]\left[\mathrm{H}_{s}^{+}\right]}$ , Z is the charge of the ion (1 in the case of H+); $\left[H_{s}^{+}\right]$ is H+ activity on the solid surface.
The electrostatic or coulombic effect can be quantified as:
From equation (10), we know that $K_{1}^{\text {int } r}$ and $\psi_{0}$ are needed in order to calculate $K_{1}^{a p p}$. The $K_{1}^{\text {int } r}$ is typically estimated using zero charge extrapolation or using the double extrapolation method as discussed in literature. Under low ionic strength conditions where $\psi_{0} \cong 0$, the intrinsic and apparent constants are equivalent.
For different minerals, the number of surface sites differs significantly, depending on their surface properties. The abundance of surface sites is important in determining the total sorption capacity. The concentration of surface sites can be calculated as follows:
Here Csite is the concentration of surface sites (mol/g mineral), ρsites is the surface density of surface hydroxyl sites (mol/m2), Aspecific is the specific surface area (SSA)(m2/g). This equation says that surface site concentration depends on the surface site density and specific surface area. Please note that if you are working with porous media, you will need to calculate the total gram of minerals for surface complexation to get the total number of available sites.
Specific surface area (SSA) and site density values can be determined experimentally from Brunauer -Emmett-Teller (BET) surface area and tritium exchange measurements, respectively. The units of site/nm is often used in literature, where 1 site/nm = 1.66x10-6 mol/m2. Typical values of specific surface area (SSA) and site densities for different types of minerals are listed in Table 1. The total number of surface sites for a particular system (mol) can be calculated by multiplying site density (mol/m2) with SSA (m2/g) and the mineral mass (g). Minerals such as clays tend to have a large surface area and have a large capacity to sorb chemicals.
Absorbent | SAA (m2/g) | Site density (mol/m2) Strong Site |
Site density (mol/m2) Weak Site |
Reference |
---|---|---|---|---|
Goethite | 14.7 | 1.76×10-6 | 3.22×10-6 | (Müller and Sigg, 1992) |
Kaolinite | 19.5 | 2.20x10-6 | 3.00×10-6 | (Lackovic et al., 2003) |
Illite | 66.8 | 1.30x10-6 | 2.27x10-6 | (Gu and Evans, 2007) |
Smectite | 56.4 | 4.77x10-8 | 9.54x10-7 | (Bradbury and Baeyens, 2005) |
Organic and inorganic chemicals are usually sorbed at hydroxyl surface functional groups that are located at the broken bonds and edge sites on minerals with excess negative charges [Baeyens and Bradbury, 1997]. We often classify two kinds of sorption sites: "strong" sites $\left(\equiv \mathrm{S}^{\mathrm{S}} \mathrm{OH}\right)$ and "weak" sites$\left(\equiv \mathrm{S}^{\mathrm{W}} \mathrm{OH}\right)$. "Strong" sites have a low capacity and a high sorption affinity and dominate the uptake of adsorbate at low concentrations. "Weak" sites have a considerably larger capacity however much lower sorption affinity. Table 2 shows reactions and equilibrium constants for U(VI) sorption on ferrihydrite, where FesOH represents strong site with orders of magnitude higher intrinsic equilibrium constants than those of the weak sites $\left(\equiv \mathrm{FE}^{\mathrm{W}} \mathrm{OH}\right)$ (Zheng et al., 2003). In this soil with the presence of ferryhdrate, the site density ratio of weak to strong site is 476:1 (i.e., 0.21% of the total surface sites, 99.79% for the weak sites).
Reactions | LogKintr |
---|---|
$2 \equiv \mathrm{Fe}^{S} \mathrm{OH}+\mathrm{UO}_{2}^{2+} \Leftrightarrow \equiv\left(\mathrm{Fe}^{S} \mathrm{O}\right)_{2} \mathrm{UO}_{2}+2 H^{+}$ | -2.35 |
$2 \equiv \mathrm{Fe}^{W} \mathrm{OH}+\mathrm{UO}_{2}^{2+} \Leftrightarrow \equiv\left(\mathrm{Fe}^{W} \mathrm{O}\right)_{2} \mathrm{UO}_{2}+2 H^{+}$ | -6.06 |
$2 \equiv \mathrm{Fe}^{S} \mathrm{OH}+\mathrm{UO}_{2}^{2+}+\mathrm{CO}_{3}^{2-} \Leftrightarrow \equiv\left(\mathrm{Fe}^{S} \mathrm{O}\right)_{2} \mathrm{UO}_{2} \mathrm{CO}_{3}^{2-}+2 \mathrm{H}^{+}$ | 4.33 |
$2 \equiv \mathrm{Fe}^{\mathrm{W}} \mathrm{OH}+\mathrm{UO}_{2}^{2+}+\mathrm{CO}_{3}^{2-} \Leftrightarrow \equiv\left(\mathrm{Fe}^{W} \mathrm{O}\right)_{2} \mathrm{UO}_{2} \mathrm{CO}_{3}^{2-}+2 H^{+}$ | -0.24 |
Surface complexation leads to surface-charged solid surfaces. Electric surface charges govern characteristic chemical and physical phenomena such as ion exchange, adsorption, swelling, colloidal stability, and flow behavior (Sposito, 1981). It is well known that the surface charges on layered silicates and insoluble oxides depend on the pH of aqueous solutions The pH of the point of zero charge (PZC), where the net total particle charge is zero, is a convenient reference for describing the pH dependence of surface charges. (Appel et al., 2003). When solution pH is above PZC, the solid surface has a negative charge and predominantly exhibits an ability to exchange cations, while the solid surface retains anions (electrostatically) if pH is below its PZC. A list of common substances and their associated PCZs is shown in Table 3.
Chemical Formula | pHPZC |
---|---|
$\text { Kaolinite }$ | 4.6 |
$\text { Montmorillonite }$ | < 2.5 |
$\text { Corundum, } \alpha-\mathrm{Al}_{2} \mathrm{O}_{3}$ | 9.1 |
$\gamma-\mathrm{Al}_{2} \mathrm{O}_{3}$ | 8.5 |
$\text { alpha- } \mathrm{Al}(\mathrm{OH})_{3}$ | 5.0 |
$\text { Hematite, } \alpha-\mathrm{Fe}_{2} \mathrm{O}_{3}$ | 8.5 |
$\text { Goethite, } \alpha \text {-FeOOH }$ | 9.3 |
$\text { Birnessite, } \delta-\mathrm{MnO}_{2}$ | 2.2 |
$\mathrm{Fe}(\mathrm{OH})_{3}$ | 8.5 |
$\text { Quartz, } \mathrm{SiO}_{2}$ | 2.9 |
$\text { Calcite, } \mathrm{CaCO}_{3}$ | 9.5 |
Example 3.1: Cr(VI) surface complexation on illite. Chromium is a common containment in natural subsurface due to its natural occurrence and wide industrial usage, including electroplating, pigmenting, and dye synthesis. Anionic Cr (VI) is highly mobile and poses a tremendous risk to human and ecosystem health. Clay minerals such as illite are important in controlling Cr(VI) sorption and natural attenuation due to its large surface area and ubiquitous presence (Wang and Li, 2015).
We have an initial solution listed in Table 4. The illite grains in the solution have specified surface area and site density of $\equiv \mathrm{SiOH}$. The surface site $\equiv \mathrm{SiOH}$ goes through several surface complexation reactions as listed in Table 4. Please calculate:
Initial conditions (total concentrations) | Value |
---|---|
Temperature | 25oC |
Solution volume | 250 mL |
pH | 8.0 |
CrO42- | 9.61x10-5mol/L |
Na+ | 0.01 mol/L |
Cl- | 0.01 mol/L |
K+ | 1.93x10-4mol/L |
Al3+ | 1.00× 10-6 mol/L |
Mg2+ | 1.00× 10-6mol/L |
SiO2(aq) | 1.00× 10-5mol/L |
Site density $\equiv \mathrm{SiOH}$ | 1.00× 10-6mol/L |
Illite specific surface area | 15.36 m2/g |
Illite volume fraction | 0.003 |
Reactions | Log Kapp |
$\equiv \mathrm{SiOH}+\mathrm{H}^{+} \Leftrightarrow \equiv \mathrm{SiOH}_{2}^{+}$ | 0.95 |
$\equiv \mathrm{SiOH} \Leftrightarrow \mathrm{SiO}^{-}+\mathrm{H}^{+}$ | -6.59 |
$\equiv \mathrm{SiOH}+\mathrm{Na}^{+} \Leftrightarrow \equiv \mathrm{SiONa}+\mathrm{H}^{+}$ | -6.60 |
$\equiv \mathrm{SiOH}+\mathrm{CrO}_{4}^{2-}+2 H^{+} \Leftrightarrow\left(\equiv \mathrm{SiOH}^{0}-\mathrm{H}_{2} \mathrm{Cr} O_{4}^{0}\right)^{0}$ | 14.50 |
Before setting up the simulations in CrunchFlow, let's think about how to represent this sytem, a well-mixed reactor, in a mathematical form, how many chemical species do we have, how many algebraic relationships that we have, and how many equations we need to solve. Please watch the following video (13:12).
Presenter: We are going to show an example of Surface Complexation reaction today. This is somewhat similar to one of the previous lessons on Aqueous Complexation. But the difference is really, now we have solid phase and complexations are being formed between aqueous species and the species on a solid phase. So what I have here is example 3.1. You also have that in the online material. So this is the example, if we think about a system that you have well-mixed again a batch reactor.
So well-mixed meaning all the concentration in the water phase will remain the same. It's uniform. So we don't really solve for concentration difference in different parts of the batch reactor. Now in this system we have Illite grains, which is a very common type of clay. And then we have the water that has this chromium 6 (CrVI) on it. And we know that this species will solve our surface complexes with species on Illite. So what we have here is these grains and then this water.
But also there's some background species, like sodium chloride, that's providing some salinity. And Illite itself will be slowly dissolving up. So there are some other species, for example, magnesium silica. We'll talk about that later. So in order for solve for system, we think about this system again. Surface complexation is usually considered also a very fast reaction, similar to aqueous complexation. So we can usually think about the thermodynamics of these reactions instead of kinetics of these reactions.
So let's just go over the chemistry of the system. So first of all, we have these reactions, right? And we think about this as there's both reactions happening in water phase and also at the interface of water and solid. So in the water phase we actually will simplify the system to only include, for example, the water dissociation to become hydrogen ion (H+) and hydroxide (OH-). This is a reaction that must be there. But also it includes a chromium related reaction.
Chromium 6 can have three different forms. You have H2CrO4 (aq) can become H+ and then this species $$HCrO_4^-$$. And this can further dissociate to have hydrogen ion (H+) and this $CrO_4^-$ form. So in the water phase we are actually, there could be a lot of other reaction happening. But for simplicity we would only include these three.
So that's for the water phase reactions. And also at the water and solid interface, we're really talking about water and Illite interfaces. We have these solid species, like surface species, right? So this, if you look at the form, we are kind of using this to represent a solid surface. And then you have the SiOH as a functional group on the solid surface.
So this surface specie can react with hydrogen ion (H+) to form this. And also dissociate, the hydrogen ion comes out to become this. $\equiv \mathrm{SiOH}_{2}^{+} / \equiv \mathrm{SiOH} +\mathrm{H}^{+}=$, But also there's, for example, when there's sodium in the water and chromium in the water, they can also form these surface complexes. So you probably notice that in the different reactions here, these reactions, we write, for example, the same type of laws of mass action like in aqueous complexation.
So we had this activity of hydrogen ion, activity OH- is timed together, equal to Kw. And similar for chromium A1, A2, right? I'm not writing everything out. Because this is all in similar form. You write activity of species in the right side of the reaction divided by activity of species in the left side of the reaction. So obviously the K's are constants. So again, we have these three reactions. But also then we have 4, 5, 6, 7. So we have three aqueous phase reactions and a another four reactions that occur at the water and solid interface.
And each of them you can have these expressions of laws of mass action, which I'm not detailing out. But also, as I mentioned, the Illite itself would dissolve slowly. So in control we actually also have this reaction in the background, except that it occurs so slowly it doesn't change much of the chemistry of the system. But when we set up the control, to input a file, we do need to have these reactions, these chemical species there. That is actually part of the Illite. But it's not really explicitly talked about in this aqueous phase.
So these are the chemical components of Illite that we have to put these there as primary species. Now, if we think about this sort of system, so we have this many reactions. And we think about how many different chemical species, we have, if we just list them out, you have H, of course, hydrogen ion, OH-.
And then you have three chromium related species, which was different hydrogen ion there, $CrO_4^-$. And then you have these solid. We also need to solve the concentration for the solid species as well, right? So you have these $\equiv \mathrm{SiOH}_{2}^{+}$, SiOH self, SiO-. And then this forms complexly to have SiOHNa or SiO. And then you also have these OH with $\mathrm{H}_{2} \mathrm{CrO}_{4}^{0}$. So these are the five possible or potential surface complexes that can be formed. Now on top of that, you also have, for example, sodium chloride Na+,Cl-. And then the chemical composition of these Illites, right? So you have Mg2+, potassium K+, aluminum, SiO2, aqueous.
OK, so let's count. We need to solve for all these different species, right? So let's count this, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. And then you have another 6. So in total we have 16 species, including all the possible aqueous species and solid species or surface species. So we have 16 species. That means that we are going to solve for 16 unknowns. And we already have 16 unknowns. Because we have 16 species.
Now we already know we have seven different relationships, 1, 2, 3, 4, 5, 6, 7, right? So these are the reaction we specify. And we know every time we specify one reaction, there's an algebraic relationship relative to that. So we have 16 minus 7 equal to 9. So we have 16 unknowns. 7 we know the relationship meaning these concentrations, all activity are dependent on each other, this relationship. That would mean we need to specify 9 additional conditions for completely solving these reactions.
So what we can do here is, for example, a lot of times we know PH. So these conditions should be already given to you. And what we have, for example, typically, let's say, we know PH. Then we should know the activity of hydrogen ion. Or we know the question should give the total concentration of chromium 6. And this would be equal to, for example, concentration of CrO42- plus concentration of $\mathrm{H}_2\mathrm{CrO}_4^{ }$ plus concentration of $\mathrm{H}_{2} \mathrm{CrO}_{4}^{0}$ aqueous, right? And this together should be equal to whatever constant they gave to you. Which, I'm not writing those.
And you also should know, it should also give you a concentration of sodium, give the concentration of chloride, and give you a concentration of potassium, aluminum, magnesium, and SiO2(aq). Another condition they should give you is also, how much total site you have. So this would be something like so this is a total site on the solid surface should be the concentration of all these five potential species adding together, right?
So you can think about this as total sites, C sites. And this would be adding all the surface complex species. For example, $\mathrm{SiOH}_{2}^{+}$ plus CSiOH plus CSiO- and then the Csi of sodium plus CsiOHH2CrO4. So these are the five different surface complexes that can be formed. And this should be equal to a constant total concentration of sites on the Illite grain. And the total concentration of sites should be equal to, for example, how much Illite grain you have, how many grams. And also the site density, times the site density times the surface area.
So these should be conditions you should have. So if you look at this, you have 7, and then 8, 9, 10, 11, 12, 13, 14, 15, 16. So this is close to form, right? You have 16 unknowns. You have 16 relationships to-- You have seven relationships, but you also know nine conditions that specify the system that you can solve for the whole system.
Now what you end up with is the concentration of each species, both aqueous and solid at equilibrium. Because the system reaches equilibrium really quickly. So essentially you have 7 relationships and then 9 conditions to completely solve the concentration of all species involved in the system.
Here are the equations and key points.
Setting up a simulation for surface complexation involves both input and database file. Relevant reading materials on surface complexation in CrunchFlow includes keywords on pages 63, 64, 69.
In the input file, the keyword block for surface complexation is the SURFACE_COMPLEXATION block. Complexation must occur on a specific mineral, so a valid mineral name (listed in the MINERALS keyword block) must be given in the MINERAL keyword block as well. An example:
SURFACE_COMPLEXATION
$\equiv \mathrm{SiOH}$ on Illite
END
Here the $\equiv \mathrm{SiOH}$ is a surface site on the mineral Illite. The mineral must be present in the database. To specify a non-electrostatic model, the mineral name should be followed by the keyword –no_edl. For example:
$\equiv \mathrm{SiOH}$ on Illite -no_edl
The term “-no_edl” means no electrical double layer.
In the database file, you need to specify the surface complexation reactions in Table 4 in the “Begin Surface Complexation” section. In addition, you need to specify charges of the surface species in “Begin Surface Complexation parameters” section.
The exercise 4 in the CrunchFlowExampleExercise is also for surface complexation.
If you try to set up in Phreeqc, Phreeqc manual includes the introduction of surface complexation calculation and the key words such as SURFACE, SURFACE_MASTER_SPECIES, SURFACE_SPECIES. Example 8 in Phreeqc is a good reference for setting up surface complexation reactions.
If instead we have two types of surface sites on illite surface with the reactions and parameters shown in Tables 5. All other conditions remain the same as in example 3.1.
Initial conditions | Value |
---|---|
Temperature | 25oC |
Solution volume | 250 mL |
pH | 8.0 |
Total CrO4-- | $9.61 \times 10^{-5} \mathrm{~mol} / \mathrm{L}$ |
Na+ | 0.01 mol/L |
Cl- | 0.01 mol/L |
K+ | $18.5 \times 10^{-5} \mathrm{~mol} / \mathrm{L}$ |
Site density $\equiv \mathrm{SiOH}$ | $1.0 \times 10^{-6} \mathrm{~mol} / \mathrm{m}^{2}$ |
Site density $\equiv A l O H$ | 0.1 10-6 mol/m2 |
Illite specific surface area | 15.36 m2/g |
Illite volume fraction | 0.003 |
Reactions | Log K |
---|---|
$\equiv\mathrm{SiOH}+\mathrm{H}^+\Leftrightarrow\ \equiv\mathrm{SiOH}_2^+$ | 0.95 |
$\equiv\mathrm{SiOH}\quad\Leftrightarrow\ \equiv\mathrm{SiO}^-+\mathrm{H}^+$ | -6.59 |
$\equiv\mathrm{SiOH}+\mathrm{Na}^+\Leftrightarrow\ \equiv\mathrm{SiONa}+\mathrm{H}^+$ | -6.60 |
$\equiv \mathrm{SiOH}+\mathrm{CrO}_{4}^{2-}+2 \mathrm{H}^{+} \Leftrightarrow\left(\equiv \mathrm{SiOH}^{0}-\mathrm{H}_{2} \mathrm{CrO}_{4}^{0}\right)^{0}$ | 14.50 |
$\equiv\mathrm{AlOH}+\mathrm{H}^+\Leftrightarrow\ \equiv\mathrm{AlOH}_2^+$ | 5.70 |
$\equiv AlOH\quad\Leftrightarrow\ \equiv AlO^-+H^+$ | -11.40 |
$\equiv AlOH+Na^+\Leftrightarrow\ \equiv AlONa+H^+$ | -9.15 |
$\equiv\mathrm{AlOH}+\mathrm{Cl}^-+\mathrm{H}^+\Leftrightarrow\ \equiv\mathrm{AlOH}_2\mathrm{Cl}$ | 7.90 |
$\equiv \mathrm{AlOH}+\mathrm{CrO}_{4}^{2-}+\mathrm{H}^{+} \Leftrightarrow\left(\equiv \mathrm{AlOH}_{2}^{+}-\mathrm{CrO}_{4}^{2-}\right)^{-}$ | 9.42 |
$\equiv \mathrm{AlOH}+\mathrm{CrO}_{4}^{2-}+2 \mathrm{H}^{+} \Leftrightarrow\left(\equiv \mathrm{AlOH}_{2}^{+}-\mathrm{HCrO}_{4}^{-}\right)^{0}$ | 16.30 |
Set up a batch reactor model for AsO43- sorption on Fe(OH)3 given initial conditions and parameters in Tables 6. Run the simulation to understand how sorbed concentrations are affected by different parameters and geochemical conditions.
Initial conditions | Value |
---|---|
Temperature | 25oC |
pH | 7.0 |
AsO4--- | 0.005 mol/L |
Na+ | 0.001 mol/L |
Fe+++ | 0 |
PO4--- | 0.001 mol/L |
Site density $\equiv \mathrm{FeOH}$ | $1.0 \times 10^{-6} \mathrm{~mol} / \mathrm{m}^{2}$ |
Specific surface area | 50 m2/g |
Fe(OH)3 volume fraction | 0.2 |
Reactions | Log K |
---|---|
$\equiv\mathrm{FeOH}+\mathrm{H}^+\Leftrightarrow\ \equiv\mathrm{FeOH}_2^+$ | 5.10 |
$\equiv \mathrm{FeOH} \Leftrightarrow \quad \equiv \mathrm{FeO}^{-}+\mathrm{H}^{+}$ | -10.70 |
$\equiv\mathrm{FeOH}+\mathrm{Na}^+\Leftrightarrow\ \equiv\mathrm{FeO}-\mathrm{Na}+\mathrm{H}^+$ | -9.00 |
$\equiv \mathrm{FeOH}+\mathrm{AsO}_{4}^{3-}+\mathrm{H}^{+} \Leftrightarrow\left(\equiv \mathrm{FeOAsO}_{3}\right)^{2-}+\mathrm{H}_{2} \mathrm{O}$ | 16.6 |
$\equiv \mathrm{FeOH}+\mathrm{PO}_{4}^{3-}+\mathrm{H}^{+} \Leftrightarrow\left(\equiv \mathrm{FeOPO}_{3}\right)^{2-}+\mathrm{H}_{2} \mathrm{O}$ | 16.9 |
Surface complexation occurs ubiquitously in natural and engineered systems. Understanding surface complexation reaction is important to understand and predict reactive transport and fate of chemicals (such as Cr, As, Cd , Cu, Pb) in natural waters, soils, and sediments. Here we introduce the mechanisms and importance, and controlling parameters of surface complexation reactions, and how to set up the model for ion exchange reactions in CrunchFlow.
Competitive Sorption and Transport of Heavy Metals in Soils and Geological Media.by H. M. Selim, 2012. Chapter 2. Equilibrium and kinetic modeling of competitive heavy metals sorption and transport in soils.
Kinetics of Water-Rock Interaction. Brantley, Susan, Kubicki, James, White, Art. 2008. Chapter 4. Kinetics of sorption-desorption.
Aquatic surface chemistry: chemical processes at the particle-water interface. Stumm, Werner. New York : Wiley 1987.
You have reached the end of Lesson 3! Double-check the to-do list on the Lesson 3 Overview page to make sure you have completed all of the activities listed there before you begin Lesson 4. [35]
Appel, C., Ma, L.Q., Dean Rhue, R. and Kennelley, E. (2003) Point of zero charge determination in soils and minerals via traditional methods and detection of electroacoustic mobility. Geoderma 113, 77-93.
Appelo, C.A.J. and Pstma, D. (2005) Geochemistry, groundwater, and pollution, 2nd ed. CRC Press, Taylor and Francis group, Boca Raton, FL.
Baeyens, B. and Bradbury, M.H. (1997) A mechanistic description of Ni and Zn sorption on Na-montmorillonite Part I: Titration and sorption measurements. Journal of Contaminant Hydrology 27, 199-222.
Bradbury, M.H. and Baeyens, B. (2005) Modelling the sorption of Mn(II), Co(II), Ni(II), Zn(II), Cd(II), Eu(III), Am(III), Sn(IV), Th(IV), Np(V) and U(VI) on montmorillonite: Linear free energy relationships and estimates of surface binding constants for some selected heavy metals and actinides. Geochimica et Cosmochimica Acta 69, 875-892.
Goldberg, S. (1991) Sensitivity of surface complexation modeling to the surface site density parameter. Journal of Colloid and Interface Science 145, 1-9.
Goldberg, S., Criscenti, L.J., Turner, D.R., Davis, J.A. and Cantrell, K.J. (2007) Adsorption - Desorption processes in subsurface reactive transport modeling. Vadose Zone Journal 6, 407-435.
Gu, X. and Evans, L.J. (2007) Modelling the adsorption of Cd(II), Cu(II), Ni(II), Pb(II), and Zn(II) onto Fithian illite. Journal of Colloid and Interface Science 307, 317-325.
Guswa, A.J., Celia, M.A. and Rodriguez-Iturbe, I. (2002) Models of soil moisture dynamics in ecohydrology: A comparative study. Water Resour. Res. 38, 1166.
Hayes, K.F., Redden, G., Ela, W. and Leckie, J.O. (1991) SURFACE COMPLEXATION MODELS - AN EVALUATION OF MODEL PARAMETER-ESTIMATION USING FITEQL AND OXIDE MINERAL TITRATION DATA. Journal of Colloid and Interface Science 142, 448-469.
Huang, Y., Ma, X.Y., Liang, G.Z., Yan, Y.X. and Wang, S.H. (2008) Adsorption behavior of Cr(VI) on organic-modified rectorite. Chem. Eng. J. 138, 187-193.
Lackovic, K., Angove, M.J., Wells, J.D. and Johnson, B.B. (2003) Modeling the adsorption of Cd(II) onto Muloorina illite and related clay minerals. Journal of Colloid and Interface Science 257, 31-40.
Li, W., Tang, Y.K., Zeng, Y.T., Tong, Z.F., Liang, D.W. and Cui, W.W. (2012) Adsorption behavior of Cr(VI) ions on tannin-immobilized activated clay. Chemical Engineering Journal 193, 88-95.
Müller, B. and Sigg, L. (1992) Adsorption of lead(II) on the goethite surface: Voltammetric evaluation of surface complexation parameters. Journal of Colloid and Interface Science 148, 517-532.
Riese, A.C. (1982) Adsorption of Radium and Thorium onto Quartz and Kaolinite: A Comparison of Solution/Surface Equilibria Model. Ph.D. Dissertation, Colorado School of Mines 292p.
Serrano, S., O'Day, P.A., Vlassopoulos, D., Garcia-Gonzalez, M.T. and Garrido, F. (2009) A surface complexation and ion exchange model of Pb and Cd competitive sorption on natural soils. Geochimica Et Cosmochimica Acta 73, 543-558.
Sposito, G. (1981) THE OPERATIONAL DEFINITION OF THE ZERO-POINT OF CHARGE IN SOILS. Soil Sci Soc Am J 45, 292-297.
Um, W., Serne, R.J., Brown, C.F. and Rod, K.A. (2008) Uranium(VI) sorption on iron oxides in Hanford Site sediment: Application of a surface complexation model. Applied Geochemistry 23, 2649-2657.
Wang, L. and Li, L. (2015) Illite Spatial Distribution Patterns Dictate Cr(VI) Sorption Macrocapacity and Macrokinetics. Environmental Science & Technology 49, 1374-1383.
Webb, P.A. (2003 ) Introduction to Chemical Adsorption Analytical Techniques and their Applications to Catalysis. MIC Technical Publications.
Zachara, J.M., Cowan, C.E., Schmidt, R.L. and Ainsworth, C.C. (1988) Chromate adsorption by kaolinite. Clay Clay Miner 36, 317-326.
Zeng, H., Fisher, B. and Giammar, D.E. (2008) Individual and competitive adsorption of arsenate and phosphate to a high-surface-area iron oxide-based sorbent. Environmental Science & Technology 42, 147-152.
Zheng, Z.P., Tokunaga, T.K. and Wan, J.M. (2003) Influence of calcium carbonate on U(VI) sorption to soils. Environmental Science & Technology 37, 5603-5608.
This lesson introduces general principles of ion exchange reactions and setup simulations for ion exchange reactions in well-mixed batch reactors using CrunchFlow.
By the end of this lesson, you should be able to:
Recommended Readings |
|
|
---|---|---|
To Do |
|
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
Ion exchange reactions occur when ions in water exchange with those electrostatically bound to the solid phase. They occur commonly in the presence of iron oxides, organic matters, and clay minerals with large surface area. Ion exchange reactions are important in determining natural composition in surface waters (e.g., rivers, lakes) and ground water. They can alter water composition and trigger other reactions including mineral dissolution and precipitation. For example, in the coastal freshwater aquifer, Na+ in seawater displaces presorbed Ca2+ from solid phase (Figure 1). In this reaction, Ca2+ is replaced by Na+ and the water changes from Na-rich to Ca-rich (Appelo and Willemsen, 1987; G., 2008; Slomp and Van Cappellen, 2004) while the solid surface change from Ca-rich to Na-rich. In northeastern United States, the use of road salt as deicer increases the salinity of fresh water and mobilize metals through ion exchange reactions (Kaushal et al., 2005).
Applications of ion exchange reactions in industry include drinking water softening (Figure 2), desalination, ultra-pure water production (Rodrigues, 1986), and chromatography. In petroleum refining processes, ion exchange reactions are often used to purify, separate, and dry natural gases (Marinsky and Marcus, 1995). Natural gas extraction in the Marcellus shale has led to the production of large quantity of wastewater with high metal concentrations. Ion exchange is commonly used to remove metals such as Ba2+ and Sr2+ (Gregory et al., 2011).
Ion exchange reactions occur when ions exchange their positions at the soid-surface – water interface. These reactions are typically assumed reversible and occur instantaneously at time scales ranging from microseconds to hours. Ion exchange is typically represented in the following form:
Here (s) and (aq) refer to solid and aqueous phases, respectively; X- denotes the negatively charged surface sites that bound cations Au+ and Bv+, with u and v being the charges of A and B, respectively. In this reaction, A and B are cations that compete for sorption sites (Appelo and Postma, 1993; Sposito et al., 1981; Vanselow, 1932). The equilibrium constant Keq of reaction (1) can be expressed as
Here the parentheses [ ] represent activities. Aqueous concentrations are easily related to activities through concentration and activity coefficients. Activities of ions on solid phases are typically expressed as a fraction of the total, either as molar fraction or as equivalent fractions. The total number can be based on the number of exchange sites or as the number of exchangeable cations.
The units meq is often used in ion exchange calculations. An meq is the number of ions that sums a specific quantity of electrical charges. For example, an meq of K+ is about 6.02 x 1020 positive charges. On the other hand, an meq of Ca2+ is also 6.02 x 1020 positive charges, however only 3.01 x 1020 ions because Ca2+ has two positive charges. For an ion Ii+ with a charge of i, the equivalent fraction bI is calculated as:
where I, J, K are exchangeable cation with charges i, j, k, respectively. A molar fraction $\beta_{I}^{M}$ is obtained by the following form:
Here TEC is the total exchangeable cations in mmol/kg sediment, not cation exchange capacity. The use of fractions should give the summation of fractions being 1, that is, .
Three common conventions used in writing ion exchange equilibrium constants are Gaines-Thomas convention, Gapon convention, and Vanselow convention. If the ion exchange is between cations of the same valence (homovalent exchange), the convention does not make a difference. If the exchange is between cations of different valences (heterovalent), the convention makes a difference. For example, the ion exchange reaction between Na and Ca can be written as follows:
With
Here if the equivalent fraction of the exchangeable cations is used for $β$ values, the Gaines-Thomas convention is followed (Graines and Thomas, 1953). If we use molar fractions for $β$ values, we follow the Vanselow convention (Vanselow, 1932). If the activities of cations on exchange sites are expressed as a fraction of the number of exchange sites (X-), we follow the Gapon convention, the reaction will be written as follows:
Where
Here the activities are expressed in terms of the mole fraction of the total number of exchangeable sites. To fully understand the difference between different conventions and calculation, please go over Example 6.3 in chapter 6 of Appelo and Postma (2005).
The capabilities of ions to compete for exchange sites are governed by their affinity to the surface of the exchangers. Selectivity coefficients have been reported in literature for common ions including Na+, K+, Ca2+, and Mg2+, however rarely for trace metals. The larger the selectivity coefficient, the higher affinity of the ion to the solid phase and the more competitive for exchange. In general, ions have high affinity to exchange sites when they have higher valence, are less solvated with water molecules, and react strongly with the surface sites. The following Table 1 lists selectivity coefficients reported in literatures for ion exchange on kaolinite (Appelo and Postma, 2005; Bundschuh and Zilberbrand, 2011).
No. | Ion exchange reaction | log K |
---|---|---|
1 | $2 \mathrm{Na}-X+\mathrm{Mg}^{2+} \Leftrightarrow \mathrm{Na}^{+}+\mathrm{Mg}-X_{2}$ | 0.60 (0.44 ~ 0.78) |
2 | $2 \mathrm{Na}-\mathrm{X}+\mathrm{Ca}^{2+} \Leftrightarrow \mathrm{Na}^{+}+\mathrm{Ca}-X_{2}$ | 0.80 (0.44 ~ 0.104) |
3 | $2 \mathrm{Na}-X+\mathrm{Ba}^{2+} \Leftrightarrow \mathrm{Na}^{+}+\mathrm{Ba}-X_{2}$ | 0.91 (0.44 ~ 0.104) |
4 | $2 \mathrm{Na}-X+\mathrm{Sr}^{2+} \Leftrightarrow \mathrm{Na}^{+}+\mathrm{Sr}-X_{2}$ | 0.91 (0.44 ~ 0.104) |
For typical ion exchangers, the sequence of affinity is as follows: Ba2+> Pb2+ > Sr2+ > Ca2+ > Ni2+ > Cd2+ > Cu2+ > Co2+> Zn2+ > Mg2+ > Ag+> Cs+ > K+ >NH4+ >Na+ >H+. When there are multiple cations co-exist in the solution, the surface exchange site composition is largely determined by the cation aqueous concentration and their affinity to the exchange sites. To calculate the exchange site composition, that is, the mole fraction of each cation on the exchange site, please follow the example 6.4 in chapter 6 of Appelo and Postma (2005).
The cation exchange capacity (CEC, meq/kg solid) is a measure of the solid phase capacity for ion exchange reactions (Meunier, 2005). The CEC of different porous media are very much associated with their clay content, organic carbon, and grain size. Different materials have different CEC values. CEC values of clay minerals such as muscovite, illite, kaolinite, and chlorite are high for their grain sizes smaller than 2 μm (Drever, 1982). In general, organic matter has the highest CEC values (1500 - 4000 meq/kg). Iron oxides play a vital role in natural processes and controls nutrient availability and heavy-metal mobility (Houben and Kaufhold, 2011). Iron oxides, such as goethite and hematite, have CEC values from 40 to 1000 meq/kg. Many solid materials have iron oxides or organic matters coated on their surface and therefore have large CEC values. For soils, CEC value is a function of solution pH depending on hydrolysis reactions of surface sites. In general, cation exchange occurs due to the broken bonds around the crystal edges, the substitutions within the lattice, and the hydrogen of exposed surface hydroxyls that may be exchanged. Higher pH values give rise to more negative charges on clay, resulting in higher CEC. CEC values increase as the grain size decreases due to the large surface area associated with smaller grains.
Minerals | Grain size (μm) | Surface area (m2/g) | CEC (meq/kg) |
---|---|---|---|
Kaolinite | 0.1-5.0 | 5-20 | 30-150 |
Illite | 0.1-2.0 | 15-40 | 150-400 |
Montmorillonite | 0.01-1.0 | 600-800 | 800-1200 |
Marcellus Shale flow back and produced waters typically contains high concentrations of metals, including Na(I), Ca(II), Mg(II), Sr(II), and Ba(II) (Chapman et al., 2012). The release of Marcellus Shale waste water into natural soils and sediment can lead to ion exchange reactions and partition of surface sites between major cations. Here we set up a batch experiment in CrunchFlow to model the major cations exchange on clay surface from Marcellus shale produced water.
Assume that we have the following Marcellus water in a 300 ml batch reactor. The water is in equilibrium with partial pressure of CO2 at 3.15×10-4 atmosphere. Please note that the concentrations are in ppm and in mol/kgw. Both can be used directly in CrunchFlow.
Component | Conc. (ppm) | Molar (mol/kgw) |
---|---|---|
pH | 7.02 | |
Co2 (aq) | Co2(g), 3.15e-4 atm | |
Na(I) | 7900.0 | 3.43e-1 |
Ca(II) | 2774.0 | 6.94e-2 |
K(I) | 82.5 | 2.11e-3 |
Mg(II) | 239.0 | 9.83e-3 |
Ba(II) | 183.0 | 1.34e-3 |
Sr(II) | 6.5 | 7.38e-5 |
C1(-I) | charge |
“Charge” means that the ion concentration is calculated from charge balance.
The water is in equilibrium with a sediment with kaolinite being the major clay with a specific surface area of 13.9 m2/g and a CEC value is 100.0 meq/kg (the units in CrunchFlow is eq/g). The kaolinite occupies 0.005 of the total volume.
The selectivity coefficient of the cations are as follows in the database (values from (Appelo and Postma, 1993; Li et al., 2010)):
No. | Ion exchange reaction | log K |
---|---|---|
1 | $\mathrm{Na}X\Leftrightarrow\mathrm{Na}^++X^-$ | 0.0 |
2 | $\mathrm{K} X \Leftrightarrow \mathrm{K}^{+}+X^{-}$ | -0.69 |
3 | $\mathrm{Ca} X_{2} \Leftrightarrow \mathrm{Ca}^{2+}+2 X^{-}$ | -0.39 |
4 | $\mathrm{Mg} X_{2} \Leftrightarrow \mathrm{Mg}^{2+}+2 X^{-}$ | -0.30 |
5 | $\mathrm{Ba} X_{2} \Leftrightarrow \mathrm{Ba}^{2+}+2 X^{-}$ | -0.45 |
6 | $\operatorname{Sr} X_{2} \Leftrightarrow \mathrm{Sr}^{2+}+2 X^{-}$ | -0.45 |
To set up CrunchFlow simulation please read CrunchFlow manual page 62 – 63 for details. Please also refer to the Exercise 5: Cs exchange on Hanford sediments.
You will need to set up the ION_EXCHANGE key word block in the input file to define the name of the exchange site and give parameters related to ion exchange in the CONDITION block. In addition, you will also need to define exchange reactions and equilibrium constants in the “Begin Exchange” block. For example, if you call your exchange site “XKao”, then you have the following block:
ION_EXCHANGE
exchange XKao- on Kaolinite
convention Vanselow
END
The first line specifies name of the exchange site and the mineral that the exchange is on. The second line specify the calculation convention that is used.
As an extension of the example 4.1, we can look at different parameters and understand how they change the cation concentrations on the exchange site. Please include all Cl- and OH- aqueous complexes as in example 4.1. Please calculate the major cation concentrations on the exchange site (K, Na, Ca, Mg, Ba, Sr) by changing only the parameter discussed in each sub-question, with all other parameters being the same as those in example 4.1. For each sub questions, please plot major cation concentrations on the exchange site as the changing parameter.
For each question, calculate the mole fraction of each species on the exchange site, and the mole fraction of each species compared to its own original total mass. Please make a table and a figure comparing different species. Which one has the largest percentage of its own mass on exchange sites? Why?
Which parameter(s) have the largest impact on ion exchange reactions? In lesson 3 we learned that pH is important for surface complexation reactions. Does it make a difference for ion exchange reactions?
Read the paper (Valocchi et al., 1981). Native groundwater in the injection test of this paper had the composition of $Na^+=\ 86.5,\mathrm{\ Mg}^{2+}=18.2,\text{ and }\mathrm{Ca}^{2+}=11.1\mathrm{\ mmol}/\mathrm{L}$. Injected water has 14.66 meq Cl-/L. Selectivity coefficient were (Gaines and Thomas convention, activity = molal concentration) $\mathrm{K}_{\text{Na\Mg}}=0.54\text{ and }\mathrm{K}_{\text{Na\Ca}}=0.41$. Sediment $\mathrm{CEC}=750\mathrm{\ meq}/\mathrm{L}$ pore water. Calculate
Ion exchange is an important chemical process in natural and engineered systems. Understanding ion exchange reaction is important for us to understand and predict the reactive transport and fate of chemicals in nature. Here we introduce the reaction thermodynamics of ion exchange, and the set up of ion exchange simulation in a well-mixed batch reactor in CrunchFlow.
(2008) St. Johns River Water Management District,. WaterWays: More about Water Below the Ground.
Appelo, C. and Willemsen, A. (1987) Geochemical calculations and observations on salt water intrusions, I. A combined geochemical/minxing cell model. Journal of Hydrology 94, 313-330.
Appelo, C.A.J. and Postma, D. (1993) Geochemistry, Groundwater and Pollution. A. A. BALKEMA PUBLISHERS
Appelo, C.A.J. and Postma, D. (2005) Geochemistry, groundwater and pollution. CRC Press.
Bundschuh, J. and Zilberbrand, M. (2011) Geochemical Modeling of Groundwater, Vadose and Geothermal Systems. CRC Press.
Chapman, E.C., Capo, R.C., Stewart, B.W., Kirby, C.S., Hammack, R.W., Schroeder, K.T. and Edenborn, H.M. (2012) Geochemical and Strontium Isotope Characterization of Produced Waters from Marcellus Shale Natural Gas Extraction. Environmental Science & Technology 46, 3545-3553.
Drever, J.I. (1982) The geochemistry of natural waters. Prentice-Hall, Englewood Cliffs, N.J.
G., C.J. (2008) Overexploitation and Contamination of Shared Groundwater Resources. Springer Netherlands.
Gregory, K.B., Vidic, R.D. and Dzombak, D.A. (2011) Water management challenges associated with the production of shale gas by hydraulic fracturing. Elements 7, 181-186.
Houben, G. and Kaufhold, S. (2011) Multi-method characterization of the ferrihydrite to goethite transformation. Clay Minerals 46, 387-395.
Kaushal, S.S., Groffman, P.M., Likens, G.E., Belt, K.T., Stack, W.P., Kelly, V.R., Band, L.E. and Fisher, G.T. (2005) Increased salinization of fresh water in the northeastern United States. P Natl Acad Sci USA 102, 13517-13520.
Li, L., Steefel, C.I., Kowalsky, M.B., Englert, A. and Hubbard, S.S. (2010) Effects of physical and geochemical heterogeneities on mineral transformation and biomass accumulation during biostimulation experiments at Rifle, Colorado. Journal of Contaminant Hydrology 112, 45-63.
Ma, C. and Eggleton, R.A. (1999) Cation exchange capacity of kaolinite. Clays and Clay Minerals 47, 174-180.
Marinsky, J.A. and Marcus, Y. (1995) Ion Exchange and Solvent Extraction: A Series of Advances. CRC Press.
Meunier, A. (2005) Clays. Springer, New York, Berlin.
Rodrigues, A.E. (1986) Ion exchange: science and technology. Springer.
Selim, H.M. (2012) Competitive Sorption and Transport of Heavy Metals in Soils and Geological Media. CRC Press.
Slomp, C.P. and Van Cappellen, P. (2004) Nutrient inputs to the coastal ocean through submarine groundwater discharge: controls and potential impact. Journal of Hydrology 295, 64-86.
Sposito, G., Holtzclaw, K.M., Johnston, C.T. and Levesquemadore, C.S. (1981) Thermodynamics of Sodium-Copper Exchange on Wyoming Bentonite at 298-Degrees-K. Soil Science Society of America Journal 45, 1079-1084.
Valocchi, A.J., Street, R.L. and Roberts, P.V. (1981) TRANSPORT OF ION-EXCHANGING SOLUTES IN GROUNDWATER - CHROMATOGRAPHIC THEORY AND FIELD SIMULATION. Water Resources Research 17, 1517-1527.
Vanselow, A.P. (1932) The utilization of the base exchange reaction for the determination of activity coefficients in mixed electrolytes. Journal of the American Chemical Society 54, 5.
In this lesson we will learn about the thermodynamics and kinetics of microbe-mediated redox reactions, the biogeochemical redox ladder in natural systems, and setting up microbe-mediated reactions in well-mixed batch reactors.
By the end of this lesson, you should be able to:
References (optional) |
|
---|---|
To Do |
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
Microorganisms are everywhere in the Earth surface environment. Their presence has significant influence in the natural and engineered environments. They play a key role in the biogeochemical cycling of elements such as carbon (C), oxygen (O) and nitrogen (N). For example, primary production involves photosynthetic microorganisms that transform CO2 in the atmosphere into organic (cellular) materials of vegetation in terrestrial systems. Planktonic algae and cyanobacteria on the other hand, are the “grass of the sea”. They account for ~ 50% of the primary production on the planet and are the primary carbon source for marine life. The fixation of N from the N2 gas also occurs through N-fixing microorganisms. Microorganisms are believed to be critical in the origin of life and in the transformation of rocks to soils that create life-accommodating environments on Earth. In the modern time when humans generated large quantities of waste and contaminants, microorganisms have overtaken the daunting tasks of cleaning up and recycling wastes in engineered systems such as wastewater treatment plants. They have also been indispensable in remediating contaminated water, soils, and aquifers.
For any redox reactions to occur, we need an electron donor and electron acceptor. The oxidation state of the electron donor increases during redox reactions, whereas that of the electron acceptor decreases. In natural environments, organic carbon often serves as the electron donor in microbe-mediated reactions and become oxidized from organic carbon to inorganic form (e.g., CO2), while multiple electron acceptors often co-exist, including oxygen, nitrate, manganese, iron oxides, among others. There are typically multiple functioning microbial groups that use different electron acceptors.
Microorganisms typically do not use multiple electron acceptors simultaneously. Instead they use electron acceptors in the order of biogeochemical redox ladder, as shown in Figure 1 for some of the naturally-occurring electron acceptors. Different redox reactions release different amount of energy. For example, aerobic oxidation reactions generate much more energy than other redox reactions. The oxidation of glucose can produce 2,880 kJ / mol of C6H12O6; sulfate reduction can produce 492 kJ /mol of C6H12O6, as shown in Figure 2, a much smaller amount. The microorgnanisms that use high-energy-output redox reactions therefore has growth advantage and can outgrow other species.
Reactions | Chemical Formula | Free Energy kJ/mol glucose |
---|---|---|
Aerobic Oxidation | $\mathrm{C}_{6} \mathrm{H}_{12} \mathrm{O}_{6}+6 \mathrm{O}_{2} \rightarrow 6 \mathrm{CO}_{2}+6 \mathrm{H}_{2} \mathrm{O}$ | -2,880 |
Denitrification | $5\ \mathrm{C}_6\mathrm{H}_{12}\mathrm{O}_6+24\mathrm{\ NO^-}_3+24\mathrm{\ H}^+\rightarrow\ 30\mathrm{\ CO}_2+42\mathrm{\ H}_2\mathrm{O}+12\mathrm{\ N}_2$ | - 2,720 |
Sulfate reduction | $6\mathrm{\ C}_6\mathrm{H}_{12}\mathrm{O}_6\ +\ 6\mathrm{\ SO}_4^{2-}+9\mathrm{\ H}^+\rightarrow\ 12\mathrm{\ CO}_2+12\ \mathrm{H}_2\mathrm{O}+3\mathrm{H}_2\mathrm{S}+3\mathrm{HS}^-$ | - 492 |
Methanogenesis | $\mathrm{C}_6\mathrm{H}_{12}\mathrm{O}_6\ \rightarrow\ 3\mathrm{\ C}\mathrm{O}_2+\ 3\mathrm{\ CH}_4$ | -428 |
Ethanol Fermentation | $\mathrm{C}_6\mathrm{H}_{12}\mathrm{O}_6\ \rightarrow\ 2\mathrm{CO}_2\ +\ 2\mathrm{CH}_3\mathrm{CH}_2\mathrm{OH}$ | -244 |
*Table modified from Rittman and McCarthy, 2001. Used with Permission
Redox reactions written in the form shown in Figure 2, do not consider microorganisms or biomass as part of the reaction products. In order to do so, we need to consider reaction energetics and metabolic pathways. The catabolic pathway breaks down the electron donor (e.g., organic carbon) into smaller molecules and generates energy. The anabolic pathway uses organic carbon and energy harnessed in the catabolic pathway to synthesize large cell molecules. The two pathways complement each other in that the energy released from one is used up by the other. The degradative process of a catabolic pathway provides the energy required to conduct a biosynthesis of an anabolic pathway. Both pathways use electrons from the electron donor. Therefore, electrons from the electron donor partition into the two pathways. To write the full reaction equations, we will need to know the fractions of electrons being used for energy production (fe) and for cell synthesis (fs). The summation of $\ f_e\ and\ f_s\ is\ 1.0$. For different redox reactions, these fractions differ. Those higher in the biogeochemical redox latter generate more energy per electron flow so they need smaller fractions of electrons for energy production (smaller fe) and can channel more electron into cell synthesis. Therefore, values of fe increase and fs decrease as we go from aerobic oxidation that is high in the redox ladder to those lower in the redox ladder. That is, less microbial cells are being produced as we go down the redox ladder with the same amount of electron donor. To write the full reaction equation, we also need chemical formula for microbial cells. A commonly used one is C5H7O2N, approximating the ratios of major elements C, N, O, H, in biomass without including trace elements such as phosphorous, sulfur, metals, etc. Other examples include C8H13O3N2, often used to represent microbial sludge produced in waste water treatment plants. Details of how to write such reaction equations are given in Chapter 2 of Rittmann and McCarthy (2001). Here we show a few examples of reaction equations with biomass as the product using acetate as the electron donor (Cheng et al., 2016; Li et al., 2010). Note that for the text below, the overall reactions R is developed using half reactions of electron donor (Rd), electron acceptor (Ra) and cell synthesis (Rc). The energy production reaction Re = Ra - Rd and cell synthesis reaction Rs = Rc - Rd, respectively. The overall reaction $\ R\ =\ f_e\ R_e+f_s\ R_s=\ f_e\ \left(R_a-R_d\right)+f_s\ \left(R_c-R_d\right)=f_e\ R_a+f_s\ R_c-R_d$. A few examples are shown below.
$R_a\ \left(Reaction\ of\ electron\ acceptor\right):\frac{1}{4}O_2+H^++e^-=\frac{1}{2}H_2O$
$\mathrm{R_d\ \left(Reaction\ of\ electron\ donor\right)}:\ \frac{1}{4}\mathrm{HCO}_3^-+\frac{9}{8}\mathrm{H}^++e^-=\frac{1}{8}\mathrm{CH}_3\mathrm{COO}^-+\frac{1}{2}\mathrm{H}_2\mathrm{O}$
$R_c\ \left(Cell\ synthesis\ reaction\right):\ \frac{1}{4}HCO_3^-+\frac{1}{20}NH_4^++\frac{6}{5}H^++e^-=\frac{1}{20}C_5H_7O_2N_{AOB}+\frac{13}{20}H_2O$
Here AOB represent aerobic oxidating bacteria.
The catabolic pathway Re = -Ra - Rd that yields: $\frac{1}{8}CH_3COO^-+\frac{1}{4}O_2=\frac{1}{4}HCO_3^-+\frac{1}{8}H^+$
The anabolic pathway RS = Rc - Rd that yields:
$\frac{1}{8}CH_3COO^-+\frac{1}{20}NH_4^++\frac{3}{40}H^+=\frac{1}{20}C_5H_7O_2N_{AOB}+\frac{1}{20}H_2O$
Combining the two pathways using $R\ =f_e\cdot\left(R_a-R_d\right)+f_s\cdot\left(R_c-R_d\right)=f_e\cdot R_a+f_s\cdot R_c-R_d,$
with $\ f_e=0.4\ and\ f_s=0.6$, we have the following:
$\begin{array}{l}R:0.100\mathrm{O}_2(aq)+0.125\mathrm{CH}_3\mathrm{COO}^-+0.030\mathrm{NH}_4^+\rightarrow\\
0.030\mathrm{C}_5\mathrm{H}_7\mathrm{O}_2\mathrm{~N}_{\mathrm{AOB}}+0.100\mathrm{HCO}_3^-\ +0.090\mathrm{H}_2\mathrm{O}+0.005\mathrm{H}^+\end{array}$
$R a: \frac{1}{5} N O_{3}^{-}+\frac{6}{5} H^{+}+e^{-}=\frac{1}{10} N_{2}+\frac{3}{5} H_{2} O$
$\mathrm{Rd}: \frac{1}{4} \mathrm{HCO}_{3}^{-}+\frac{9}{8} \mathrm{H}^{+}+e^{-}=\frac{1}{8} \mathrm{CH}_{3} \mathrm{COO}^{-}+\frac{1}{2} \mathrm{H}_{2} \mathrm{O}$
$R c: \frac{1}{4} H C O_{3}^{-}+\frac{1}{20} N H_{4}^{+}+\frac{6}{5} H^{+}+e^{-}=\frac{1}{20} C_{5} H_{7} O_{2} N_{N R B}+\frac{13}{20} H_{2} O$
$\text{ Following }R=f_e\cdot R_a+f_s\cdot R_c-R_d,\text{ using a higher }f_e=0.45,\ f_s=0.55,$ we have:
$\begin{array}{l}R:\ 0.090\mathrm{NO}_3^-+0.125\mathrm{CH}_3\mathrm{COO}^-+0.0275\mathrm{NH}_4^++0.075\mathrm{H}^+\rightarrow\\
0.0275\mathrm{C}_5\mathrm{H}_7\mathrm{O}_2\mathrm{N}_{NRB}+0.1125\mathrm{HCO}_3^-+0.1275\mathrm{H}_2\mathrm{O}+0.015\mathrm{N}_2(aq)\end{array}$
$\mathrm{Ra}:\mathrm{\ FeOOH}(s)+3\mathrm{H}^++e^-=\mathrm{Fe}2^+2\mathrm{H}_2\mathrm{O}$
$\mathrm{Rd}:\ \frac{1}{4}\mathrm{HCO}_3^-+\frac{9}{8}\mathrm{H}^++e^-=\frac{1}{8}\mathrm{CH}_3\mathrm{COO}^-+\frac{1}{2}\mathrm{H}_2\mathrm{O}$
$Rc:\frac{\ 1}{4}HCO_3^-+\frac{1}{20}NH_4^++\frac{6}{5}H^++e^-=\frac{1}{20}C_5H_7O_2N_{\mathrm{Fe}RB}+\frac{13}{20}H_2O$
$\text { Following } R=f_{e} \cdot R_{a}+f_{s} \cdot R_{c}-R_{d}, \text { where } f_{e}=0.60, f_{s}=0.40$
$\begin{array}{l}\mathrm{R}:\mathrm{\ FeOOH}(\mathrm{s})+0.208\mathrm{CH}_3\mathrm{COO}^-+0.033\mathrm{NH}_4^++1.925\mathrm{H}^+\rightarrow\\
\mathrm{R}:\mathrm{\ FeOOH}(\mathrm{s})+0.208\mathrm{CH}_3\mathrm{COO}^-+0.033\mathrm{NH}_4^++1.925\mathrm{H}^+\end{array}$
$\mathrm{Ra}:\ \frac{1}{8}\mathrm{SO}_4^{2-}+\frac{9}{8}\mathrm{H}^++e^-=\frac{1}{8}\mathrm{HS}^-+\frac{1}{2}\mathrm{H}_2\mathrm{O}$
$\mathrm{Rd}:\ \frac{1}{4}\mathrm{HCO}_3^-+\frac{9}{8}\mathrm{H}^++e^-=\frac{1}{8}\mathrm{CH}_3\mathrm{COO}^-+\frac{1}{2}\mathrm{H}_2\mathrm{O}$
$Rc:\ \frac{1}{4}HCO_3^-+\frac{1}{20}NH_4^++\frac{6}{5}H^++e^-=\frac{1}{20}C_5H_7O_2N_{\mathrm{Fe}RB}+\frac{13}{20}H_2O$
$\text{ Following }R=f_e\cdot R_a+f_s\cdot R_c-R_d,\text{ where }f_e=0.92,\ \ f_s=0.08,$
$\begin{array}{l}R:\ 0.125\mathrm{SO}_4^{2-}+0.13525\mathrm{CH}_3\mathrm{COO}^-+0.004375\mathrm{NH}_4^++0.0065\mathrm{H}^+\rightarrow\\
R:\ 0.125\mathrm{SO}_4^{2-}+0.13525\mathrm{CH}_3\mathrm{COO}^-+0.004375\mathrm{NH}_4^++0.0065\mathrm{H}^+\end{array}$
In each of these reactions, the ratio of the biomass carbon number (C) in the biomass formula (C5H7O2N) produced per C of organic C is called yield coefficient. For example, for aerobic oxidation, the yield coefficient is 0.030*5/(0.125*2) = 0.60; for denitrification, the yield coefficient is 0.0275*5/(0.125*2) = 0.55; for the FeOOH reduction reaction, the yield coefficient is 0.033*5/(2*0.208) = 0.40; for the sulfate reduction reaction, the yield coefficient is 0.004375*5/(2*0.13525) = 0.08 C biomass /C organic. In general, large fs values leads to high yield coefficients. Values of fs are typically in the range of 0.6 – 0.75 for aerobic oxidation, 0.55 – 0.7 for denitrification, and 0.08 – 0.30 for sulfate reduction. For typical values of different types of redox reactions, please refer to Rittmann and McCarthy (2001).
The kinetics of microbe-mediated reactions are often described by Monod rate laws in the following form:
Here $\mu$ is the rate constant (mol/s/microbe cell), $C_{C_{5} H_{7} O_{2} N}$ is the concentration of microbe cells (cells/m3), Cd and Ca are the concentrations of electron donor and acceptor of the reaction (mol/m3). The Km,D and Km,A are the half-saturation coefficients of the electron donor and acceptors (mol/m3), respectively. These coefficients are the concentrations at which half of the maximum rates are reached for the electron donor and acceptor, respectively. If the electron donor is not limiting, it means that $C_{D} ? K_{m, D}$, so that the term $\frac{C_{D}}{K_{m, D}+C_{D}}$ is essentially 1.
In order to represent the biogeochemical redox ladder, we will need to introduce an additional term in the dual Monod rate law, the inhibition term. A Monod rate law with an inhibition term looks as follows:
Here the KI,H is the inhibition coefficient for the inhibiting chemical H. In contrast to the Monod terms, the inhibition terms become 1 (not inhibiting) only when KI,H?CH.
As an example, in a system where oxygen and nitrate coexist, which is very common in agricultural soils, aerobic oxidation will occur first before denitrification occurs. The sequence of that can be represented by the following:
These two rate laws, with proper parameterization, will ensure that denitrification reaction occurs only when O2 is depleted to a certain extent that the term $\frac{K_{I, O_{2}}}{K_{I, O_{2}}+C_{O_{2}}}$ approaches 1.0. If other electron acceptors that are lower than nitrate in the redox ladder also occur, then multiple inhibition terms are needed. For example, if there is also iron oxide in the system, we will need the following for the iron reduction rate law:
Where the additional litrate inhibition term will allow iron reduction to occur at sufficiently significant rates only when nitrate level is low compared to the inhibition constant value.
The system: I have a 200 ml bottle in my kitchen. I accidentally drop some recently-fertilized soil grains from my backyard into the bottle. I close the bottle after the accident. The soil has some microbial cells, along with some organic carbon (assuming the formula of acetate), and nitrate. The initial acetate and NO3 concentration in the bottle are 5.0 and 1.0 mmol/L, respectively. Assume initially there are some dissolved O2 in the water at the concentration of 0.28 mmol/L while all other chemicals in the soil are not reactive. The initial biomass concentration of O2-reducing and nitrate-reducing biomass are 2.0×10-6 and 1.5×10-6 mol-biomass/L, respectively. In the bottle, I put in some magic powders (to be invented in the future by my daughter Melinda Wu!) that would automatically show the concentrations of chemical species without the need of measurements (that would be my dream come true!). The microbe-mediated reactions advance as the follows:
Aerobic oxidation $\left(f_{s}=0.6 \text { and } f_{e}=0.4\right)$:
$\begin{array}{l} 0.500 \mathrm{O}_{2(a q)}+0.417 \mathrm{CH}_{3} \mathrm{COO}^{-}+0.067 \mathrm{NH}_{4}^{+} \rightarrow \\ 0.067 \mathrm{C}_{5} \mathrm{H}_{7} \mathrm{O}_{2} \mathrm{~N}_{\mathrm{AOB}}+0.500 \mathrm{HCO}_{3}^{-}+0.200 \mathrm{H}_{2} \mathrm{O}+0.150 \mathrm{H}^{+} \end{array}$
Denitrification $\left(f_{s}=0.55 \text { and } f_{e}=0.45\right)$:
$\begin{array}{l} 0.090 \mathrm{NO}_{3}^{-}+0.125 \mathrm{CH}_{3} \mathrm{COO}^{-}+0.0275 \mathrm{NH}_{4}^{+}+0.075 \mathrm{H}^{+} \rightarrow \\ 0.0275 \mathrm{C}_{5} \mathrm{H}_{7} \mathrm{O}_{2} \mathrm{~N}_{\mathrm{NRB}}+0.1125 \mathrm{HCO}_{3}^{-}+0.1275 \mathrm{H}_{2} \mathrm{O}+0.045 \mathrm{~N}_{2(a q)} \end{array}$
Reactions | $\mu_{\max }(\text { mol/mol-biomass/yr) }$ | $\begin{array}{c} K_{m, \text { acceptor }} \\ \text { (mol/kgw) } \end{array}$ | $\begin{array}{c} K_{m, \text { donor }} \\ \text { (mol/kgw) } \end{array}$ | $\begin{array}{r} K_{I, O 2(a q)} \\ (\mathrm{mol} / \mathrm{kg} \mathrm{w}) \end{array}$ |
---|---|---|---|---|
Aerobic | 50000 | 1.00x10-4 | 1.00x10-3 | |
Denitrification | 20000 | 1.00x10-3 | 1.00x10-3 | 1.00x10-6 |
*Range of relevant parameters is from Cheng et al., 2016; Li et al., 2010.
Questions:
Please watch the following video: Microbe Mediated Reactions (24:50)
Li Li: So let's get started. This is lesson 5, so we'll be talking about microbe-mediated reactions. I put in the reading materials the importance of microbe-mediated reactions-- how ubiquitous they are in natural environment. And we also talked about the biogeochemical redox ladder in the reading material. So things happen in sequence.
And then we have this example of-- we're talking about soils, so I grabbed some soil in the back yard, putting it in a water bottle in the kitchen. So this is the cup, I would say. So that's the system we have. In that example, we put some soil there. Let's close it. We have some water there. So we have some soil. Let's just put these as grains.
And the water originally has some oxygen in it. And I told you that the soil that I have we just recently fertilized. That means you will have nitrate in the system. So essentially, in the system, you have oxygen and nitrate as electron acceptors. And typically, for example, the soil usually have bacteria already there. Microbes is everywhere. They do a lot of work.
And also, usually, you would expect the soil carbon have-- the soil have some organic material. So essentially, it also have organic carbon, which we use acetate to represent that. So essentially, we're talking about a system that has water in it. It's well mixed. It has bacteria in it. So let's use that to represent bacteria. And it has two electron receptors, oxygen and nitrate. Perfect.
We have a system. We know it's going to have something going on, right? And again, this is well mixed. We are not talking about the advective transport, dispersion-diffusion yet. So this is a well-mixed system, meaning there's no concentration gradient in the system.
Now, then, the other system you have-- so you have microbes We use microbes The formula for bacteria is C5 H7 ON. Something like that. And it depends on what type of bacteria it is. If it's using oxygen, then we put oxygen on the side, just to indicate this is oxygen-reducing bacteria. If it's nitrate-reducing bacteria, we put nitrate there, just to distinguish between the two different bacteria. So essentially, you have two types of bacteria there in the system.
That's the system you have. And what's there? All the players, right? Now, we think about the reactions, what are the reactions? So first of all, we talk about, in the biogeochemical redox ladder, oxygen is the one that is going to be used first, because usually it has fast reaction rate. The microbes also get more energy out of it. So that means it has growth advantage, because it can grow more with the same amount of organic carbon they have. So the first one we call aerobic oxidation, meaning you have the transformation from CH3-- we use acetate to represent organic carbon. I'm not putting all these numbers into stoichiometrics. But essentially, just writing what are the major products.
So organic carbon will become oxidized to become bicarbonate or other forms of carbon. If it's relatively acidic condition, you might have CO2 gas coming out. And you have water or hydrogen under whatever condition it allows. And then you would grow bacteria.
So bacteria is also one of the products that we talk about. It'll has N. And then producing the oxygen. It's oxygen-reducing or aerobic bacteria or microbes. So that's the oxygen related reaction.
So we talk about, OK, that reaction will happen first. And remember, this is a closed system, meaning we have a certain amount of oxygen dissolved at first, and then it's going to be depleted over time because of reaction. So the microbes will be using this up. And then, the next steps, the next reaction that's going to happen will be denitrification.
So You will have, again-- that's assuming we will have plenty of organic carbon there. So then you will have nitrate. Again, there will be inorganic carbon become produced. It's the same process, If you think about it, as the human system. We breathe oxygen. We eat food, like grains, bread, rice, different type of, essentially, organic carbon. And then we breathe out CO2, which essentially is another form of bicarbonate, right? And then the water. Let's not write that.
And then you have another type of bacteria coming-- microbe coming out, which will be the denitrification microbe. So these are the two drivers of the system. But then, typically, you would also have-- I'm just writing at the side-- there should be a somewhat fast reaction compared to these microbe-mediated reaction. Would be the transition between CO2 aq, bicarbonate, and carbon. These are the inorganic carbon species that we know is going to happen. We talk about that in aqueous speciation or aqueous complexation reactions that we discussed before. So these are very fast reactions. So the major drivers are really these. And if we have other species, for example, calcium or magnesium, if there's plenty of them, there might be precipitation come, because these reactions will produce a lot of bicarbonate. And bicarbonate can-- when the condition allows, it could precipitate calcite, for example.
We're not going to focus on these. Our major goal right here, right now, is talking about these drivers. And if there's applications that you will need to consider these as something that would be specific about particular applications.
So we have these reactions. These are the processes we focus on. And then we'll be thinking about the species. What other species are there? So, for sure, we would have all the electron acceptors there, which is oxygen-- it dissolved oxygen. You have nitrate. You have its CH3 COO. These are the variables that we are solving for in these equations.
And then we also have-- let's put bicarbonate. But we know this is going to be also carbonate, CO2 aq, or H2 CO3. These three are almost like equivalent with each other. Hydrogen ion, for sure, and then it goes with OH-. Electron acceptor, electron donor, products. Nitrogen is-- did I write? OK. Here, I should add another, which is the product of denitrification. One product of denitrification will be N2. There could be a lot of-- a few other intermediate reaction products, like NO2-. People have seen that. Or ammonia, even. This could happen too. We're just using this as an example. So N2. So these are the major species.
But also, when we solve for these, when we are writing these four reactions in terms of these reactions. Microbe as products. We also should have C5 H7 ON0 O2. And then C5. These are the two major products too. It's not just the abiotic species. We would consider the microbes are the products of this reaction too. So these are the species that you are solving.
So we are not talking about other complexation reactions. So really, we're focusing on these species. And you probably quickly know, by this time, these are the-- these could be the second species, second species, second species. These are the things that we consider. And then the primary species would be oxygen, nitrate, this, this, hydrogen, N2, N2, and this. These are the primary species. Let me just write as primary to indicate the color. So that's what you have in terms of variables that we are solving for. Now, again, we are talking about the well-mixed system. So when we think about the reactions, the equations, it would be really the reaction in terms-- for example, for oxygen, you have something like this. It would be dc dt for oxygen, for example.
Oxygen is being consumed. So the volume-- so it's really volume of the bottle, the mass, and then times this dc dt. But then you also have, essentially, the mu. The mu-- let's call this mu as the stoichiometric coefficient. Fish It's consuming. It's in the left-hand side of the reaction, meaning it's being consumed. So this is stoichiometric. If, for example, this is 0.5, this value will be 0.5. And then we will be-- we think about the consumption. The rate of this will follow Monod, dual Monod kinetics. So you would have this mu max for oxygen term. And then we have this-- how much bacteria there. So it'll be the C5 H7 O2 N.
I should have an O2 here. That's why I feel something is missing. It's weird. O2. It should all have O2 there. All right. So you have these species, O2, and this is oxygen-reducing bacteria. So $\mu$ max, concentration of biomass, and then you should have two Monod terms. One is for concentration of oxygen, and then the Km oxygen plus concentration of O2. That's an electron acceptor. And then you have the C-- . Sorry, I am writing a bit of-- so this would time, right? So Km, and then the plus C acetate. CH3. I can't write anymore, but you know what I mean. So this is electron. Donor term electron. Acceptor term biomass. Maximum reaction rates. So that's for oxygen.
Similarly, you would write this for V dc. It's also being consumed. So essentially you would have the same, but you will have mu of CH3 COO- and the same expression. Right? It's the same reaction. We just consume different chemical species in different proportions. And you would write again, also, for bicarbonate. That's one of the products. So you really should have positive sign before this. So you have to write V dc dt. This would be like total carbon. And then $\mu$ of bicarbonate. So all this goes on. So you will be writing for all these.
Now, when you are write for-- so if you only have oxygen, you would have oxygen, this acetate, bicarbonate, and other species. That's good enough. But if you have nitrogen, then your equation for nitrogen would be a little bit different, in the sense that you would have the would same. But they would be replaced by nitrate acetate. But you will have the inhibition term by oxygen, because oxygen would act as inhibitor for denitrification to occur. When you have plenty of oxygen in the system, denitrification is not going to happen. Only after this concentration of oxygen has been decreased to levels that this term is more or less much higher than this, then this value would be close to 1. Then denitrification really kicks in. So that's for nitrate. Also notice, when you have both electron acceptors, this reaction will be happening at the same time just, at first, oxygen what will be the dominant one, because this term will have very small value, will inhibit denitrification to happen. But then this, the rate of this will become smaller and smaller because oxygen becomes smaller concentrations, this term becomes bigger and bigger. So there's a switch between different electron acceptors. The other thing is when you have multiple electron acceptors, when you write these acetate equations or when you write these for the produced inorganic carbon, you need multiple terms too, because the denitrification process also contributes to the consumption of acetate and production of bicarbonate. These are common products or common reactants. So you have multiple terms if you have the nitrate-related term.
I'm not going to detail everything. But essentially, that's what you think would have. So solving these equations. Let's say you have all the species. You'll be writing how many independent equations for primary species. And you have three fast reactions that you can form algebraic relationship to solve these other secondary species. You solve everything as a function of time. So this is dt. This is not partial t. Now, in that case, what do you think will happen when we have these products? So I talk about we have some magic powder putting in the bottle will help us. We read our numbers. What do you think will be the trend of, for example, oxygen? And what's going to be the trend of temporal evolution of nitrate. Using nitrate for this one. How do you think they are going to look like in terms of curve? So that will be actually the output of the model. You can see about how it was going to happen. So aerobic oxidation going to happen. First, oxygen concentration is going to be consumed first, then decreased. So let's say it's maybe oxygen starts from here. You see a decreasing trend of oxygen. Let's say it's something like that.
And then, what about nitrate? So nitrate is-- at first, because it was in the presence of oxygen, you wouldn't see much of nitrate decreasing. So it probably, at the beginning, it's going to, let's say, start from somewhere here. It'll be relatively flat. But it also, when the nitrogen-- I'm sorry. When the oxygen starts to decrease, it just start to show a sign of-- it depends on rates. How fast, how long it'll take. It depends. So when it's like almost completely depleted, then you would have nitrate kind of going full swing. And then it's going to decrease really fast. But then, when you think about the products, this produced N2. So N2, let's say, at the beginning, it's not much. It's almost zero. Now, when nitrate starts to be denitrified, then N2 starts to increase. Increase. It's almost like in the mid-- it almost mirrors nitrate. So this would be N2. That's the product of denitrification. Oxygen decrease, followed by nitrate decrease, and oxygen increasing.
What if you also have sulfate? Let's use another color. Let's use blue again. What if you have nitrate? So let's assume there's no ion there. So when you have sulfate-- so there's another electron acceptor, which is sulfate, which is also pretty active in subsurface environment. Probably not in very shallow soil, but in bit of deeper, when you don't have a lot of oxygen, you tend to have sulfate.
So sulfate would be probably, let's say, something like-- let's start here. It would start to become decreased when this denitrification happen. Then the concentration decreases. Then it will also follow maybe one-- depends on how much electron donors you have and acceptors. This is representing sulfate. And the product of that is sulfides. So sulfides probably would be small at the beginning. But then it'd be increasing over time. Something like that. Call it HS-. This is like randomly ending. Another electron, except it's not easy equation. But if you need to write for sulfate too, then you would have three major microbe-mediated reactions. And in addition to that, you would also have an equation related to sulfate and sulfide production. And notice, also, sulfate reduction will be inhibited by both oxygen and nitrates. So you have two inhibition terms. So in that case, only when oxygen and nitrate is depleted to pretty small amount, then the sulfate reduction can start occurring. So these inhibition terms really would ensure the biogeochemical redox ladder in the system.
OK. So that's what we have for this unit, for this lesson. And you can use this. Will set up stage for you to do the homework. Think about different reactions that are going to happen in sequence in natural systems-- in soil, aquifers, and then how microbes evolve and the different concentrations involved with time.
Actually, we didn't talk much about microbes. But the microbial concentration will increase over time, following where the-- different types will be following their own electron acceptor. Electron acceptor evolution. Oxygen decrease and aerobic microbial will be increased. Things like that. OK. So let me stop here, and you are ready to do the homework. Thank you.
Please watch the following video (29:51):
Li Li: OK. Let's look at this example 5.1, which we are going to set up control rounds for microbe-mediated reaction. So system is still a closed system with a bottle, which means we don't have flow comes in, flow comes out. We're really focused on the reaction itself.
And you imagine there's a bottle, 200 milliliter in volume. We have soil which contains some organic carbon. We also have nitrate and oxygen in the system.
And then we also have these magic powder, et cetera, which we're very excited about. Everything is ready to do this. Now, these two reactions that I put here, aerobic oxidation, denitrification reaction, these two reaction is already written as an outcome of a half reactions that we discussed earlier with these specified energy petitioning proficient in term how much goes into catabolic pathways, how much goes into the anabolic pathways. OK. So here we are looking at these two interaction, imagining that these are happening. Aerobic oxidation should occur first before denitrification. I also put the rate of kinetic parameters here.
So we have $\mu$ max, Km, which are half saturation constants, and inhibition coefficient, which is Ki for denitrification reaction because it's going to be inhibited by the presence of oxygen. OK. And then we are charged to answer the question, which reaction would occur first, which I already told you. But we are going to look at what is the magic powder, of course. We are going to look at just the evolution of the concentrations for different species. How does acetate concentration, which is organic carbon, evolve over time? In addition, we are looking at oxygen nitrate, totally organic carbon, biomass, et cetera.
So let's look at our folder. Now here when you look at this there are several reactions. I'm sorry. There are several files. The same one is the Input file. There's one Database file, which is Old Rifle Database. That one is a regular database that you have seen before similar to like the database datcom .dbs, which has a lot of fraction equivalent constants, reaction stoichiometry, and primary species, second species, and all that. That's the same.
One major thing that is different is now here you have two more database file, which is AqueousMicrobe and SolidMicrobe. AqueousMicrobe list all the microbe-mediated reaction stoichiometry that occur in aqueous phase. For example, for this example 5.1, you have aerobic oxidation. You have nitrate. You have denitrification. All these reaction happen in the aqueous phase. So this file is for aqueous phase microbe-mediated reactions. There is another file which is SolidMicrobe, which has catabolic pathways for reactions that occur that involve solid phase. For example, if you have reduction of iron oxides, these reaction stoichiometry would appear in this file instead of the previous aqueous file. So in this particular example, we actually don't use this SolidMicrobe.
This is putting here just for those of you who will be running microbe-mediated reaction involving solid phase. It's going to be handy for you later when you need these reaction's stoichiometry. So what's in this aqueous microbe?
The reaction's stoichiometry is actually more complicated from you just see in the example file. So it actual written as these two half reactions. One is catabolic.
So it's a catabolic. What catabolic mean is breaking down things. So essentially, it's the oxidation of acetate here becoming inorganic carbon, which is bicarbonate. So essentially, you have these reactions that happen, right? You think about this keq value. This should be capital K. This keq value is more like a dummy number. It's big, which means it's not going to functioning as limiting the reactions.
The reaction almost always new dox reactions. Reaction almost occur always in the direction from left to right, meaning acetate become bicarbonate as long as you have the reactants that are still there. So that's catabolic. But you also notice there's another anabolic reaction which includes, for example, the build up of bacteria. So you have C5H7O2 Anabolic mean building up things.
So essentially, this reaction is indicating using the ammonium and acetate to build up microbial material C5H7P2. The bacteria is using oxygen. So this is how you figure, OK, this actually is for the oxygen reducing bacteria or micro-organism. So you have these two. And in the online material, we talk about how much the energy produced from the reduction of oxygen goes to the catabolic pathway and anabolic pathway. And this is fs and fe values, essentially, for these purposes.
So if you look back in these, so here we're saying fs and fe is the number for aerobic. And it's different for denitrification, which is determined by reaction thermodynamics. OK. So you will be looking at this. And so what's in the Input file? Let me open it. And the Input file, so I mentioned the key difference is you have these two additional database. Essentially, in this case, it would be really one additional database. So everything else would be very similar. You have this runtime. You're reading database. You have output keyword spatial_profile. Putting these different times, times the areas.
And you also have these, for example, discretization which you only have one grid. You should have the total volume of 200 milliliter, which is 200 centimeter cubed, which is same as 200 milliliter. You don't have boundary condition, because it's closed system. You have initial conditions. So let's look at aqueous kinetics. So if we think about aqueous kinetics, right, you have this anaerobic oxidation. And you're essentially picking up the name.
This block is different in the sense, OK, you will need to specify aerobic oxidation reaction-- the pathway you have for aerobic oxidation half. So this name is there to be the same as this name, which is the catabolic pathway. So once it recognizes the name, it recognizes this reaction stoichiometry. And then also in the Input file you have another the pathway which is with this name. Which the code will be picking up this pathway, looking and picking up this name and process. OK. So here, you have aqueous kinetics block.
And this aerobic oxidation H5 is defined as false aerobic oxidation reaction. we actually don't need this hw5. This is really not necessary. Actually, all the code is looking for is this name and the other name here. So this is for the catabolic pathways. This is for the anabolic pathways that build up, so C5H7O2 which is oxygen reducing microbe. So if you go back into, for example, the aqueous microbe, the corresponding name is here. So it would take up for catabolic you would have this reaction stoichiometrically or reducing by the code. And when it sees this, it will have this set of reaction's stoichiometry reading. Now, we also talk about those $f$s and $f$e values, right? So in the Input file, you are actually putting the 0.4 the fs and fe value here and here. So this is essentially the $f$e values, 0.4 and 0.6.
So this is are how you specify these $f$s, $f$e values. So essentially, the code is already programmed to reading these half reactions and then use these $f$s and $f$d values to come up with overall reaction as we indicated here. So you don't need to worry about changes, like this detailed calculation. As long as you have the $f$s and $f$e values, you're able to construct that. And then this rate is a $\mu$ max. So actually in this increase microbe, you also have-- so the top part is the reaction stoichiometry. And there's aqueous kinetics in the bottom part.
And so the bottom part, essentially, is telling you the reaction to the form-- OK, rate expressions. So form of the rate law-- or actually, it says this term needs-- if you have H5 here-- hw5 is what I'm talking about. This name needs to be the same name as what you have here. So if the database says have this dash or underscore h homework five, then you do need it here then to be consistent as a kinetic block. Anyway, so look at this. We talk about all the reaction for microbial kinetics. So this rate, essentially $\mu$ max, is the same as the Input file.
And then you have monod_term. Which mean it depends on total acetate aqueous oxygen concentration. These value are the key m value for acetate. This is a key m value for oxygen, which is the same as what you have here. So all these reaction kinetics terms are essentially in the later part of the bottom part in the database aqueous microbe.
Same for denitrification reaction, which I'm not going to detail. So essentially what you have here is, OK, aerobic oxidation H5. The coder will be looking for both reaction's stoichiometry of the two half reaction, build up the reaction, and then look for other kinetic information-- same for nitrate reduction inhibition homework five. So you use corresponding rocks for this reaction and also for catabolic and anabolic. Notice that here your $f$s and $f$e numbers of different now, right? This would be the $f$e actually. This is $f$s. So they're different, because of the reactions I've done for microbe-mediated reactions.
So these are the aqueous kinetics. And pay attention that with all these reactions when you look at the reaction in the database, it's important as you put these-- for example, primary species acetate, ammonium, nitrate, oxygen-- in the primary species. Bicarbonate, for example. The other things that are really important is that if you think about how we will present the biomass-- right, so biomass a lot time once they are coming out of the consumption of organic carbon, they can either stay in the water phase or a lot of times they attach to a solid phase and form biofilm. And it's not as mobile as in the water phase.
And actually a lot of literature has been showing that this is a big tendency, forming biomass or biofilms. They are not very mobile. So essentially, what we do in the code-- if you look at these reaction stoichiometry, they are really producing aqueous microbe, right? These are aqueous faces. But then we are saying once they have produced, they will quickly become solid phase attach on the surface. So we're actually using a kind of precipitation reaction to approximate that process. So there is a reaction related to this. And there's a solid phase call a mineral. Mean Yeah. So there's a solid phase species. So you have this C5H7O2NO2, C5H7O2NNO3, right? So these are O2NO3.
These are the species that actually goes into representing the biofilm. So we put a relatively small kind of-- it still follow this, still follows TST rate law as an approximation. But you can put smaller keq or larger keq values to control how much the biomass can be in aqueous phase. Is The majority of them usually come into the solid phase or biofilm phase quickly. And use these reaction rates to represent how fast they become biofilm. If you put these very small values, like 10 to minus 8, 10 to the minus 10, then a lot of biomass will be flushing out. Because the resin time is fast. The resin time is short compared to the precipitation reaction to happen. So you want to put these rates relatively high if you want to see a lot of biomass actually go through solid phase.
And quartz is representing the majority of the soil here to say, OK, essentially soil, like soil is a majority. A lot of them are quartz. That's why it's here. But it's actually not necessary to be here. Because it's not involving any other reactions. Let me see. So primary species, secondary species-- these are all necessary reaction. And then you have these conditions, right? So initial condition you want to be consistent with the example.
There are the different chemical species. The condition doesn't give you this. We usually put this for small numbers. Bacteria would build up that quickly. And let me see. So these conditions and then others-- so you do need to put a soil microbe to represent the biomass. And we do put in some volume fraction here to represent initial amount of these biomass here. A lot of time, they are like volume fraction essentially. If you put very small numbers, the system will need to take a long time to build up enough biomass. So you'll see concentrate decrease at some point. So if you put small numbers, there will be a lag phase in terms of how much biomass need to grow first before they actually functioning to allow you to see changes in concentrations.
So this is for the Input file. For a solid microbe, actually there's a set up in the system that you only have reactions stoichiometric there. So reaction kinetics, that actually is doing the Old Rifle Database. Because these reactions, original was already there. So aqueous microbe was added later on. So it's a bit different there. But anyway, in one of these or two of these file, we find the reaction's stoichiometry and reaction kinetics. So let's run this.
So everything is put in there. Let's run this. And we can kind of quickly walk through the Output file. So if we look at Output file, each of these have four time, 1, 2, 3, 4. And that's because in this file, we have put in a spatial profile. We have 1, 10, 20, 30-- four time points. Really for the batch reactor it sees spatial profile, not necessarily you can comment them out. If you would like to do that, then we can actually-- let me just do this, so you Output file is a bit cleaner.
You really don't need that other file. So I'm deleting all this. And I'm going to change-- I just comment out, so then it wouldn't have spatial profile coming out. And let's rerun it. Oh. OK. Sorry. So I do need that. Because I do need this spatial profile to give the time points, when it's going to end. So let's put 30, which means it's going to be end in a month. So you do need this. This is necessary for a way to look at a time evolution. You can just put one, so you would only have one spatial profile coming out. Let's rerun it. OK. Now, you have all this. And you notice every one of these species only have one, right? Because I only give one time point for the code to spit out a spatial profile.
But what we are really looking at should be time series. Which, essentially, it's kind of-- later one, when we introduce column expander. It's kind of like bricks through curves. And the grid block is the grid block number. We have one grid block number. This will give the time evolution of these chemical species in this particular a grid block. So let's look at this first. So this will be the file that it will rely on to plot out the concentration of different species.
See, OK. It's here. OK. So you can look at time, pH, oxygen, NO3, N2, acetate, bicarbonate. And this is a bit harder to see. But the easiest way would be you expose them either in Excel or read them in MATLAB, and then directly plot. And actually so let me just close this. So output in for this example is really this. So we have already plotted out You can see this. So Reaction I, of course, would occur first. And you can look at the concentration as a function of time. And organic carbon, right? So what this tell you is this is based that carbon chemical output, right?
So we talk about oxygen will be consumed first. So O2 is coming down and not coming back again, because it's being consumed. And then at first, nitrate is still remain like its flat, meaning it's not much reduced. But then it wasn't reduced right after oxygen is consumed. The reason why-- if you look at the Input file, you realize that nitrate have a relatively small concentration at the beginning. So that's what I'm talking about. First of all, it has lower rate than aerobic. For example, mu max is much smaller. And also initially, doesn't have as much biomass there. So essentially, the denitrification biomass need to build up over some time. So that's why you see very flat even after oxygen is consumed, remain flat for a certain period of time, about 10 days, and then start to get down. Because at that time, there's enough bacteria to do that.
And it's actually you can see that in this bacteria biomass figure, right? Biomass is a function of time. And you have the nitrate reducing bacteria build up. Early on, it's relatively small. It takes some time to build up, right? Take a bit of time to build up. And then the oxygen reducing bacteria increase quickly, because it has faster rates. And once they get to the point about beyond five, it remain constant. Because it does not grow anymore. Oxygen is all gone.
Now, correspondingly, you see for acetate, so both reaction consume acetate. And they are producing DIC, which is bicarbonate. So acetate is being reduced.
Why you see this kind of changing almost like in slope here? Because x oxygen reduction is much faster just as this is same slope right here. And then this one is almost mirror the nitrate concentration. So essentially, this is corresponding to the consumption of acetate by nitrate. And these two essentially mirror each other with this almost symmetric reactions stoichiometry, right? So DIC increase quickly first, but then slow down until denitrification with certain rates that are higher. So these are things that you can actually look at and think about. So is this one, you do see the redox later. Oxygen occur first, nitrate occur second. But it doesn't occur right away after oxygen, because it has a bit of lag time.Well, this is example 5.1. And homework have some extension questions that are related. But there's also another reaction that I asked you to look at, which was sulfate. You will be looking at the different files and the sulfate related reactions.
In addition to the chemicals and biomass specified in example 5.1, assume there is also sulfate in the water at the concentration of 3.0 mmol/L. The initial sulfate reducing bacteria has the concentration of 1.0×10-6 mol-biomass/L.
Sulfate reduction (Fs = 0.08 and Fe = 0.92) goes as follows:
$\begin{array}{l}
0.125 \mathrm{SO}_{4}^{2-}+0.13525 \mathrm{CH}_{3} \mathrm{COO}^{-}+0.004375 \mathrm{NH}_{4}^{+}+0.0065 \mathrm{H}^{+} \rightarrow \\
0.004375 \mathrm{C}_{5} \mathrm{H}_{7} \mathrm{O}_{2} \mathrm{~N}_{S R B}+0.250 \mathrm{HCO}_{3}^{-}+0.013 \mathrm{H}_{2} \mathrm{O}+0.125
\end{array}$
Sulfate reduction rate parameters are specified in Table 2. Note that it has two inhibition terms because it has both O2 and NO3 above in the redox ladder.
Reactions | $\mu_{\max }(\mathrm{mol} / \mathrm{mol}-\mathrm{biomass/yr})$ | $\begin{array}{l} K_{m, \text { acceptor }} \\ \text { (mol/kgw) } \end{array}$ | $\begin{array}{c} K_{m, \text { donor }} \\ \text { (molkgw) } \end{array}$ | $\begin{array}{c} K_{I, O 2(a q)} \\ \text { (molkgw) } \end{array}$ | $\begin{array}{c} K_{I, \mathrm{NO} 3(a q)} \\ (\mathrm{molkgw}) \end{array}$ |
---|---|---|---|---|---|
Sulfate reduction | 35000 | 1.25x10-3 | 1.25x10-3 | 1.00x10-6 | 1.00x10-3 |
*Range of relevant parameters is from Cheng et al., 2016; Li et al., 2010.
Questions:Extension: In Problem 1, we mainly discuss how O2(aq) affects the microbe-mediated reactions. Other thermodynamics and kinetic parameters, including maximum biomass grow rate, half saturation of the electron donor and acceptor, and the concentrations of electron donor and acceptor, also affect the biomass reactions. Please use the input and database files from Problem 1 as the base case to do the following analysis comparing N2(aq) and sulfide evolution:
In this lesson, we have learned to write microbe-mediated reactions, the Monod rate law, the set up of reactions in CrunchFlow, and the importance of different parameters in controlling the rates and onset of different reactions using the biogeochemical redox ladder.
Borch, T., Kretzschmar, R., Kappler, A., Cappellen, P.V., Ginder-Vogel, M., Voegelin, A. and Campbell, K. (2010) Biogeochemical Redox Processes and their Impact on Contaminant Dynamics. Environmental Science & Technology 44, 15-23.
Cheng, Y., Hubbard, C.G., Li, L., Bouskill, N., Molins, S., Zheng, L., Sonnenthal, E., Conrad, M.E., Engelbrektson, A., Coates, J.D. and Ajo-Franklin, J.B. (2016) Reactive Transport Model of Sulfur Cycling as Impacted by Perchlorate and Nitrate Treatments. Environmental Science & Technology 50, 7010–7018.
Li, L., Steefel, C.I., Kowalsky, M.B., Englert, A. and Hubbard, S.S. (2010) Effects of physical and geochemical heterogeneities on mineral transformation and biomass accumulation during biostimulation experiments at Rifle, Colorado. J. Contam. Hydrol. 112, 45-63.
Rittmann, B.E. and McCarty, P.L. (2001) Environmental Biotechnology: Principles and Applications. McGraw-Hill, New York.
This lesson introduces physical processes, including advection, diffusion, and dispersion processes, in 1D system. This is the first time in this course that we introduce the space dimension. An example will be shown about how to set up a one-dimentional (1D) flow and transport simulation in a homogeneous column in CrunchFlow.
By the end of this lesson, you should be able to:
Required Readings |
|
---|---|
Optional Reading (if phreeqc) |
|
To Do |
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
Understanding flow and transport processes in the natural subsurface are important for a wide range of applications and disciplines. Flow and transport processes play a critical role in ground water and surface water management and environmental protection, energy extraction from deep hydrocarbon reservoirs, chemical weathering, and soil management. Here we primarily discuss fundamental flow and transport processes in natural subsurface systems.
For a conservative tracer that does not participate in reactions, advection, dispersion and diffusion are the major transport processes that control its transport. Advection, also called convection in some disciplines, determines how fast a tracer moves with the fluid flow; dispersion and diffusion processes are driven by concentration gradients and/or the spatial variation of flow velocities and therefore determine the extent of spreading.
Advection is the transport process where solutes flow with the bulk fluid phase. This is like you let go of yourself when you swim so you have the same velocity of the flowing fluid. The advective flux, $J_{a d v}\left(\mathrm{~mol} / \mathrm{m}^{2} / \mathrm{s}\right)$ of the solute, can be expressed as
where $\phi$ is the porosity of porous media; v is the linear fluid velocity in poroud media (m/s); and C is the solute concentration (mol/m3). Flow through a porous medium is described with Darcy’s Law:
where $u$ is the Darcy flux ($\left(m_{\text {fluid }}^{3} / m_{\text {medium }}^{2} / s\right)$) that is proportional to the gradient in the hydraulic head $\nabla h(\mathrm{~m})$; K is the hydraulic conductivity (m/s); One can also write Darcy’s Law in terms of hydraulic head by defining the hydraulic head as
where z is the depth (m), P is the fluid pressure (Pa), is the fluid density (kg/m3), and g is the acceleration of gravity (9.8 N/m2). The hydraulic conductivity (m/s) can be expressed as
where $\kappa$ is the permeability of the porous media (m2) and is independent of fluid property, $\mu$ is the fluid hydraulic viscosity (Pa·s). Representative values of hydraulic conductivity and permeability are listed in Table 1 for various subsurface materials.
K (m/s) | 100 | 10-1 | 10-2 | 10-3 | 10-4 | 10-5 | 10-6 | 10-7 | 10-8 | 10-9 | 10-10 | 10-11 | 10-12 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
10-7 | 10-8 | 10-9 | 10-10 | 10-11 | 10-12 | 10-13 | 10-14 | 10-15 | 10-16 | 10-17 | 10-18 | 10-19 | |
Unconsolidated Sand & Gravel | Clean Gravel | Clean Sand or Sand & Gravel | Very Fine Sand, Silt, Loess, Loam | ||||||||||
Unconsolidated Clay & Organic | Peat | Stratified Clay | Unweathered Clay | ||||||||||
Consolidated Rocks | Highly Fractured Rocks | Oil Rocks Reservoir | Sandstone | Limestone | Granite |
By combing Eqn. (2)-(4), Darcy’s Law can also be written in terms of the fluid pressure, permeability, and the viscosity
Here, $\nabla P$ is the fluid pressure gradient. If the gravity term is negligible compared to the pressure gradient, Eqn. (5) can be simplified to
The characteristic time of the advection is the residence time, i.e., how long does a fluid particle stays within a given system. The residence time, $\tau_{a}$, can be calculated as follows:
Here Vpore is the pore volume (m3), Q is the flow rate (m3/s), L is the length, A is the cross-section of the porous media in the direction perpendicular to the flow.
Darcy’s Law is applicable at the continuum scale where a representative elementary volume (REV) is significantly larger than the average grain size. The range for the validity of the Darcy’s Law can be checked using the Reynolds number Re:
where the volumetric flow rate Q (m3/s) and p is the perimeter of a channel or grain size (m). The upper limit of the validity of the Darcy’s Law is when Re is between 1 and 10 [Bear, 2013].
Molecular diffusion is driven by concentration gradient and is described by the Fick’s First Law:
Here Jdiff is the diffusive mass flux per unit area (mol/m2/s); D0 is the molecular diffusion coefficient in water (m2/s); and x is the distance (m).
Diffusion coefficients in water are typically in the order of 10-9 m2/s and depend on chemical species, temperature and fluid viscosity. Table 2 lists the diffusion coefficients in water at room temperature for some species. Given the diffusion coefficient at a specific temperature, the diffusion coefficient at another temperature can be calculated as follows:
where $D_{0, T_{0}}$ is the diffusion coefficient (m2/s) at a reference temperature T0 (K), and $\mu_{T}$ and $\mu_{T_0}$ are the dynamic viscosities $\left(Pa\cdot s\right)$ at temperatures T and , T0 respectively.
Cations | D0 (10-9 m2/s) | Anions | D0 (10-9 m2/s) |
---|---|---|---|
H+ |
9.311 |
OH- |
5.273 |
Li+ |
1.029 |
F- |
1.475 |
Na+ |
1.334 |
Cl- |
2.032 |
K+ |
1.957 |
I- |
2.045 |
NH4+ |
1.957 |
NO3- |
1.902 |
Mg2+ |
0.706 |
HCO3- |
1.185 |
Ca2+ |
0.792 |
HSO4- |
1.331 |
Al3+ |
0.541 |
H2PO4- |
0.879 |
Fe2+ |
0.719 |
SO42- |
1.065 |
Fe3+ |
0.604 |
CO32- |
0.923 |
Fick’s second law combines the mass conservation and Fick’s first law:
where t is the time (s). Eqn. (9) can be analytically solved with appropriate boundary conditions.
Diffusion in porous media differs from diffusion in homogeneous aqueous solutions. Diffusion in porous media occurs through tortuous and irregularly shaped pores, as shown in Figure 1, and is therefore slower than that in homogeneous solutions. The effective diffusion coefficient describes diffusion in porous media and can be estimated using several forms:
Here F is the formation resistivity factor (dimensionless), an intrinsic property of porous media; the cementation factor (dimensionless) m quantifies the resistivity caused by the network of pores. Reported cementation factors vary between 1.3 and 5.0 [Brunet et al., 2013; Dullien, 1991]. For many subsurface materials, m is equal to 2 [Oelkers, 1996]. TL is the tortuosity (dimensionless) defined as the ratio of the path length in solution relative to the tortuous path length in porous media.
Dispersion describes the mixing of a solute due to fluctuations around the average velocity. This is caused by three factors: 1) microscopic heterogeneity, which make the fluid moves faster at the center of the pore and slower at the water grain boundary; 2) variations in pore sizes, in which cases fluid particles will move through larger pores faster; 3) variations in path length, causing some fluid particles going longer paths than others.
Mechanical dispersion is a result of variations in flow velocities. Dispersion coefficients in porous media is typically defined as the product of the average fluid velocity and dispersivity $\alpha$:
where u is the average flow velocity (m/s) and $\alpha$ refers to the dispersivity (m). In systems with more than one direction, the longitudinal dispersivity DL in the principle flow direction is typically higher than DT in the direction transverse to the main flow.
Dispersion is a scale-dependent process with larger dispersivity values observed at larger spatial scales. At the column scale, a typical dispersivity may be on the order of centimeters. At the field scale, apparent dispersivities can vary from 10 to 100 m, as shown in Figure 2.
The spreading of the solute mass as a result of diffusion and dispersion is similar to diffusion. This has led to the use of Fick’s First Law to describe the dispersion process as follows:
where Dh is the hydrodynamic dispersion coefficient defined as the sum of effective diffusion coefficient De and mechanical dispersion coefficient Dm:
As such, the hydrodynamic dispersion includes both diffusion and mechanical dispersion processes.
By combining the transport processes outlined above, we can derive an expression for the mass conservation of a non-reactive solute as follows:
Substitution of Eqn. (1) and (13) into Eqn. (16) yields
Equation (17) is the classical Advection-Dispersion equation (ADE). For one-dimensional systems, Equation (17) is simplified into
Analytical solution of ADE is available for homogeneous porous media [Zheng and Bennett, 2002]:
with the initial and boundary conditions:
where C0 is the injecting concentration of the tracer, and erfc(B) is complementary error function:
Péclet number (Pe) is often used to describe the relative importance of advection and dispersion/diffusion in terms of their respective time scales $\tau_{a} \text { and } \tau_{d}$:
where L is the length of the domain of interest (m), u is the average Darcy flow velocity in the direction of interest (m/s), $D_{h}$ is the dynamic dispersion coefficient (m2/s). There are also some mathematical equations to define the time scales of these processes with similar concepts, mostly depending on the selected characteristic length [Elkhoury et al., 2013; Huysmans and Dassargues, 2005; Steefel and Maher, 2009; Szymczak and Ladd, 2009]. For example, L can also be the grid spacing (m) or correlation length (m) [Huysmans and Dassargues, 2005]. As shown in Figure 3, increasing Pe values indicate increasing dominance of advective transport and sharper front in breakthrough curves.
Please watch the following video: Advection-Dispersion Equation (ADE) for non-reactive tracers (16:42)
So in is this video we'll be talking about advection-dispersion equation. This is first time actually we're not talking about chemistry but talking about the physical processes. So what do we have? Thinking about this is, let's say you have chemicals that are nonreactive. So the system you have we could use it as example is a column. So think about, for example, you are doing an experiment. You pack up a column, a column with like sand grains. And then you inject a chemical, let's say bromide, from the left into the right. So essentially you would have these chemical species you inject for a short period of time. You can imagine this is chemical will be moving along with the flow. So over time it will be eventually moving out of system with certain velocities. So v is the velocity here. And you are talking about here for example is the last of this column is L. So we want to know, as you can imagine, the constitution of this chemical would change with time and with space. Depends on what time you're looking at the snapshot. This chemical species might be here, might be here and other places there's no of this chemical species because we consider maybe, for example, we start with clean water.
So how do we mathematically solve this kind of equation and get the solution so we know the concentration of this chemical species in different time and different locations? And this is what we are going to talk about which is the advection-dispersion equation. We call it ADE with the other reaction terms. So ADE, when we think about it, I'm not going to derive in detail of where is this equation coming from. In general these equations are coming from these mass conservation principles. So if we look at these terms, different terms of the ADE, so the first home that we called the mass accumulation term. And it should have the units of accumulation, like mass per time. It's how fast things are changing. And this C is the concentration of the chemical species in the water phase. So C here is concentration of this-- let's call that tracer in water So everything we are solving is for how much they have in the water because that's what we really care about. And this should have the units of, for example, mass per volume.
So the first term is the mass accumulation. The second term here we call advective transport. So this is a process where there's a chemical almost you think about swimming. The chemical is a tracer. It actually flows together with the water at the same speed as the water flow. So that's called advective transport. And the last term is a dispersive transport. It's somewhat similar to diffusion process. It's driven by, first of all, concentration gradient, but also in this type of medias. You have grains. They tend to have different grain size, like different size of pore space and some spatial [INAUDIBLE] that leads to different flow velocity, mixing processes that lead to the concentration difference in different places. So these are three major terms in this equation. And for phi-- so this is phi is we call porosity, is how much space you have, how much pore space you have in a given volume. So the porous media has both pore space and solid space. So this is how much pore space in terms of percentage, per unit of total volume.
This v is velocity, flow velocity. It's a linear velocity in the pore space. So it's different from the u. We usually we call that's the velocity. And the relationship between the two, typically we knows that-- velocity-- we know that u will be equal to 5 porosity times the linear velocity. So there's a relationship between the two. And this Dh-- Dh is a very important term for this dispersive transport. This would be equal to-- I talk about it could be coming both from diffusion and the mechanical dispersion in porous media. So we have this kind of two terms adding together. This is coming from diffusion in porous media-- diffusion coefficient in porous media. And then the second term is essentially taking into account the mechanical dispersion. And you can think about alpha is called dispersivity, which is a parameter to describe how fast mechanical dispersion happens. And it's related to velocity. The whole term is also related to the velocity. So the flow, you actually have larger term of this.
So we have this equation. All the terms are here. And usually when this equation, there's an analytical solution for this equation if you give the right initial boundary condition. When we're numerically solving this, we will be discrete sizes of this equation in time and space, and then you get a solution for that. But before we do that, typical we need two conditions, right? This is the first time we introduced a space diminishing, x here. Before in the mineral dispersion precipitation, last time we actually induced the time. So you'll notice here that this equation, we call this is a partial differential equation, meaning it has two independent variables. One is the time. The other is space. So this is first time we introduce a space dimension here.
Now in order to solve this equation we need to know, like at t equal to 0, what are the concentrations? So this is initial condition. When t equal to 0, what are the concentrations? Now usually this is given for a given system. And we also need to know at the two boundaries, x equal to 0, x equal to L, what are the conditions that is specified? Is it, for example, a no-flow boundary or pure advective, or pure dispersive? These are different types of boundary conditions you can specify. But you will need to, in both the initial condition boundary conditions, to solve this. And different types of condition will give you different solutions because it matters what is the concentration in the t equal to 0. If your already start with something high you will see a different concentration versus, for example, at t equal to 0 you have clean water.
So terms of a solution, let's assume we have done all this work to solve the equation. What do we expect to see after we solve this equation? So I'm talking about a system. I'll be using an example system. So you will be injecting kind of shot pulse of this chemical bromide into a system. So if you conceptually think about it, at t equal to 0 it's somewhere here. And let's say everywhere it started with clean water. So you will have, at some point here, let's say you have kind of a little bit smeared at other time. And as time goes on you will have more and more smeared, but going in the direction of flow and moving along and maybe become more. So total mass will remain the same, but at the end you will see kind of wider and wider over time and over a longer distance. So this is conceptually how you would think about this solution you expect to see. Now when we think about from mathematical term, let's draw this. After we solve it, let's see. We'll look at the concentration as a function of distance. And what do you expect to see at different time?
So first of all, let's say at initial time you probably would see something like this. And here you see a pulse of this chemical. So is t at about 0, maybe a little bit past 0. But you think about, OK, over time, this pulse will move along. And then you should see different distributions. So this would be, let's say, at t equal to v times some small consolation, tu with C at another place. But this will become a little bit wider, smeared out. This is t1. And over time as you're going further distance, this becomes like more and more smooth. And I'm not drawing accurate. So total mass will remain the same. And I don't think I'm drawing nicely in terms of mass conservation. But in any case you will see over time total mass remain the same. But then the center of this will move along. So it's a v. So if it's a speed of this plume, how fast it moves over a certain time will be determined by-- so speed, v determines the speed of the center, of the rate of the center moving to downstream. But also another case will be-- we notice this Dh, which is a very important parameter as well. So if you're think about two different situations, let's say they both have the same v, but they will have different Dh values. Let's say we have another case with much higher Dh. Your will probably see something like this, a wider distribution. Still the same total mass, but it's much more spread. So this will be representing your large-- I'm sorry, a small Dh. The blue one would be representing a much larger Dh.
So now I draw the Dh. It will be more spread. So essentially you can think about Dh determines the width of the plume. So this is t1, t2. And over time, for example, if there's a t3 we care about, it should be even more spread out. This is t3. So if we wait long enough, this is going to spread a lot. So that's the type of solution you would expect to see when you have pulse of injection of a tracer. Now I'm just very briefly mentioning the characteristic time. There are several times we think is important. Once is residence time. This is directly related to how fast the flow goes and how long it actually stays in this column, how long this tracer actually stays in this column. So this we call tau a is equal to porosity times the length over mu. Or you can just say L over v. This is residence time. Another time is how long it takes for the dispersion process to uniform the whole concentration field. So this is what we call tau d. It's related to dispersion. So it's L squared divided by Dh. You can almost think about this as like if it's not open system, it would be how long it takes for diffusion to uniform, to homogenize the whole concentration field. And then a lot of times we define this Peclet number is tau d over tau a. So it's a relative between these two time scales or the rate of flow versus the rate of dispersion. And this should have a unit of L u phi Dh. This is a dimensionless number. So these times will determine like how fast-- you will realize this are kind of grouping dimensionless numbers. The Peclet number will determine the relative importance of dispersion versus advection, which in homework I asked you to do some exercises under different conditions with different Peclet numbers.
Click here for Example 6.1 CrunchFlow file package [39]
This example introduces setting up numerical simulation of the flow and transport processes for a non-reactive tracer in a 1D column of 10 cm long. A tracer Br- is injected into the column at the concentration of 1.2×10-4 mol/L. The permeability of the column is 1.75×10-13 m2 and the porosity is 0.40. A constant differential pressure is imposed at the x direction and results in a Darcy flow velocity of 4.20×10-6 m/s. The molecular diffusion coefficient in water D0 for the tracer bromide is 1.8×10-9 m2/s (between that of F- and Cl- in Table 1). The cementation exponent m is 1.0. The dispersivity α is 0.07 cm. In order to keep consistent with the value of water viscosity provided in the CrunchTope original code, the water viscosity we applied here is 1.00×10-3 Pas at 20 0C.
Click here for input and database file package [39]
For the 1-D flow, Darcy flow equation can be simplified to:
$u=-\frac{\kappa}{\mu} \nabla P$CrunchFlow setup brief: This is the first time that we set up spatial dimension in CrunchFlow. Please read CrunchFlow manual the DISCRETIZATION key word block (page 46 – 47), INITIAL_CONDITION and BOUNDARY_CONDITION (page 72) as well as the TRANSPORT keyword block (page 72-75). The CrunchFlow exercise example 6 also teaches how to set up non-reactive tracer test.
Solution 2:
In this lesson, we covered definition and principles of transport processes including advection, diffusion, and dispersion processes. We also discussed the Advection-Dispersion-Equation (ADE), its analytical solution, and how to set up solving ADE in CrunchFlow.
Flow and Transport in Porous Formations (1989) by Gedeon Dagan;
Dynamics of fluids in porous media (2013) by Jacob Bear, Courier Dover Publications;
Applied contaminant transport modeling (2002), 2nd edition, by Chunmiao Zheng and Gordon D. Bennett, John Wiley and Sons, Inc.
Bear, J. (2013), Dynamics of fluids in porous media, Courier Dover Publications.
Brunet, J. P. L., L. Li, Z. T. Karpyn, B. G. Kutchko, B. Strazisar, and G. Bromhal (2013), Dynamic Evolution of Cement Composition and Transport Properties under Conditions Relevant to Geological Carbon Sequestration, Energy & Fuels, 27(8), 4208-4220.
Dullien, F. A. (1991), Porous media: fluid transport and pore structure, Academic press.
Elkhoury, J. E., P. Ameli, and R. L. Detwiler (2013), Dissolution and deformation in fractured carbonates caused by flow of CO< sub> 2</sub>-rich brine under reservoir conditions, International Journal of Greenhouse Gas Control, 16, S203-S215.
Gelhar, L. W. (1986), Stochastic subsurface hydrology from theory to applications, Water Resour Res, 22(9S), 135S-145S.
Haynes, W. M. (2012), CRC handbook of chemistry and physics, CRC press.
Huysmans, M., and A. Dassargues (2005), Review of the use of Péclet numbers to determine the relative importance of advection and diffusion in low permeability environments, Hydrogeol J, 13(5-6), 895-904.
Oelkers, E. H. (1996), Physical and chemical properties of rocks and fluids for chemical mass transport calculations, Reviews in Mineralogy and Geochemistry, 34(1), 131-191.
Steefel, C. I., and K. Maher (2009), Fluid-Rock Interaction: A Reactive Transport Approach, Rev Mineral Geochem, 70, 485-532.
Steefel, C. I., D. J. DePaolo, and P. C. Lichtner (2005), Reactive transport modeling: An essential tool and a new research approach for the Earth sciences, Earth Planet Sc Lett, 240(3-4), 539-558.
Szymczak, P., and A. J. C. Ladd (2009), Wormhole formation in dissolving fractures, J Geophys Res-Sol Ea, 114.
Zheng, C., and G. D. Bennett (2002), Applied contaminant transport modeling: theory and practice, 2nd ed., 621 pp., John Wiley and Sons, Inc., New York.
This unit introduces general concepts of spatial heterogeneities, principles of physical flow and transport processes in heterogeneous porous media, as well as how to set up simulations for flow and non-reactive transport in a 2D heterogeneous domain in CrunchFlow.
By the end of this lesson, you should be able to:
Recommended Readings |
|
---|---|
To Do |
|
If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.
Flow and transport processes, including advection, dispersion and diffusion are described in the lesson Flow and Transport Processes in 1D System. In natural system we often need to consider these processes in multiple dimensions and in heterogeneous systems where the physical and geochemical properties of the subsurface are not evenly distributed. As an example, in Figure 1, we show that the distribution of permeability dictates the spatial distribution of injected tracer from injection wells during a biostimulation experiments at Old Rifle, Colorado (Li et al., 2010).
The general Advection-Dispersion Equation (ADE) in multiple dimensions is as follows:
Here C is the tracer concentration (mol/m3 pore volume), t is time (s), D is the combined dispersion–diffusion tensor (m2/s), v (m/s) is the flow velocity vector and can be decomposed into vL and vT in the directions longitudinal and transverse to the main flow in a 2D system. The dispersion-diffusion tensor D is defined as the sum of the mechanical dispersion coefficient and the effective diffusion coefficient in porous media De(m2/s). At any particular location (grid block), their corresponding diffusion / dispersion coefficients DL (m2/s) and DT (m2/s) are calculated as follows:
Here $\alpha_L$ and $\alpha_T$ are the longitudinal and transverse dispersivity (m). The dispersion coefficients vary spatially due to the non-uniform permeability distribution. Values of $\alpha_T$ are typically at least one order of magnitude smaller than $\alpha_L$ (Gelhar et al., 1992; Olsson and Grathwohl, 2007).
The terms homogeneity and heterogeneity are used to describe the uniformity and regularity in spatial distribution of geomaterial properties in natural subsurface systems. Homogeneity means spatially uniform-distributed properties. In natural systems, geological media are almost always heterogeneous. That is, their physical and (bio)geochemical properties vary spatially. Spatial heterogeneity can refer to both physical and geochemical properties. Physical properties include mineral grain size, porosity, and permeability / conductivity. Geochemical properties include, for example, mineral types, lithology, mineral surface area, and cation exchange capacity. Figure 2 shows a picture of an outcrop with layers of different types of geomaterials from the Macrodispersion site (MADE) in Columbia, Mississippi (Zheng and Gorelick, 2003).
Spatial heterogeneities can lead to significantly different chemical reactions and physical transport processes from their homogeneous counterparts. As a result, heterogeneities play an important role in determining processes and applications in subsurface systems. Examples include contaminant reactive transport (Li et al., 2011), oil and gas production (Chen et al., 2014; Hewett, 1986), and environmental bioremediation (Murphy et al., 1997; Song et al., 2014). For example, as shown in Figure 3, the different geomaterials in Figure 2 leads to orders of magnitude in the spatial variation of hydraulic conductivity (Figure 3A), and abnormal solute transport (Figure 3B) (Zheng et al., 2011).
Geostatistics has been well developed to quantitatively describe the characteristics of spatial physical heterogeneities, especially conductivity / permeability, in the natural subsurface. Here we briefly introduce a few geostatistical measures. Interested readers are referred to geostatistical books and software tools for more detailed information (Deutsch and Journel, 1992; Remy and Wu., 2009).
It measures the average water-conducting capability of porous media. There are different definitions widely used in subsurface hydrology:
Arithmetic mean $\kappa_T$ is defined as:
Here n is the total number of zones or subsystems; $\kappa_i$ is the local permeability of subsystem.
Harmonic mean $\kappa_h$ is defined as:
Geometric mean $\kappa_g$ is defined as:
Effective permeability$\kappa_e$ is derived from the Darcy’s Law:
where u is the average flow velocity (m/s) of the field; $\mu$ is the flow viscosity (Pa·s); $\Delta L$ is the length (m) along the main flow direction; $\Delta P$ is the pressure differential (Pa) along the main flow direction. The effective permeability describes the real water-conducting capability of the porous media. The upper bound of the effective permeability is often the arithmetic mean, and the lower bound is the harmonic mean. They correspond to flow through a perfectly layered systems, parallel and perpendicular to the layering, respectively (Heidari and Li, 2014; Zinn and Harvey, 2003).
The extent of deviation from the mean is calculated as follows:
where $\operatorname{Var}(\kappa)$ is the permeability variance (m4). Larger variance indicates larger extent of spatial heterogeneity. Permeability of geological media differs significantly in natural subsurface systems. Conductivity (K) is also commonly used as a measure of conducting capability of porous media. The Macrodispersion Experiment (MADE) site in Mississippi has both small silt and sand regions. It has a centimeter scale ln(K) variance as large as 20, while its variance at the meter-scale measured by borehole flowmeters is approximately 4.0, which is considerably smaller because small-scale variability is averaged out (Feehley et al., 2000; Harvey and Gorelick, 2000; Zinn and Harvey, 2003). The Wilcox aquifer in Texas also has high ln(K) variance of about 10, as does the Livermore site (ln(K) variance >20) (Fogg, 1986; LaBolle and Fogg, 2002). The Bordon site at Borden, Ontario, however, is relatively homogeneous with ln(K) variance around 0.20 (Mackay et al., 1986; Sudicky, 1986). Note that ln(K) mentioned here is the natural logarithm of hydraulic conductivity K (m/s).
Correlation length is a measure of correlation in spatial variation. Two points that are separated by a distance larger than the correlation length have fluctuations that are relatively independent, or uncorrelated. In contrast, properties of two points that are within correlation length are correlated. For detailed calculation and application of correlation length in subsurface fields, readers are referred to (Bruderer-Weng et al., 2004; Gelhar, 1986; Heidari and Li, 2014; Salamon et al., 2007). Correlation lengths differ significantly in natural subsurface systems, as shown in Table 1. Figure 4 shows the patterns of permeability fields (m2) with different correlation length generated from Gaussian Sequential Simulation Method.
Source | Geological Media | Correlation Length (m): Horizontal | Correlation Length (m): Vertical | Overall Length (m): Horizontal | Overall Length (m): Vertical |
---|---|---|---|---|---|
Socorro, New Mexico | fluvial sand | >3 | 0.1 | 14 | 5 |
Rio Grande Valley, New Mexico | fluvial sand | 7.6 | 760 | ||
Cape Cod, Massachusetts | glacial outwash sand | 5 | 0.26 | 20 | 5 |
East central Illinois | sandstone aquifer | 4.5×104 | -- | 5×105 | -- |
Las Vruces, New Mexico | alluvial silty-clay loam soil | 0.1 | -- | 6 | -- |
A series of hydraulic connectivity have been defined as integrative measures of spatial heterogeneity characteristics. Static connectivity can be calculated by the connective structure of the hydraulic conductivity or geological facies. Dynamic connectivity directly relates to physical flow or transport processes. For interested readers, we refer them to (Knudby and Carrera, 2005; Knudby and Carrera, 2006; Renard and Allard, 2013; Siirila-Woodburn and Maxwell, 2014; Willmann et al., 2008).
Setting up 2D heterogeneous domain for flow and transport calculation generally includes three primary steps. First, the 2D domain needs to be defined with the targeted size, dimension, and resolution. The relevant keyword blocks for setting up the domain include DISCRETIZATION, INITIAL_CONDITION, and BOUNDARY_CONDITION. The second step involves the set up for the calculation of flow, which involves the use of various keywords in the FLOW key block, including constant_flow, calculate_flow, read_permeabilityFile, and pressure. For 2D heterogeneous systems we always use “calculate_flow” to calculate flow by specifying a pressure gradient using “pressure” at the two main boundaries. The permeability distribution can be either defined in INITIAL_CONDITIONS with relevant permeability values in specific grid blocks defined under GEOCHEMICAL_CONDITIONS, or by reading in permeability distribution using the keyword read_permeabilityFile. Detailed format of the permeability file is provided in the manual. After that, the transport keywords can be specified in the TRANSPORT keyword block, similar to the ones discussed in the 1D homogeneous setup.
Click here [41] for CrunchFlow Files for Example 8.1 and other related files for Li et al. 2014
The following gives an example of setting up flow and transport calculations in 2D homogeneous and heterogeneous domains. Here we use the physical set up of 2D cross-sections of the Mixed and One-zone column in (Li et al., 2014). The authors packed 4 columns of 2.5 cm in diameter and 10.0 cm in length with relatively similar amount of magnesite and quartz distributed in different spatial patterns. Here we focus on two columns that represent two extreme cases: the Mixed column with uniform distributed magnesite and quartz, as well as the One-zone column with magnesite all gathered in one cylindrical zone of diameter of 1.0 centimeter. To keep it simple, here we will do the calculation for 2D cross-sections, instead of following the steps in section 3.2.4 in the paper to convert 2D to 3D. So our numbers might be different from the paper. We are also assuming that in the middle one zone magnesite is 100% of the solid phase. The 2D system has a size of 25 mm by 100 mm. A constant differential pressure is imposed at the boundaries in the z (vertical) direction, leading the primary water flow direction in the z direction from the bottom to the top. No flux boundaries are specified in the X direction (Li et al., 2014).
Columns | Mg zone |
Qtz zone | $^2a_L$ (cm) |
$^4a_T$ (cm) |
$\mathrm{\kappa}_e$ (×10-13m2) |
$\phi_{\mathrm{avg}}$ |
---|---|---|---|---|---|---|
Mixed | 0.05 | 0.005 | 8.26 | 0.44 | ||
One-zone | Width: 1.0 cm
$\Phi_{Mg:}\ 0.54$ |
Width: 1.5 cm $\Phi _{Qtz:}\ 0.38$ |
0.07 | 0.004 | 10.74 | 0.44 |
*The permeability of the pure sand columns of the same grain size was measured to be $8.7 \times 10^{-13} \mathrm{~m}^{2}$.
Please do the following:
Run the simulation at the grid block resolution of 1 mm by 1 mm.
Note: in plotting 2D figures, you can use softwares such as Tecplot or Matlab. Penn State provides free Matlab web access on WebLabs at Penn State [42]. You can access matlab from there. In addition, you can also get video tutorials from LinkedInLearning [43] at Penn State about using matlab.
In Heidari and Li (2014), 2D flow experiments and modeling were used to understand non-reactive solute transport in heterogeneous porous media with different spatial patterns. There are three two-dimensional (2D) sandboxes (21.9 cm × 20.6 cm) packed with the same 20% (v/v) fine and 80% (v/v) coarse sands in three patterns that differ in correlation length: Mixed, Four-zone, and One-zone. The Mixed cases contain uniformly distributed fine and coarse grains. The Four-zone and One-zone cases have four and one square fine zones, respectively (Figure 8).
Read the paper carefully and follow the example in this lesson to generate the 2D flow fields and tracer breakthrough curves for the three sandboxes of high permeability contrast (HC). All parameters and properties of the three flow cells are in Table 1 and Table 2 in the paper. You can compare your modeling output of the flow field and breakthrough curves with Figure 5 (Flow field and local breakthrough curves) and Figure 6 (overall breakthrough curves from 2D model) in the paper. Note that one pore volume is the same as one residence time.
Apologies for no original CrunchFlow files for this hw. you can follow Example 8.1 and similarly set it up for this homework.
In this lesson we studied general concepts and geostatistical measures of spatial heterogeneities, as well as how to set them up in 2D domains.
You have reached the end of Lesson 8! Please make sure you have completed all of the lesson activities before you begin Lesson 7.
Boulanger, R. and Duncan, J. (2000) Geotechnical Engineering Photo Album: University of California at Davis, accessed March 15, 2006.
Bruderer-Weng, C., Cowie, P., Bernabé, Y. and Main, I. (2004) Relating flow channelling to tracer dispersion in heterogeneous networks. Adv. Water Resour. 27, 843-855.
Chen, C., Balhoff, M.T. and Mohanty, K.K. (2014) Effect of Reservoir Heterogeneity on Primary Recovery and CO2 Huff'n'Puff Recovery in Shale-Oil Reservoirs. SPE Reservoir Evaluation & Engineering.
Deutsch, C.V. and Journel, A.G. (1992) Geostatistical software library and user’s guide. New York.
Feehley, C.E., Zheng, C. and Molz, F.J. (2000) A dual‐domain mass transfer approach for modeling solute transport in heterogeneous aquifers: Application to the Macrodispersion Experiment (MADE) site. Water Resour Res 36, 2501-2515.
Fogg, G.E. (1986) Groundwater Flow and Sand Body Interconnectedness in a Thick, Multiple‐Aquifer System. Water Resour Res 22, 679-694.
Gelhar, L.W. (1986) Stochastic subsurface hydrology from theory to applications. Water Resour Res 22, 135S-145S.
Gelhar, L.W., Welty, C. and Rehfeldt, K.R. (1992) A critical review of data on field‐scale dispersion in aquifers. Water Resour Res 28, 1955-1974.
Harvey, C. and Gorelick, S.M. (2000) Rate‐limited mass transfer or macrodispersion: Which dominates plume evolution at the macrodispersion experiment (MADE) site? Water Resour Res 36, 637-650.
Heidari, P. and Li, L. (2014) Solute transport in low‐heterogeneity sandboxes: The role of correlation length and permeability variance. Water Resour Res 50, 8240-8264.
Hewett, T. (1986) Fractal distributions of reservoir heterogeneity and their influence on fluid transport, SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
Knudby, C. and Carrera, J. (2005) On the relationship between indicators of geostatistical, flow and transport connectivity. Adv Water Resour 28, 405-421.
Knudby, C. and Carrera, J. (2006) On the use of apparent hydraulic diffusivity as an indicator of connectivity. J Hydrol 329, 377-389.
Koltermann, C.E. and Gorelick, S.M. (1996) Heterogeneity in sedimentary deposits: A review of structure-imitating, process-imitating, and descriptive approaches. Water Resour Res 32, 2617-2658.
LaBolle, E.M. and Fogg, G.E. (2002) Role of molecular diffusion in contaminant migration and recovery in an alluvial aquifer system, Dispersion in Heterogeneous Geological Formations. Springer, pp. 155-179.
Li, L., Gawande, N., Kowalsky, M.B., Steefel, C.I. and Hubbard, S.S. (2011) Physicochemical heterogeneity controls on uranium bioreduction rates at the field scale. Environ Sci Technol 45, 9959-9966.
Li, L., Salehikhoo, F., Brantley, S.L. and Heidari, P. (2014) Spatial zonation limits magnesite dissolution in porous media. Geochimica Et Cosmochimica Acta 126, 555-573.
Li, L., Steefel, C., Kowalsky, M., Englert, A. and Hubbard, S. (2010) Effects of physical and geochemical heterogeneities on mineral transformation and biomass accumulation during a biostimulation experiment at Rifle, Colorado. J Contamin Hydrol 112, 45 - 63.
Mackay, D.M., Freyberg, D., Roberts, P. and Cherry, J. (1986) A natural gradient experiment on solute transport in a sand aquifer: 1. Approach and overview of plume movement. Water Resour Res 22, 2017-2029.
Murphy, E.M., Ginn, T.R., Chilakapati, A., Resch, C.T., Phillips, J.L., Wietsma, T.W. and Spadoni, C.M. (1997) The influence of physical heterogeneity on microbial degradation and distribution in porous media. Water Resour Res 33, 1087-1103.
Olsson, Å. and Grathwohl, P. (2007) Transverse dispersion of non-reactive tracers in porous media: a new nonlinear relationship to predict dispersion coefficients. J Contam Hydrol 92, 149-161.
Remy, N., Alexandre Boucher, and Wu., J. (2009) Applied Geostatistics with SGeMS. Cambridge University Press.
Renard, P. and Allard, D. (2013) Connectivity metrics for subsurface flow and transport. Adv Water Resour 51, 168-196.
Salamon, P., Fernàndez‐Garcia, D. and Gómez‐Hernández, J. (2007) Modeling tracer transport at the MADE site: the importance of heterogeneity. Water Resour Res 43.
Siirila-Woodburn, E.R. and Maxwell, R.M. (2014) A heterogeneity model comparison of highly resolved statistically anisotropic aquifers. Adv Water Resour.
Song, X., Hong, E. and Seagren, E.A. (2014) Laboratory-scale in situ bioremediation in heterogeneous porous media: Biokinetics-limited scenario. J Contam Hydrol 158, 78-92.
Sudicky, E.A. (1986) A natural gradient experiment on solute transport in a sand aquifer: Spatial variability of hydraulic conductivity and its role in the dispersion process. Water Resour Res 22, 2069-2082.
Willmann, M., Carrera, J. and Sánchez‐Vila, X. (2008) Transport upscaling in heterogeneous aquifers: What physical parameters control memory functions? Water Resour Res 44.
Zheng, C., Bianchi, M. and Gorelick, S.M. (2011) Lessons Learned from 25 Years of Research at the MADE Site. Ground Water 49, 649-662.
Zheng, C.M. and Gorelick, S.M. (2003) Analysis of solute transport in flow fields influenced by preferential flowpaths at the decimeter scale. Ground Water 41, 142-155.
Zinn, B. and Harvey, C.F. (2003) When good statistical models of aquifer heterogeneity go bad: A comparison of flow, dispersion, and mass transfer in connected and multivariate Gaussian hydraulic conductivity fields. Water Resour Res 39.
In this lesson we use a mineral reaction example to illustrate 1) how physical processes and geochemical processes are coupled in the short and long time scales, and 2) how to set up physical and geochemical processes in a 1D column in Crunchflow. The set up combines mineral dissolution (and precipitation) and physical flow processes in previous lessons. Please review these two lessons if needed.
By the end of this lesson, you should be able to:
References (optional) | Review reading materials and examples on transport in 1D and mineral dissolution and precipitation |
---|---|
To Do |
|
Mineral dissolution and precipitation occur ubiquitously in natural systems. As minerals dissolve, chemicals in rocks are released into water, reducing solid mass and increasing aqueous concentrations. Mineral dissolution also opens up porosity for water storage. The opposite occurs when mineral precipitates and clog pore spaces. In the deep subsurface where we extract energy (e.g., oil, gas, geothermal energy) from, mineral dissolution and precipitation change the geochemical surface and water-conducting capacity, therefore regulating the energy extraction processes. Over geological times, the chemical and physical property alteration transforms rocks into soils and is called chemical weathering. Chemical weathering, together with physical weathering (erosion), shapes the Earth’s surface. Chemical weathering consumes carbon dioxide and locks them into carbonate minerals, which regulates the atmospheric CO2 level and sustain relatively clement earth conditions (Kump et al., 2000).
A bedrock column is 10.0 meter deep with a porosity of 0.20. Mineral composition is simple with 80% v/v of quartz and 20% v/v of K-Feldspar. As the rainwater infiltrates through the rock, these primary minerals dissolve whereas a secondary mineral kaolinite precipitates, as listed in Table 1.
The annual net runoff (precipitation – evapotranspiration) is 0.5 meter/year. Although we know that not all runoff infiltrates through the bedrock, for simplification here we assume the runoff flows through the bedrock and exits from the bottom of the bedrock. Assume the diffusion coefficients of all species are 1.0E-10 m2/s, and dispersivity is 0.1 m. The cementation factors are the default values of 1.0.
Let’s assume the initial pore water composition (all concentrations in molal): pH 7.0, Cl- (charge balance), K+ (4.0E-5), Na+ (4.0E-5), Al3+ (4.0E-6), HCO3- (1.0E-5), and SiO2(aq) (8.0E-5). One of the major drivers of chemical weathering is acidity. The acidity in soils comes primarily from the soil CO2 gas from root respiration and organic acids from root exudates. Here we set the total HCO3- concentration (equivalent to total inorganic carbon (TIC)) in the rainwater to 0.3 mM, which is equilibrated to soil gas CO2 of 9000 ppmv.
Assume chemical weathering occurs for 0.1 million years. Please set up the simulation in Crunchflow. Plot the profiles of pH, Saturation index $(\mathrm{SI}=[\log10(\mathrm{IAP}/\mathrm{K_{eq}})])$, $\phi$ ( volume of mineral/ volume of total porous media) of minerals and $\tau$ values of K and Si (relative to Zr in this soil column, we will explain later) at t = 0, 5,000, 10,000, 20,000, 40,000, and 100,000 years. How does the chemical composition and physical properties (e.g., porosity, permeability) evolve over time?
Setting up this lesson will involve what you learned from the lesson on mineral dissolution and precipitation, and the lesson on 1D transport. If needed please revisit these two lessons.
Primary mineral dissolution | log10Keq |
log10k (mol/m2/s) |
Specific surface Area (m2/g) |
---|---|---|---|
$\mathrm{SiO}_{2(\mathrm{s})}\text{(quartz)}\ \leftrightarrow\mathrm{\ SiO}_2\text{ (aq) }$ |
-4.00 | -13.39 | 0.01 |
$\begin{array}{l} \mathrm{K} \mathrm{AlSi}_{3} \mathrm{O}_{8} \text { (K-Feldspar) }+4 \mathrm{H}^{+} \leftrightarrow \\ \mathrm{Al}^{3+}+\mathrm{K}^{+}+2 \mathrm{H}_{2} \mathrm{O}+3 \mathrm{SiO}_{2}(\mathrm{aq}) \end{array}$ |
-0.28 |
-13.00 -9.80aH+0.5 -10.15aOH-0.54 |
0.01 |
Secondary mineral precipitation | |||
$\begin{array}{l} \mathrm{Al}_{2} \mathrm{Si}_{2} \mathrm{O}_{5}(\mathrm{OH})_{4} \text { (Kaolinite) }+6 \mathrm{H}^{+} \leftrightarrow \\ 2 \mathrm{Al}^{3+}+2 \mathrm{SiO}_{2} \text { (aq) }+5 \mathrm{H}_{2} \mathrm{O} \end{array}$ |
6.81 | -13.00 | 0.01 |
*: Thermodynamic and kinetic parameters from EQ3/6 database [Wolery, 1992]. The three different rate constant values indicate the rate constants for neutral, acidic and basic kinetic pathways, respectively.
A general governing equation for the mass conservation of a chemical species for a primary species i is the following:
where $\phi$ is porosity, Ci is the concentration (mol/m3), Dh is the hydrodynamic dispersion (m2/s), v is the flow velocity (m/s), ri,tot is the summation of rates of multiple reactions that the species i is involved (mol/ m3/s), and Np is the total number of primary species. Here we notice that this equation is very similar to the advection-dispersion equation that we discussed in previous lessons except for two major differences. One is that it has an additional reaction term ri,tot to take into account sources and sinks coming from mineral dissolution and precipitation. The other difference is that this is written for the species i, which is one of the Np primary species. Note that if the species is K+, the ri,tot in Equation (3) will be the summation of the K-feldspar reaction rate rK-feldspar and Kaolinite reaction rate rKaolinite. This is because K+ is involved in both of the kinetic reactions, the K-feldspar and Kaolinite reactions (either dissolution or precipitation). The mineral reaction rate laws follow the Transition State Theory (TST):
Where kH2O, kH, and kOH are the reaction rate constants with values of 10-13.0 mol/m2/s, 10-9.80 mol/m2/s under acidic conditions, and is the reaction rate constant of 10-10.15 mol/m2/s under alkaline conditions as shown, respectively (Table 1). The exponents of the H+ and OH- activities indicate the extent of rate dependence on H+ and OH-. A is the reactive mineral surface area (m2/m3 porous media). For kaolinite we have
Where kH2O is the reaction rate constant of 10-13.0 mol/m2/s as shown in Table 1.
In the chemical weathering context, K-feldspar dissolves and releases Al3+, K+, SiO2(aq), whereas Kaolinite precipitates and is the sink for Al3+ and SiO2(aq). The rate term in the parenthesis (1-IAP/Keq) is positive when minerals dissolve and is negative when minerals precipitate so we don’t need to change the sign of the rate term in the equation.
Species. The chemical species from the minerals are Al3+, K+, SiO2(aq). The species from the rainwater include H+, OH-, CO2(aq), HCO3-, CO32-, Na+, and Cl-. Here the major anions are the carbonate species so we can assume the aqueous complexes are KHCO3 (aq) and KCl (aq) and we have a total of 12 species. We can pick Al3+, K+, SiO2(aq), H+, HCO3-, Na+, and Cl- as primary species (7). This means we have secondary species including KHCO3(aq), KCl(aq), OH-, HCO3-, and CO32-. We need to solve for 12 (total number of all species) – 5 (total number of secondary species) = 7 primary species (governing equations).
As the rainfall contacts the rock, cations and anions are released from the rock. The processes change properties of both solid and aqueous phases. However, the solid rock is dense in terms of amount of chemical mass compared to the aqueous phase so that the solid phase properties change at rates that are orders of magnitude slower than the aqueous phase. The system therefore often reaches steady state where the concentrations of species relatively constant over time.
The relative magnitude of advection, diffusion/dispersion and mineral dissolution are often compared using dimensionless numbers. The dimensionless numbers are quantitative measures of the relative magnitude of the characteristic time scales of different processes. In this case we have three characteristic time scales: the residence time related to advection $\tau\ \text{advection}=\frac{l}{v}$, the time for diffusion/dispersion $\tau\ \text {diffusion/dispersion}=\frac{l^{2}}{D_{h}}$, and the reaction time to reach equilibrium $\tau\ \text {reaction }=\frac{C_{e q}}{k A}$. Here we introduce two dimensionless numbers comparing the different time scales: DaI (Damköhler number for advection), DaII (Damköhler number for dispersion):
where l is the characteristics length of a domain (m), k is the rate constant (moles m−2 s−1), A is the reactive surface area (m2m−3), Dh is the diffusion coefficient in porous media (m2 s−1), and Ceq is the mineral solubility (moles m−3). When Damköhler number >>1, reaction times are much faster than transport times so that the system is transport controlled. The relative importance of advective versus dispersive transport is compared by the Péclet number that we have introduced before.
when Pe>>1, advective transport dominates and dispersive and diffusive transport are negligible.
If you have paid attention to a roadside cut, you will notice that the color, texture or structure along the direction of flow varies. This is because rocks at different depth have been subject to different extent of chemical weathering, resulting in gradients in mineral composition, reaction tendency, and porosity. Soils are generated over thousands to millions of years of chemical and physical alteration of their parent bedrock, which typically have much lower porosity and permeability than soils. Why and how are soils at different depth weathered differently? How do different parameters control chemical weathering?
Using mass conservation principle, the code calculates the change in mass and volume in each mineral phase, which updates individual mineral volume fractions using the following equation:
Where $\phi_{j m, t}$ is the volume fraction of each mineral phase jm, Vjm is the molar volume of mineral jm, $r_{j m, t+\Delta t}$ is the the reaction rate of the mineral jm at time $t+\Delta t$. This can be done for each grid block in the domain.
The porosity at any time $t$ can be updated as follows:
Where $\phi_t$ is the porosity at time $t$ and mtot is the number of solid phases.
The reactive surface area for mineral jm at time $t$, $A_{j m, t}$ is calculated based on the change in porosity and mineral volume fraction compared to their values at initial time 0 (Steefel et al., 2015):
The local permeability in individual grid blocks is calculated using local porosity values from equation and the Kozeny-Carman equation (Steefel et al,. 2015):
Where $K_{t+\Delta t}$ is the permeability at time $t+D t$updated from time $. With updated permeability, flow velocities can be updated using Darcy's law. In this example, however, we use constant flow velocity for simplicity.
Within a soil profile, $\tau$ is used to assess the effect of the chemical and physical alteration. Note that this is different from the time scale $\tau$ that we used previously on dimensionless numbers. We have to stick to this notation because the chemical weathering community uses $\tau$ for mass transfer coefficient, which is calculated by:
where $\tau_{i, j}$ is the mass transfer coefficient of element j relative to reference element i. cj,w and cj,p are the concentration of element j in weathered soil and in parent rock, respectively. ci,w and ci,p are the concentration of element i in weathered soil and in parent rock, respectively. The element i is considered as an immobile reference element to exclude the effects of the physical forces such as expansion/compaction. Titanium and zirconium are usually considered as immobile elements in calculating mass transfer coefficient. A negative $\tau_{i, j}$ value indicates mass loss while a positive $\tau_{i, j}$ value indicates mass enrichment for element j relative to element i. When $\tau_{i,j}=-1$, element j has been completely depleted as a result of the chemical weathering process. In calculating $\tau$ of K and Si relative to Zr, it is typically assumed that the soil column always has a constant concentration of immobile reference element Zr, which means $\frac{c_{Z r, p}}{c_{Z r, w}}=1$ all the time.
Assume the chemical weathering occurred for 0.1 million years. Please set up the simulation in Crunchflow. Plot the profiles of pH, Saturation index (SI, equal to [log10(IAP/Keq)]), $\phi$ ( volume of mineral/ volume of total porous media) of minerals and $\tau$ of K and Si relative to Zr in this soil column at the initial time, 5 thousand years, 10 thousand years, 20 thousand years, 40 thousand years and 100 thousand years. In calculating $\tau$ of K and Si relative to Zr, please assume the soil column always has a constant concentration of immobile reference element Zr, which means $\mathrm{C}_{\mathrm{Zr}, \mathrm{p}} / \mathrm{c}_{\mathrm{Zr}, \mathrm{w}}=1$ all the time. Discuss what you observe.
Setting up this lesson will involve what you learned from the lesson on mineral dissolution and precipitation, as well as the lesson on 1D transport. If needed please revisit these two lessons. In this lesson we do not provide videos and template anymore.
Please watch the following video: Mineral Dissolution and Precipitation 1D Columns (20:09)
Click here for CrunchFlow files [44]
Li Li: OK, so this lesson we're going to talk about mineral dissolution, precipitation in one dimension, like, 1D columns. So this is different from the previous lessons that we had before on mineral dissolution and precipitation, which was in batch reactor, well-mixed systems. So it is really no dimension. Or you can call it a zero dimension because, in every spatial point, they are the same. So here, we are talking about 1D. What that means is that, essentially, you have a column. So we're talking about processes like, for example, rainfall comes in, hitting parent rock, and things dissolving out. And then, over a long period of time, these chemical species release out. And rock becomes soil.
So I'm going to use this as one example to talk about mineral dissolution and precipitation in 1D columns and how we set up these equations. But it's not really only about chemical weathering. It's really about these mineral dissolution precipitation reactions in general. In the short-terms, these reactions actually change, for example, water chemistry. If you think about a ground-water system, it would actually tend to-- we can see the change if there is mineral dissolution precipitation -- happen pretty fast. But in order to see change in the solid phase, we do need to wait a bit longer for enough change to happen. So today, I'm going to use this system. So I'll talk about this. If we're talking about rock transformation to become soil, we are setting up a system that has rainfall come in at a rate, maybe, for example, some way that is about annual precipitation minus e-t, something like that, at this flow velocity. And first thing, usually, we would have an unsaturated system in these topsoils, but for simplicity, we would just think about everything at the fully saturated system.
So it's easy at the beginning. So we have these systems that have essentially two minerals as parent rock. One is quartz -- and dissolves pretty slowly. Everyone knows about that. The other is K-feldspar, which is a very common silicate in rocks. K-feldspar has a formula of KAlSi and then O3-8. So essentially, you can think about this system-- if you think about rainfall, it comes in, dissolving out, and then releases these chemical species, that you can think about several processes as happening at the same time. One is advection. So this is the same advection we talked about in one unit set with a chaser, essentially, as the flow brings out the chemicals, but also there's a dispersion, diffusion. This again, is the same as what we talked about last time in the AD equations. But then the last one that is different from the previous, one other reaction, which is mineral dissolution precipitation reaction. So essentially this is like combining the AD equation unit with the mineral dissolution precipitation reaction unit. So if you need to go back, you're welcome to review this material. But essentially, because these three different processes are happening at the same time-- but then we also have multiple minerals -- we know the qualities are pretty slowly. And we can more or less think about this as almost non-reactive. Because compared to K-feldspar, it's a relatively very slow reaction. And when we think about the reaction that's happening, the K-feldspar will be dissolving out. So you would have this K-feldspar dissolving out to become released calcium, which is very important nutrients, metals, cations. And then it would be also released aluminum plus SiO2, so that's aqueous. These are species that are going to be coming out.
But at the same time, you also think about, some of these chemical species, when they reach a certain level, actually they can precipitate out as solid phase again. For example, aluminum plus SOi2 actually will become another mirror, which we call kaolinite. So essentially it will have the formula of Ai-Al2-Si2-O5-OH4. This is kaolinite. This is one type of very common clay in our system. Or if you go digging some soil, you will very easily see this type of clay. So essentially, you can think about the whole process as, like, you have parent bedrock. Rainfall comes in, releases some of the chemical species. But then some of these species also re-precipitate out. So it's almost like a redistribution of minerals. But at the same time, because some of species dissolving out over a long period of time, it would change the properties of the system, for example, porosity, permeability. Make this solid phase more permeable and have more pore space for water to flow through. So there's a change in physical properties as well.
So one question we often ask is how does this water and rock change over time and also over space? We want to know how fast things change and how much they change and what they have become. And a lot of time, we want to be able to predict that. So in order to do that, we have to come set up the equations, reactions, and all of that, putting everything together, and think about how we solve them.
So again, with the system, every time we need to think about the chemical system and how they evolve over time and space, we think-- the first thing we need to think about-- the chemical species first. So you have several primary species again. If you review back some of the material we covered before, we had primary species, then existent system-- because we have aluminum, silica-- that they have to be there. This has to be part of the building block of the system. You have potassium too. And when the rainfall comes in, usually it brings some acidity. So H-plus should be there. You would also have CO2-aq. In addition, a lot of times, these rainfalls have a little bit of salt in them. So to be representative, we are also putting a bit of sodium chloride in that example-- sodium chloride as part of the rainfall. That means, also, in the second species, we need to decide, OK, what complexities would be formed, how much can complex reaction happen, and things like this. So to keep it simple, I'm going to just have potassium and then HCO3-aq and maybe KCl-aq as two complexities. And then, of course, we also need other species, for example, H-minus. We need the other form of bicarbonate and carbonate. So it's quickly become a long list of species.
So if we are talking these chemical species, we'll be essentially solving for 1, 2, 3, 4, 5, 6, 7. We'll have seven primary species. So that essentially means we would have seven independent equations to solve for. So what I'm writing here is a general type of governing equation, for one species, I. And again, this first term-- we call, mass accumulation term, that we discussed before. Is a summation of overall change. And then you have the advection term, diffusion term, dispersion term. These doesn't change, because of the system. But for example, the primary would change, but the terms wouldn't change. You could have different porosity phi or velocity or diffusion-dispersion coefficient for a different chemical species. But the term itself wouldn't change. And then the last term is new. It's a reaction term. So it would take into account, essentially, different types or the rates of different types of reactions that this species, I, is involved in. And you think about how fast these different reactions would change in the continuation of species I. And of course, this I needs to be written for 1 to Np, which is the number of primary species. So you need to write, for this system, seven different equations to solve. But then, on top of that, again, you have this number of secondary species you need to solve for in the form of algebraic relationships.
So how does this r-i-j look? These are terms different from the previous AD equations, that without this term. Now here, this r is a [INAUDIBLE] term. So you have that. So for this, r's, it needs to be-- everything about this particular r-- let's make this i is potassium. So for potassium, means this system-- you can think about two reactions involved in. One is that the K-feldspar dissolving out to release out calcium. And then it comes out. And it goes. So actually, one reaction-- so it would be just K-feldspar dissolving out. So if it's K, then you would just have this rate. And then you follow TST rate all right. So you have, probably, this complicated K-H-plus, dependence on H-plus, whatever term, empirical numbers, K-H2O and then a water plus KOH-minus and then Activity OH-minus-- some empirical exponent here as well. These are the regular part, early part, that depends on pH. But then, again, you have surface area. And the last term would be 1 minus IAP over K-aq. Again, it's how fast this is close to equilibrium or far away from equilibrium. So that's the general form of reaction rate law. So when we solve equals to this term, we'll be putting there-- so this is going to be a complicated term. But also, if you think about, for example, species like aluminum or silica, then they're actually involving two different reactions. One is K-feldspar dissolving out to release these. But a kaolinite also precipitated out. So there's one source reaction term and r-i term. And the other would be a sink term for kaolinite precipitating out. So it's a sink for aluminum or silica.
So we solve all of these equations. And what we come up with is all these Ci as a function of time and space. So you have X versus-- essentially, once we know all the parameters, initial conditions, bounded conditions-- what water comes in, what water comes out the beginning, what's the pore water composition in the rock and all that, this will give us this distribution of concentration from time and space. And over geological time, when this changes line-- and when I say geological time, I mean at least thousands of years, thousands or 10 thousands or close to 1 million years, things like that. So over this time scale, you would have these things keep on changing. And then the rock becomes soil, property of the porous media change. You have porosity increases typically. And you have, also, permeability change, because the water-- the resistance of this porous media to water becomes lower and lower. So you have permeability increase as well. So this would be the chemical weathering over long-term. Over short-term, is all of these reaction changes or water chemistry.
But the last part of it, I want to mention a little bit about the dimensionless number, from the mass point-of-view, how these would change, how we analyze these systems. So last time, when we talked about this equation of AD with other reaction term, we talked about the Pe or Peclet number, which is the tau of the advection. It's a tau diffusion-dispersion over tau advection. Tau is a time scale, right? So the time scale of diffusion and dispersion, which is L squared over D-- so it's length of this column-- over diffusion-- because that's the time scale for diffusion-dispersion-- divided by-- or times, if I do it this way-- times 1 over tau-a, which is the opposite of lance over advection. So you would have v-L. So that cancels out. It becomes L-v over D. And so that's a Peclet number, that we talked of before. So this term essentially compares the relative magnitude of advection versus diffusion-dispersion. Which point is dominant under phosphoro or slow conditions? But then, because we have the reaction term, we have two more different dimensionless numbers, we call Damkohler numbers. So the first step-- because think about, we have two transport processes. One is diffusion-dispersion. Other is the advection. So the first one is comparing the time scale of reaction-- the relative time scale of advection to the relative time scale of reaction. And in that case, again, tau-a would be the same as L over v. But then, you'd be-- this is over tau reaction. For reaction, we are thinking about how much things can dissolve and, eventually, reach equilibrium. So you would have total bulk volume times porosity, which is how much pore space you have that can hold water. So that's a lot of volume times the maximum concentrations and you can dissolving out, divided by the rate constant as a maximum, K times A. So this would be, essentially, the reaction time, the characteristic time for reaction. And in notes, I have simplified this. And you can have that. You can look at these. For the second Damkohler number, you would be comparing the time scale of diffusion-dispersion over tau-r. And again, this would be going back to L squared over D, divided by the same thing. Because the reserve reaction term would be the same.
So if you think about-- we typically talk about, under fast reaction system, fast compared to transfer versus when you have very slow reaction systems. So when you have very fast reaction systems, the systems tend to reach equilibrium quickly. And we tend to be transporter-controlled. It would be determined by the bottleneck of all processes, which mean it's going to be transport-controlled. If the reaction is actually slower-- so you will see, by this analysis, you will be able to tell how the system would behave, whether it's transport-controlled or surface reaction-controlled. And how these will be changing over distance and time will give us some kind of grouping analysis that we can do.
So I'm going to end here. I'm going to go back and look at notes. And this, I went through pretty fast. So look at notes we have. Do some homework. And you will realize, these analyses will really help you to understand how things are different under different conditions.
Figure 2 Evolutions of (A) pH (B) porosity (C) Saturation Index $\left[\log 10\left(I A P / K_{e q}\right)\right]$ of K-Quartz, (D) $\phi$(Quartz) (volume of Quartz / volume of total porous media) (E) Saturation Index $\left[\log 10\left(I A P / K_{e q}\right)\right]$ of K-Feldspar (F) $\phi$ (K-Feldspar) (volume of K-Feldspar / volume of total porous media) (G) Saturation Index $\left[\log 10\left(I A P / K_{e q}\right)\right]$ of Kaolinite, (H) $\phi$(Kaolinite) (volume of Kaolinite / volume of total porous media), $\text { (I) } \tau(K)$ and $\text { (J) } \tau(S i)$ as a function of soil depth at different times.
Initially, the pH in the rock is 6 and porosity is 0.2. After the simulation starts, rainwater flows into the soil column from the top with a lower pH so the pH in the upper column (close to the surface) decreased. In the lower part of the column, incoming hydrogen ion is consumed during mineral dissolution so pH increases. As mineral dissolves, porosity is enlarged. Porosity increases more where mineral dissolution is faster. Because rainwater has lower pH and lower solute concentrations, Quartz is under saturated and dissolves from the top. However, because the lower part of column has a higher pH, Quartz is re-precipitating in the lower part of column. K-Feldspar is under saturated throughout the column so it always dissolves. Kaolinite precipitates as a secondary mineral as K-Feldspar and quartz dissolve and elevate the aqueous concentrations.
From the discussion above it is apparent that soil at different depth weathers differently primarily due to the availability of reactants for the dissolution reactions at different depth. The delivery through flow and transport and the consumption of reactants compete with each other in the soil profile, which determines the location and shape of the weathering front. In fact, the extent and speed of chemical weathering can often be interpreted inversely based on understandings of the vertical profile through numerical simulations [Brantley and White, 2009].
assess the role of dissolution rates, dissolution kinetics, rainwater chemistry, and rainwater abundance.
Continuing along example 1, please analyze the role of different parameters/conditions in determining chemical weathering rates and profile. In each of the questions below, please compare the base case in the example to two more cases with different parameter values. In each question, please keep all other parameters and conditions the same so we can fully focus on the effects of the changing parameter. In each question, draw the depth profiles of porosity, volume fractions of Quartz, K-Feldspar, and Kaolinite, and $\tau$ figures at 100,000 years:
Click here for CrunchFlow Solutions. [45]
In a homogeneous column at the length of 10.0 cm and diameter of 2.5 cm, calcite and sand grains are uniformly distributed in the column. Detailed physical and geochemical properties are in Table 2.
Parameters | Values |
---|---|
Calcite (gram) | 14.62 |
Quartz (gram) | 76.10 |
Grain Size Calcite ($\mu m$) | 225-350 |
Grain Size Quartz ($\mu m$) | 225-350 |
BET SSA of Calcite (m2/g) | 0.115 |
BET SSA of Quartz (m2/g) | 0.41 |
AT Calcite (m2) | 1.68 |
$\phi \text { ave }$ | 0.40 |
$k_{e f f}\left(x 10^{-13} m^{2}\right)$ | 8.20 |
Local longitudinal dispersivity $a_{L}(\mathrm{~cm})$ | 0.20 |
The initial and inlet solution conditions are listed in Table 3. Flushing through experiments were carried out at flow velocities of 0.1, 0.3, 3.6, 7.2 and 18.5 m/d at pH 4.0. Prior to dissolution experiments, each of the columns was flushed with 10-3 mol/L NaCl solution at pH 8.8 at 18.0 m/d to flush out pre-dissolved Ca(II) for a relatively similar starting point.
Species | Initial Concentrations (mol/L, except pH) | Inlet Concentrations (mol/L, except pH) |
---|---|---|
pH | 8.8 | 4.0 for all columns and 6.7 only for Kratio,Ca/Qtz at 0.5 (in dissolution experiment) 8.8 (in tracer experiment) |
Total Inorganic Carbon (TIC) | 1.0x10-3 (Approximate, close to equilibrium with calcite) | 1.0x10-10 to 1.0x10-5 (depending on experimental conditions, some contain CO2 bubbles) |
Ca(II) | Varies between 1.0x10-5 to 1.5x10-4 depending on experimental conditions | 1.0x10-20 |
Na(I) | 1.0x10-3 | 1.1x10-3 |
Cl(-I) | 1.0x10-3 | 1.0x10-3 |
Br(-I) | 1.0x10-20 | 1.2x10-4 |
Number | Kinetic reactions | Log Keq | kCa (mol/m2/s) |
---|---|---|---|
(1) | $\mathrm{CaCO}_{3}(s)+\mathrm{H}^{+} \Leftrightarrow \mathrm{Ca}^{2+}+\mathrm{HCO}_{3}^{-}$ | 1.85 | 1.0x10-2 |
(2) | $\mathrm{CaCO}_{3}(s)+\mathrm{H}_{2} \mathrm{CO}_{3}^{0} \Leftrightarrow \mathrm{Ca}^{2+}+2 \mathrm{HCO}_{3}^{-}$ | -- | 5.60x10-6 |
(3) | $\mathrm{CaCO}_{3}(s) \Leftrightarrow \mathrm{Ca}^{2+}+\mathrm{CO}_{3}^{2-}$ | -- | 7.24x10-9 |
Reaction (1) dominates when pH<6.0; reaction (3) dominates at higher pH conditions: reaction (2) are important under CO2- rich conditions.
Please do the following:
2.1) Calculate the pore volume (total volume of pore space), residence time $\tau_{a}$, and Peclet number (Pe) at each flow velocity;
2.2) Set up simulation for each flow velocity, and plot the breakthrough curves (BTCs, concentrations as a function of $\tau_{a}$) under each flow condition for Ca(II), IAP/Keq, and pH; plot one figure for each (Br, Ca(II), IAP/Keq, and pH) so you have all BTCs under different flow conditions in the same figure; comment on the role of flow in determining calcite dissolution and why; Here breakthrough curve is defined as concentrations vs time at the last grid block of the column.
2.3) Calculate the steady-state column-scale rates under each flow velocity by R = Q(Ceffluent – Cinfluent). Note that Q (volume/time) = u * Ac, where u is the Darcy flow velocity, Ac is the cross-sectional area of the column; Also note that "steady-state conditions" means the conditions under which concentrations in each grid do not change with time any more.
2.4) Calculate the DaI and DaII (again, under steady state) for each flow velocity;
2.5) Make a table of v, R, $\tau$, Pe, and Da numbers at each flow velocity. Plot R as a function of Pe and DaI and DaII.
In this lesson, we introduced the basic principles of transport and mineral reaction coupling in the context of chemical weathering process and how soils are formed.
Brantley, S. L., and A. F. White (2009), Approaches to Modeling Weathered Regolith, Thermodynamics and Kinetics of Water-Rock Interaction, 70, 435-484.
Kump, L. R., Brantley, S.L. and Arthur, M. A. (2000), Chemical weathering, atmospheric CO2, and climate, Annual Review of Earth and Planet Science, 28, 611 - 667.
Steefel, C.I., Appelo, C.A.J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., Lagneau, V., Lichtner, P.C., Mayer, K.U., Meeussen, J.C.L., Molins, S., Moulton, D., Shao, H., Šimůnek, J., Spycher, N., Yabusaki, S.B. and Yeh, G.T. (2015) Reactive transport codes for subsurface environmental simulation. Computational Geosciences 19, 445-478.
Wolery, T. J. (1992), EQ3/6: A software package for geochemical modeling of aqueous systems: Package overview and installation guide (version 7.0), Lawrence Livermore National Laboratory Livermore, CA.
In this lesson we will set up 2D reactive transport systems in both homogeneous and heterogeneous domains. You will learn 1) how physical and geochemical processes are coupled, and 2) how spatial heterogeneity controls reaction processes. In essence, this lesson combines 2D transport and mineral dissolution and precipitation lesson. Please review these two lessons if necessary.
By the end of this lesson, you should be able to:
In this lesson we are not going to introduce new concepts. The flow and transport processes in homogeneous and heterogeneous systems (1D and 2D) have been introduced before. Similarly, mineral dissolution has also been discussed in a previous lesson. These processes are typically coupled in heterogeneous porous media, meaning they can affect each other. That is, flow influences reactions, and reactions can also have impacts on flow.
In a system where both physical flow and transport and reactions coexist, a general governing equation for any dimension system is as follows:
where Φ is porosity, Ci is the concentration of a primary aqueous species i (mol/m3), D is the hydrodynamic dispersion tensor (m2/s), u is the flow velocity (m/s), ri,tot is the summation of rates of multiple reactions that the species i is involved (mol/ m3/s) , and Np is the total number of primary species. Solving Np equations simultaneously gives the concentrations of primary species at different time and locations, which can then be used to solve for the concentrations of secondary species through laws of mass action (equilibrium relationship of fast reactions).
The equation (1) is general for systems with different number of dimensions (one, two, or three dimensions). Note that this is similar to the governing equation that we discuss in the 1D chemical weathering lesson, except that here we do not specify the dimension. If we specify this for 1D system, it will then be the same from as the governing equation in the 1D chemical weathering lesson.
Also note that in this equation, the dispersion coefficient and the velocities are written as tensors or matrix (in bold symbols) and can have different values in different locations. In fact, all porous medium properties, including porosity, permeability, and surface area can have spatial variations, as we will see in the example.
Here we will build upon example 8.1 in lesson 8 to illustrate a magnesite dissolution example in homogeneous and heterogeneous 2D systems based on the experiments in Li et al. (2014). I am reiterating the example 8.1 here.
The following gives an example of setting up flow and transport calculations in 2D homogeneous and heterogeneous domains. Here we use the physical set up of 2D cross-sections of the Mixed and One-zone column in (Li, Salehikhoo et al. 2014). The authors packed 4 columns of 2.5 cm in diameter and 10.0 cm in length with relatively similar amount of magnesite and quartz distributed in different spatial patterns. Here we focus on two columns that represent two extreme cases: the Mixed column with uniform distributed magnesite and quartz, as well as the One-zone column with magnesite all clustered in one cylindrical zone of diameter of 1.0 centimeter. To keep it simple, here we will do the calculation for 2D cross-sections, instead of following the steps in section 3.2.4 in the paper to convert 2D to 3D. So our numbers might be different from the paper. We are also assuming that in the middle one zone magnesite is 100% of the solid phase and the Mixed column has the same total amount of magnesite as in the One-zone case. The 2D system has a size of 25 mm by 100 mm. A constant differential pressure is imposed at the boundaries in the z (vertical) direction, leading the primary water flow direction in the z direction from the bottom to the top. No flux boundaries are specified in the X direction (Li, Salehikhoo et al. 2014).
Columns | Mg zone | Qtz zone | 2aL (cm) |
4aT (cm) |
ke (× 10-13m2) |
Φ avg |
---|---|---|---|---|---|---|
Mixed | - | - | 0.05 | 0.005 | 8.26 | 0.44 |
One-zone | Width: 1.0 cm ΦMg: 0.54 |
Width: 1.5 cm ΦQtz: 0.38 |
0.07 | 0.004 | 10.74 | 0.44 |
*The permeability of the pure sand columns of the same grain size was measured to be 8.7×10-13 m2.
We have shown in the example and video how we calculate the flow field, solute transport, and breakthrough of these two columns. In this lesson we will add the reaction component of this experiment, which is magnesite dissolution. The reaction network and the inlet and boundary conditions are listed in Table 1 and Table 2, respectively. The experiment injected acidic water (inlet) into the two columns, dissolving magnesite. As you probably notice here, the only kinetically-controlled reactions are magnesite dissolution.
Log Keq | k (mol/m2/s) | SSA (m2/g) BET, measured |
|
---|---|---|---|
Aqueous speciation (at equilibrium) | |||
-14.00 | - | - | |
-6.35 | - | - | |
-10.33 | - | - | |
-1.04 | - | - | |
-2.98 | - | - | |
Kinetic reactions (logK value is logKsp value) | |||
- | 6.20×10-5 | 1.87 | |
- | 5.25×10-6 | 1.87 | |
-7.83 | 1.00×10-10 | 1.87 |
Species | Initial Concentrations (mol/L, except pH) |
Inlet Concentrations (mol//L, except pH) |
---|---|---|
pH | 8.8 | 4.0 |
Total Inorganic Carbon (TIC) | 3.43E-3 (Approximate, close to equilibrium with magnesite) | 1.07E-5 (in equilibrium with CO2 gas) |
Mg(II) | Varies between 0.52E-5 to 1.20E-5, depending on experimental conditions | 0.0 |
Na(I) | 1.00E-3 | 1.00E-3 (in dissolution experiment) 1.12E-3 (in tracer experiments) |
Cl(-I) | 1.00E-3 | 1.00E-3 |
Br(-I) | 0.0 | 0.0 (in dissolution experiments) 1.20E-4 (in tracer experiments) |
Si(VI) | 1.00E-5 | 1.00E-4 |
* The measured quantities include pH, total aqueous Mg(II), Na(I), Cl(-I), Br(-I), and Si(VI). Concentrations of all individual species were calculated using the speciation calculation in CrunchFlow based on thermodynamic data in Eq3/6 (Wolery et al., 1990)
In this example, the species from the inlet water include H+, OH-, CO2(aq), HCO3-, CO32-, Na+, Cl-, and Mg2+. Here the major anion is carbonate species so we assume the major aqueous complexes areMgHCO3+ and MgCO3 (aq), as shown in Table 1. This means that we need to solve for 8 (total number of all species) – 5 (total number of secondary species) = 3 governing equations. Note that the ri,tot here should be the total magnesite dissolution rate through the three reaction pathways listed in Table 1:
Where all rate constants and surface area as shown in Table 1. Please note that in such columns, you can calculate the overall magnesite dissolution rates at the column scale using mass balance of the column:
Here RMgCO3 is the column-scale magnesite dissolution rates, CMg(II),out is the average Mg(II) concentration coming out of the column outlet, and CMg(II),in is the average Mg(II) concentration coming into of the column. Here please note that we use Mg(II) to represent the total concentrations of Mg(II), instead of individual species such as Mg2+ and MgHCO3+. , where CMg(II),i and qi are the concentration and flow rate from the grid block i in the outlet cross-section. With the same inlet solution across the inlet cross-section, you can directly use the inlet concentration for the calculation.
Given the above geochemical conditions above and the physical flow and transport given in lesson 8, please set up the 2D Mixed and One-zone magnesite dissolution columns. Please answer the following questions:
In this lesson, we learned how to couple physical flow and transport with mineral dissolution in 2D systems with both homogeneous and heterogeneous spatial distribution of minerals. We can do this similarly for other types of reactions, including surface complexation, ion exchange, and microbe-mediated reactions. Flow has fundamental impacts on how fast reactions occur. In long time scales where the extent of reaction is sufficient to change porous medium properties, the reactions also affect the flow processes.
Li, L., F. Salehikhoo, S. L. Brantley and P. Heidari (2014). "Spatial zonation limits magnesite dissolution in porous media." Geochimica Et Cosmochimica Acta 126(0): 555-573.
Links
[1] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/ExtraExercisesFiles/CrunchFlow.zip
[2] http://esd.lbl.gov/about/staff/carlsteefel/, https://scholar.google.com/citations?user=Yv-KZG0AAAAJ&hl=en
[3] http://www.csteefel.com/CrunchFlowIntroduction.html
[4] https://creativecommons.org/licenses/by-nc-sa/4.0/
[5] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/CrunchFlowManual.pdf
[6] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/0LessonCrunchFlowOrientation.pdf
[7] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/CrunchTope-32bit.exe
[8] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/libiomp5md.dll
[9] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/CrunchTope-64bit.exe
[10] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/libiomp5md_0.dll
[11] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/CrunchTope
[12] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/CrunchTope-Mac
[13] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/transport1D.in
[14] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/datacom.dbs
[15] https://www.e-education.psu.edu/png550/698
[16] http://webbook.nist.gov/chemistry/
[17] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/images/lesson01/Example1.1.zip
[18] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/images/lesson01/Lesson1TakeHome1.1Solution.pdf
[19] http://en.wikipedia.org/wiki/Water_(molecule)
[20] http://en.wikipedia.org/wiki/Chloride
[21] http://en.wikipedia.org/wiki/Sodium
[22] http://en.wikipedia.org/wiki/Magnesium
[23] http://en.wikipedia.org/wiki/Sulfate
[24] http://en.wikipedia.org/wiki/Calcium
[25] http://en.wikipedia.org/wiki/Potassium
[26] http://en.wikipedia.org/wiki/Total_inorganic_carbon
[27] http://en.wikipedia.org/wiki/Bromide
[28] http://en.wikipedia.org/wiki/Total_boron
[29] http://en.wikipedia.org/wiki/Strontium
[30] http://en.wikipedia.org/wiki/Fluoride
[31] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson1/HW1.zip
[32] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson2/HW2.zip
[33] http://en.wikipedia.org/wiki/Double_layer_(interfacial)
[34] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson3/HW3%20solution.zip
[35] https://www.e-education.psu.edu/png550/517
[36] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson4/HW4%20Solution.zip
[37] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/images/Lesson_5/example5.1CrunchFlowNotes.pdf
[38] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/images/Lesson_5/Example%205.1Solution.pdf
[39] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson6/example6.1.zip
[40] https://www.e-education.psu.edu/png550/827
[41] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson8/Hw8%20Solution.zip
[42] https://www.webapps.psu.edu
[43] http://lynda.psu.edu
[44] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson9/Ex9.1%20Solutions.zip
[45] https://www.e-education.psu.edu/png550/sites/www.e-education.psu.edu.png550/files/FileUploads/lesson9/Hw9%20solution.zip