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). 

Figure 0.1. A schematic representation of the oxidation–reduction zones that may develop in an aquifer downstream from an organic-rich landfill. A zone of methanogenesis develops and is followed down gradient progressively by zones of sulfate reduction, dissimilatory iron reduction, denitrification, and aerobic respiration that develop as the plume becomes oxidized through the influx of oxygenated water. Within the iron reduction zone, a pore scale image shows that the influx of dissolved organics provides electrons for iron reduction mediated by a biofilm, with a suite of other reaction products including Fe²⁺, HCO₃-, and OH^-, as well as solid phase precipitates, such as siderite or calcite, which end up reducing the porosity and permeability of the material. Sorption of Fe²⁺ may also occur on clays, displacing other cations originally present on the mineral surface (from Steefel et al., 2005 with permission). 

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)

The following is a representative mass conservation equation for a primary species I in the aqueous phase:

Here C_i is the total concentration of species i (mol/m³ pore volume), t is the time (s), n is the number of primary species, D is the combined dispersion–diffusion tensor (m²/s), u (m/s) is the Darcy flow velocity vector and can be decomposed into u_x and u_z in the directions parallel and transverse to the main flow direction. Nr is the total number of kinetic aqueous reactions that involve species i, v_ir is the stoichiometric coefficient of the species i associated with the reaction r, Rr is the rate of aqueous reaction r.

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^*(m²/s). At any particular location (grid block) with flow velocities in longitudinal (L) and transverse (T) directions, their corresponding diffusion / dispersion coefficients D_L (m²/s) and D_T (m²/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:

1. CrunchFlow -  a reactive transport code
2. For windows, you need to have these files run Crunchflow: executable file (.exe), a dll file (.dll), an input file (.in), a database file (.dbs)
3. Input file: contains information the code need to set up the systems
4. Database file: includes information related to reaction stoichiometry, reaction equilibrium constants, and reaction rate parameters

In input file:

1. Comments should start with “!” and start from the beginning of the line, not middle of the line.
2. Do not use Tab. Use space if you need to lay out things neatly. 
3. Output files from early runs will be overridden by later runs. Keep output files from different runs in different folders if you would like to have a record. Please use the folder name that self explains so that you know the content of different runs easily. 
4. Search CrunchFlow manual for the syntax of keywords.
5. The exercises in the exercise folder are complementary to examples that we go through. 
    

Please watch the video (32:46 minutes) Lesson 0 Crunch Flow Orientation that introduces CrunchFlow and its database and input file structure.

Summary

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.

Reminder - Complete all of the Lesson 0 tasks!

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. 

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 (cm³/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. 

1. What are the species?
   
   List of species (5): H⁺, OH^-, H₂CO₃⁰, HCO₃^-, CO₃^2-, (H₂O is typically included implicitly).
    
2. How many algebraic relationships do we have that define the dependence between activities of different 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. 

3. How many primary species do we have and what are they?
   
   Here we have 5 species in total and 3 dependencies (equations). So we should have 5-3 = 2 primary species. The primary species should be defined so that all other secondary species can be written in terms of primary species.

• Can we choose H⁺ and OH^- as primary species? No, because we cannot write the carbonate species as combinations of H⁺ and OH^-,
• Can we choose H⁺ and CO₃^2-? Yes. See the following table. The species in the top row are primary species. The first left column includes all species.

  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 CO₃^2-:

$$\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.

4. Is the choice of primary species unique? No. As long as the secondary species can be written as combinations of primary species, the list is legitimate. For example, we can also use OH^- and $\mathrm{HCO}_{3}$ as primary species in this example.

We can then write all species in terms of primary species:

Similarly, (H⁺, HCO₃^-), (H⁺, H₂CO₃), (OH^-, CO₃^2-) are also legitimate choices for primary species. However (H₂CO₃, CO₃^2-), (HCO₃^-, CO₃^2-), are not. You can practice using these species to write the expression of secondary species.

In this subsection we will discuss setting up aqueous complexation calculations in CrunchFlow.

Example 1.1 RTM setup

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

Take home practice 1.1 RTM Set up

Example 1.1 Files [17]

Practice 1.1 Solution [18]

Question 1

Metal complexation in seawater. The seawater composition is as follows:

1. Please set up CrunchFlow to do the aqueous complexation calculation for the seawater system and identify the dominant species (the top 2 species with the highest concentrations) for the major metals (Na, Ca, Mg, Sr). Is charge balanced? If not, what species should you put to balance charge?
2. Try the keyword “database_sweep” (look it up in the manual). What are the dominant secondary species formed in seawater? Pick the dominant secondary species for each cation. What is the difference by including and not including the dominant secondary species?
3. A Pb-containing solution is accidentally added to a tank of seawater, resulting in a total concentration of all Pb-containing species being 10^-2 mol/kg. All other concentrations remain roughly the same. In calculation, please make sure charge is balanced.
   • Redo the calculation and identify the dominant species for all major metals, including Pb.
   • Note that the total concentrations of Ca and Pb are similar. Do they differ in their dominant species? If there are differences, which parameter leads to such differences? 

Question 2 open carbonate system

Here we assume we have an open carbonate system instead of a closed system. This means H₂CO₃, or CO₂(aq), is in equilibrium with atmospheric CO₂ concentration at P_CO2 of 10^-3.5 atm. Henry’s law constant for CO₂ dissolution is K_H = 10^-1.47 mol/L/atm. Suppose we have a solution with a given pH and H₂CO₃ 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).

1. At pH=7, what are the concentrations of individual species (H₂CO₃, HCO₃^-, CO₃^2-)? 
2. Calculate the concentrations of TIC, H₂CO₃, HCO₃^-, CO₃^2-, H⁺, and OH^-, under pH of 2, 4, 6, 8, and 10.
3. Plot the concentrations of all individual species as a function of pH.
4. Observing from the plot, what are the dominant species under different pH conditions.
5. How does this open system behave differently from the closed system in example 1.1?

Question 3 buffering effect of carbonate system

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]

