Explaining behavior of Kaolinite dissolution/plotting rates

[smolen]

Hey all, new here. Please excuse any silly questions or misunderstandings on my part.

I am attempting to model the time it takes for a batch reactor containing kaolinite and pure water to come to solubility equilibrium. I am using a RATES block courtesy of Zhu et al. (2019), and the rest I have taken from bits and pieces of other code that I have found. Here is the code:

Code: [Select]

# Kinetic dissolution of kaolinite

RATES

Kaolinite

# from Marty et al 2015
# pre-exponent coefficient A is calculated from logk using equation A=k/exp(-Ea/RT)
# experimental condition range T=22-80C, pH=0.5-12

-start
1 rem unit should be mol,Liter-1 and second-1
2 rem parm(1) is surface area in the unit of m2/L
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals
5 rem M0 is the initial moles of minerals
6 rem parm(2) is a scaling factor
10 rem acid solution parameters
11 a1=2.56E-04
12 E1=43000
13 n1=0.51
20 rem neutral solution parameters
21 a2=5.0E-08
22 E2=38000
30 rem base solution parameters
31 a3=2.87E-03
32 E3=46000
33 n2=0.58
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Kaolinite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= Rate*Time
110 rem do not dissolve more minerals than present
115 if (moles>M) then moles=M
200 save moles
-end


KINETICS
Kaolinite
 -m0  0.387                        # initial moles of kaolinite
 -parms 2380  1
 -time_step 0.01 day in 400           # X days in Y steps
INCREMENTAL_REACTIONS true

SOLUTION 1
 -water 0.5
 -temp 200

USER_GRAPH
 -chart_title "Kaolinite dissolution"
 -axis_titles Hours "SR" "KIN_DELTA"
 #-axis_titles ACT("H+") "SR"

 -initial_solutions  true
 -start
 10 graph_x total_time / 3600  # time in hours on x-axis
 20 graph_y SR("Kaolinite")
 #30 graph_sy -la("H+")   
 #30 graph_sy KIN_DELTA("Kaolinite")
 30 graph_sy Rate1
 #30 PLOT_XY ACT("H+"),SR("Kaolinite")   
 -end
END

I am trying to explain the trend I am seeing from the graph results, wherein the progress towards solubility equilibrium is very slow up to ~0.04 hr, then rapidly increases (run the code for a plot, I am unable to add attachments - apologies for that). Dissolution should be fastest initially, when the pure water is free of solutes.

I believe it has something to do with formation of other species that are considered within phreeqc, and/or the pH (H+ activity). But I don't know how to investigate this within the code, besides plotting things on a secondary axis along with the SR. It appears that KIN_DELTA ("Kaolinite"), which shows the amount of Kaolinite dissolving, shows a change in behavior at that same time ~0.04 hr.

I would like to plot Rate1, Rate2, and Rate 3 on the secondary axis to investigate further, however I am not sure how to do so. When I put:

Code: [Select]

graph_sy Rate1
it does not work and just shows as 0 for the whole time frame.

Can anyone help me to explain this behavior, or assist me in how to plot the Rates? Thanks in advance to anyone taking the time to look at this - your help is greatly appreciated!

- Jon

[dlparkhurst]

The sign convention is that a positive number for the SAVEd value decreases the number of moles of kinetic reaction; that is, KIN_DELTA will be negative (as is the case in your script).

I think the plot of SR is probably confusing you because the concentrations are changing orders of magnitude while the SR remains close to zero. A plot of SI (saturation index) or the increase in silica concentration will conform more to your expectations, with a monotonic, concave-downward increase to equilibrium.

The RATES Basic program runs independently from the USER_GRAPH Basic program. All Basic programs have independent name spaces; they do not share the same variables. There is a way to pass values among Basic programs with PUT and GET. You can PUT a value in RATES and GET the value in USER_GRAPH as I have done in the example. This is not exactly right because the RATES Basic program is executed many times and the final execution before USER_GRAPH may not be at the end of the integration time interval. In most cases, the values from the last RATES execution are close to the values that would be calculated at the end of the interval. To do the calculate the rates at the end of the integration interval exactly, the code from RATES should be included in USER_GRAPH.

