Brian Farrell

Hi Masoud, That plot is not taken directly from Gtplot, but is actually two plots (minerals and aqueous species) overlain together in a program like PowerPoint or Illustrator. The key is to use consistent units (mmoles here) and axes when overlaying one diagram on another. Of course, for the slopes of the lines method, the x and y axes must match each other. The point of the plot is to see what changes over the course of the reaction path. If a particular species is present at very small concentrations, or its concentration does not change at all, it can be removed from the plot since it is not involved in the reaction. If you plot the aqueous species, you'll see that HCO3-, CO2, and K+ are the only species which change concentration appreciably. Sometimes a delta plot can be helpful for figuring out which species or minerals are important in a reaction path. Hope this helps, Brian

Dear GWB Users, Just a gentle reminder: June 25 is the last day to early register for the pre-Goldschmidt Reactive Transport Modeling workshop. Register by Tuesday and save $100! The same deadline applies for early registration for the Goldschmidt conference, if you're attending. Regards, Brian

Dear GWB Users, Please join us October 25-26 in Denver, CO for a workshop on geochemical modeling using the GWB. Then stay on in Denver for the GSA Annual Meeting. Regards, Brian Farrell Aqueous Solutions LLC

Dear GWB users, We are pleased to announce our latest maintenance release, GWB 9.0.7. The 9.0.7 update improves functioning under Windows 8 and fixes all known glitches. Update from 9.0 through 9.0.6 at no charge to ensure you have all the newest features and bug fixes. Existing installations should automatically update to this release, unless auto-update is disabled. In that case, users should update their installations from the Help pulldown on any GWB app. Regards, Brian Farrell Aqueous Solutions LLC

Hi Frank, I think you'll want to write a control script for this. Try taking a look at the Control Scripts and Report Command appendices in the GWB Reference Manual for details on setting up a control script. Once you familiarize yourself with those, the Multiple Analyses appendix contains an example script in which SpecE8 runs analyses for a number of samples, then produces a file with the desired output. As an alternative, you could you the Plug-in Feature of React. You would write a program which would use React to handle the chemistry calculations, then retrieve the results to your program, which could write the results to a file of your exact specification. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hello Mustafa, I'll try to focus on how The Geochemist’s Workbench can make solving this problem easier for you. Perhaps someone else can check your hand calculations. A very brief background: GWB uses thermodynamic datasets in its calculations. Each dataset contains a list of Basis species (H2O, O2(aq) for O, H+ for H, HCO3- for C, Ca++ for Ca, etc.). These Basis species are the building blocks for constructing all secondary species (H2CO3 or CO2(aq), CO3--, etc.), redox species (Basis species in alternate redox states, like CH4(aq), CH3COO-, H2(aq), etc.), minerals, or gases (CH4(g), H2(g)) that you want to consider in your system. The thermo dataset is a compilation of reactions for these various secondary species, redox species, minerals, or gases, written in terms of the Basis species. Each reaction has an associated equilibrium constant at several discrete temperatures (generally ranging from 0 to 300 C). I recommend looking for the Geochemical and Biogeochemical Reaction Modeling textbook by Craig Bethke. Chapter 7, Redox disequilibrium, contains an example for the program SpecE8 which describes what you’re looking to do. Basically, you’ll want to take your water analysis and use it to add species to your Basis, then constrain your entries by setting pH, concentration, fugacity, etc.. You should decouple the redox pairs between species in different oxidation states (HCO3- and CH4, for example), then constrain the concentration of each independently. Since you have a partial pressure for CH4(g), you would replace CH4(aq) in your Basis with CH4(g) by performing a “Basis swap.” Then you would specify the fugacity of CH4, which is like the partial pressure in atm (be sure to use the total, not guage pressure, to calculate the partial pressure of CH4). You can do the same thing for CO2(g) by swapping it in for HCO3-. You can set the temperature of your system directly. The program will use the grid of equilibrium constants at different temperatures to solve for the equilibrium state at your conditions of interest. That way, you won’t have to mess around with calculating integration constants or using the Van ‘t Hoff equation. When you run your model, the program will perform a speciation analysis, listing the actual concentrations (and thermodynamic activity) of all the species in your system. It will also calculate the redox potential (Eh) for each electron donating and accepting half-cell reaction. By combining different donating and accepting reactions, you can calculate the free energy available from various redox reactions. Since you can easily modify your input (pH, fugacity, concentration, T) you can quickly run a number of analyses without being slowed down by the repetitive hand calculations. I'll contact you about a trial license of the software. Please let me know if you have any questions. Regards, Brian Farrell Aqueous Solutions LLC

