Lesson 0: Introduction to Reactive Transport Modeling

Overview

This lesson introduces Reactive Transport Models (RTMs), primarily focusing on a brief history of RTM development, governing equations, and key concepts. It includes a lightboard video about the governing RT equations, a video that introduces the code that we will use in this class, CrunchFlow, and the reading materials here. The idea here is to give you an overview of RTM. Many concepts introduced here will be detailed in later lessons so it is OK if you do not fully grasp them in this lesson.

Please note that this is not a course that teaches how to numerically solve for reactive transport equations, which deserve a separate course by itself. Instead, this is a course that teaches fundamental reactive transport concepts and how to use an existing software CrunchFlow to solve and answer specific questions. So this is a model application course, not a numerical method course.

Learning Outcomes

By the end of this lesson, you should:

Know the RTM development history and recognize the names of general available reactive transport codes;
Understand what RTMs do in a broad sense and how they are different from other types of models ;
Gather what you need to run CrunchFlow simulations;
Run CrunchFlow simulations using the provided example input and database files.

There are example files and hw files in each lesson (almost). If you would like extra exercise files, click here [1].

Lesson Roadmap


(Optional) Reading	Steefel, C.I. 2007. Geochemical kinetics and transport. in Kinetics of Water-Rock Interaction. S. Brantley, J.K., A. White (ed), pp. 545-589, Springer New York. Steefel et al., 2005. "Reactive transport modeling: An essential tool and a new research approach for the Earth Sciences." Earth and Planetary Science Letters 240: 539-558. Li et al., 2017. "Expanding the role of reactive transport models in critical zone processes." Earth Science Reviews. 165: 280-301 CrunchFlow manual Note: these are your references. You can skip through quickly in this lesson to get some overview idea. You will use these again and again later.
To Do	Read the online lesson material; Watch the lightboard video (~ 7 minutes) about reactive transport equations; Watch the orientation video 0lessonCrunchFlowOrientation; Run the simulation with the files provided; let the professor and / or teaching assistant know of any problems that arise.

Questions?

If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.

0.1 History of reactive transport modeling

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

Schematic representation of the oxidation-reduction zones that may develop in an aquifer- see caption

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

0.2 Reactive transport equations and key concepts

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)

Reactive Transport Reactions
Click for a transcript of the Reactive Transport Reaction video

Reactive Transport Reactions

PRESENTER: What I'm going to do today is really to introduce the general idea of reactive transport equation and the overview of it. This equation that I put here is we call mass conservation equation for chemical species in aqueous phase, for one of the representative species. I don't expect you to know every term, or the details of every term in this equation. Because we will be talking more about each term later. But this is really to give you a general idea and overview before we start about everything.

And the importance of this is that we run these reactive transport codes. And there are a lot of built up architecture behind the code, in terms of what they solve. And it's always a good idea to know some of these, like what they solve, and what are the things behind these codes. So what I will do today is really talk through each of these terms and talk about the [? fake ?] meaning of each term, so that you get a general idea what are they really for.

So the first term we call the mass conservation term. It's called mass accumulation rate. And it should have the units of mol per length cubed, which is volume, of porous media per time, whatever time you pick. But all the terms have to be consistent, have the same length and time units. So this equation really says mass accumulation rate depends on several different processes, right?

So the first term is the rate itself. And the overall rate, the second term, what we call dispersive and diffusive transport. So you think about how a chemical species in water has a change over time, which is this. So some of these rates coming from, for example, the chemical species gold to get different concentrations in different locations. For example, you think about a dye put in a cup of water. And over time, they tend to have the same color everywhere. So this is one of the driving forces, in terms of what we call diffusive or dispersive transport. And they should have the same unit as the first term. Every term should have the same unit. So that's a second term.

And then, the third term is what we call the advective transport. And this is a process where, for example, you think about rivers, right? And the chemical species will flow together with the water. And so essentially, the water brings the chemical species to different places. So this advective transport, in this term, you have the hue, which is we call Darcy velocity. And then, the concentration actually I probably should explain here, the Ci will be the concentration of one chemical species of species i, a representative species, i. And the Ci everywhere is the same. So this is advective transport.

But also in a lot of systems, you have reactions, right? So the last term, four, is for the total reaction rates. And again overall, it has the same units as the first term. But essentially, this could it be a summation of several different-- let's call this equal to summation of i, which is for the chemical species, i. But this i could be involved in, let's say, ik, different numbers of reactions. So essentially, you will be adding all these reactions that this can chemical species, i, is participating in. And this ik would be a total number of reactions.