More generally, the order that the Basic programs are evaluated will affect GET and PUT. I think the RATES are evaluated in the order given in KINETICS, followed by USER_GRAPHs at the end of the integration in the order of their user numbers, and USER_PRINT is evaluated after USER_GRAPHs (but I am not completely sure).

The average rate over the time interval may be more representative, and it can be calculated as KIN_DELTA("Kaolinite")/KIN_TIME.

Code: [Select]

# Kinetic dissolution of kaolinite

RATES

Kaolinite

# from Marty et al 2015
# pre-exponent coefficient A is calculated from logk using equation A=k/exp(-Ea/RT)
# experimental condition range T=22-80C, pH=0.5-12

-start
1 rem unit should be mol,Liter-1 and second-1
2 rem parm(1) is surface area in the unit of m2/L
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals
5 rem M0 is the initial moles of minerals
6 rem parm(2) is a scaling factor
10 rem acid solution parameters
11 a1=2.56E-04
12 E1=43000
13 n1=0.51
20 rem neutral solution parameters
21 a2=5.0E-08
22 E2=38000
30 rem base solution parameters
31 a3=2.87E-03
32 E3=46000
33 n2=0.58
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Kaolinite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= Rate*Time
110 rem do not dissolve more minerals than present
115 if (moles>M) then moles=M
200 save moles
210 PUT(rate1 *(1-Sr_mineral)*SA*parm(2), 1)
220 PUT(rate2 *(1-Sr_mineral)*SA*parm(2), 2)
230 PUT(rate3 *(1-Sr_mineral)*SA*parm(2), 3)
-end

KINETICS
Kaolinite
 -m0  0.387                        # initial moles of kaolinite
 -parms 2380  1
 -time_step 0.01 day in 400           # X days in Y steps
INCREMENTAL_REACTIONS true

SOLUTION 1
 -water 0.5
 -temp 200

USER_GRAPH
 -chart_title "Kaolinite dissolution"
 -headings hours Si Rate1 Rate2 Rate3 Rate
 -axis_titles Hours "Si, mol/kgw" "Rate, mol/s"
 -axis_scale sy_axis auto auto auto auto log
 -initial_solutions  true
 -start
 10 graph_x total_time / 3600  # time in hours on x-axis
 20 GRAPH_Y TOT("Si")
 30 GRAPH_SY GET(1), GET(2), GET(3), -KIN_DELTA("Kaolinite") / KIN_TIME
 -end
END

[smolen]

D.L.,

The code you provided is exactly what I was looking for, and I understand that portion of the syntax now. Thank you so very much for your help.