Reference:

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.

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 (CaCO₃) dissolves into carbonate species (H₂CO₃, HCO₃^-, CO₃^2-) 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 CO₂ leakage from CO₂ 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 CO₂ consumption during chemical weathering (through mineral dissolution) helps maintain the clement conditions for life on earth. 

Calcite Dissolution Kinetics

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 CO₂- 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 (m²) k₁, k₂, and k₃ are reaction rate constants (mol/m²/s) for the three parallel reactions, respectively; $a_{H^{+}}$ and $a_{H_{2}C{O}_3^0}$are the activities of hydrogen ion and carbonic acid: IAP and K_eq 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 H₂CO₃⁰, and how far away the reactions are from equilibrium (IAP/K_eq).

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/m²/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. 

Figure 1. Calcite dissolution rates measured at 298 K and under various pH and CO₂ partial pressure conditions (Brantley, 2008). At low pH, rates are transport controlled and are not dependent upon CO₂ partial pressures, and the logarithm of rates linearly depend on pH. At pH>3.5, dissolution rates are nearly constant at relatively high CO2 partial pressure (0-10³ atm) whereas decrease with pH at low CO₂ partial pressure (e.g., 10^-6 atm).

Generalized Transition State Theory (TST) Rate Law

A general TST rate law is as follows:

Here nk is the total number of parallel reactions, k_j is the rate constant of the parallel reaction $j(mol/(m^2·s))$, the term $a_{H^{+}}^{n_{H+}}$describes the rate dependence on pH, and the term $a_i{}^{n_i}$ describes the rate dependence on other aqueous species that potentially accelerate or limit reactions. In Equation (4) for calcite dissolution rates, $a_{H_{2}C{O}_3^0}$ also accelerates calcite dissolution. The affinity term $IA{P}_j/K_{eq,j}$ quantifies the distance of solution from the equilibrium state of the mineral reaction j. The exponents m_1,j and m_2,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.

Figure 2. Dissolution rates of silicates as a function of pH (25 deg C).

This is Figure 8.16 from Appelo and Pstma, 2005, use with permission

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. 

Rate Constant Dependence on Temperature

The Arrhenius equation is used to quantify the reaction rate dependence on temperature.

Here A_F is the pre-exponential factor and is usually considered as a constant that is independent of temperature: the activation energy $E_a\left(kcal/mol\right)$ is always positive; here R is the ideal gas constant $(1.987cal/(mol·K))$, T is the absolute temperature (K). If we take the logarithm of this equation, we obtain:

This means if we draw lnk versus $\text{1/T}$, we get a straight line, as shown in Figure 2. For a particular reaction, if ${\text{E}}_{\text{a}}$ and ${\text{k}}_{\text{1}}$ at the temperature T₁ is known, ${\text{k}}_{\text{2}}$ at another temperature T₂ can be calculated by

Figure 3 shows the dependence of calcite rate constants on temperature under different CO₂ partial pressure.  

Figure 3. The effect of temperature (298-373 K) on calcite dissolution rates at constant stirring speed of 425 rpm, pH=4.0 in the solution of 0.1 M NaCl for 2, 10, 30 and 50 atm pCO2. The calcite dissolution rates increase with the increasing temperature.

(Pokrovsky et al., 2009, used with permission).

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. 

Example 2.1: Calcite Dissolution

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 $3.15\times {10}^{-4}mol/kgw$. The kinetic rate constant for three parallel reactions are $1.00\times {10}^{-2},5.62\times {10}^{-6}$, and $7.24\times {10}^{-9}mol/m^2/s$ respectively. The volume fraction of the calcite grains in the batch reactor is 0.5%, with a specific surface area (SSA) of 0.24 m²/g. Please set up the simulation and plot the following quantities as a function of time:

1. Total Ca(II) and total inorganic carbon;
2. Individual Ca-containing species and carbonate species;
3. pH and the saturation index $IAP/K_{eq}$

Here is a light board video (12:44) of what equations the code solves for this system.

Click for a transcript of the mineral dissolution and precipitation video.

Mineral Dissolution & Precipitation

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.

Credit: Li Li @ Penn State University is licensed under CC BY-NC-SA 4.0 [4]

Here is a video demo in CrunchFlow on how to set up the simulation in CunchFlow (43:51 minutes).

Click here for a transcript of the Lesson 2 Ex CrunchFlow video.

Li Li: OK, so let's start. This is lesson 2 on mineral dissolution and precipitation. Now, with example 1, it's to set up as if you are doing a batch reactor experiment with some calcite grains in the solution. And you essentially will be running numerical experiments, imagining that your calcite dissolution rates are in these given numbers-- you have how much calcite grain in the system, et cetera.

So here-- let's see. For this lesson, the difference between this one and previous one in lesson 1 is, lesson 1 only have aqueous complexation. So it doesn't have any kinetics, it only has thermodynamics. The other one is zero-space dimensioning, and also, no-time dimensioning.

Now here, by introducing the reaction kinetics, we are essentially now looking at it as a time dimension. We'll still have a well-mixed system, which is a zero-dimensional system. We're not looking at things change in the functional space, so you have uniform contusion in the batch reactor.

So that's a major difference. Let's read the example question, and then I would explain a little bit, and we can try to set it up in CrunchFlow. OK, so example 1 is about a calcite dissolution. Calcite dissolution proceeds through parallel pathways, as shown in reaction 1 to 3. That's what we have-- we're seeing lesson 2.

