Brian Farrell

Hi Yuan, A numerical model isn’t guaranteed to handle everything that’s thrown at it. Your simulation with the built-in rate law runs to completion, but it still struggles quite a bit because the problem isn’t especially well posed. The dolomite reacts orders of magnitude faster than any of the other kinetic minerals. This forces really small time steps. And because the mineral dissolves completely, in a moving front as the model marches through time, you have incredibly high reaction rates adjacent to zero rates. Considering these substantial differences, there’s really no point in using kinetics for the dolomite. An equilibrium model should suffice. Your simulation should run faster, too, because the other kinetic minerals don’t require such small time steps. Or, if you really want to continue using kinetics, you should decrease the rate slightly. Decreasing the pre-exponential factor by even an order of magnitude lets the model run to completion. There’s essentially no difference in the output, considering the current course grid, which can accentuate differences. You could even set in your custom rate law a reasonable maximum reaction rate, and take the minimum of that and your calculated rate. The point is that there’s little difference in the output whether you use an amazingly fast or just a really fast rate for dolomite dissolution, apart from the numerical instabilities you introduce with the amazingly fast rate. It might seem wrong to disregard or tinker with the kinetics here, but keep in mind that there are always a lot of variables that go into reactive transport models. Certain conditions make the models very sensitive to some parameters and insensitive to others. The precise value you set for the diffusion coefficient, for example, won’t change your output compared to using the default value because the term is insignificant wherever groundwater flows at an appreciable rate. Sometimes you need to take a step back and determine what is most important. By the way, the values you’ve set for Courant and Xstable are far too large. Their default values, 1, are the largest they should ever be. You should only set smaller values for these parameters as needed. In this simulation the kinetics of dolomite dissolution end up limiting time steps anyway, but in general you shouldn’t ever increase those values. For more information, please see the GWB Command Reference. Hope this helps, Brian

Hi BryanW, The Piper diagram units cannot be changed. Because they’re percentages, though, it doesn’t matter whether you use meq/kg or meq/l. You can see this for yourself by plotting a ternary diagram (these accommodates any unit choice, such as mg/kg or mmol/l, not just electrical equivalent units) or making a quick Excel calculation. If you’d like, you can edit the text in PowerPoint or Illustrator. Just go to Edit > Copy/Copy As and paste special/paste the .emf or .ai file into your desired program. Can you please attach a plot with the text overlapping the legend border? Or your GSS file? It might be that certain characters in the sample name aren’t being accounted for correctly when the legend border is drawn. For now, you can also edit this in PowerPoint or Illustrator. Thanks, Brian Farrell Aqueous Solutions

Hi Yuan, I'm looking into it. Thanks for your patience. In the meantime, does the built-in rate law work for your purposes? Regards, Brian

Hi Yuan, I don't have your custom thermo dataset, but I was able to check out your script using the installed thermo.com.V8.R6+.tdat. I think the difference between your custom rate law and the buildt-in rate law for Dolomite involves the surface area term. The surface area carried in the built-in rate law is the product of a specific surface area in cm2/g (set with the "surface" keyword) and the current mass in g of the mineral (you have set the "volume%" keyword, and the program converts this to g). In your custom rate law, though, you include only the specific surface area (a value of 530). You need to multiply that by the current mass of the mineral in your custom rate law. You can use the helper functions mw() and mass() to return the mole weight and mass in moles, respectively, of any species, so your rate law might include: *530*mass("Dolomite")*mw("Dolomite")*. Alternatively, you can use the internal parameters mw and mass within a rate law to refer to that specific mineral: *530*mass*mw*. If you want to include the parameter for the specific surface area, it might look like *sparea*mass*mw*. Or, if you want to jump straight to the actual surface area, just use the "surface" internal parameter: *surface*. Unfortunately, "surface" has different meanings in the context of the built-in rate law's keywords and the internal parameters for the custom rate laws. For more information, please see Table 5.1 and Table 5.2 in the GWB Reaction Modeling Guide. Regards, Brian