Still, my system appears to exhibit fundamentally different behavior (not just quantity) from the example with quartz that I began with (https://hydrochemistry.eu/exmpls/kin_qu.html), in terms of SR(mineral) vs. time.

Follow-up question (for anyone that may be reading):
Why does the rate/time to SR=1 change when I multiply initial moles (M0) and solution volume (water) by the same factor? For instance, if I divide each by 10, it slows down, and if I divide each by 100 or 1000, it gives unreasonable plot results that I cannot explain. In other words, why does changing the individual values while keeping the M0/water ratio constant effect the time to solubility equilibrium?

[dlparkhurst]

Basically, your rate expression is the same, regardless of M0 or water volume. You are not changing the surface area proportionally with M0.

With Phreeqc version 3, we realized that it is better to write the parameter in KINETICS as area/mol rather than strictly area. Using area/mol allows a kinetics definition to scale correctly and more intuitively. In PHAST, KINETICS definitions may be scaled depending on water volume or rock volume in a cell, so it is important for your little experiment to work as you expected. So, I would reconsider how to change the RATES definition to use PARM(1) in terms of area per mole.

[smolen]

D.L,

It took me quite a while, but I believe I understand what you are saying (or at least have come to an explanation of some sort).

I am changing the surface area proportionally by changing M0 and solution volume by the same factor, but the rate expression remains equal (at least for the first time step). Thus, for the first time step, the rate of dissolution will be equal. Then, M/M0 will be less than after the first time step of the system containing more mineral (though same SA m2/L ratio), which factors into the rate. Overall, it takes LESS time to reach solubility equilibrium with a system of smaller size, even if the m2/L ratio is the same.

This same rate at the first time step is the reason why dividing by too large a factor and not decreasing the time step range provided meaningless results.

This is physically inaccurate, as the absolute rate should increase with increasing absolute m2 of mineral... thus, you recommend I put PARM(1) in terms of m2/mol rather than m2/L? Does this make sense?

[dlparkhurst]

Sorry, I don't think I was clear. I suggest that you change PARM(1) to be m^2/mol and change the rate  expression accordingly. So, here is the KINETICS definition, where 6150 = 2380 / 0.387 m^2/mol.

Code: [Select]

Kaolinite
 -m0  0.387                        # initial moles of kaolinite
 -parms 6150 1 #2380  1

The surface area is then calculated as

Code: [Select]

43 if (M0<=0) then SA=PARM(1) else SA=[PARM(1)* M0] * (M/M0)^0.67

So, when M0 changes, the surface area changes proportionally. If you double the solution volume and double M0, the surface area will be double, your rate will double, and the approach to equilibrium will be the same for the doubled and undoubled system.

[Probably too much information, but PHAST uses the data block KINETICS_MIX (http://water.usgs.gov/water-resources/software/PHREEQC/Phreeqc_ReleaseNotes.tx) to scale reactants. When a KINETICS definition is multiplied by a coefficient, M and M0 are multiplied by the coefficient, but the PARMs are not.]

[smolen]

D.L.,

Got it - implemented that and it works exactly as expected.

Now, isn't this the physically accurate model? I am just trying to understand why Zhu et al. (2019) deal with PARM(1) in m^2/L, because this does not scale properly and seems to be physically incorrect to me...

P.S. I have not dealt with PHAST and I'm not really sure what it is, but I certainly will. Thanks for the insight. Unfortunately, the link you provided does not appear to work properly.

[dlparkhurst]

If you use 1 L and do not use different M0, then it does not matter whether you write it with M^2, M^2/L, or M^2/mol. Constant conditions apply to many many batch reaction calculations, it is only when you changed M0 and solution volume that a problem occurred.

Sorry, the complete link was corrupted: http://water.usgs.gov/water-resources/software/PHREEQC/Phreeqc_ReleaseNotes.txt

PHAST is a 3D reactive-transport model based on Phreeqc and PhreeqcRM https://www.usgs.gov/software/phast-computer-program-simulating-groundwater-flow-solute-transport-and-multicomponent. I probably should not even have mentioned PHAST, except it provides a situation where M and M0 may change depending on the media properties of cells (porosity, cell volume, and others). In this case, the KINETIC reactions need to scale properly (possibly differently for different cells).

[smolen]

D.L.,

Makes sense. Thank you so much again, I think I am in a good place for now with that solubility model. I will look into PHAST, as the features you mentioned may be of interest to me as this project progresses.

It appears that, of the phases considered in the model, Gibbsite should be precipitating (SR=12-18 at the end of my simulations). How can I see if it does, in fact, precipitate? TOTMOLE argument seems to work only with elements, not solid species.

Also, any chance you've heard of Phreeqc being used for modeling Ostwald ripening? This is likely the next piece of my project.

[dlparkhurst]

If you want gibbsite to precipitate if supersaturated, then you need to add "EQUILIBRIUM_PHASES; Gibbsite 0 0" to your script; the second 0 indicates gibbsite is not present initially.

You can print the amount of Gibbsite that precipitates with "SELECTED_OUTPUT; -equi Gibbsite" or by using the Basic function EQUI("Gibbsite") in USER_PUNCH.

I don't know much about Ostwalt ripening. A colleague Dennis Eberl did experiments and wrote several papers on the subject. A Google Scholar search for "D Eberl ripening" will find them.

[smolen]

Brilliant.

Including gibbsite (which I believe is more accurate) completely changed the time to saturation for kaolinite, by 24x (depending on T)! This only amounts to ~0.009 weight% gibbsite, which is likely why I have not detected it in my XRD analysis. Yet, clearly, it significantly effects the kinetics of kaolinite dissolution. Something very important to consider.

And I have looked into Eberl's papers, very relevant to me. I can't thank you enough for all of the insight.