In addition to these kinetic reactions, there are also a series of instantaneous aqueous complexation reactions, that occur at the same time. The three fast reactions should be the same as the ones in the example 2 of lesson 1. Also, that's for these fast reactions on the aqueous speciation reaction. So this essentially specify what aqueous speciation reaction you should include in this example.

Now, in a batch reactor, the initial solution is at a pH 5.0, with close to 0 salinity, meaning it is not a constituted solution. The initial total carbon concentration is 3.15 times 10 to the minus 4 mole per kilogram water. Kinetic rate constant for three parallel reactions are 1.00 times to 10 to the minus 2, 5.62 times 10 to the minus 6. And then the other number-- this is for the corresponding three reactions in 1, 2, 3.

And then the volume fraction of the calcite grains in the batch reactor is 0.5%. So if you have 100 milliliter, then you have 0.5 milliliter of calcite. You can convert that to how much mass is this with a specific surface area of 0.25 meters squared per gram. Please set up the simulation and plot the following quantities as a function of time.

So now we're looking at the time dimension, right? We're starting from the initial condition. So initially, you have some total carbon concentration there-- 3.15 times 10 to the minus 4 mole per kilogram. But essentially there, you don't have other species.

Now, the question it's asking you, first of all-- so total calcium concentration and total inorganic carbon. So you need to plot that as a function of time. Second was the individual calcium-containing species. And carbon is species. And so here, you want to see, essentially asking, which one the dominant species and everything. And then also, the pH and IAP over Keq, which is a saturation index. Maybe that-- make it more explicit.

OK, so let's do this. Let's go to that folder again. OK, again-- so I have four files here waiting for you. Let's open the input file. Oh, I hate it when-- let's do it in Notepad. I don't like all these red line and everything-- underlining all my words. Let's open it with Notepad. OK. Much better.

All right, so this is a template, essentially. And we need to fill in, OK? So this did not say TITLE will be "lesson 2-- calcite dissolution in a batch reactor." These runtime specifications are still there. Database-- make sure it's datacom.dbs and everything. OK.

So here, make sure you don't do speciate only, because when you do speciate only, the time is not going to be stamped, and you cannot do this calculation. Because this is for kinetic, and you want your time stamping Output file-- the output procure block. We need to specify what we want the system to output.

So first of all, we need to specify, what units are we using? Let's do seconds. We know calcite dissolution very fast-- or maybe minutes. Let's do minutes. And let's call it time. We need to specify time series. And let's call that calcite dissolution .out.

And for batch reactor, we should only have one group block. So let's just specify this as a group block 1. Time series print. So what time series print do is, you specify what species you want to print. So you can put a specific freedom-- calcium.

And this need to be the-- I'm just going to make it explain. Needs to be primary species. It cannot output times there for secondary species. So we can specify it's the calcium pH or bicarbonate. But it's very important that maybe we'll do all-- "all" meaning all primary species would be printed out.

Discretization-- we should only have one species, right? So first, we need to specify the units' distance-- units. Let's say it's meters. And then, when you need to specify kind of zones-- x-zones, y-zones. So you have x-zones should be 1.1. So you only have one group block.

Now, for the other dimensions-- if you only specify 1, so that essentially you have 1 meter, 1 meter, 1 meter. If you don't specify other, it will assume it's meter-- the same units-- and it's 1. Quantity is 1. So essentially, you have 1 meter cubed. Maybe we'll make it centimeters, just to be more realistic.

And you have 10.0. It shouldn't matter, because we are specifying volume fractioning, not a specific number. And so if I'm specifying-- maybe I will do 100. Because a typical batch or a flux or something, that you would do batch reacting as hundreds of milliliters. So I'm specifying here-- this is 200 centimeter cubed, essentially. The other two dimension-- y-direction and z-direction into it-- by default, it will be 1 centimeter, 1 centimeter. So it should be 200 centimeter cubed-- essentially 200 milliliter. That's the discretization.

Now, we actually also need-- in output, we need to specify how long we want this to run-- spatial profile. So maybe we'll do, let's say, 1 minute. And then, maybe 20 minutes? 10 minutes. OK. Maybe 10 minutes. It is out pretty fast. Let's see. It's 5.0, 10.0, 50.0. All right. So we run it for 100 minutes. So it's the end point time. So this submission would be run until 100 minutes.

Now, we don't have spatial discretization. We don't have. We only have one group block. It's a well-mixed system, so we don't need to specify-- we don't have something coming from the boundary. We should be fine. So let's specify-- initial conditions should be whatever condition we specified it. So let's do it later.

We don't have transport. We can safely delete that. We also don't have flow, we don't have porosity. We can now delete all of that. So the primary species you want to put in, we talked about before. You should have calcium. First of all, all we need-- H plus, actually. What else you need? I'm going to look at the question again.

Let's put a bit of sodium chloride there, just in case we need-- to make charge balanced, we can use sodium chloride. So usually, you use species that are not reactive to do charge balance. Otherwise, you would be changing the condition we put in the species.

What else? And we need bicarbonate speciates-- still use bicarbonate. That should be it for primary species. Second species-- again, OH minus would be the top one. Or the other carbonate species-- A, Q-- 0 3. Potentially with these, we could form COH plus-- that's another aqueous species.

You can check by doing speciation only-- by data sweeping or something-- to see which species will be dominating, if you would like. But that's your calcium. Carbon we also form-- aqueous calcium carbonate. And if you would have carbide, possibly you would also have that. OK, I think that's good for our second speciate, I believe.