Hi Erik, Yes, a precipitated mineral in X1t or X2t can later dissolve if the system changes and it becomes undersaturated. Regarding removing precipitated minerals, the flow-through configuration in React does this at each step of the calculation so that they don’t dissolve back into the fluid. The conceptual model is that you’re following the vantage point of a single packet of water flowing through an aquifer. It’s a “quasi” reactive transport approach. In X1t or X2t, though, you actually have a spatially discretized domain, so removing minerals isn’t necessary. A single packet of water won’t react with the minerals in an upstream node after it’s been transported away. You have a continuous addition of new packets of water, though, that should be able to react with the minerals they flow though. If you’d like, you can specify a custom rate law that would allow a kinetic mineral to precipitate according to one equation but dissolve according to another (perhaps set the rate to 0). That way the mineral doesn’t react, even though it’s still present. Yes, porosity reflects any changes in mineral volume due to dissolution or precipitation. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Yuan, You can report the factor limiting X2t's time steps by checking "explain steps" in the Config > Stepping dialog, and then "Follow Output" on the Results pane. This might give you a clue about why your run is failing. For more information, please see 9.40 explain_step in the GWB Command Reference. If that doesn't offer any good clues, you can attach your script and we'll try to take a look. Regards, Brian Farrell Aqueous Solutions

Hi Johan, We’ve looked into your problem. Unfortunately, there are some limitations imposed by the Windows OS. As I mentioned, you can bypass the GUI using Gwb.exe to include any number of reactants. Alternatively, you can try editing your registry to allow more user objects. The default Windows allows is 10000, but you can increase it to 18000. There’s still ultimately a limit, but you may be able to add almost twice as many reactants. The specific registry entry is: HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows NT\CurrentVersion\Windows\USERProcessHandleQuota To edit your registry, type regedit into your Seach bar, then click to Run. Navigate from HKEY_LOCAL_MACHINE down to the various subfolders until you get to Windows, then right-click on USERProcessHandleQuota and select Modify. With the “Base” set to Decimal, enter 18000 for the “Value data:”, then hit OK. You’ll need to reboot your computer. After this, please try working with your script again. You should be able to add some more reactants. There will still ultimately be a limit, though. We’re exploring options to make this work better in future releases. Please let us know if you have any more questions. Regards, Brian

Hmm, hopefully both Thermoddem datasets are equivalent. Perhaps it's only a matter of the Phreeqc dataset using O2(g) in certain places instead of O2(aq), or vice versa. I'm not sure I have more to contribute regarding the Phreeqc dataset, but if you have more questions about Act2 or the GWB in general, let me know.

Hi Jerome, You're welcome. I just downloaded the dataset from here: http://thermoddem.brgm.fr/databases/geochemists-workbench. The log K values I see don't match what you describe. Are you working with a different version of the dataset? FYI, since this is an older format dataset (the second line of the dataset has the oct94 tag, and it has a .dat extension), it opens up by default in a text editor like Notepad. If you view it there, you'll see it includes the Eh or half-cell reaction in the header tables. That reaction is 2 H2O = 4e- + O2(g) + 4H+. The reaction and more detailed information are included in the Thermo Datasets chapter of the GWB Reference Manual. The more recent dataset formats (oct 13 and jul17) end with .tdat extensions, and they'll open up by default in TEdit, the graphical thermo data editor. With either format, though, you can right-click on the file and open with another program. I say this because newer format dataset don't contain the Eh reaction. Instead, it's rebalanced as the free electron reaction. That reaction can be written in several ways, but the GWB apps including TEdit will take that old Eh reaction and balance it as e- + .25 O2(aq) + H+ = .5 H2O. The log K values from the header table are divided by -4, then. Just want to clear up any confusion there. With that complication out of the way, the dataset I downloaded has (aq) or (g) appended to the names of the aqueous and gaseous O2 and H2 species, respectively. There is a little flexibility in redox coupling reactions, such as using O2(aq), O2(g), or e- when O2(aq) is in the basis species, or H2(aq), H2(g), or e- when H2(aq) is in the basis. But in any case, the reactions should be explicit whether a gaseous or aqueous form is used. If your dataset is different, and you'd like me to take another look at your problem, please attach the dataset along with your Act2 script. Hope this helps, Brian