OK, so this is one representative equation. We write this equation form for the species, i. But I put the i from 1 to n here, meaning you can have any arbitrary number of what we call a primary species. So if you have, let's say, 10 different chemical species, then the primary species you have n equal to 10. And you will be writing 10 of these equations.

And these equations, essentially, if we solve these equations, you get the concentration or different chemical species as a function of time and space. So essentially, the outcome of this is the temporal and spatial distribution of chemical species. So essentially, you can tracing after each chemical species and look at how they change as a function of time and space and how, in different parts, they have different rates, and all that. Species i and i from 1 to n.

OK, so that's what you hope you guys will be exposed to in using this code. You will be solving this equations for particular, specific questions, problem applications. And then, you're usually given a set of initial conditions, like where is the concentration of different species at time 0 at different locations, and then over time, how these countries have different species evolve over time. And you will see you will learn a lot about these, using this code. And we'll be talking more about this in each of these terms, what are they, over time in different lessons that will follow.

Source: The Pennsylvania State University

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

\frac{\partial (C_{i})}{\partial t} + \nabla \cdot (- D \nabla C_{i} + v C_{i}) + \sum_{r = 1}^{N_{r}} v_{i r} R_{r} + \sum_{m = 1}^{N_{m}} v_{i m} R_{m} = 0, i = 1, n

(1)

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:

D_{L} = D^{*} +_{a L} v_{Z}

(2)

D_{T} = D^{*} + a_{T} v_{T}

(3)

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.

0.3 Course Structure

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.

0.4 CrunchFlow

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:

CrunchFlow - a reactive transport code
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)
Input file: contains information the code need to set up the systems
Database file: includes information related to reaction stoichiometry, reaction equilibrium constants, and reaction rate parameters

In input file:

Comments should start with “!” and start from the beginning of the line, not middle of the line.
Do not use Tab. Use space if you need to lay out things neatly.
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.
Search CrunchFlow manual for the syntax of keywords.
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.

Click here for a transcript of the Crunch Flow Orientation video.

PRESENTER: Let's talk about the CrunchFlow. So this short video is supposed to be orientation for CrunchFlow. I introduced in the lesson, in the orientation lessons, that there are a lot of other reactive transport codes. Besides CrunchFlow, you have TOUGHREACT, PFLOTRAN, PHREEQC, [INAUDIBLE]. Different reactive transport has a bit of different flavors, but they all solve similar mass conservation equations for different species.

So the backbone of different reactions code are the same. But they have different solvers. They could have different computational speed because of use of different solvers. Also, using different solvers, it could have different convergence capabilities. So the way it's numerically solved, it could be very different. And that directly determines how fast things are going to run in your system.

I choose to teach CrunchFlow for several reasons among all of these other reactive transport codes. Mostly because I think it's very flexible in terms of setting up At the beginning, the learning curve might be a little bit-- you might think it's a bit sharp because of it's not a GUI interface. It's not graphic. And you might not be used to text file and opening objects exactly. But once you get used to it, this allows you a lot of flexibility to run things and put in things that are much faster than a lot of other reactive transport code does.

And also, it uses advanced solver like this [INAUDIBLE] package which has a lot of fast solvers. package was developed in a national lab. I believe it's Argonne, if I remember correctly. So it runs fast. And it tends to be robust. It usually doesn't get into non-convergence problem unless your system is ill conditioned or there's a bit of somewhere, there's a bug there.

But then it also is very flexible in terms of setting up heterogeneous systems. Like if you have different low permission, high permission, some mineral in somewhere that is more abundant than in other locations and less abundant, you can act expressly set it up in CrunchFlow, which I think it's the biggest advantage for CrunchFlow.

So in order to be able to run CrunchFlow, you need-- actually, so in this class, I put all the different executables on the course website. So there's also the manual files. There's a folder of CrunchFlow exercises which has a lot of different files there that [INAUDIBLE], when he taught CrunchFlow short courses before, he tended to use these exercise files.

So you can see these exercises in addition to what we teach here. And if at the end, we are running your project, you will be looking for things that you can use. And those exercises could be a very nice complementary in terms of adding additional examples and template input files for you to use.

So I think in these classes, when you try to learn CrunchFlow, you should use the menu frequently. You can search for the keywords and look at the explanation. You don't have to read everything at the beginning because I think that doesn't-- if you do that, the information doesn't register as much because you haven't gotten much of a handle on understanding the code yet. So I always advertise, learn it when you use it. I think that's where the information registers the best.