OK. So for this one, we need initial condition. Because the system started with-- this is where t equal to 0. So let's call that condition initial. And you will need all the prime-- you need to specify all the primary species, so copy the list of primary species here.

So question says, initial pH is 5.0. So let's specify pH would be 5.0 first, instead of saying H plus. Let's do 5.0. So initial calcium concentration should be almost 0. Now usually, I try not to put directly zeroes there, because sometimes computers have problem with 0. And I always put a very small number there, just in case. So for the very small-- essentially none-- I'm putting very-- these 10 to minus 10 concentrations. All right.

Bicarbonate is a good-- the question said, initially, it's 3.15 times 10 to minus 4. OK, they say what you put in there. That's the total concentration saturate, not specific. And question says, total carbonate concentration. That's consistent.

Now here, we also have either-- in solution, we should also have calcite, because calcite need to dissolve out. So we need to have the calcite mineral. This is very important, because otherwise, there would be no representation calcite going in there. So let's put calcite. And it says, volume fraction is 0.05%.

So the first number should be the volume fractioning. And then, it also says, specific surface-- the key word is specific-- surface area is 0.25. So let me just comment on that. The first number, volume fractioning, is 0.05. This SSA is 0.24 meter squared per gram. So that's the unit.

You can also specify total surface area, I believe. Let me just search the manual. I believe you can do-- this should be similar. Oh, you can specify BICG surface area or specific surface area. If it's BICG surface, you'll be putting the absolute value. If it's specific surface area, it should be in the units meter squared per gram.

So it's the typical-- the way you're putting is, name of solid face, volume fractioning, surface area option, and value. BICG surface area. So you can look through these if you would like to have more details. Go into page 68, 69 in the manual. All right.

Now, we cannot put in the grain size. But the specific surface is more or less refracting the grain size you have. When you have small grains, you tend to have larger specific surface area. And when you have larger grain, you have small specific surface area. But also, of course, for different mineral, they're also different. OK. So the initial condition.

But one thing we're keeping-- we need to keep in mind is that the solid face need to be recognized by the mineral-- key word, block. This is where you're putting the kinetic information. Without that, you will either use the kinetic information grid block, or in the database-- which I usually prefer to putting input file, because then you have everything in the same place.

So here, you will need to put in the rate constants in the log 10 units. These rates should be in log 10 rates in units of mole per meter squared per second. All right? So the way you're putting is calcite, and then you have the label, if it splits it using default.

Then your rate is-- OK. Let me just check the three reaction pathways. We are saying, for calcite, you have these three reactions. So you have k1, k2, and k3. They all have the same-- so one is, k1 depends on each process. k2 depends on H2O CO3. k3-- there's no dependence. So let's just also check the database.

Let's save it, and let's look at database. Make sure it's a default in everything. Make sure we know what are the format for each label. Search for calcite. OK. So this said, T-S-T rates are-- default is-- OK. H plus 1 is one that depends on each process. That's the k1, right? So let's see, label. So you need to look back and forth a bit. So let's do H plus first.

And I said, the rate constant is what? 1 times 10 to minus 2, for the 3 or 4 k one. OK, for the H plus dependence once. So this one is-- in units of mole per meter squared per second. So this one should be minus 2. And then, you have the-- so second one is-- let's call it a CO2, so it's less confusing. Let's call it a CO2 dependence. And it's dependent on this with 1.0.

So this dependence need to specified in the database, but not the one-- so I probably should say H2O3. And then, the default is 1 without any dependence. OK, so if we specify the number here-- OK, so that will be CO2. And then rate.

So the second rate constant is 5.62 times 10 to the minus 6 in log unit. Let me just-- 5.62. OK, 10 to the minus 6-- minus 6. So it's 5.25. All right. And then, the default ones without dependence is 7.24 times 10 to the minus-- so 7 point. Minus 9-- minus 8.14. For the three parallel reactions-- 1, 2, and 3. OK.

So now, what-- so these rate constants should be override to whatever that is in the database. OK, you have condition-- let's just check it. You have a condition-- you have mineral specified with kinetics. You have primary and second species. OK, you need to specify what are the initial condition.

So you should have say, initial condition-- initial, we call it-- yeah. We called that initial. So this condition will be in the first code block. We only have one group block. This is specify the location and-- you can dig into manual for them-- search for this key word block, and the manual explain in detail what you can put in.

So here-- this initial condition is interesting. OK, I'm using this initial condition for this group block 1. Let's see if they run. I think you should. And so I made a mistake somewhere. So that's 2 lesson, example 1. Oh-- rate label. Cannot have enough on-- in database. Looking for calcite rate label CO2. Let's see. Hm. Interesting. Maybe I didn't save it.

OK, let's do it again and just save it 2 lesson example. Now it runs. There's something else coming up. So the other time I changed it, the to CO2, but I didn't save it. So the code doesn't recognize the CO2 label. So I save it to the database with the CO2. Then the code not recognize that label.

And then, now is speciation-- OK. So it-- speciation is finished. No aqueous kinetic block found. Retardation parameters not found. Both the number of cells and the grid block spacing should be-- oh, OK. Only one value provided. Yeah, I've seen the problem with discretization. So I need to provide-- let's go to-- just make sure we-- discretization key word block-- x-zone, y-zone, z-zone.

OK. At discretization, you can specify x-zone, z-zone, y-zone. So the syntax is-- it's supposed to say x-zone. And it's how many number of cells and spacing. OK, good. And number of cells should be integer, and spacing is a real number. So we are saying, we have one group block. And the spacing is-- it's 200 centimeter. I think that should make it work. Let's try it again. So I just look it up and do it again.