Hi, We’re looking into this issue. There may be a limitation to the GUI imposed by Windows, but the calculation engine itself can handle any number of reactants. In the meantime, if you run X2t by entering commands into Gwb.exe (go to the support pane of the GWB dashboard), you shouldn’t experience the problem. Regards, Brian Farrell Aqueous Solutions LLC

Hi, I think the lower stability limit is being drawn correctly. A nice feature of Act2 is that it can draw the water stability limits at a pressure of your choice, rather than assuming 1atm all the time. By default, Act2 uses pressure taken from the header tables of the thermo dataset in calculating the position of the water limits. This pressure in the Thermoddem tables is ever-so-slightly different from 1 atm. As a result, the log f H2(g) term in the equilibrium equation for the lower limit doesn’t disappear completely, as you assume. Your perceived problem is because the log K for the reaction H2(g) = 2H+ + 2e- is ~ 0 at your temperature of interest. Here are the log K values for the reaction H2(g) = 2H+ + 2e- calculated at the principal temperatures, using the Thermoddem dataset: Log K's: 0 C 0.0000 150 C 0.0008 25 C 0.0000 200 C 0.0030 60 C 0.0000 250 C 0.0071 100 C -0.0000 300 C 0.0143 You can verify this yourself by loading the dataset into Rxn and typing react H2(g) swap e- for O2(aq) or by adding the log Ks for the following reactions that appear directly in the dataset: .5 x [2 H2O = 4e- + O2(g) + 4H+, the half-cell reaction in the header section of the dataset] .5 x [O2(g) = O2(aq)] 1 x [H2(g) + .5O2(aq) = H2O] When you run a calculation at a non-principal temperature, the program fits all the data to a polynomial. At 90 C, your log K works out to -1.785E-05, a very small number. So we have a small log K term added to another really small log f H2(g), so the y intercept in Act2’s output isn’t simply the log K/2. You don’t notice the effect of the pressure term in the upper stability limit because the log K is so much larger than the log f O2(g) term. You don’t see a problem with the lower limit in the default dataset for the same reason - the log Ks for the H2(g) reaction are far enough from 0. In any case, if you set pressure directly to 1 atm in your script, the term will go away and you should get a y-intercept exactly equal to the log K/2 in the text output, though it won’t affect your plot in any noticeable way. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Peter, Thanks for bringing this problem to our attention. The fix will be available in the next maintenance release. Regards, Brian Farrell Aqueous Solutions LLC

Hi Manny, This looks like a problem from an earlier GWB12 release that has already been fixed. You can check your release from any app by going to Help > About "app...", and check for updates by going to Help > Check for Updates... The current release is 12.0.4. Please let us know if you're still having issues after updating. Regards, Brian Farrell Aqueous Solutions LLC

By the way, it sounds like you're using ExtendRun() to daisy-chain a series of steady runs (each with a different flow rate) to simulate transient flow. You can, however, set up transient flow without the ExtendRun() function. For steady flow, call the “FlowRate()” member function once, before entering the time marching loop. For transient flow, alternatively, call “FlowRate()” at the head of each pass through the loop.

Hi Johan, I was just getting ready to get back to you. The ExtendRun() member function is intended to continue a simulation after a time marching loop ends so that you can chain together additional paths and have everything within a single plot file. The system chemistry is retained from the end of the previous loop, as are the masses of any equilibrium minerals, but all reactants are kept at their original values by default. A fixed CO2 buffer, for example, will be maintained. Similarly, the amount to be added of a simple reactant or kinetic mineral reactant is reset to its original value. You might have acidic water flow into a reservoir and want to treat it by titrating limestone into it, then have another inflow of water with a second treatment round. Perhaps you’d want to use the same quantity of the same reagent, or perhaps you’d alter either the amount added or the type of reagent entirely. ExtendRun() assumes your reactants remain the same by default, but an important feature is the ability to add or remove reactants before extending the run. In the User’s Guide, there’s an example in which a simple reactant NaOH from the first stage is removed, and HCl is added for the second stage. In your case, you can retrieve the endpoint mass of Calcite and feed that into a Config() command before using ExtendRun(). In a general sense, your client might look like this: # First time marching loop # Retrieve current mass of kinetic reactant cp.Report1("mass_remaining Calcite", "kg") # Reconfigure reactants to use endpoint mass, then extend the run cp.Config("react %f kg Calcite" % cp.Report1("mass_remaining Calcite", "kg")) cp.Extend(20, “yr”) # Second time marching loop If the kinetic reactant completely dissolves away, you can alternatively reconfigure and extend with something less general, like: # Reconfigure reactants by setting mass to 0, then extend the run cp.Config("react 0 kg Calcite")) cp.Extend(20, “yr”) or # Reconfigure reactants by removing kinetic mineral, then extend the run: cp.Config("remove reactant Calcite")) cp.Extend(20, “yr”) The feature is different from extending a completed run in X1t or X2t, as described in 2.22 Running a model in the Reactive Transport Modeling Guide. In that case you’re simply continuing the run without making any adjustments except for the end time. Here, you’re allowed to make adjustments to any reactants before you trigger the ExtendRun() member function. I understand the confusion, and will discuss with the development team whether we want to make any adjustments to the documentation, or less likely, how the member function actually works. I hope you enjoy using the software. Regards, Brian