Hi Maki, If you'd like to independently constrain the initial concentrations of species in various redox states, you'll need to decouple the appropriate redox pair(s). For example, if both lactic acid and HCO3- are present in your system, but not at equilibrium with the "Master Eh," then you'll need to decouple the Lactic_acid(aq)/HCO3- redox pair. Setting a "Master Eh" (O2(aq), Eh, pe, etc.) will have no effect on the distribution of mass between lactic acid and HCO3-. Only a kinetic rate law (redox or microbial) will control the manner in which a chemical species is transformed from one redox state to another. Because the Master Eh will only affect those redox couples which remain coupled, it doesn't really have meaning for decoupled redox pairs. You should probably track the Eh of the redox pairs that you are interested in. Regards, Brian

Hi Kyle, In the flush configuration, a reactant fluid (solvent water plus solutes) displaces the existing fluid from the system. Kinetic minerals added from the Reactants pane are not gradually titrated into the system (as if they were a simple reactant) or removed from the system (as in a flow-through model). Rather, their entire mass coexists with the fluid from the start, and they are allowed to dissolve into (or precipitate from) the fluid at a rate determined by the kinetic rate law you specify. The number you specify for "reactants times" will be multiplied by the mass you specify for a kinetic mineral, so you might divide the mass you specify by the "reactants times" to get the actual mass that you want. As an example, say you have a kg of solvent water with a few buffer minerals at equilibrium, plus an additional mineral which is not in equilibrium with that water-rock assemblage. Set up the equilibrium system in the Basis pane, set a flush configuration, add 1 kg water in the reactants pane, and specify a number for reactants times - let's say 10. If you want 100 g of some non-equilibrium mineral in your system, then add a kinetic mineral to the Reactants pane and specify a mass of 10 g. 100 g of the non-equilibrium mineral will be exposed to the equilibrium system. Hope this helps, Brian

Hi Masoud, Sections 13.2 and 15.1 of the Geochemical and Biogeochemical Reaction Modeling text describe the "slopes of the lines" method for determining the overall reaction that occurs in a system. You can also go to Config - Output to set the amount of information that you would like to appear in the text output file. By enabling the reactions loaded, you can view each reaction that is considered in your model, along with its equilibrium constant. Keep in mind this will only tell you the possible reactions, not which are actually occurring to any significant extent. Hope this helps, Brian

Hi Maki, In Gtplot, "Eh" generally refers to the "Master Eh" of the system. In an equilibrium model, you might set an overall oxidation state for your initial system (i.e. DO, Eh, O2(g) fugacity, aSO4--/aHS-, ...). This will affect all coupled redox pairs. When you start decoupling redox pairs, however, these are considered separately from the "Master Eh". Since you've decoupled Lactate, CH4, and H2, there are no more coupled redox reactions in your system, and Eh becomes identical to the Eh(O2(aq)/H2O). Keep in mind that Eh is relative to pH (try comparing the slanted lines of an Eh-pH diagram with the flat lines of a log f O2-pH diagram). Your system is becoming more acidic with reaction progress, so even though oxygen fugacity is decreasing, the "Eh" actually goes up a little. If you look at O2 fugacity, you'll see the decrease in oxidation state, though not by the amount you expect. Although most of the oxidant O2(aq) has been removed from your system, you don't have enough reductant to drive the redox state down to very low values. If you look at FD and TPF for the methanogen, you'll see that these are very close to zero (the reaction is limited both kinetically and thermodynamically). This is why your methanogen is not active. You need more H2(aq) in your system. Although reaction is possible at H2(aq) = 1e-5 molal, it won't proceed very far since that is still very little H2(aq). The overall Eh (or Eh(O2(aq)/H2O) in this model) has no effect in the model on the H2(aq)/H+ redox state, so that's why it is possible for the methanogen to be active even in the presence of oxygen. If the methanogen should be inhibited by oxygen physiologically, you could probably use a rate law script. For example, something like: IF molality"O2(aq)" > 0.001, THEN 20 ELSE 40 20: rate = 0 40: rate = rate_con * ... You are correct that H+ is produced by the lactate oxidation. The methanogen could use this H+, along with HCO3- and H2(aq) to drive its metabolism. However, 1e-20 molal H2(aq) is simply not enough for this to be possible. The way you've written the reaction for the methanogen looks correct. You just need a source of H2(aq) in your system. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Kyle, Nice to hear from you again after the Melbourne course. After you run the first model, take a look at the text file and compare the "Aqueous species" and "Original basis" sections. The former section lists the concentration and activity of all the actual species in the system (H+, CO2(aq), HCO3-, NaHCO3, CaHCO3, CO3--, etc.) while the latter lists the thermodynamic components (H+, HCO3-, etc.). When React picks up the your calculation results, it uses the components, not the species concentrations. It looks like you're taking values from the actual species, but you don't need to do this. As for the number of buffer minerals allowed, that is a limitation arising from thermodynamics, not the GWB. In a true equilibrium state, you can't have a fluid in equilibrium with "extra" phases. What you can do, however, is to add extra minerals to the Reactants pane. Section 26.5 of the Geochemical and Biogeochemical Reaction Modeling text (or Section 4.2 in the GWB Reaction Modeling Guide) contains an example that might be useful to you. In the example, a fluid at equilibrium with Kaolinite, Muscovite, and Quartz (but not feldspar) is exposed to a reactant mineral Albite, which dissolves into the fluid at a rate calculated from a kinetic rate law. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Johan, Just wanted to let you know that Maintenance release 9.0.6 is out. It includes a fix for reading pH in the rate law parser, so your original basic script should work now. Please let us know if you have any questions. Regards, Brian

