Jump to content
Geochemist's Workbench Support Forum

custom rate laws for siderite dissolution and calcite ppt

Recommended Posts



I am trying to simulate high-pH leachate (~12) flowing through a reactive porous medium composed of siderite (starting out with 20% by volume). I have rate laws that are in the form of mol/sec:


siderite dissolution:


(rate_con + 5.55e-22 * activity("H+")) * (mass("siderite")) * 115.86) * (1 - Q/K)


with a rate_con of 9.77e-12 mol/g/sec


(note - the equation I have specifies the middle term as mass of siderite in grams - so I manually altered it to multiply the mass("siderite") term by the gram formula weight since this command returns mass in terms of moles)




calcite dissolution:


rate_con * (Q/K - 1)


with a rate_con of 6.75e-10 mol/sec


I am trying to send the high Ca and pH leachate through a 30 m permeable reactive barrier and achieve pH control by dissolution of siderite, precip of calcite, and precip of iron hydroxides, which I will assume as equilibrium controlled because this is a very fast reaction.


My main problems are my rate constants are not in the form specified by GWB (mol/cm2 sec), and I don't have surface areas in my equations. However, GWB requires inputs for both these variables. I am getting a dissolution rate of 0 for siderite and precipitation rate of -6e-11 for calcite when I run the model, which is not surprising because the form I am using does not appear to be compatible with GWB.


I am attaching my x1t script (GWB v. 9).









Link to comment
Share on other sites

Hi Zeno,

You’re required by the program to enter a value for specific surface area with kinetic minerals, but if you’re using a custom rate law that doesn’t make use of surface area then the value you specify will be ignored. Take a look at example 3.7, quartz precipitation in a vein, in the GWB Reactive Transport Modeling Guide for an example that should prove useful.

Similarly, the programs expect you to enter a rate constant, but you don't need to use it (or the internal parameter "rate_con" that corresponds to it). If both were in the same units this would make sense, but since they're not you can just define a variable, MyRateConstant, in your custom rate law script. Or, in a less elegant solution, you could hard-code the value of the rate constant into the rate law.

By the way, you can multiply siderite's mass in moles by the internal parameter mw (or the mw("Siderite") helper function) instead of typing the number in yourself. Either way it should work, though. Take a look at Tables 5.1 and 5.2 in the GWB Reaction Modeling Guide for a list of internal parameters and helper functions you can utilize if you wish.

Hope this helps,

Brian Farrell

Aqueous Solutions LLC

Link to comment
Share on other sites

Hi Brian,


Thank you, this was very helpful.


I want to confirm that the custom rate law equation in GWB calculates a dissolution rate in the units of mol/sec, is that correct?


I think the problem with my siderite rate law was that I was using the "mass" helper function, instead of the "moles" function. I rewrote my rate law in this form:


(9.77e-12) * (moles("Siderite") * mw("Siderite")) * (1 - QoverK)


and it seems to be working...


However, my calcite rate law does not seem to be working


(6.75e-10) * (QoverK - 1)


I am still getting a negative and extremely small precipitation rate. When I check "System Parameters" in xtplot I get low decimals fractions for Q/K, Calcite; but when I check "mineral saturation" for calcite I have Q/K's in the 1000's, which makes more sense since I am running very high Ca and alkalinity waters through the model domain, and siderite is dissolving.


Any idea why this might be? I am attaching my script.





Link to comment
Share on other sites

Hi Zeno,


Glad to hear that it helped. Yes, you’re calculating in mol/sec. I think either the mass or moles internal parameters should work, as should the mass(“ ”) helper function. FYI, (moles(“ ”) isn’t listed in the documentation but it works too. Also, the internal parameters are specific to each kinetic reactant. In a rate law for Siderite, for example, the mass internal parameter always refers to the number of moles of Siderite while in a rate law for Calcite that same parameter would refer to Calcite. With a helper function, though, you can refer to any species of interest (using the (“ ”)) from any rate law, for example mass(“Siderite”), mass (“Calcite”), or mass(H+”).


The Q/K values for all minerals in the Mineral saturation list correspond to how they’re written in the thermo dataset (as destruction reactions, i.e. dissolution, dissociation, etc.). The Q/K values for kinetic reactants in the Kinetic parameters (or System parameters, pre-GWB11) correspond to how the rate law is set up. Using the default method of forward reactions (destruction reactions like dissolution and dissociation), the Q/K will appear like those in Mineral saturation, but for reverse reactions (formation reactions, like precipitation and association) you’ll see the inverse. Take a look and confirm this for yourself.


I played around a little with different ways of writing the rate law. Using the normal forward (dissolution) option, I found reasonable-looking, but not necessarily correct results (siderite dissolved and calcite precipitated near the inlet) using


rate = (-6.75e-10) * (QoverK - 1)


or rate = (6.75e-10) * (1 - QoverK)


I’m not sure how the rate law in your reference was derived or whether the assumptions about parameters and reaction direction match up with what GWB expects. I recommend checking out the original paper and running a simple test in a batch system to make sure you have everything straight.


You should also take a look at section 4.3, Forward and reverse reactions, in the GWB Reaction Modeling Guide for more information.


Hope this helps,


Link to comment
Share on other sites

Hi Brian,


Thank you- this has been hugely helpful. I will look into how the rate laws I used were derived.


I have one final problem that is holding me up on this project. I wanted to incorporate DO and have Fe(OH)3 precipitate instead of Fe(OH)2. I added DO and Fe+++ (I made Fe+++ starting concentrations negligible) to the basis, and suppressed minerals I thought to be unlikely to form under these conditions. However, my model fails to converge when DO is included.


When I run the model Fe(OH)3 (ppd) is supersaturated and is swapped for Fe(OH)4-, and the model fails to converge soon after.


This should (hopefully) be the last query on this thread, and I thank you for your help thus far.




Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

  • Create New...