Hi Johan, React has a few special mass-transfer configurations, such as the flush, flash, and flow-through models, which are included in ChemPlugin for completeness. They might be useful in a ChemPlugin client designed to emulate React, or to behave similarly to it. One example is the mReact program described in the ChemPlugin User's Guide. In constructing a reactive transport model using ChemPlugin instances, however, these configurations would most likely not be useful. As you've observed, you explicitly specify flow rates into or out of ChemPlugin instances. We're looking into the ExtendRun() member function. Can you confirm that's what you were originally trying to use? And that you attempted to use pickup after your results looked incorrect to you? Thanks, Brian

Hello, I apologize, the pore volume report command in ChemPlugin is actually a stub of a planned future feature. It was not meant to be documented in the User's Guide. It will be available in the next release for sure. We will look into making it available in the current release as well. Once again, I apologize for the inconvenience. Regards, Brian Farrell Aqueous Solutions LLC

Hi Maria, What release of the software are you using? Can you attach your SpecE8 script so that I can take a look? Regards, Brian Farrell Aqueous Solutions LLC

Hi Coralie, Well, it’s easiest if you have a complete fluid analysis. In that case you’d just run the fluid through SpecE8 or React to calculate saturation indices for every possible mineral. That way you could figure out what you need to swap into your basis. When you’re working with limited data, as you’re doing, it’s much more difficult. You can try different assemblages of minerals you’ve observed, and at least try to rule out unlikely assemblages when you can’t get the initial system to converge. You might also be able to use some experience or consult some stability diagrams to rule out combinations of minerals that aren’t likely to be in equilibrium together under particular conditions. You can add minerals as simple reactants, in which case they’re continuously “titrated” into the system at a constant rate, or as kinetic reactants, in which case the rate of dissolution is set by a kinetic rate law. Reactant minerals in transport models are probably set as kinetic reactants more commonly, but you should make that decision yourself. For more information, please see sections 3 Tracing Reaction Paths and 3.1 Titration paths, as well as 4 Kinetic Reaction Paths, 4.1 Setting kinetic reactions, and 4.2 Kinetics of precipitation and dissolution in the GWB Reaction Modeling Guide. Sure, you can approximate reactive transport in a very simplistic way with React’s flush and flow-through configurations. The flush configuration follows the point of view of a section of porous media through which water flows. You might be interested in alteration products as a reactant fluid flows in and reacts with the existing minerals. The flow-through configuration, conversely, follows a packet of fluid as it migrates. Precipitated minerals aren’t allowed to re-dissolve into the fluid because the fluid packet travels downstream, away from the immobile minerals. You might try either of them before moving to X1t. But the very first step should be to find an initial system that you’re happy with, using Go Initial in React. Hope this helps, Brian

