PNG 550
Reactive Transport in the Subsurface

0.4 CrunchFlow

PrintPrint

CrunchFlow was developed by Carl I. Steefel 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. 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.

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
Downloads:
CrunchFlow manual; CrunchFlow running notes
Windows32bit: executable, dll file; Windows64bit: executable, dll file
Example files to try: input File, Database File