So what I have here is a folder that has four files-- so executable CrunchTope 64-bit for the Windows system I have. There's a DLL. There's the input file and DBS file. So these four files are the necessary files in order to run the CrunchFlow in a folder. You can not just set up exactly in the system to make it run default without copying to different folder. But for now, let's just keep it simple and keep it in this folder.

So we have this example, which we knew is a family of the CrunchFlow executable. It's called now CrunchTope because it actually incorporated isotope geochemistry. So the code should be able to run isotope information and everything.

This DLL file is a library file. How the code does is when it's solving the equations, it actually calls some dynamic library. So this library file is necessary in order for the code to finish. And then you have the input file and the database file.

So we are going to introduce these two files a bit just to introduce the structure of it. I'm not going to talk about it in detail, but just give you kind of a sense of what's going on. And then in the later lessons, most of the keywords will be introduced. And you learn by working on exercises and examples.

So the input file is mostly divided into keyword blocks. For example, here you have title block, which you can have a name and introduce what is this input file about it. It's good record keeping. Every keyword block, typically tried to differentiate keyword block name and keyword. So you use a capital to indicate the keyword block names. So it's easier to see. And keyword block should end with an END, like what we have here.

So typically, this is flow and transport in 1D system. So there's multiple keyword blocks here. There's RUNTIME, there's OUTPUT, DISCRETIZATION, and all that. Today, I probably will just introduce the RUNTIME mostly.

A lot of these keywords, it's related to a numerical solution process. Typically, I don't suggest changing early on because you might-- whatever default that you use it. One thing that's sometime useful is this specifies a solver, just name a solver. So those several different solvers if you look up in the menu. Several different solvers, you potentially can change one to the other if you have, for example, non-convergence problem at some point.

There's a pc, pclevel. They're all explained in CrunchFlow. So typically when I do it, I tend to kind of open the manual and I learn it. I tend to look at all that and try to learn what it means, why I'm using that, for example, pclevel, pc, pclevel. So the manual explains what it means. So you can dig into that and do things.

And one thing that is important to keep in mind is for you not to put in comments. For example, there's a line you want to explain what the number is from and what it means and what are units and everything. You can use a command line starting with this, exclamation mark.

You cannot do exclamation mark in the middle of line, for example, like doing this this. The code is not going to recognize. So you always need to start a new line with starting with the exclamation mark to do comments. So the code will be actually read as a text not as a keyword.

One thing it's important to keep in mind is your name of database file need to the exact same as the name of a database file in the same folder. For example, here we have datacom.dps. You should put that there. If you put another name there and then you don't have that database in this folder, then it's not going to run. Or if you have multiple databases in the same folder, you need to specify which one you are going to use.

These are the different ways of solving. So screen output is the number of steps that you skip before the screen actually has an output. If it's one, it's putting up every timestamp. If it's 10, it's putting up every 10 timestamps.

And then there's OUTPUT keyword block. You can putting time unit like minutes, seconds, hours, et cetera. So this is the time. You can put several times you want to see the snapshot of the system. And the last one is how long the system is going to run to. So in this system, it's going to run to 250 minutes.

Time series is where you want to have, for example, in a particular grid block, you want to have a times series from time beginning to end of simulation. You want to see how the concentration, for example, evolve over time. So that's the name of the output file that is going to have that information. And this specifies which grid block you want to have that time series. You can do multiple times series with different names and different grid blocks.

And this time set interval determines how fast-- how many-- like every once or every time step the code is going to write on the breakthrough curve. If you want that output for every 10-- every 10 times timestamps, that's fine too. Then it makes the breakthrough file a bit smaller because you have less output. Anyway, so this is output.

I'm not going to introduce everything else in the system because we will be introduced in detail all the keyword blocks in one of the lessons on 1D physical process, with diffusion dispersion process. So I'm going to leave these to there for the introduction of these blocks.

So all we teach of this is more or less get familiar with the keywords for each different process. What do you need to do in order to setting up this meaning for it and making sense.

Now before we close, let's talk about the database. So this datacom. This is a long file. And it's borrowed from EQ3 database, essentially. It started with a line called temperature points. So these are different. There are eight different temperature points.

If you want to understand what are the details, for example, there are eight different temperature points. These are for different degree Celsius. So there's log K when we talk about it. This is a reaction, different reactions at eight different K values corresponding to the eight different temperatures, which is already done in the EQ3/6.