Hi Coralie, Go Initial is primarily for testing out your initial chemical system (what you specify on the Basis pane of SpecE8/ React or the Initial pane of X1t). It calculates the initial state before any reactants are added or boundary fluids flow in. It’s a convenient way to understand the starting point for your calculation without having to worry about the complexity that will be added during the reaction path. Since there are no minerals in the Basis pane of the React script you attached, and the fluid isn’t supersaturated with respect to any minerals, doing a Go Initial will leave you with only a fluid. The reactant minerals haven’t been added yet. There are various ways to conceptualize models like this. If you think that one or more minerals are in equilibrium with the porewater at the very start of your transport model, you should swap those into the Basis. Not all observed minerals are necessarily in equilibrium with the pore fluid, though, so only those that are truly in equilibrium with your fluid should be swapped into the basis. That’s why I suggested trying out different basis configurations and using Go Initial - so you don’t have to worry about the transport or other complexities. You could add the non-equilibrium minerals as reactants later on, after the initial system is sorted out. If, on the other hand, you think the starting point for your transport model is the water-rock system resulting from the titration of your basalt minerals into a brine (the React script you attached), then run the React script to completion (Go, not Go Initial). There should be minerals in the resulting equilibrium system, which you might transfer to the Initial pane of X1t as described previously. Hope this helps, Brian

Hi Sam, It’s not uncommon to detect a change in fluid chemistry due to mineral dissolution without observing a significant change in the mineral’s mass. You might be interested in the concept of relaxation times for groundwater and aquifer minerals, which can differ by several orders of magnitude. For more information on the subject, please see section 27.1 of Craig Bethke's Geochemical and Biogeochemical Reaction Modeling textbook. Perhaps the bigger issue in your example is how you plot the data. In your example, the kinetic mineral exists in the first few nodes only. In a plot of mineral mass along the aquifer, the default scale will range from 0 to just over the maximum mineral mass. That will tend to hide any small variations in mineral mass. You could instead make a plot of mineral mass vs. time at one of the first few nodes to see the variation in mass more easily. A delta scale in a plot of mass vs. time will be especially useful for determining how much mineral dissolved. Additionally, you might wish to plot the rate of the kinetic mineral’s dissolution vs. distance or time. Change the variable type to Reactant properties and choose “Dissolution rate, Hydroxyapatite”. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Erik, I've discussed this idea with the development team and we're looking into it as a possibility for the next release. Cheers, Brian

Hi Coralie, You’re welcome. You can drag and drop compositions from one app to another. Right-drag (right-click and hold) from the Basis pane of React to the Initial pane of X1t to transfer the initial settings. If you’d like to transfer the results of a React calculation (perhaps a supersaturated mineral precipitated) you can right-drag from the Results pane of React into the Initial pane of X1t. For more on drag and drop, please see 1.7 Drag and drop feature in the GWB Essentials Guide. Regards, Brian

Hi Erik, In X2t you can set up recirculation wells, in which fluid produced from one well can be injected into another. You cannot currently circulate fluid from X1t’s outlet to its inlet, however. For more information, please see 4.2 Wells in the GWB Reactive Transport Modeling Guide. Regards, Brian Farrell Aqueous Solutions LLC

Hi Betty, As far as I am aware specific surface area isn’t used for ion exchange calculations. That's true in the GWB, at least. You specify the cation or anion exchange capacity for an ion exchange surface in SpecE8’s Sorbing Surfaces dialog. You can specify the CEC or AEC in terms of total electrical equivalents. Or, you can specify them in terms of electrical equivalents per mass of rock, in which case the program multiplies this value by the total mass of equilibrium and kinetic minerals, as well as inert mineral mass. It’s hard to say definitively. Your website mentions first solving problems assuming no aqueous complexation of the cations (i.e. complete dissociation), and then making a second calculation including aqueous complexes (e.g. CaCl+). You should make sure you’re comparing simulations with the same assumptions, like those above. To ignore complexes, you’ll want to use the convenient Suppress feature in SpecE8. For more information please see the Suppress command in the SpecE8 chapter of the GWB Command Reference. You should also verify that the underlying thermodynamic and surface data matches the original source. At a quick glance, I think the selectivity coefficients for the ion exchange reactions you’ve specified are incorrect. The website has a table of log K values. The GWB, however, requires that selectivity coefficients for the ion exchange model be expressed in the surface dataset as linear, not logarithmic values. Please see the Preamble and Bibliography of the IonEx.sdat example installed in the GWB’s Gtdata folder for more information. Please see as well section 9.2.5.4 Ion Exchange in the Using TEdit chapter of the GWB Essentials Guide. Hope this helps, Brian Farrell Aqueous Solutions LLC