Hi Jay, Just wanted to let you know that we've come out with maintenance release 9.0.6, which fixes a problem in which React takes too small time steps when running flow-through models containing both kinetic reactants and buffer minerals. I think that should take care of the problem you mentioned. Regards, Brian

Hi Catriona, Just wanted to let you know that maintenance release 9.0.6 is out. Included is a fix that might be of interest to you: Fix a problem in which React sets too small time steps when running flow-through models containing both kinetic reactants and buffer minerals. I think there is still an issue with kinetic reactions of widely varying rates in your script. For example, the rates of pyrite and calcite dissolution in your script are approximately 7 orders of magnitude faster than the slowest kinetic reactions. The fast reactions cause React to take extremely small time steps, bogging down the program. A good strategy is to divide mineral reactions into three groups: reactions that proceed quickly over the time scale of the calculation - assume equilibrium minerals that react negligibly over the calculation - suppress these reactions that proceed at a slow but measureable rate over your time scale of interest - use kinetics for these Hope this helps. By the way, have you made progress with an X1t model? Regards, Brian

Hi maki, Just wanted to let you know that maintenance release 9.0.6 is out. It includes the above-mentioned fix to the node-by-node editor for specifying x and y-discharge in X2t. Regards, Brian

Hello awebb, Just wanted to let you know that our latest maintenance release, 9.0.6, has the following change that should fix your problem: Disable control script mode when reading an input file from within a control script (unless the input file itself turns control script mode back on) and restore control script mode when finished. Please let us know if you have any questions. Regards, Brian Farrell Aqueous Solutions LLC

Dear GWB users, We are pleased to announce our latest maintenance release, GWB 9.0.6. The 9.0.6 update includes the following changes: Allow search for missing input files Add extra column and row, respectively, to x-discharge and y-discharge node-by-node editors when specifying the flow field in X2t Get correct value for pH when parsing rate laws Fix a problem in which React sets too small time steps when running flow-through models containing both kinetic reactants and buffer minerals Fixes for all known issues. Update from 9.0 through 9.0.5 at no charge to ensure you have all the newest features and bug fixes. Existing installations should automatically update to this release, unless auto-update is disabled. In that case, users should update their installations from the Help pulldown on any GWB app. Regards, Brian Farrell Aqueous Solutions LLC

Hi Kaizen, There is no one "correct" way to constrain an activity diagram. It's a very simplified, visual interpretation of a chemical system. What you've done, swapping in multiple different SO4-- complexes (SO4--, plus CaSO4 for Ca++, MgSO4 for Mg++, FeSO4+ for Fe++), is fine. Because mass action equations express the relationship between various species, you could just as easily constrain your Basis without any swaps (SO4--, Ca++, Mg++, Fe++) to get the same result. Just make sure you match up the correct activity for each species. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Howard, That warning message is certainly something you'll want to look into, but the charge imbalance on the reactants pane is actually not your main issue. Your initial system (which you've defined on the Basis pane) fails to converge before the Reactants pane is even considered. Try using SpecE8 to equilibrate your initial system before making a more complicated reaction path model. You should ask yourself whether all of those minerals are really in equilibrium with your system. As a simple test, I used Rxn to write the following reaction: Mg4Al2O7^10H2O + 8 H+ + 2 Quartz = 12 H2O + 4 Mg++ + Kaolinite I then set pH to 5.4, T to 25, and activity of H2O to 1. The resulting equilibrium equation states that the activity of the Mg++ ion in this assemblage should be 10^3.7888 - quite a large value. This leads me to believe there is something wrong with the way you defined your system. For example, a higher pH would result in a more reasonable Mg++ activity. Perhaps you should check whether the thermo data is correct. Another item to consider is whether it makes sense to have so much oxygen in a system that is in equilibrium with Pyrrhotite. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Maki, We're looking into your question regarding copying values pasted into Excel from the color map. As for your second question, molal concentration refers to moles of solute per kg solvent (water only), whereas moles per kg refers to kg of solution (water plus all dissolved solutes). Regards, Brian

