Brian Farrell

Hi Silvain, As long as you don’t specify “free” or choose a unit that implies a free constraint (e.g. activity or pH), your constraint is for the bulk concentration of a component. To test this out, fire up SpecE8 and enter HCO3- = .001 molal pH = 6 Na+ = 1 molal balance on Cl- go You can look in the text output file to see the concentrations of several carbon species: CO2(aq) has the highest concentration at .0005060 molal, then in decreasing abundance HCO3-, NaHCO3, CO3--, and NaCO3-. Add them up and you’ll find the molal concentrations will sum to .001 molal, which is the bulk constraint you supplied. Scroll down to “Original basis total moles” in the text file and you’ll see HCO3- = .001 moles, which matches your bulk constraint. You can of course look in the plot file as well. The “species concentration” variable type refers to all the individual aqueous species, and the “components in fluid” refers to bulk composition. You could alternatively use the input below and you’d get the same results. swap CO2(aq) for HCO3- CO2(aq) = .001 molal pH = 6 Na+ = 1 molal balance on Cl- go Of course, if you use units like mg/kg you’ll have to account for the mole weight of whatever species is in your basis, or use the “as” setting to specify concentration in terms of mg Carbon, mg CO2, etc. Since we calculated the complete distribution of species, we can take our results from before and test out the “free” constraint setting: swap CO2(aq) for HCO3- CO2(aq) = .0005060 free molal pH = 6 Na+ = 1 molal balance on Cl- go Running this latest version should give equivalent results. For your second question: Act2 calculates the simple type of diagram that geochemists have traditionally drawn by hand. By design, a number of simplifications make the calculation straightforward (but still laborious). For example, you can only have logarithmic axes (log activity, pH, pe or Eh, etc.). There’s no mass balance. You have to work in terms of activity, rather than concentration. The diagrams are in many cases fairly similar, though. If you want a general picture of how chemistry of a particular system works, a traditional calculation might be fine. Certain applications might demand a more rigorous solution. Beyond recreating activity or redox-pH diagrams, though, Phase2 calculates other diagrams that Act2 cannot even contemplate. You can include surface complexes in your calculations, for example. Or, you can diagram how various properties (solubilities, saturation indices, gas pressures, and so on) vary across the diagram using color maps or contours. You can account for isotope fractionation or kinetic reactions. You can plot assemblages (combinations) of stable minerals under different geochemical conditions. Basically, it has all the capabilities of React, so it's almost endlessly configurable. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Polly, For the two-layer surface complexation dataset, you need to specify a site density or densities for each sorbing mineral. If you don’t, the program can’t account for the existence of the surface. If Ferrihydrite contains both of the sites you defined, >(w1)SOH and >(w2)SOH, your entry might look like this: Ferrihydrite surface area= 600.0000 m2/g 2 sorption sites >(w1)SOH site density= .0050 mol/mol mineral >(w2)SOH site density= .2000 mol/mol mineral You should, of course, supply values appropriate for your sorbing mineral. In GWB14, by the way, you can alternatively specify site density in sites/nm2, as you sometimes see in the literature. You should also remove the entry for >(w) in the surface species section. I’m not sure if it was intentionally put in the dataset like this, or it was a stub of an entry that you didn’t finish, but there’s no reaction, stability, or mole weight and that’s causing a problem. Finally, the surface dataset has a field where the user specifies the thermo dataset to use with it. The aqueous species that are included in the surface reactions, as well as the sorbing minerals, are drawn from the thermo dataset you specify. If you plan to use this surface dataset with your custom thermo dataset, it might be best to specify that custom thermo dataset within the surface dataset. Note that the sorbing mineral you’ve chosen, Ferrihydite, is not in your custom thermo dataset, so you may need to make further modifications to either your thermo or surface dataset. You’ll probably have to make these modifications in a text editor, like Notepad, since TEdit can’t open the unproperly formatted dataset. As for your thermo dataset, I don’t think you should include SiO2, H4SiO4, and Si(OH)4(aq) as separate basis entries. They all represent essentially the same thing. Pick one and write all reactions in terms of that species. Hope this helps, Brian Farrell Aqueous Solutions