All right. So the runs is finished. That's good. OK. And suddenly, you discover there are a lot of output files. Each of these output file-- for each of these-- these are all out, right? So you have each of these every time. You have every 1, every 2, every 3, 4, 5. Each of these is corresponding to one of these times that you specified in the spatial profile. Remember, in spatial profile, I put 1 minute.

So it's spit out the output at 1 minutes, 5 minutes, 10 minutes, 50 minutes, and 100 minutes. So the 1 will be representing 1 minutes. 2 will be 5 in this. So each of these have 5. Right? I said 2, perhaps there is 5. Concentration have 5. Guess mineral percentage, because you have mineral. pH, porosity-- a lot of information.

OK, these are standard operations. Some of them are not useful for this example. But I want to-- let's go through. So you don't have to answer questions at these. You want to print-- plot everything in the function of time. All these other areas-- like, area 1 to 5-- because we only have one grid block. So the spatial profile is not really useful for us.

If we have 100 grid block, we want to see how this varies spatial rate. Now you only have one grid block. So it's not really interesting. It's more interesting to look at that as a function of time. And then we specified the breakthrough curve should be in-- OK. Calcite dissolution data out. So let's look at that file. This will be what you should look for for the breakthrough.

Calcite. This-- OK, this good. So as we said, essentially output everything, right? OK. So again, first it's time stepping varied with very tiny time steps. And then, later, it becomes more larger times then, because it's approaching to equilibrium.

Now look, and so you have pH, you have H plus, calcium, sodium chloride, always a function of time. We said "all," so it also spit out sodium chloride. But we actually don't really need that. But why-- OK, we also should have-- OK. Hm, interesting. Bicarbonate, H plus, CO2, AQ, 03, calcium. So it actually does spit out all the aqueous speciation, as well.

And these are inactive variables. These are in, actually, log units. For the primary species, it's an absolute value. It's not log. For the second species there, you log in log. So this, it would be 10 to the minus 8.99, 10 to the minus 3.52, 10 to minus 9.5, for example. So they should be increasing.

So you actually can directly read, if you want, to plotting. So you can either do Matlab, Excel, whichever ways are the easiest for you. Let me just show you how to open that in Excel, so if that's easiest for you. I'm opening up Excel. And I will say, Open. Hm. Why don't you-- oh. It's says no book. Sorry.

OK, let's do this one-- let's see how. OK. So you would say, so the code, I can actually fix the widths. It will space out for this. OK, it probably doesn't matter. OK, so you have all the different species-- pH, H plus, pH. So you can plot-- whichever that I asked for you can plot. Problem is species. Pay attention, there's this thing.

OK, so this also give the calcium in the log units. So these are individual species, and these are the total for the primary species. So if you want to plot a primary species, like the individual species, you do want to use these columns, starting from H plus here, when it has a log-- after it's already been log.

So this is how you can do, for example, getting all the time series and plot out the files. So that's an example. You can see all the change in the function of time, how the saturation changed the function of time. So here, I'm giving you the solution already, but you want to be able to generate that by yourself. I want to be able to see that.

And you shouldn't just-- before you do anything, you shouldn't look at this solution yet. But I'm still telling you, here-- these solutions, OK. For the individual concentrations, OK, these are total calcium, total bicarbonate. So it should be from the first several columns on primary species.

And you can see, the carbon is always higher than the calcium, because initially, there is some carbon in species there. So if you look carefully, the difference is actually exactly the initial concentration of carbonate. Because essentially, the calcium-- when the calcite dissolves, it should be dissolving calcium and total carbon in one-to-one stoichmetric coefficient.

When the calcite dissolves, pH increases over time, until it reach around 8, so it's concentration in which certain values stopped, because it reach equilibrium. IAP over Keq get to about 1, which is at equilibrium at about 120 or 130 minutes, which is similar to where concentration reached a plateau.

And looking at the individual concentrations-- so this is a bit hard to see. Total calcium should be the highest, and blue is-- so the free calcium is a dominant species. And then this light blue, you see it-- bicarbonate. So calcium bicarbonate is the second one. Red one is Ca CO3, and then you have COH plus, which is reasonable, because you have a lot of carbonate species in the system.

So that should give you a sense of how you would do with the mineral dissolution precipitation. And as an extension of the example 1, in the homework assignment, I asked you to do the carbon dissolution, changing several key parameters, essentially.

Example 1 did all this in one particular condition. But we all know, in the lessons, we discussed where the constant is important, surface area is important. Initial pH, salinity-- all these things are important. So how do they change the dissolution process, if we change these numbers? So in every question of these, it asks you to change your numbers, and compare three different cases, in term of how things are different with different sets of parameters.

And then, the question 1 is about feldspar dissolution. So in example 1, give you everything about calcite. In question 2, you are supposed to calculate to set a feldspar dissolution using the figure 7 in Blum and Stillings, with the rate dependence on H plus, OH minus, or the neutral one without any pH dependence. And you should be able to study that up to run these simulations.

All right, I think I'm going to stop here. And I hope you have fun playing with CrunchFlow in setting up mineral dissolution precipitation. Again, this is in batch reactor. So imagine that you are doing a numerical experiment of the batch reactor, and you're putting calcite grains-- different grain size.

So you have larger, smaller specific surface area, and all that. OK, I'm going to stop by, and hope you have fun here. And we can discuss later, if you have questions. OK.

Click here for the solution.

Figure 4 below shows the evolution of aqueous chemistry during calcite dissolution. Total Ca(II) and carbonate increases over time and should follow the 1:1 stoichiometric ratio. However, total carbonate concentration is higher due to the initial presence of total carbonate at the concentration of $3.15\times {10}^{-4}mol/kgw$. The individual speciation of Ca(II)-containing species are in D shows that Ca²⁺ is the dominant species. During the dissolution, pH and $IAP/Keq$ increase over time, until $IAP/Keq$ reaching close to 1 at about 100 minutes, after which the system stabilizes and no more dissolution occurred after that.