So if you are specifying temperature is 150, the code is going to read the log K value for the different reactions from the fifth number, not from the default 25 degree Celsius. That's what it's for. So that's the temperature point.

These are the Debye-Huckel coefficients. Because the code uses the Debye-Huckel equation to calculate activity coefficients. And then starting after the Debye-Huckel block, you have from starting from water should be end of primary. So the answer from up until here, it's a list of primary species.

And each primary species has three parameters there. The first one is the Debye-Huckel parameters that you can use for that species. This is charge. And this is molecular weight for all the different primary species. So you can actually add in your own primary species if you use artificial makeup and elements or something. So sometimes, you need to manipulate code to do something that you want. There should be a presence in the primary species because the primary species is essentially the building block of the system.

So after primary species, you would have, starting from the line-- following the end of primary, you will start to have the secondary species. So all these secondary species are written, for example, like this. Let's look at HS.

So for the relevant primary species for the sulfide species is a sulfate. It didn't exist in this database. You actually can put in another sulfide if you want. But here, the primary species here is sulfate. So there's no sulfides there. You can do the sulfide species when you give it another name. So you can specify the sulfide species as written in terms of redox reaction with oxygen. So a sulfide becomes reduced to become sulfide. Or it's just plus O2 gas to become sulfate being oxidized.

Let's look at another maybe bicarbonate or something. For carbonates, look at carbonate. For carbon species, what you specify in the system is-- here, this CO2. There should be a CO2. Let's see. Let's search for. Oh, yeah. There's a bicarbonate there. So the bicarbonate is a primary species.

And then everything else, for example, CO2 aq or bicarbonate should be written as bicarbonate. So here when it's-- so let's look at this line specifically. This is essentially writing CO3 bicarbonate species. It involved two different species. One is hydrogen ion. The other is bicarbonate.

So if you write it, it'd be CO3 minus minus plus. When it's minus meaning it's in left hand side of equation. Plus equals two bicarbonate. So the primary species need to be in the right hand side.

And what are these? So these should be the equilibrium constants. These are not the right equivalent constant. It's the equivalent constant for following the way the reaction is written. So it's saying, bicarbonate species can be-- if you write it that way, it will be 9.6 log K.

This is the log K value. So this is not 10.03, which means it's at a different temperature. And it's using the same value. So it's not corresponding to the temperature.

So if you want to look at a temperature effect, every reaction, later on when you're writing, you will need to specify. If you are actually running high temperature, you want to look up this number to make sure this number in database is correct, just to confirm because sometimes these databases are changed by people. So you want to make sure the numbers are right. Because these are very important numbers.

If we look at CO2 aq, it's essentially saying CO2 aq plus water equals H plus and bicarbonate. So there's three species involved. And for this reaction, you have the log K equal at zero is minus 6.58. And at high temperature, it keeps on decreasing. So this set of numbers is making sense.

But you still want to check. Every number you use and you end up using in the modeling process, you want to check these numbers. So 3.0 is, again, is the Debye-Huckel. So these three numbers is essentially like the same three numbers as used in the primary species. It was Debye-Huckel parameters, charge, and molecular weight. So that's the list of secondary primary.

And I do Control-F, I will see end of secondary. So it will say end of secondary. This is towards the end of the secondary species list. Everything in secondary species is written in terms of primary species. And then you see a block of everything with a g. Process g means a gas phase. So all the gas phases are listed here. And then there should be an end of gases. These are end of gases.

And then after that, after gases, you would have all these mineral phases. So you can see the different reactions, different-- let's search for calcite. Calcite, so again, for the calcite, you have molecular weight. The reaction written in terms of two species. This is a stoichiometric coefficient, one calcium and one carbonate.

And again, you have these are the log K. Again here, the value is not necessarily right because you have all the same. But I believe this number should be the log K equals at 25 degrees Celsius.

For the mineral phase, you only have molecular weight. You don't have Debye-Huckel and charge because these are not in the water for aqueous phase. So then this should be-- so each mineral is written this way. So if you go to end minerals, it's the end of mineral phase.

And then it will begin the surface complexation session. It will specify the way surface complexation is written as well as, again, this is a reaction. And then this log K value in different conditions. Every time you see all the 500.000, these are numbers that are kind of dummy numbers that you put in because you don't know what number to use.