Hello Bron, Running your example in GWB 8 (and earlier) I get a charge imbalance error of -70%. In GWB 9, which uses a different convention, I get -35%, which seems to agree with the PHREEQC and hand calculations. I'll look into this further, but I think the explanation offered in the thread you mentioned mixed up the two conventions. Based on a quick spreadsheet I made up, GWB 9 appears to use the same convention that you typed out here, while GWB 8 and earlier versions divide the denominator by 2 (which is the same as multiplying the whole thing by 2). If you'd like to check your own calculations using GWB 9, I think we can arrange to send you a trial version. Just send an email to support@gwb.com. Regards, Brian Farrell Aqueous Solutions LLC

Hi Maki, The X2t file looks fine. Keep in mind, you've set the total concentration of Cl- in the Initial pane of X2t (also known as a bulk constraint, or the Cl- component). This value, 1 mol/L, is the sum of the various Cl species in solution: Cl- and NaCl(aq). Due to mass balance on the Cl- component, if you look at the concentration of the Cl- and NaCl(aq) species in Xtplot, you'll find that these add up to 1 mol/L. In the printed output, the 1.009 value you mention is the concentration of just the free Cl- species, but the units are molal, not mol/L. If I look at the concentration of the free Cl- species in molal units in Xtplot, I get the same answer: 1.009 molal. I did not see anywhere a value equal to 1.0000000385. I'm currently out of the office teaching a short course, but I'll try to answer the rest of your question when I return next week. Regards, Brian Farrell Aqueous Solutions LLC

Hi, Do you have a script that you can post? And a thermo database with H2O2 added? Thanks, Brian Farrell Aqueous Solutions LLC

Hi Masoud, In the Basis pane of React/ SpecE8, or the Initial pane of X1t, X2t, you are specifying the composition of a fluid, or optionally a fluid-rock system, at equilibrium. So Calcite and Ca++ will in fact be related. The presence of Calcite, in equilibrium with your fluid, can constrain the Ca++ in your fluid, or perhaps the HCO3- or pH. As for the brine, you can constrain it separately in the Inlet pane. As the model moves forward, the Inlet fluid will displace the original pore fluid. The Inlet fluid will react with the original pore fluid as well as the minerals making up the rock matrix. Furthermore, if your fluid is in equilibrium with a mineral like K-feldspar, an aluminosilicate, then you will in fact have SiO2 and Al+++ in your system (even if it's a very small amount). You may not have measurements, though, and guessing unreasonable numbers might be why your model is not converging. You can always simplify your system, but you cannot leave out the Basis species which must be present in the system you've configured. To illustrate this, try opening up Rxn. Go to the command pane and type "react K-feldspar" (or select K-Feldspar from the balance reaction for section). When you move to the Results pane to balance the reaction, you'll see the Basis species you must have in your system in order for K-feldspar to exist in your system. The same principals apply in React, X1t, and X2t. Hope this helps, Brian

Hi Dirk, I'm answering your email query here. If you want ln(Q/K) to be evaluated in your rate law script, not log(Q/K), then you won't be able to use the internal parameter "logQoverK", which is the saturation index, or log10(Q/K). Instead, you'll need to use some of the "helper functions" from Table 5.2 (GWB Reaction Modeling Guide), and the rate law script syntax described in Table 5.3. Specifically, I think you'll need the helper function QovK("Siderite") and the mathematical operator log() instead of log10(). This is kind of confusing, mostly due to the different conventions for logs commonly used by geochemists and mathematicians. So instead of writing "logQoverK" in your script, to get the natural log you'll need log(QovK("Siderite")). You mentioned the script not working when you removed both negative signs from line 90 in your rate law script. I think you only need to remove the first one (I get the model to run when I do this). This is the same as adding a minus sign to each of your rate_acid, rate_neutral, and rate_base parameters, or just multiplying rate_prec by negative 1 (I think this is where it makes most sense, based on the way you set up your model). Please give this a shot and think on it. If that doesn't take care of the problem satisfactorily for you, there are a few things you can try in general to troubleshoot kinetic reactions. Try setting the variable for dxplot to 0 (Config - Output) so that every step in the calculation shows up in Gtplot. You can then take a close look at system parameters like the Siderite dissolution rate, Q/K, etc. You can also set the variable delQ (Config - Stepping) to a number smaller than 0.1, it's default (try 0.01). This limits the change in Q, the activity product, from one step to another, thus limiting the length of a time step and helping with stability. Hope this helps, Brian