BCF Crystal Growth Theory

[dirk]

I am interested in writing a simple script that incorporates the BCF crystal growth theory equation but when I run React I get fluctuating precipitation/dissolution for the precipitating mineral. The result is fluctuating positive then negative mineral mass for siderite (with initial mass = 0). This occurs when the mineral mass is very low 1e-186 but this low number is the nature of the nonlinear BCF. When left to run, the mass keeps going from positive to negative.

Is there some way to write in the script something to prevent the negative numbers? It appears to be a complete transition from positive to negative to positive, ie. the next logical number in the series during the reaction can end up being the opposite sign.

 

Any help is much appreciated.

 

Attached is the script for siderite.

 

Dirk

Siderite1.rtf

[Brian Farrell]

Hi Dirk,

 

Do you have a React file to use with this rate-law script?

 

Thanks,

 

Brian Farrell

Aqueous Solutions LLC

[dirk]

Hi Brian, see the attached. You will have to add the 2 minerals in the 'add to thermo' file to your thermo.com.V8.R6.dat file.

Hutton CO2 test 4.rea

add to thermo.txt

[Brian Farrell]

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

[dirk]

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

 

The BCF growth theory is the Burton-Cabrera-Frank (Burton et al. 1951) is a non-linear crystal growth model that is commonly accepted (see various publications by Luttge, Arvidson, Aargaard, Hellvang amoungst others). The second term of the Q/K>1 portion is the precipitation and nucleation term derived from classical nucleation theory. It is a small number at low levels of supersaturation but increases significantly especially above lnQ/K of 4 or 5 (depending on the gamma_n term).

The equation comes from Hellvang and Aagaard (2013) IJGGC, V. 15, p. 3-15. They ran it using PHREEQC and provide the PHREEQC input file.

 

Dirk

[Brian Farrell]

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

[dirk]

Hi Brian, yup, look like an issue with the sign, I read dissolution as positive and precipitation as negative rates based on the GWB output so I followed that convention. Getting rid of both the - sign at the front of term 1 and the minus before term 2 fixed the problem.

 

Cheers

Dirk

[Brian Farrell]

Hi Dirk,

 

I'm answering your email query here. If you want ln(Q/K) to be evaluated in your rate law script, not log(Q/K), then you won't be able to use the internal parameter "logQoverK", which is the saturation index, or log10(Q/K). Instead, you'll need to use some of the "helper functions" from Table 5.2 (GWB Reaction Modeling Guide), and the rate law script syntax described in Table 5.3. Specifically, I think you'll need the helper function QovK("Siderite") and the mathematical operator log() instead of log10(). This is kind of confusing, mostly due to the different conventions for logs commonly used by geochemists and mathematicians. So instead of writing "logQoverK" in your script, to get the natural log you'll need log(QovK("Siderite")).

 

You mentioned the script not working when you removed both negative signs from line 90 in your rate law script. I think you only need to remove the first one (I get the model to run when I do this). This is the same as adding a minus sign to each of your rate_acid, rate_neutral, and rate_base parameters, or just multiplying rate_prec by negative 1 (I think this is where it makes most sense, based on the way you set up your model). Please give this a shot and think on it.

 

If that doesn't take care of the problem satisfactorily for you, there are a few things you can try in general to troubleshoot kinetic reactions. Try setting the variable for dxplot to 0 (Config - Output) so that every step in the calculation shows up in Gtplot. You can then take a close look at system parameters like the Siderite dissolution rate, Q/K, etc. You can also set the variable delQ (Config - Stepping) to a number smaller than 0.1, it's default (try 0.01). This limits the change in Q, the activity product, from one step to another, thus limiting the length of a time step and helping with stability.

 

Hope this helps,

Brian