Brian Farrell

Hi Maki, Currently the only way to specify the flow field is with discharge across node boundaries, not discharge within the nodal blocks themselves. Entering values for the nodes themselves will give you incorrect results. The best way to input your discharge data is to prepare two tab delimited spreadsheets (for x and y directions), which you can edit in Excel. As for the node-by-node editor, it should include an extra column for x-discharge and an extra row for y-discharge. We'll get this fixed in the next maintenance release. Regards, Brian Farrell Aqueous Solutions LLC

Hi Masoud, Try reading section 7 (Using SpecE8) in the GWB Essentials Guide for info about setting up your equilibrium system and performing Basis swaps. I think using the GUI will be very helpful for you to understand why the additional constraints you supply are unnecessary. Perhaps you could even paste in one line at a time using the command pane, then take a look at the GUI to see the effect of that command. Once you get used to setting up your chemical system in SpecE8, try taking a look at the steam flood example from section 3.11 of the GWB Reactive Transport Modeling Guide. In that example several minerals are swapped into the Basis in the Initial system, but not in the Inlet fluid. Hopefully this helps you to understand the process a little more. Regards, Brian

Hi Masoud, The problem is that you swap minerals into the Basis, then unswap them for the original Basis entries a few lines down. Although you're not using the "unswap" command, by setting the concentration of the original Basis entries this is implied. The input below is for the Initial system. scope = initial swap Quartz for SiO2(aq) swap Calcite for Ca++ swap Pyrite for H+ swap K-feldspar for K+ swap Smectite-high-Fe-Mg for Na+ 27 vol% Quartz 2 vol% K-feldspar 2 vol% Pyrite 40 vol% Calcite 2 vol% Smectite-high-Fe-Mg pH = 6.35 HCO3- = 976 mg/kg Ca++ = 208 mg/kg Mg++ = 53 mg/kg Na+ = 82394 mg/kg K+ = 832 mg/kg Cl- = 128488 mg/kg SO4-- = 39 mg/kg Fe++ = 3 mg/kg Quartz is preserved in your initial system because you specify its mass, but do not specify the concentration of SiO2(aq) later on in the Initial system. Calcite, which you swap in for Ca++, is a different story. First you specify the mass of Calcite, then a few lines down the concentration of dissolved Ca++, which overrides your previous entry for Calcite. You can only provide one constraint for each Basis entry. The constraints for Ca++, pH, K+, and Na+ in the Initial system need to be removed. Once you do this, you'll receive an error message and find that you'll need to add a few species to the Basis. For example, since you're adding Pyrite, you'll need to add O2(aq) to the Basis (you can of course swap it out after) and the K-feldspar will require you to add Al+++. Does this make sense? Regards, Brian

Hi Dirk, I don't have access to those publications, but I made a spreadsheet (sent via email) which will calculate the various rate terms (rate_acid,...) based on the chemistry (QoverK, H+ activity, surface area) at a particular timepoint. I then evaluated this at a few timesteps (0.37 years, for example) after running the model as is to the point of crashing (try setting end time to 0.5 years). According to your rate law, the "rate_prec" is negative when siderite is supersaturated (QoverK ~13.94). Since this is a rate law for dissolution, a negative rate (precipitation) makes sense given the supersaturation. Evaluating the total rate for QoverK >=1, I find that the rate is positive. Do you have a negative sign where you shouldn't, perhaps? Not sure if this is correct, but I removed the negative sign from the line below and afterwards the model ran to 100 years without a problem. In my test, I changed 90: rate = -(surface * lamda * rate_prec *(QoverK - 1)^2) - k_n * exp(-(gamma_n * expterm)) to 90: rate = (surface * lamda * rate_prec *(QoverK - 1)^2) - k_n * exp(-(gamma_n * expterm)) With this change to your rate law script, the run looks fairly similar, at least qualitatively, to one in which siderite is allowed to precipitate as an equilibrium mineral. Over 100 years that may not be unreasonable. Take a look and see whether this seems correct. Regards, Brian

Hi Dirk, I'm not familiar with BCF theory, but based on your rate law, it will be tough for any mineral to form without initially being present (you set the initial mass of Siderite to 0). Both your dissolution (QoverK <=1) and precipitation (QoverK >=1) rates are multiplied by the "surface" command - the total surface area of Siderite. The precipitation rate includes an additional term, - k_n * exp(-(gamma_n * expterm)), which is added to the portion multiplied by Siderite surface area. It appears this second term will always be an extremely small number, so I'm not sure how you can get any significant precipitation. Are you sure the rate law script is set up correctly? Have you tested a simpler example, or replicated a published experiment or implementation of the BCF theory? Regards, Brian