So this number means if you are looking, for example, at this here, we know that 25 degree Celsius log K for this reaction is 8.82. But for all other temperatures, we don't know. So if you need to use at these other temperatures, you need to look it up in the literature. And then this ending end of surface complexation.

Now the other thing I want to point out is there's also surface complexation at end. So there's the beginnings of the complexation parameters that give you the charges of the different surface species. So these are actually the surface charges of the surface reactions directly related to surface complexation. So this is surface complexation.

Then you have aqueous kinetic block and then mineral kinetic block. Again, I always suggest that you will be looking up literature for real numbers. For example, let's say here you are specifying, let's look at calcite again as an example.

So for example for calcite here, it's specified several different brick blocks for different [INAUDIBLE]. So every one has a label. In this block, the label default and you are actually looking at this reaction and with dependence on H plus. Here, it didn't specify dependence, so it should be-- it says no dependents here.

Like here, it will be dependence active to H plus raised to 1.0. And this would read, constant activation energy, the type of [INAUDIBLE] used, and all that. And for the other one, for this one, you're saying is dependent on CO2 aq to the 1.0 for the calcite dissolution reaction. So some figure them somewhere else, you also write pyrite oxidation reaction and you have all these.

So you can also in the input file put the rate constant like this. So if you're putting input files your rate constant, then it's going to ignore what you have in the database and use the one in the input file. Anyway, so that's for the mineral kinetics.

And if the mineral that you are interested in is not already there, you can add another block specifying all these numbers. And as long as you specify, for example, the mineral is part of the mineral list, if you have a new mineral that is not in the list, you can add that in in the mineral list which specifies a reaction's log K values for that mineral. And then you add another mineral block for the kinetics here. So that's for the mineral kinetics.

And then you would have the exchange block, you can specify a multiple site ion exchange or a single site ion exchange. These are specified, essentially. Again, the reactions for the half life ion exchange reaction the exchange coefficient.

So after this, at the end, you have this. After this, so later part, you have some numbers These tend to be [INAUDIBLE]. You want to move around a little bit, you can just temporary start there. So you can more or less ignore this line if you want. I could just remove all that.

So let's do just a run. So what you do is when you try to run control, you click on the executable file. And then you will be putting the name. It will ask for the name of the input file transport. 1D.ing. You always have your input file in .ing so that the code recognizes it's the input file.

Now, after this you see a lot of more files. Because with all these, you should look at extension. All these output file are the output. So because I specify in the special profile one time, it will be output at Z1 time for all the different velocity, for the total concentration which we'll be talking that later, total concentration speciation, pressure, gas, concentration. And there's also an isotope file, and also a breakthrough curve because I want this time series. So you will see a lot of these output files.

If you change your input file at this point to a, let's say it's such and such a number, and rerun again here, the new output file will replace the old output file. So if you want to keep a record of both, you need to have another folder to run the simulation in order to make it be able to kind of compatible to still have the record of each. Otherwise, it's going to be overwritten.

Let's look at transport 1D out. It's essentially kind of echoes what has been [INAUDIBLE]. It's essentially going through the file to see the input file and kind of keep a record of what has been [INAUDIBLE]. And then at end, it says, OK, I'm going to stop starting, initialization completed, starting timestamping.

So sometime if the code has stopped running in the middle or something, you can look at this file to see if everything has been written in correctly. If not, that means there's something wrong in your input file. Another thing to keep in mind is in your input file, you're not supposed to the use tab. You have to space if you want to have things lay out nicely and everything. Don't use tab. The code doesn't recognize tab.

Let's look at the breakthrough. So here what we specify is the species bromide. So this is essentially for the [INAUDIBLE] cell, which is at outlet. You have time. You have bromide concentration. So early on, they'd run very-- because I put one [INAUDIBLE] for time interval. So it run every timestamping. So early on, it's like starting from 10 to minus 10 minutes. It's essentially initial condition. You don't have any bromide. And then in the input file, in that you inject input files, as the bromide [INAUDIBLE]. So you start to have higher and higher concentration at some point later.

So the inner concentration is this. So at somewhere, for example, around maybe 100 minutes or something, you start to have something come out. And then eventually this concentration stabilized to the end. So this is a system of just transport.

And these should be the log concentration. This is the concentration of bromide and log concentration. So this runs. And if you need to end this concentration, total concentration, for this non-reactive species, it's not really important.

So I'm going to stop here. And hopefully you have a bit of a sense of what's going on. Several things that you should remember is you always need the executable file, database, input file, and the library file to run the simulation. And the CrunchFlow manual will be your good friend in this class, besides the reading materials that you will learn about the general principles.