Figure 4. Predicted concentration evolution in the batch reactor. (A) total Ca(II) concentration and total CO2(aq) concentration, (B) logarithmic concentrations of Ca-containing species, (C) logarithmic concentrations of carbonate species, (D) pH, and (E) IAP/Keq

As can be seen from these figures, the dissolution reactions reach equilibrium at about 200 minutes, where pH increases to about 8. Total Ca(II) include Ca²⁺, CaOH⁺, CaCo₃(aq), and CaHCO₃⁺, with Ca2⁺ being the dominant species (shown in figure 4B). Total inorganic carbon includes CO₂ (aq), HCO₃^-, CO₃^--, CaHCO₃^-, and CaCO₃(aq). HCO_3^- is the dominant species at most of the time except very early time when pH is still at about 5, where CO2(aq) dominates. 

Credit: Li Li @ Penn State University is licensed under CC BY-NC-SA 4.0 [4]

1. Carbonate Dissolution (Example 1 extension) (Total 70 points, each sub question is 10 points):

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 m²/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 m²/g. Do changing k and A have the same impact on dissolution rates?

3) Equilibrium constant K_eq 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 K_eq 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⁺, CaCO₃(aq), and CaHCO₃⁺. 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?

2. Feldspar Dissolution (optional, bonus points: 30, for students who develop clear solution and reasoning):

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.

Click Here for HW2 solution package [32]. 

Summary

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.

Reminder - Complete all of the Lesson 2 tasks!

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.

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 (Cd²⁺, Cr³⁺, Co²⁺, Cu²⁺, Fe³⁺, Pb²⁺ or Zn²⁺) 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.

Figure 1. Schematic illustration of the formation of inner-sphere and outer-sphere complexes at the solid–solution interface

Credit: [Goldberg et al., 2007]. Used with permission.

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.

Surface Configuration of the Solid–Solution Interface

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.

Figure 2. A schematic diagram of solid surface, double layer, and bulk liquid with a negatively charged solid surface.

http://en.wikipedia.org/wiki/Double_layer_(interfacial) [33]

Surface Complexation Reactions

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 K^app 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 L^l- represents an anionic ligand of charge l⁻.

Apparent (K^app) and Intrinsic (K^intr) Equilibrium Constants

The apparent equilibrium constant of surface complexation reactions, K^app, is an important parameter because it determines the ion partition between aqueous and solid phases. Large K^app values indicate high affinity of the ions to the solid surface. The relationship between total Gibbs free energy $\Delta G_{\text {tot }}$ and K^app 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.

Surface sorption site

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 C_site is the concentration of surface sites (mol/g mineral), ρ_sites is the surface density of surface hydroxyl sites (mol/m²), A_specific is the specific surface area (SSA)(m²/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/m². 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/m²) with SSA (m²/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.

Different types of surface sites

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 Fe^sOH 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).

Surface charge and point of zero charge (PZC)

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.

Click here for HW3 solution package. [34]

Summary

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.

Other reference books

Reminder - Complete all of the Lesson 3 tasks!

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]

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 Ca²⁺ from solid phase (Figure 1). In this reaction, Ca²⁺ 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 Ba²⁺ and Sr²⁺ (Gregory et al., 2011).

Figure 1. Intrusion of seawater into the coastal aquifer. Divalent Ca is partially displaced by monovalent Na from seawater through ion exchange reactions. The water changes from Na-rich to Ca-rich, while the solid surface changes from Ca-rich to Na-rich.

Credit: Figure modified from St. Johns River Water Management District (2008)

Figure 2. Process of water softening. When hard water with high concentration of Ca and Mg flows through the softener resin tank, Ca and Mg exchange with Na, resulting in Ca and Mg concentration decrease in water (soft water).

Cation exchange reactions

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 A^u+ and B^v+, 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 K_eq 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 10²⁰ positive charges. On the other hand, an meq of Ca²⁺ is also 6.02 x 10²⁰ positive charges, however only 3.01 x 10²⁰ ions because Ca²⁺ has two positive charges. For an ion Iⁱ⁺ with a charge of i, the equivalent fraction b_I 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, ${\displaystyle \sum _{}\beta =1}$ .

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).

Selectivity coefficient

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⁺, Ca²⁺, and Mg²⁺, 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).

For typical ion exchangers, the sequence of affinity is as follows: Ba²⁺> Pb²⁺ > Sr²⁺ > Ca²⁺ > Ni²⁺ > Cd²⁺ > Cu²⁺ > Co²⁺> Zn²⁺ > Mg²⁺ > Ag⁺> Cs⁺ > K⁺ >NH⁴⁺ >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).

Cation exchange capacity

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.

Example 4.1.

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 CO₂ 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.

“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 m²/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)):

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.

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 N₂ 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., CO₂), 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.

Biogeochemical redox ladder

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 C₆H₁₂O₆; sulfate reduction can produce 492 kJ /mol of C₆H₁₂O₆, 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.

Figure 1: Example biogeochemical redox couples for environmentally relevant species. Values were calculated at pH 7 with concentrations of 1 M for all dissolved constituents except for $\mathrm{Fe}^{2}+\left(1 \times 10^{-5} \mathrm{M}\right)$, $U6+\left(1 \times 10^{-6} \mathrm{M}\right)$, and $\mathrm{CO_3}^{2-}\left(3\times10^{-3}\right)$

source: Borch et al., 2010. Used with permission.

*Table modified from Rittman and McCarthy, 2001. Used with Permission