Hi Thomas, I hope you're doing well. I happened to come across this old post and thought you might be interested to know about Phase2, an app introduced with GWB12. The program essentially traces a stacked series of reaction paths, as you'd run in React, to traverse two geochemical variables of interest. You can set up a diagram with sliding log f O2(g) and fixed pH along the y axis, then sliding pH with fixed f O2(g) along the x axis. The basis fluid is defined in terms of total concentrations, as in React, and mass is conserved throughout the calculation. You can also titrate a species into a fluid initially devoid of it to consider a range of total concentrations. By titrating SO4-- in log steps you can make a diagram much like you've envisioned, with the sulfur speciation depending on the y axis variable, log f O2(g). You can render the calculation results in various types of 2D diagrams or in horizontal or vertical cross-sections through the diagram. For the 2D diagrams, you can plot "true predominance" for any basis species or element (the species accounting for the most mass predominates, not the species with the highest activity), mineral assemblage diagrams (which show every stable mineral or combinations of minerals), and render any variable as a color map, mask, or contour. In a log f O2-pH diagram, for example, you can diagram the stable iron minerals under various conditions and contour the concentration of dissolved Fe. Please visit GWB.com/phase2.php to learn more. I'm happy to send a demo if you'd like to try it out. Cheers, Brian Farrell Aqueous Solutions LLC

Hi Polly, Thanks for providing the thermo dataset. I'm taking a look to see if I can offer any suggestions. Regards, Brian

Johan, Some time ago you tried to report the pore volumes displaced from a ChemPlugin instance in GWB12. I’m writing to let you know that GWB14 is now available, and ChemPlugin instances now plot pore volumes displaced and have the value available in the report command. Additionally, React now plots pore volumes displaced from flush and flash models. I hope you enjoy using the software. Cheers, Brian

Dear GWB users, GWB14 is here! A new compute engine makes it a surface chemistry powerhouse. Plus, our latest release features user equations in GSS datasheets, high-temperature Pitzer coefficients, cluster computing, and flexible input for thermo data, just to start. Visit our GWB14 page to learn more. GWB subscribers upgrade automatically. Can’t wait another minute? Click “Check for updates” on your GWB dashboard. Not subscribing yet? Take advantage of one of our flexible plans—fixed or floating, 3 months to 3 years! Or, convert that old GWB license to a modern subscription and a free kick-off period is on us. Contact us today for a quote, or go shopping at our online store and you’ll be up and running in minutes. Sincerely, Brian Farrell Aqueous Solutions LLC

Hi Sanjoy, I’m writing to let you know about a new feature in the GWB14 release. You can control species loading in a calculation by specifying a temperature range over which thermodynamic data should be available. In your case, many species have thermo data available at 25 C, but not other temperatures, so your calculations at 25 C and 26 C loaded a different set of species, which caused your results to differ more than expected. By specifying a temperature range, say 10-150 C, you can ensure that any calculations within that range load the same exact set of species. It’s basically like automatically suppressing a set of species with limited thermo data, simply by setting a temperature range. You can read about the span command in the GWB Command Reference. The new release has many more new features as well. Please let us know if you’re interested in trying GWB14. I’m happy to send you a demo. Cheers, Brian

Hi Bill, I’m writing to let you know about a new feature available in GWB14. You can control species loading in a calculation by specifying a temperature range over which thermodynamic data should be available. You can additionally force the program to always get log K values by evaluating a polynomial, even when the calculation is isothermal at a principal temperature. In this case, the program would normally use the value directly from the thermo dataset. The difference is small, but it’s why the endpoint of your polythermal path and your isothermal calculation post-pickup were slightly different. You can read about the feature in section 6.92 span in the GWB Command Reference. The new release has many more new features as well. Please let us know if you’re interested in trying GWB14. I’m happy to send you a demo. Cheers, Brian