Hi Masoud, The procedure you've described sounds correct. Could you please attach your script so that I can take a look? Thanks, Brian

Hi Jay, The problem is that you're using kinetics to model an incredibly fast reaction over a long period of time. This forces the program to take tiny timesteps, eventually causing it to crash. Try doing a "Go Initial" for this example, and take a look at the Calcite dissolution rate in the first node - it's huge. If your formation is mostly Calcite, set it as an equilibrium mineral (swap for Ca++ or HCO3-, for example) in the Initial pane. You should look into what's called the Damkohler number (Da), which represents the rate at which a component reacts chemically relative to the rate at which it is transported by advection. In your example you've specified a very fast reaction and a very slow flowrate (0.0001 m head drop across a 1000 m domain). This results in a very high Da, which indicates that you should stick with a local equilibrium approach. Try checking out Chapter 21 and 27 in the Geochemical and Biogeochemical Reaction Modeling textbook (and the references therein) for more info on Da and equilibrium vs. kinetics. Regards, Brian

Hi Johan, We've found a small error in the rate-law parser when reading pH. This will be fixed in the next maintenance release. In the meantime, I think you can use the rate-law scrip below in its place. Be sure to verify that it is behaving as you would expect it to (perhaps run a batch model). IF activity("H+") < 1e-3 THEN 20 ELSE 40 20: rate = 6.4565e-13 * surface * (1 - Q/K) * (activity("H+")^-0.11) * (activity("O2(aq)")^0.5) GOTO 60 40: rate = rate_con * surface * (1 - Q/K) * (activity("Fe+++")^0.93) * (activity("Fe++")^-0.4) 60: RETURN rate Apologies for the inconvenience. Regards, Brian

Hi Dirk, Do you have a React file to use with this rate-law script? Thanks, Brian Farrell Aqueous Solutions LLC

Hi, I forgot to add H+ to the list of Basis species you'll need. You'll need that for other species or minerals to be considered (to show up in the "Species and minerals in main system" section) when Illite(IMt2) is chosen as the main species. Whether they actually plot on your diagram depends on the chemical conditions you specify. Changing the temperature, or the activity of various Basis species, can alter your diagram. Hope this helps, Brian

Hi Jay, Could you try attaching the model again? Thanks, Brian Farrell Aqueous Solutions LLC

Hi mchesh9973, Have you tried adding additional Basis species to your system using the "in the presence of field"? I tried adding Na+, Mg++, Fe++, and Fe+++ (after decoupling the Fe++/Fe+++ redox pair, or you could add Fe+++ and O2(aq) instead) and was able to generate a stability diagram with more than just Illite(IMt2). Try taking a look at the printed output (click "View Results") to get a better idea of how Act2 works. Pay attention to the "Basis species" section (which will include your main species, water, species on the axes, and the ligand species) and the "Species and minerals in main system" section. This shows the different species or minerals that the main species (Illite(IMt2) in your case) can become based on the Basis you've specified. Hope this helps, Brian Farrell Aqueous Solutions LLC

Hi Catriona, Not sure yet when the maintenance release will come out. Since the flow-through model is basically an approximation/ simplification of a reactive transport problem, you might try out your problem in X1t. Regards, Brian

Hi Johan, I've made some progress testing the rate law scripts in simpler examples, but I'm not sure yet why your X1t example isn't working. Regards, Brian

Hi Catriona, One thing I notice right away is that you're using kinetics for a very large number of minerals. Typically these all won't dissolve or precipitate at the same rate. A carbonate mineral like calcite, for example, is likely to dissolve much more quickly than some of the silicate minerals. When you set kinetics for extremely fast reactions (relative to the timescale of interest) you force extremely small timesteps on the program, which can slow down your run at best, or prevent it from reaching completion in certain cases if the time step keeps getting cut between iterations. It's best to start with a simpler example - an equilibrium model - and gradually add in kinetics as necessary. We did notice an issue with kinetic reaction paths and flow-through models, as described in this forum post. We're looking into this and it should be cleared up in the next maintenance release. Regards, Brian Farrell Aqueous Solutions LLC

