Jump to content
Geochemist's Workbench Support Forum

Brian Farrell

  • Content count

  • Joined

  • Last visited

  • Days Won


Brian Farrell last won the day on August 16

Brian Farrell had the most liked content!

Community Reputation

10 Good

About Brian Farrell

  • Rank
    Advanced Member

Profile Information

  • Gender

Recent Profile Visitors

The recent visitors block is disabled and is not being shown to other users.

  1. Brian Farrell

    Pourbaix diagram using Phase2

    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
  2. Brian Farrell

    Kinetic Mineral Custom Rate Law

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

    Input organics

    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
  4. Brian Farrell

    P2Plot data output

    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
  5. Brian Farrell

    Kinetic Mineral Custom Rate Law

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

    P2Plot data output

    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
  7. Brian Farrell

    P2Plot data output

    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
  8. 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
  9. 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
  10. 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
  11. 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
  12. Hi Abdulaziz, In a reactive transport model, it's best not to constrain the amount of a mineral using absolute units, like g or cm3, because changing your discretization will change the system. Relative units are much better. Since calcite is the only mineral, and you know the porosity (.16), set the mineral abundance to .84 vol. fract. You'll definitely need a lot more than three nodes to get an accurate solution. Start with 20 or 50, perhaps, to get a basic understanding of your system without costing too much in terms of computational effort, then increase the resolution to 100 or several hundred, perhaps. For a discussion, please see the Numerical Dispersion lesson in the GWB Online Academy. It doesn't sounds like the limestone was in contact with the initial dilute fluid for very long, so I'm not sure it's appropriate to assume the initial pore fluid is in equilibrium with calcite. In other words, I wouldn't swap calcite into the basis on the Initial pane. Your Aziz model.x1t is ok in this regard, but I'd put smaller values than .1 mg/kg for the various solutes, unless those values are actually known. I would add Calcite as a reactant, specifying it's abundance as described above. I think that takes care of your questions. You should not set a value for "reactants times". Looking at your inlet fluid, I would set smaller, but nonzero, values for the Ca++ and HCO3- in solution. As for the Cl- component, I think it's correct to set it as the charge balancing ion, but the value you set doesn't actually matter because the program adjusts the Cl- component to balance out the other ions in solution. The "as HCl" setting is unnecessary here. Finally, I think you'll want to double-check your width and height settings to ensure you have the same cross-section area as your cylindrical column. You also might want to double-check your flow rate. Regards, Brian
  13. Brian Farrell

    ChemPlugin Time Marching Loop

    Hi Erik, If you were writing a simple model that considered only advective transport of a non-reacting solute, you might figure the limiting time step from only the Courant condition. And if the velocity didn’t change with time, you could hard-code that value into your time marching loop. When you construct a reactive transport model that accounts for a variety of other processes, though, such as diffusion, heat transfer, and kinetic reactions, you need to account for the stability of the solution to each of their governing equations. ChemPlugin’s ReportTimeStep() member function in fact does this. The X1t and X2t programs use similar logic, which is described in section 2.20 Time marching in the GWB Reactive Transport Modeling Guide. In your program, it might be ok to set your own time step without querying ChemPlugin, as long as it’s smaller than the largest possible time step that ChemPlugin would allow. But really, you should always have ChemPlugin instances report the stable time step. You can compare with your own desired time step, if you’d like, then use the lowest of the values. Based on a previous conversation, it sounds like you might be worried that kinetic reactions are making the simulation take too long. Simply ignoring the reported stability limits would be a bad idea. Instead, you should consider whether you really need so many kinetic reactions, especially really, really fast kinetic reactions. Use kinetics for the more slowly reacting minerals, and equilibrium for the others. Regards, Brian Farrell Aqueous Solutions
  14. Hi Jeonghwan Hwang Aqueous species, such as the free H+ or OH- ions, must have positive concentrations. Thermodynamic components (the original basis that you see in the text output file), however, can have negative concentrations. For a discussion, please see section 3.2.2 Components with negative masses in Craig Bethke’s Geochemical and Biogeochemical Reaction Modeling text. FYI, the pH is listed at the top of the block of output for each step in the calculation. Keep in mind that the theory of chemical equilibrium does not include any information about time. You’ve set up an equilibrium model, so the time span you’ve set is meaningless. If you set time-dependent processes, such as a rate of simple reactant addition, an internal heat source, a dual porosity model with diffusion into a stagnant zone, or a kinetic rate law, however, the time span will actually have meaning and be used in the calculation. I’m not sure how important it is to your calculation, but if you want to consider how a mineral like montmorillonite dissolves with time, you need to set a kinetic rate law. You might be able to find kinetic rate laws and parameters in the literature, but you should be careful to ensure that they are appropriate for the system you’re modeling. For more information, please see Chapter 16 Kinetics of dissolution and precipitation in the GBRM text, as well as Chapter 4 Kinetic Reaction Paths in the GWB Reaction Modeling Guide. Hope this helps, Brian Farrell Aqueous Solutions
  15. Brian Farrell


    You can start with a delta plot showing all the species in the fluid, then one-by-one hide the species with the largest concentration changes that aren’t due to H+ transfer. For example, the concentration of CaCl+ decreases more than anything else, and the free Ca++ and Cl- ions increase almost stoichiometrically. It's clear there's no H+ transfer involved, so right-click on the curve corresponding to each species and select "Hide this variable". After a few more similar steps, you’ll see that the species we’ve highlighted have the largest change in concentration and thus best explain the pH change in our chemical model. I’m glad to hear you figured out how to overlay the plots outside our software. For anyone else reading, there are a couple different options for using MS Office products to superimpose results from separate calculations together. One method is to export the numerical values from your plots into Excel (Edit > Copy As > Spreadsheet), then plot them together there. A second method is to copy each plot as an enhanced metafile (Edit > Copy As > Enhanced Metafile) into its own slide in PowerPoint. You should make sure that the axis range spans the entire data range of interest (in your case, from 0-12 mg/kg), so that the two plots fit together correctly. Next, ungroup each image, then copy one of the curves and paste it into the other slide. For more information, please see the “How do I retrieve numerical data from my plots?” and “How do I overlay my diagrams?” slideshows in the Using GWB section of https://www.gwb.com/tutorials.php. Regards, Brian