I'm going to stop here. So from Lesson One, we'll introduce four individual processes, what are the important keywords and parameters, and we'll go through exercises, and then run some simulations, do homework to kind of re-emphasize what you have learned in class. All right. That's what do we have for--

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

Downloads:

CrunchFlow manual [5]; CrunchFlow running notes [6],

Windows32bit: executable [7], dll file [8]; Windows64bit: executable [9], dll file [10]

MacExe: exeCompiledWithIntel11.1 [11]; exeCompiledWithIntel14.03 [12]

Example files to try: input File [13], Database File [14]

References

Apps, J. A., L. Zheng, Y. Zhang, T. Xu and J. T. Birkholzer (2010). "Evaluation of Potential Changes in Groundwater Quality in Response to CO2 Leakage from Deep Geologic Storage." Transport in Porous Media 82(1): 215-246.

Atchley, A. L., A. K. Navarre-Sitchler and R. M. Maxwell (2014). "The effects of physical and geochemical heterogeneities on hydro-geochemical transport and effective reaction rates." Journal of Contaminant Hydrology 165(0): 53-64.

Bao, C., H. Wu, L. Li, D. Newcomer, P. E. Long and K. H. Williams (2014). "Uranium Bioreduction Rates across Scales: Biogeochemical Hot Moments and Hot Spots during a Biostimulation Experiment at Rifle, Colorado." Environmental Science & Technology 48(17): 10116-10127.

Bao, C., Li, L., Shi, Y. and Duffy, C. (2017) Understanding watershed hydrogeochemistry: 1. Development of RT-Flux-PIHM. Water Resources Research 53, 2328-2345.

Beaulieu, E., Y. Goddéris, D. Labat, C. Roelandt, D. Calmels and J. Gaillardet (2011). "Modeling of water-rock interaction in the Mackenzie basin: Competition between sulfuric and carbonic acids." Chemical Geology 289(1–2): 114-123.

Bethke, C. (1996). Geochemical reaction modeling : concepts and applications. New York, Oxford University Press.

Bolton, E. W., A. C. Lasaga and D. M. Rye (1996). "A model for the kinetic control of quartz dissolution and precipitation in porous media flow with spatially variable permeability: Formulation and examples of thermal convection." Journal of Geophysical Research: Solid Earth 101(B10): 22157-22187.

Boudreau, B. P. (1996). "A method-of-lines code for carbon and nutrient diagenesis in aquatic sediments." Computers & Geosciences 22(5): 479-496.

Brown, B. and D. Rolston (1980). "Transport and transformation of methyl bromide in soils." Soil science 130(2): 68-75.

Brunet, J. L., L. Li, Z. T. Karpyn, B. Kutchko, B. Strazisar and G. Bromhal (2013). "Dynamic evolution of compositional and transport properties under conditions relevant to geological carbon sequestration." Energy & Fuels.

Chapman, B. M. (1982). "Numerical simulation of the transport and speciation of nonconservative chemical reactants in rivers." Water Resources Research 18(1): 155-167.

Chapman, B. M., R. O. James, R. F. Jung and H. G. Washington (1982). "Modelling the transport of reacting chemical contaminants in natural streams." Marine and Freshwater Research 33(4): 617-628.

Druhan, J. L., C. I. Steefel, M. E. Conrad and D. J. DePaolo (2014). "A large column analog experiment of stable isotope variations during reactive transport: I. A comprehensive model of sulfur cycling and delta S-34 fractionation." Geochimica Et Cosmochimica Acta 124: 366-393.

Fang, Y., T. D. Scheibe, R. Mahadevan, S. Garg, P. E. Long and D. R. Lovley (2011). "Direct coupling of a genome-scale microbial in silico model and a groundwater reactive transport model." Journal of Contaminant Hydrology 122(1-4): 96-103.

Godderis, Y., L. M. Francois, A. Probst, J. Schott, D. Moncoulon, D. Labat and D. Viville (2006). "Modelling weathering processes at the catchment scale: The WITCH numerical model." Geochimica et Cosmochimica Acta 70: 1128-1147.

Helgeson, R. C. (1968). "Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions: I. Thermodynamic relations." Geochim. Cosmochim. Acta 32: 853–877.

Helgeson, R. C., R.M. Garrels and F. T. MacKenzie (1969). "Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions: II. Applications,." Geochim. Cosmochim. Acta 33: 455– 482.