Writing microbe-mediated reactions

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 (f_e) and for cell synthesis (f_s). 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 f_e) and can channel more electron into cell synthesis. Therefore, values of f_e increase and f_s 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 C₅H₇O₂N, 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 C₈H₁₃O₃N₂, 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 (R_d), electron acceptor (R_a) and cell synthesis (R_c). The energy production reaction R_e = R_a - R_d and cell synthesis reaction R_s = R_c - R_d, 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.  

Aerobic oxidation:

$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 = -R_a - R_d that yields:   $\frac{1}{8}CH_3COO^-+\frac{1}{4}O_2=\frac{1}{4}HCO_3^-+\frac{1}{8}H^+$

The anabolic pathway R_S = R_c - R_d 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}$

Denitrification

$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}$

Iron reduction reaction

$\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}$

Sulfate Reduction

$\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 (C₅H₇O₂N) 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).

Dual Monod kinetics

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/m³), C_d and C_a are the concentrations of electron donor and acceptor of the reaction (mol/m³). The K_m,D and K_m,A are the half-saturation coefficients of the electron donor and acceptors (mol/m³), 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.

Dual Monod kinetics with inhibition term

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 K_I,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 K_I,H?C_H.

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 O₂ 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 NO₃ concentration in the bottle are 5.0 and 1.0 mmol/L, respectively. Assume initially there are some dissolved O₂ 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 O₂-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}$

*Range of relevant parameters is from Cheng et al., 2016; Li et al., 2010.

Questions:

1. Which reaction would occur first?
2. What do you expect to see the evolution of the concentrations for the species: acetate, O₂(aq), nitrate, total inorganic carbon, and biomass for denitrifiers. Plot their concentrations as a function of time in the same figure.

Please watch the following video: Microbe Mediated Reactions (24:50) 

How to simulate microbe-mediated reactions in CruchFlow (based on the simulations in previous lessons):

Please watch the following video (29:51):

Notes on setting up in CrunchFlow [37]

Solution to example 5.1 [38]

Problem 5.1

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 (F_s = 0.08 and F_e = 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.

*Range of relevant parameters is from Cheng et al., 2016; Li et al., 2010.

1. What are the types of the reactions and the order of reactions in this system?
2. Plot out the additional sulfate concentration evolution, together with other chemicals in the example;
3. Increase both K_I,NO3 and K_I,O2(aq) values by a factor of 2 for the sulfate reduction reaction, what difference do you see in the time evolution figures.
4. Initial O_2(aq): Simulate one more case with no O_2(aq) concentration. What are the differences in the concentration evolution of different species? And why are the differences?

Problem 5.2

Extension: In Problem 1, we mainly discuss how O^2(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 N_2(aq) and sulfide evolution:

1. How does the growth rate of NRB (maximum grow rate of nitrate reducing biomass _μmax,NO3) change the evolution? The _μmax,NO3 value in the base case is 20,000 mol/mol-biomass/yr. Simulate two more cases with _μmax,NO3 of 40,000 and 60,000 mol/mol-biomass/yr.
2. How does the half saturation of the electron donor for the nitrate-reducing reaction (e.g., K_m,acetate) influence the evolution? The K_m,acetate value in the base case is 1.0 mmol/kgw. Simulate two more cases using K_m,acetate being 0.2 and 5.0 mmol/kgw.
3. How does the abundance of electron donor (acetate) affect the evolution? The initial concentration of acetate C_{initial,acetate} in the base case is 5.0 mmol/L. Simulate two more cases by increasing and decreasing C_{initial,acetate} by an order of magnitude.

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/m³). 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), $\mathrm{\$\backslash rho\$}$ is the fluid density (kg/m³), and g is the acceleration of gravity (9.8 N/m²). The hydraulic conductivity (m/s) can be expressed as