Hi Mauricio, I’m writing to let you know that GWB 14 is now available. The new release includes support for several different polydentate surface complexation formalisms. You had to convert the log K for a bidentate reaction in your surface dataset to be consistent with the method in GWB9, but now you can choose one of four methods consistent with your log Ks and specify that in your surface dataset. React will read the dataset and evaluate mass action laws according to the convention you specified. In your case, you could set the stoichiometric approach in the dataset’s header and use the original literature reference’s log K value of 4.6. The calculation will reproduce figure 3A in your plot. However, this approach is not satisfactory, as described in the reference above, and it will be in error at other concentrations of the sorbing mineral. The better approach is to convert the log K, as you did before, to use the mole fraction approach, which is entitled Hiemstra-VanRiemsdjik in GWB14. This is the default approach in GWB14, so you can leave the method unspecified in your surface dataset, but to be clear it’s best to set the method in the header. With this approach, changing the amount of the sorbing mineral should still yield reasonable results. Please let us know if you’re interested in trying GWB14. I’m happy to send you a demo. You can read about the feature in section 2.5.8 Polydentate sorption in the GWB Essentials Guide. Regards, Brian

Hi again Dave, Upon closer inspection, you can include a user-defined analyte in a radial plot as long as you put it in the "Components in fluid" category, since this is where the normal basis species plot in Gtplot. If you instead choose a category like "Chemical parameters" and pick Concentration for the dimension, it won't be recognized by the water chemistry plots. Still, if you want to do any thermodynamic calculations including Be, you'll need to use a dataset that includes the relevant reactions. Hope this helps, Brian