Hi Vivek, Besides the new kinetic reaction types available in GWB 9 (formation/ dissociation of aqueous and surface complexes, gas transfer), a number of improvements were made to the manner in which time steps are taken when kinetics are involved. As you mentioned, this helps with the solution of many difficult problems. In general, however, there is no way to say how many kinetic reactions are allowed in any given example. Keep in mind, adding complexity to a model will require a greater computing effort and could introduce difficulties in the numerical solution. You'll need to start with a simple model that makes sense, and then gradually refine it if necessary to accomplish your modeling goals. You shouldn't set kinetic rate laws for reactions of widely differing rates. Reactions that will proceed quickly over the time span of the calculation should assume equilibrium. Those which will proceed negligibly over the calculation should be suppressed. Only for the reactions which fit into neither category should you use kinetics - those which will proceed slowly but measurably over the time span of interest. If you want to set different rate parameters for a mineral's precipitation and dissolution reactions, you can't just add a second entry using the standard ("Built-in") rate law. You'll need to use a custom rate law in this case: either an equation, a script, or a function that you write and compile. Take a look at the Custom Rate Laws section in the GWB Reaction Modeling Guide. In section 5.2 you'll find a simple rate law script that is similar to the problem you describe. Regards, Brian

Hi Jay, Your X1t script runs to 1000 years as far as I can tell. Look for an Xi = 1 at the bottom of the Results pane. There should also be a dialog box which comes up telling you "Run complete: Simulation ended successfully at year 1000" and giving you the option to plot your results. In Xtplot, look for the Time levels tab to select the time step to plot. There is an issue with your React script however, due to some changes to the time step limitations implemented in GWB 9 (earlier versions do not have this problem). Basically the combination of flow-through path, in which mineral mass is removed at each step, and kinetic precipitation/ dissolution reactions, in which the mass of a mineral factors into the reaction rate (and the limiting time step as a result) causes React to take increasingly small time steps. This will be fixed in the next maintenance release. Regards, Brian

Hi Jay, What version of GWB are you running? Can you tell me when the simulation ends for you? Is there any problem solving the model? Do you autoscale turned off in Gtplot, by any chance? Thanks, Brian

Hi Johan, Thanks for posting this. Hopefully next week I'll have something for you. Regards, Brian

Hi Johan, Can you explain your problem a little more? Do you receive a specific error message? Could you attach your React script and any additional files I'll need to run it? Thanks, Brian Farrell Aqueous Solutions LLC

Hi Frank, Glad to hear it's useful in your work, and to see THEREDA making new thermo datasets freely available. Cheers, Brian

Hi Maki, As you mention, TDS circles are only used with a few special plots: Piper, Durov, and Ternary diagrams. Although the diagrams themselves can be constructed using output from SpecE8 or React, they were only designed to incorporate TDS circles when using scatter data stored in GSS spreadsheets. If you're using GWB 9, you can simply drag the contents of your SpecE8 Basis pane into GSS (as a new sample), then launch the diagram from within GSS. Since you can store any number of samples in GSS (and assign each a maker symbol, color, and size), you'll be able to compare all of your samples at once in a single diagram. If it would be useful, we'll look into adding TDS circle capabilities for plots generated from SpecE8 or React output. Regards, Brian Farrell Aqueous Solutions LLC

Hi Jay, The TDS isn't really a constraint in your model. Rather, GWB will use the value you specify to convert your input constraints into molal units, where necessary. The value you see in the output is actually calculated by adding the masses of all solutes in the final speciated fluid. Assuming you've accounted for everything in the fluid, the TDS measured from the lab should be pretty close to the output. In your case, I think the calculated value is different from the value you specify primarily due to the CO2 in your system. The high CO2 fugacity adds a lot of dissolved CO2 to your system. If you look at the concentrations of species in the printed output, you'll find CO2(aq) is the most abundant species by far. It turns out that fluid degassing between sampling and analysis can lower the TDS content of a fluid. If your fluid was sampled at pressure (f CO2 = 100) it surely lost a lot of CO2 prior to analysis. To see how this works, you can add the CO2 back into your fluid using a "sliding fugacity path". Start with your initial system at atmospheric conditions (I'm guessing you'll have a pH value closer to neutral, otherwise set log f CO2 to -3.5), then move to the reactants pane and click "add" - Sliding - Gas... CO2(g), then specify a target of 100. You'll see how dissolving CO2 into your fluid sends the TDS up. Hope this helps, Brian

Hi Jay, Are you sure this is the script that is giving you problems? The TDS isn't saved in your input anywhere (React will calculate TDS as a result), but it's still close to 16 g/l at the beginning of your run (Xi, or reaction progress = 0). Perhaps you need to check your axes in Gtplot and make sure you're actually looking at the starting point? Regards, Brian

Hi Jay, Oops, I guess the problem is with your React example, which is saved with the .rea extension. Your X1t example will be saved as a .x1t file. In either case, you just set up your model then go to File - Save as to save your run. When you add a reply to this topic, just look for the "Attachments" section and browse for your file, then click "Attach This File." If you post both your React and X1t examples that would be helpful. Thanks, Brian