where $\kappa$ is the permeability of the porous media (m²) 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.

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 V_pore is the pore volume (m³), Q is the flow rate (m³/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 (m³/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 J_diff is the diffusive mass flux per unit area (mol/m²/s); D₀ is the molecular diffusion coefficient in water (m²/s); and x is the distance (m).

Diffusion coefficients in water are typically in the order of 10^-9 m²/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 (m²/s) at a reference temperature T₀ (K), and $\mu_{T}$ and $\mu_{T_0}$ are the dynamic viscosities $\left(Pa\cdot s\right)$ at temperatures T and , T₀ respectively.

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 Coefficient D_e in Porous Media

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]. T_L is the tortuosity (dimensionless) defined as the ratio of the path length in solution relative to the tortuous path length in porous media.                

Figure 1. Tortuous diffusion paths in porous media. The grey color represents mineral grains. The white means pore spaces. Chemicals in the aqueous phase go through tortuous path while traveling through porous media.

Steefel and Maher, 2009

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

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 D_L in the principle flow direction is typically higher than D_T 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.

Figure 2. Longitudinal dispersivity versus scale of observation 

Gelhar, 1986

Hydrodynamic dispersion

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 D_h is the hydrodynamic dispersion coefficient defined as the sum of effective diffusion coefficient D_e and mechanical dispersion coefficient D_m:

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 C₀ 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 (m²/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.

Figure 3. Effects of Pe on the shape of concentration profiles at steady state. Large Pe indicates that the dominant advective transport compared to dispersion or diffusion.

Steefel et al., 2005

Please watch the following video: Advection-Dispersion Equation (ADE) for non-reactive tracers (16:42)

Example 6.1: Setting up 1D Flow and Transport Simulation in CrunchFlow

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 m² 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 D₀ for the tracer bromide is 1.8×10^-9 m²/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 PaŸs at 20 ⁰C.

1. Before the numerical experiment:
   1. Calculate the pressure gradient for the Darcy flow velocity.
   2. Calculate the effective diffusion coefficient De, mechanical dispersion coefficient D_m, hydrodynamic dispersion coefficient D_h. Which term dominates in D_h?
   3. Calculate the characteristic time for advection (residence time) and for dispersion, and Pe number. Based on the Pe number, predict which term dominates during the transport process.
2. Set up in CrunchFlow using 100 grid blocks and run the simulation:
   4. Plot the bromide breakthrough curve at the end of the column (C/C₀, with C₀ being the injection concentration) versus residence time. 

Click for Solution

Click here for input and database file package [39]

Solution 1:

1. For the 1-D flow, Darcy flow equation can be simplified to:

   $u=-\frac{\kappa}{\mu} \nabla P$
   
   Thus, $\nabla p=\frac{4.20 \times 10^{-6} \mathrm{~m} / \mathrm{s} \times 1.00 \times 10^{-4} \mathrm{~Pa} \cdot \mathrm{s}}{1.75 \times 10^{-13} \mathrm{~m}^{2}}=2.40 \times 10^{4} \mathrm{~Pa} / \mathrm{m}$
   This means that for a 10 cm (0.1 m) column, the pressure difference between one end to the other end of the column should be 2.40E3 Pa. Note that in Crunchflow, the pressure units is Pascal.
2. $D_{e}=\frac{D_{0}}{F}=\phi^{m} D_{0}=T_{L} D_{0}$
   
   So $D_{e}=0.4^{1.0} \times 1.8 \times 10^{-9} \mathrm{~m}^{2} / \mathrm{s}=7.2 \times 10^{-10} \mathrm{~m}^{2} / \mathrm{s}$
   
   $\mathbf{D}_{m}=\alpha \cdot \mathbf{u}$
   
   So $D_{m}=0.07 \times 10^{-2} \mathrm{~m} \times 4.20 \times 10^{-6} \mathrm{~m} / \mathrm{s}=2.94 \times 10^{-9} \mathrm{~m}^{2} / \mathrm{s}$
   
   $D_{h}=D_{e}+D_{m}=\phi^{m} D_{0}+D_{m}$
   
   So $D_{h}=3.66 \times 10^{-9} \mathrm{~m}^{2} / \mathrm{s} \text { and } D_{m} \text { dominates } D_{h}(80.3 \%)$.
3. Residence time:
   $\tau_{a}=\frac{V_{\text {pore }}}{Q}=\frac{\phi L A}{u A}=\frac{\phi L}{u}=\phi \frac{L}{u}$
   
   So $\tau_{a}=0.1 \mathrm{~m} / 4.20 \times 10^{-6} \mathrm{~m} / \mathrm{s} \times 0.4=9.524 \times 10^{3} s==0.110 d$
   
   Characteristic time for hydrodynamic dispersion:
   $\tau_{d}=\frac{L^{2}}{D_{h}}$
   
   So $\tau_{d}=0.1^{2} \mathrm{~m}^{2} /\left(3.66 \times 10^{-9} \mathrm{~m}^{2} / \mathrm{s}\right)=759 h=31.6 d$
   
   $P e=\frac{\tau_{d}}{\tau_{a}}=\frac{L u}{\phi D_{h}}$
   
   So $P e=31.6 d / 0.110 d=287>>1$, advection process dominates transport.

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:

Figure 4. Breakthrough outlet concentration of Br^- as a function of number of residence time. C₀ is the initial injecting concentration of the tracer Br^- while C is the tracer concentration at the outlet.

Text version of Figure 4 [40]

Summary

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.

Other Relevant Textbooks:

1. Flow and Transport in Porous Formations (1989) by Gedeon Dagan;

2. Dynamics of fluids in porous media (2013) by Jacob Bear, Courier Dover Publications;

3. Applied contaminant transport modeling (2002), 2^nd edition, by Chunmiao Zheng and Gordon D. Bennett, John Wiley and Sons, Inc.

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).

Figure 1. A. The spatial distribution of permeability in an biostimulation field experiment plot at Old Rifle, Colorado. Circles are injection wells and squares are monitoring well. B. The spatial distribution of injected tracer (bromide). The injected tracer tends to bypass the low permeability zone (modified from (Li et al., 2010))

The general Advection-Dispersion Equation (ADE) in multiple dimensions is as follows:

Here C is the tracer concentration (mol/m³ pore volume), t is time (s), D is the combined dispersion–diffusion tensor (m²/s), v (m/s) is the flow velocity vector and can be decomposed into v_L and v_T 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 D^e(m²/s). At any particular location (grid block), their corresponding diffusion / dispersion coefficients D_L (m²/s) and D_T (m²/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).

Figure 2. Illustration of spatial heterogeneity in the ground surface where a 1- to 2- meter thick, silty-clay layer is underlain by sand, sandy gravel, and gravelly sand layers.

Credit: Zheng and Gorelick, 2003. Use with permission.

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).

Figure 3. A. Three-dimensional visualization of the horizontal hydraulic conductivity distribution based on ordinary kriging of the flowmeter data for the MADE site (Zheng et al., 2011). B. observed MADE-2 tritium plume at 328 days after injection

Zheng and Gorelick 2003).

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).

Mean of permeability

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).

Variance of permeability

The extent of deviation from the mean is calculated as follows:

where $\operatorname{Var}(\kappa)$ is the permeability variance (m⁴). 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

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 (m²) with different correlation length generated from Gaussian Sequential Simulation Method.

Figure 4. Representative permeability fields (m²) generated using Gaussian Sequential Simulation Method with different correlation lengths.

Credit: Deutsch and Journel, 1992