Hi Dave, The radial diagram can plot uncharged species when you choose units like mg/l, as you've done. However, the special plots, including the radial plot, only diagram basis and redox species in your spreadsheet. They do not include user-defined analytes. User analytes can only be plotted in xy plots. You might try using a different thermo dataset, like thermo.com.V8.R6+.tdat. That has Be++ as a basis species. Or, you can modify the thermo dataset you're using by adding the element Be and any species that you need, then load that dataset into GSS. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Polly, You attached a .lnk file (a shortcut), not your actual GWB thermo dataset, which should be a .tdat file. For one thing, you should check to ensure that your dataset and your spreadsheet are using the same exact set of reactions and log K values. I looked at how you set up your script, but I couldn’t run it because I don’t have the thermo dataset. In any case, it doesn’t look like you took my advice. Can you explain why you’ve set your calculation up with trace quantities of all your basis entries, then titrated large amounts of all your minerals of interest? As I mentioned before, a common way to find solubility is to swap your mineral of interest into the basis, then find the total dissolved concentration of the constituent of interest. In my example to find the solubility of Kaolinite, the mineral was swapped into the basis and the concentration of the Al+++ component in the fluid was determined. The sliding pH part is a convenient method to make a type of solubility diagram, to see how the solubility changes with pH, but if you don’t care about that part, you can certainly omit it to look at solubility under a single condition – your groundwater of interest. I would repeat the calculation for each mineral of interest. Swap a mineral like Pd(OH)2(s) into the basis in place of Pd++, but use the measured concentrations for all of your other ions. Find the total concentration of Pd++ in the fluid in equilibrium with your Pd(OH2(s). Then move on to SnO2(am), swapping that for Sn++++ but leaving all the other ions, like Pd++, at their measured values. It’s likely that your spreadsheet is performing a series of calculations, so I think you need to do the same. As I also mentioned, you need to look at the output and ensure that your mineral of interest is still present in equilibrium with the fluid. If you swap amorphous silica into the basis to find its solubility, React will precipitate quartz, which is more thermodynamically stable, unless you suppress it or disable precipitation entirely. If you don’t, you’d actually be finding the solubility of quartz, which is lower than the amorphous silica, so it’s important to check your output carefully and adjust your input as necessary. Regards, Brian Farrell Aqueous Solutions LLC

Hi, Thanks for catching that. This issue has been fixed. It will be available in the next maintenance release. Regards, Brian Farrell Aqueous Solutions

Hi Silvain, The reaction for Pyrite in the LLNL dataset, thermo.tdat, is shown below: Pyrite + H2O = Fe++ + 1.75 HS- + .25 SO4-- + .25 H+ You’ll notice that the reaction is written in terms of the basis species SO4-- and the redox species HS- because the oxidation state of sulfur in pyrite is between that of SO4-- and HS-. When the reaction gets loaded into an app like Rxn, SpecE8, or GSS, though, it’s rebalanced to be in terms of the basis species only: Pyrite + H2O + 3.5 O2(aq) = Fe++ + 2 SO4-- + 2 H+ This is because each of the apps assume that redox coupling reactions are enabled by default. In any case, to be able to calculate a saturation index for a mineral like pyrite, you need a way to constrain every species in the reaction. So in addition to adding Fe++ and SO4-- to your spreadsheet, you need to add O2(aq) (or Eh or pe), as well as pH. In this case, the concentration you set for the SO4-- component should represent all the sulfur in the system. The program uses the oxidation state you supply to find the equilibrium distribution of mass between sulfate and sulfide species. You can alternatively decouple HS- from SO4-- to set up a disequilibrium model. This is done from the Config > Redox Couples dialog in Rxn and SpecE8, and from Data > Redox Couples in GSS. Alternatively, simply adding the redox species HS- to your spreadsheet triggers HS- to be decoupled from SO4--. In this case, you don’t need to supply a measure of oxidation state, but you do need to provide separate measurements for the SO4-- and HS- entries. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Frank and Andrew, Act2 uses a simple analytical method to calculate equilibrium lines and assemble them into Pourbaix diagrams, so there's no difficulty in creating a diagram that spans large Eh and pH ranges. If you wanted, you could set the Eh and pH to span the range -100 to 100, and it would draw the diagram. With Phase2, however, you're setting up numerical reaction path models. In a model of the aqueous phase, working far outside the stability limits of water (i.e., the bottom left and top right corners of an Eh-pH diagram) can be difficult. To investigate your calculations, try overlaying some contour plots on your 2D diagram, or better yet, just look at a simple xy plot of the diagram’s left edge. You can select the left-most vertical cross-section through the entire diagram, use the “Go Y” option in Phase2 to calculate only the left edge of the diagram, or set up a React script to reproduce the left edge of the calculation. Whatever you choose, start out by plotting the Mass of Solution or Fluid volume vs. Eh. You'll see they're initially enormous. If you make a plot of species concentration, you’ll see it’s due to the extremely high concentration of H2(aq) that would exist under your initial conditions, which are far outside the stability range of water. Next, try plotting the concentration of your metal component of interest (U++++ or Al+++ in the system) and you should see that your initial conditions are honored. In the U example, U++++ was set to 1e-10 mol/l, but under the initial conditions that’s equivalent to .19 molal! That’s why so many minerals are stable throughout the diagram, compared to the Act2 calculation. As for the Al example, Al+++ was set to 1e-5 mol, and there’s still 1 kg of solvent in the system, so the initial Al+++ molal concentration isn’t skewed the way the U++++ was. You didn’t attach an Act2 script, Andrew, but presumably it’s not as different as the U calculations were. As for the failures you’re encountering, unfortunately it’s likely to happen when you’re so far outside water’s stability region, as in the corners of an Eh-pH diagram with a traditional range. In a log f O2(g)-pH diagram, by contrast, you don’t have to include such extreme conditions. Since a large portion of the area in an Eh-pH diagram is likely to be masked, log f O2(g)-pH diagrams can be a nice alternative in that they fill more of the plot area with useful information, rather than blank space. If you need to know the Eh, you can always plot contours on the log f O2(g) diagram. Whatever diagram you choose, it’s common to use a narrower range in these types of diagrams. I think when you do that you’ll find that you can reproduce an Act2 diagram much more closely. There will still be differences, of course. Boundary lines can be somewhat curved, rather than straight, reflecting the fact that Phase2 is solving a complete multicomponent calculation at each point in the grid instead of drawing equilibrium lines. It includes mass balance and activity calculations, unlike Act2. And whereas an Act2 diagram shows the stability of minerals and predominance of aqueous species (in terms of highest activity), Phase2/P2plot’s predominance map strictly shows the species accounting for the most mass at each point in the diagram. In other words, a mineral’s stability field (which you can see with an assemblage map) can be slightly larger than the area in which that mineral predominates all other species. Sorry for the delay in responding. It looks like we missed the original post. Hope this helps, Brian Farrell Aqueous Solutions

Hi Erik, I'm happy to hear you figured it out after more troubleshooting. Regards, Brian

Hello, Sure, there are many calculations you can do involving dissolved organics: You can set the total concentration of a dissolved organic and determine its speciation in solution. For example, you can look at a ligand like EDTA or acetate to see how much is present in free form, in various protonated or deprotonated forms, and complexed with various metals, and see how the distribution of mass changes as a function of pH, temperature, concentration, etc.. You can create Eh-pH or Pourbaix diagrams to see under what conditions a dissolved organic is stable and when it would tend to mineralize to inorganic carbon, or form various other organic compounds. You can also set up disequilibrium simulations in which a dissolved organic is oxidized or reduced with time, at a rate determined by a kinetic rate law you specify. You can account for abiotic reactions in solution, catalysis on mineral surfaces, biodegradation (a simple model of enzymatic catalysis), or microbial catabolism coupled to growth and decay. The program does not currently consider non-aqueous phase organics, however. You can’t find the solubility of an organic in water the same way you would a mineral, unless you create a fictive “organic mineral” with a log K to control the solubility. To set up a simulations involving organics, you’ll first need to load a thermodynamic dataset that contains your organics of interest. The GWB’s default dataset, thermo.tdat, does not have any data for hexane. However, the thermo.com.V8.R6+.tdat dataset, a later release from Lawrence Livermore National Lab, does contain data for one hexane isomer, n-Hexane(aq), along with many other organic species. Many organic species are found in the redox coupling reactions section of the dataset, in a reaction with inorganic carbon. Please note that there is only one reaction for n-hexane, though. I’m not familiar with the chemistry of hexane, but if you expect it to form complexes, you might need to modify one of the existing thermo dataset. You can do so in a text editor like Notepad, or in the GWB’s graphical thermo data editor, the TEdit app. Next, you’ll have to decide in the GWB app (e.g. React) what redox reactions should be in equilibrium. To look at speciation of acetate as a function of oxidation state and pH, for example, assuming equilibrium with inorganic carbon, you’d load thermo.com.V8.R6+.tdat into React, add basis entries HCO3-, O2(aq), H+, and constrain those values. In a redox equilibrium model, the HCO3- component includes inorganic carbon species as well as various organics. The oxidation state you specify determines how much is present at equilibrium as inorganic carbon, how much as acetate, methane, ethanol, benzene, and so on (all reactions are considered by default). You can disable one or more redox coupling reactions if you’d like. If you disable all coupling reactions involving carbon except for the one involving acetate, you can look at the distribution of mass between inorganic carbon and acetate. And you decouple all coupling reactions. In that case, you can add acetate directly, without needing to specify the oxidation state or amount of inorganic carbon. Then, you can look at speciation as a function of T, pH, concentration, etc., as mentioned above. For more information, please see sections 2.1 Configuring a calculation, 2.4 Redox couples, 7.2 Equilibrium models, and 7.3 Redox disequilibrium in the GWB Essentials Guide. Please see as well section 4.6 Kinetics of redox reactions in the GWB Reaction Modeling Guide, and the Thermo datasets chapter in the GWB Reference Manual. Hope this helps, Brian Farrell Aqueous Solutions

Hi Drew, I'm happy to hear that you have Phase2 output up and running with Igor Pro. Thanks for providing your code. Most likely someone in the geochemistry community will find it useful. Cheers, Brian

Hi Erik, Can you please attach your script? Thanks, Brian Farrell Aqueous Solutions

Hi Drew, With sliding temperature, pH, fugacity, etc., the step size is consistent. Concentration is a little more complex because you’re setting up a titration path to add aliquots of mass at each step to change the concentration. Still, in most cases it’s pretty consistent. Note for either axis you can choose whether to take linear steps or log steps. If you realize your diagram would best be plotted with a log axis, like those you’ve shown, you should take log steps so that the grid is more or less equally sized with the log scale. Otherwise, the points would be equally spaced when you plot Fe(II) concentration on a linear scale, but not the log scale. Please see section 7.4 Linear and log stepping in the GWB Reaction Modeling Guide, as well as example 7.9 Mineral solubility. You can make a cross-section plot in P2plot to see the spacing. In your case, choose “Across scanning paths” for Orientation, then plot Fe++ in system vs. Rxn progress(y), using log axes for both. You can show the data markers (Format menu > Quick Toggle… > check data markers) to see whether they're equally spaced. Sure, you could look in P2plot’s exported table for the first 0 value for Q/K of the mineral, or the first non-zero value for mineral mass. Or, you could make an assemblage or predominance diagram in P2plot to see how it should look, then use React to reproduce segments of the plot along the boundaries. For example, set the fluid in equilibrium with Fe(OH)2(s) and slide the pH to see how the solubility changes (plot composition of Fe++ component in the fluid). You can export that plot to a table to see the x and y values. Presumably this is what you did with Visual Minteq? Regards, Brian

Hi Drew, With the 2D diagrams in P2plot (and map view plots in Xtplot), you can assign a variable to be rendered across the 2D grid as a color map or contour plot. The “Copy as text/spreadsheet” options will output a table of numerical values for the variables chosen. The most common usage is to perform some additional calculation with the data in Excel, but you can also pass the data to other programs for plotting. The predominance area diagram simply displays the predominant species at each point, and the mineral assemblage diagram the names of one or more minerals that exist at each point, so I’m not sure exporting the name or names at each point would be useful. You could, however, make a color map of your mineral of interest. In that case, you’d be able to produce a table with empty spaces where the mineral doesn’t exist, and some finite value where it does. Would that be useful? As for outputting the position of bounding lines, we don’t currently have any plans to do so. On the other hand, for simple modifications like adding a 0 before .1, you could copy the plot as an enhanced metafile and paste into PowerPoint, where you can ungroup the image and edit the text boxes as you like. Or, you can copy as an .ai file and work with the plot in Adobe Illustrator the same way. I’m not sure if Igor Pro accepts those file types, but you can also save the image (File > Save Image…) in a few other formats, such as .svg or .ps. Is this helpful at all? Regards, Brian Farrell Aqueous Solutions

Dear GWB users, Due to popular demand, we're returning to Chicago to present our famous Reactive Transport Modeling workshop. Register by October 7 to save $100! Reactive Transport Modeling November 4-5, 2019 Chicago, Illinois USA Early registration ends October 7 Two days of hands-on training will make you a GWB Pro. And while you're in town, take in the Second City's architecture, lakefront, parks, theater, art, shopping, dining, and so much more! Visit our Chicago workshop page to register now. Regards, Brian Farrell Aqueous Solutions

Hi Abdulaziz, Typically you would determine a dispersivity value by adding a non-reactive tracer and fitting your model to a breakthrough curve. A dispersivity of 3 cm in a 5 cm long column seems awfully big. I think it’s possible that the montmorillonite is dissolving, too. If your objective is to create a model that fits your pH and Ca++ measurements at the outlet, you might need to include montmorillonite in addition to the calcite in your model. Please note it may be difficult to fit such a model with only pH and Ca++ concentrations. Regards, Brian

Hi Abdulaziz, Of the kinetic parameters, I’d focus on adjusting the rate constant. I’m not sure you have any reason to change the specific surface area so much from your original value. There’s no point in adjusting the diffusion constant in your model. Perhaps you need to control the CO2 fugacity more tightly. If it’s possible that CO2 is in contact with the fluid in the column, you might need to set a fixed fugacity path, if exchange is rapid, or a kinetic gas transfer model if the exchange is slower. Or maybe there are other minerals that are present in small amounts yet are highly reactive? Regards, Brian

Hi Abdulaziz, You might want to check the cross-sectional area of your column. You said it has a radius of 2.5 cm, right? So, the area of the inlet face should be 2.5 * 2.5 * 3.14 = 19.625 cm2. If that’s the case, I’d put the delta y and delta z values at sqrt 19.625, or 4.43 cm. After that, you might recheck your flowrate and the number of pore volumes displaced to see if they make sense. You’ve set the log permeability to be constant by setting A (the porosity multiplier) to 0. Please note, however, that the permeability isn’t even used when you set specific discharge directly. It’s only used when you set a head drop or potential drop across the domain. There’s no guarantee that your mineral has the same specific surface area as the limestone in the study, but it’s probably a good estimate in the absence of your own measurements. As I mentioned previously, your fluid is traveling extremely quickly, so the effects of diffusion are negligible. Since you wish to find a value for the rate constant that fits your data, try adjusting that value so you can match your Ca++ and pH measurements. Regards, Brian