Kang, Q. J., P. C. Lichtner and D. X. Zhang (2006). "Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media." Journal of Geophysical Research-Solid Earth 111.

Knutson, C. E., C. J. Werth and A. J. Valocchi (2005). "Pore-scale simulation of biomass growth along the transverse mixing zone of a model two-dimensional porous medium." Water Resources Research 41(7).

Li, L., N. Gawande, M. B. Kowalsky, C. I. Steefel and S. S. Hubbard (2011). "Physicochemical heterogeneity controls on uranium bioreduction rates at the field scale." Environmental Science and Technology 45(23): 9959-9966.

Li, L., C. A. Peters and M. A. Celia (2006). "Upscaling geochemical reaction rates using pore-scale network modeling." Advances in Water Resources 29(9): 1351-1370.

Li, L., C. I. Steefel and L. Yang (2008). "Scale dependence of mineral dissolution rates within single pores and fractures." Geochimica Et Cosmochimica Acta 72(2): 360-377.

Li, L., Bao, C., Sullivan, P.L., Brantley, S., Shi, Y. and Duffy, C. (2017a) Understanding watershed hydrogeochemistry: 2. Synchronized hydrological and geochemical processes drive stream chemostatic behavior. Water Resources Research 53, 2346-2367.

Li, L., Maher, K., Navarre-Sitchler, A., Druhan, J., Meile, C., Lawrence, C., Moore, J., Perdrial, J., Sullivan, P., Thompson, A., Jin, L., Bolton, E.W., Brantley, S.L., Dietrich, W.E., Mayer, K.U., Steefel, C.I., Valocchi, A., Zachara, J., Kocar, B., McIntosh, J., Tutolo, B.M., Kumar, M., Sonnenthal, E., Bao, C. and Beisman, J. (2017b) Expanding the role of reactive transport models in critical zone processes. Earth-Sci. Rev. 165, 280-301.

Lichtner, P. C. (1985). "Continuum Model for Simultaneous Chemical-Reactions and Mass- Transport in Hydrothermal Systems." Geochimica et Cosmochimica Acta 49(3): 779-800.

Lichtner, P. C. (1988). "The Quasi-Stationary State Approximation to Coupled Mass- Transport and Fluid-Rock Interaction in a Porous-Medium." Geochim. Cosmochim. Acta. 52(1): 143-165.

Lichtner, P. C., C. I. Steefel and E. H. Oelkers, Eds. (1996). Reactive transport in porous media. Reviews in Mineralogy, MIneralogical society of America.

Liu, C., J. Shang, H. Shan and J. M. Zachara (2013). "Effect of Subgrid Heterogeneity on Scaling Geochemical and Biogeochemical Reactions: A Case of U(VI) Desorption." Environmental Science & Technology 48(3): 1745-1752.

MacQuarrie, K. T. B. and K. U. Mayer (2005). "Reactive transport modeling in fractured rock: A state-of-the-science review." Earth-Science Reviews 72(3-4): 189-227.

Mayer, K. U., S. G. Benner, E. O. Frind, S. F. Thornton and D. N. Lerner (2001). "Reactive transport modeling of processes controlling the distribution and natural attenuation of phenolic compounds in a deep sandstone aquifer." Journal of Contaminant Hydrology 53(3-4): 341-368.

Mayer, K. U., E. O. Frind and D. W. Blowes (2002). "Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions." Water Resources Research 38(9).

Molins, S., D. Trebotich, C. I. Steefel and C. Shen (2012). "An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation." Water Resources Research 48.

Navarre-Sitchler, A. K., R. M. Maxwell, E. R. Siirila, G. E. Hammond and P. C. Lichtner (2013). "Elucidating geochemical response of shallow heterogeneous aquifers to CO2 leakage using high-performance computing: Implications for monitoring of CO2 sequestration." Advances in Water Resources 53(0): 45-55.

Ortoleva, P., E. Merino, C. Moore and J. Chadam (1987). "Geochemical self-organization, I. Feedbacks, quantitative modeling." Am J Sci 287: 979–1007.

Parkhurst, D. L. and C. A. J. Appelo (1999). User’s guide to PHREEQC (Version 2)—A computer program for speciation, batch reaction, 1 dimensional transport, and inverse geochemical calculation: 310.

Regnier, P., S. Arndt, N. Goossens, C. Volta, G. G. Laruelle, R. Lauerwald and J. Hartmann (2013). "Modelling Estuarine Biogeochemical Dynamics: From the Local to the Global Scale." Aquatic geochemistry 19(5-6): 591-626.

Roelandt, C., Y. Goddéris, M. P. Bonnet and F. Sondag (2010). "Coupled modeling of biospheric and chemical weathering processes at the continental scale." Global Biogeochemical Cycles 24(2): GB2004.

Saunders, J. A. and L. E. Toran (1995). "Modeling of radionuclide and heavy metal sorption around low- and high-{pH} waste disposal sites at {Oak Ridge, Tennessee}." Appl. Geochem. 10(6): 673-684.

Scheibe, T. D., Y. Fang, C. J. Murray, E. E. Roden, J. Chen, Y.-J. Chien, S. C. Brooks and S. S. Hubbard (2006). "Transport and biogeochemical reaction of metals in a physically and chemically heterogeneous aquifer." Geosphere 2(4): doi: 10.1130/GES00029.00021.

Scheibe, T. D., R. Mahadevan, Y. Fang, S. Garg, P. E. Long and D. R. Lovley (2009). "Coupling a genome-scale metabolic model with a reactive transport model to describe in situ uranium bioremediation." Microbial Biotechnology 2(2): 274-286.

Soetaert, K., P. M. J. Herman and J. J. Middelburg (1996). "A model of early diagenetic processes from the shelf to abyssal depths." Geochimica et Cosmochimica Acta 60(6): 1019-1040.

Soler, J. M. and U. K. Mader (2005). "Interaction between hyperalkaline fluids and rocks hosting repositories for radioactive waste: Reactive transport simulations." Nuclear Science and Engineering 151(1): 128-133.

Steefel, C. I., C. A. J. Appelo, B. Arora, D. Jacques, T. Kalbacher, O. Kolditz, V. Lagneau, P. C. Lichtner, K. U. Mayer, J. C. L. Meeussen, S. Molins, D. Moulton, H. Shao, J. Šimůnek, N. Spycher, S. B. Yabusaki and G. T. Yeh (2015). "Reactive transport codes for subsurface environmental simulation." Computational Geosciences 19(3): 445-478.

Steefel, C. I., D. J. DePaolo and P. C. Lichtner (2005). "Reactive transport modeling: An essential tool and a new research approach for the Earth Sciences." Earth and Planetary Science Letters 240: 539-558.

Steefel, C. I. and A. C. Lasaga (1994). "A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems." Am. J. Sci. 294: 529-592.

Van Cappellen, P. and Y. F. Wang (1996). "Cycling of iron and manganese in surface sediments: a general theory for the coupled transport and reaction of carbon, oxygen, nitrogen, sulfur, iron, and manganese." Am J Sci 296: 197-243.

White, M. D. and M. Oostrom (2000). STOMP: Subsurface Transport Over Multiple Phases (Version 2.0): Theory Guide. Richland, WA.

Wolery, T. J., K. J. Jackson, W. L. Bourcier, C. J. Bruton, B. E. Viani, K. G. Knauss and J. M. Delany (1990). "CURRENT STATUS OF THE EQ3/6 SOFTWARE PACKAGE FOR GEOCHEMICAL MODELING." Acs Symposium Series 416: 104-116.

Xu, T. F., S. P. White, K. Pruess and G. H. Brimhall (2000). "Modeling of pyrite oxidation in saturated and unsaturated subsurface flow systems." Transport Porous Med 39(1): 25-56.

Yabusaki, S. B., Y. Fang, K. H. Williams, C. J. Murray, A. L. Ward, R. D. Dayvault, S. R. Waichler, D. R. Newcomer, F. A. Spane and P. E. Long (2011). "Variably saturated flow and multicomponent biogeochemical reactive transport modeling of a uranium bioremediation field experiment." Journal of Contaminant Hydrology 126(3-4): 271-290.

Yeh, G.-T. and V. S. Tripathi (1991). "A Model for Simulating Transport of Reactive Multispecies Components: Model Development and Demonstration." Water Resources Research 27(12): 3075-3094.

Yeh, G. T. and V. S. Tripathi (1989). "A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components." Water Resources Research 25(1): 93-108.

Yoon, H., A. J. Valocchi, C. J. Werth and T. Dewers (2012). "Pore-scale simulation of mixing-induced calcium carbonate precipitation and dissolution in a microfluidic pore network." Water Resources Research 48.

Summary and Final Tasks

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.