Beaker Dissolution

[Bulbasaur]

Hello PHREEQC Users,

I am interested to make a model which looks at the dissolution of a mineral over time. When I run the code I get no change in forsterite and no change in aqueous concentrations for Ca or Mg. I'm wondering if I am not using the rates block or kinetics block correctly to simulate dissolution.

Thanks in advance,

Code: [Select]

TITLE Beaker Dissolution of Ultramafic Rock Sample in Simulated Seawater
#Using Example 5 from PHREEQC Version 3 as a base/template
#https://phreeqcusers.org/index.php/topic,2316.0.html -- this is a good post
Solution 1 ATSM Artificial Seawater
units    mmol/kgw
    pH       8.1
    temp     25
    Na       485.55
    K        10.2
    Ca       10.28
    Mg       53.0
    Cl       546.0
    S(6)     28.24  as SO4
    C(4)     2.3    as CO3
    B        0.416  as B(OH)3
    Sr       0.091
    Si       0.066  as H4SiO4
    Alkalinity 2.4  as HCO3
KINETICS 1 dissolution over ~1 day
Forsterite
-m 0.0
-m0 0.010
-parms 200000 0.3
-steps 100 1000 10000 100000
RATES
Forsterite
-start
1 rem M = current number of moles of calcite
2 rem M0 = number of moles of calcite initially present
3    rem PARM(1) = A/V, cm^2/L
4    rem PARM(2) = exponent for M/M0
10 si_cc = SI("Forsterite")
20 if (M <= 0 and si_cc < 0) then goto 200
30    k1 = 10^(0.198 - 444.0 / TK )
40    k2 = 10^(2.84 - 2177.0 / TK)
50    if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
60    if TC > 25 then k3  = 10^(-1.1 - 1737.0 / TK )
70    t = 1
80    if M0 > 0 then t = M/M0
90    if t = 0 then t = 1
100   area = PARM(1) * (t)^PARM(2)
110   rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120   rem 1e-3 converts mmol to mol
130   rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
140   moles = rate * TIME
200 SAVE moles
-end
SELECTED_OUTPUT
-file beakerdissV2.sel
-kinetic_reactants Forsterite
-molalities Ca+2 Mg+2
END



[dlparkhurst]

Use PRINT statements in your RATES definition to debug your rate equations. Output from the PRINT statements will be in the output file.

[Bulbasaur]

Hi dlparkhurst, thank you for the reply. I think my code is working now.

I have been trying to make sense of the RATES block example for calcite in the manual so I can create my own based on publications about forsterite. I think based on what I have read that my rates block will be more dependent on pH than temperature for forsterite.

I have read through the original publication Plummer et al. 1978 and I was wondering if I can ask you a couple questions about the calcite rate block:

1) line 4  - Parm (2) is referring to. Is it just the fraction of mineral remaining or what is meant by exponent.
2) line 130 - this step results in mol/s dissolution rate I think. I'm not sure what the last step (1-10^(2/3*si_calcite) is for.

Code: [Select]

TITLE Beaker Dissolution of Ultramafic Rock Sample in Simulated Seawater
#Using Example 5 from PHREEQC Version 3 as a base/template
#https://phreeqcusers.org/index.php/topic,2316.0.html -- this is a good post
Solution 1 ATSM Artificial Seawater
units    mmol/kgw
    pH       8.1
    temp     25
    Na       485.55
    K        10.2
    Ca       10.28
    Mg       53.0
    Cl       546.0
    S(6)     28.24  as SO4
    C(4)     2.3    as CO3
    B        0.416  as B(OH)3
    Sr       0.091
    Si       0.066  as H4SiO4
    Alkalinity 2.4  as HCO3
KINETICS 1 dissolution over ~3 hours
Forsterite
-m0 0.1
-parms 500 0.3
-steps 10000 in 10
RATES
Forsterite
-start
1 rem M = current number of moles of forsterite
2 rem M0 = number of moles of forsterite initially present
3    rem PARM(1) = A/V, cm^2/L
4    rem PARM(2) = exponent for M/M0 #I am unsure what PARM(2) is for
10 si_cc = SI("Calcite") #pulls saturation indicies for calcite
20 if (M <= 0 and si_forsterite < 0) then goto 200 #if current moles of calcite is 0 and SI of calcite < 0 then save current moles
30    k1 = 10^(0.198 - 444.0 / TK ) #k1 is a rate constant dependant on temperature (TK is temperature in Kelvin)
40    k2 = 10^(2.84 - 2177.0 / TK) #k2 is a rate constant dependant on temperature (TK is temeprature in Kelvin)
50    if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK ) #k3 is a rate constant dependant on temperature (TC is temperature in celcius, this if for TC at or below 25)
60    if TC > 25 then k3  = 10^(-1.1 - 1737.0 / TK ) #k3 is a rate constant dependant on temeprature (for TC above 25 degrees)
70    t = 1
80    if M0 > 0 then t = M/M0
90    if t = 0 then t = 1
100   area = PARM(1) * (t)^PARM(2) #(surface area to volume ratio)*(fraction of forsterite remaining^(M0 exponent))
110   rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120   rem 1e-3 converts mmol to mol
130   rate = area * 1e-3 * rf * (1 - 10^(2/3*si_forsterite))
140   moles = rate * TIME
200 SAVE moles
210 PRINT
-end


SELECTED_OUTPUT
-file beakerdissV3.sel
-kinetic_reactants Forsterite
-saturation_indices Forsterite
-molalities Ca+2 Mg+2
END


[dlparkhurst]

The factor is (M/M0)^parm(2), where M/M0 is the fraction of kinetic reactant remaining, and  "^" raising to the power given by parm(2).

(1-10^(2/3*si_calcite) is a factor that ensures the reaction stops at equilibrium. The sign is such that when SI is greater than 0 precipitation occurs, and when SI is less than 0, dissolution occurs.

[Bulbasaur]

Hi dlparkhurst,

Thank you for the explanation. I found another post which created a rates block from data compiled in Palandri and Kharaka 2004. I used this as a template to create rate blocks for forsterite and lizardite. The rate block was originally in this post https://phreeqcusers.org/index.php?topic=2316.0

I was wondering if you knew how they determined PARM(2) which says it is a correction factor?

Initially I get my code to run when I use just the block for forsterite:

Code: [Select]

TITLE Beaker Dissolution of Ultramafic Rock Sample in Simulated Seawater
#Using Example 5 from PHREEQC Version 3 as a base/template
#https://phreeqcusers.org/index.php/topic,2316.0.html -- this is a good post
SOLUTION 1 Seawater
    units   mol/kgw
    temp    25.0
    pH      7.0
KINETICS 1 dissolution over ~40 day
Forsterite
-m0 0.01
-parms 6.3 0.3
-steps 3456000 in 30
RATES
Forsterite
#from _________________________________________________
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters #k1 rate for olivine dissolution under acidic condition
11 a1=8.38E+04 #a1 temperature independan pre-exponential factor
12 E1=67206 #E1 is the molar activation energy
13 n1=0.470 #n1 is the reaction order for acidic conditions
20 rem neutral solution parameters #k2 rate constant
21 a2=1.58E+03 #a2 temperature independant pre-exponential factor
22 E2=79000 #E2 is the molar activation energy
30 rem base solution parameters #k3 rate constant
31 a3=1.00E-07 #a3 temperature independant pre-exponential factor
32 E3=56637 #E3 is the molar activation energy
33 n2=-0.600 #n2 is the reaction order for basic conditions
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Forsterite") #saturation ratio (i.e. IAP/K)
41 if (M<0) then goto 200 #if moles of forsterite less than 0 then save moles (prevent negative?)
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 #2/3 or 0.67 is a geometric scaling value assuming spherical grains are decerasing
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("H+")^n2  #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-SR_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
210 PRINT
-end



SELECTED_OUTPUT
-file beakerdissV3.sel
-kinetic_reactants Forsterite
-saturation_indices Forsterite
-molalities Ca+2 Mg+2
-totals Mg
END


but when I include an additional block for Lizardite the simulation gets stuck on step 1

Code: [Select]

TITLE Beaker Dissolution of Ultramafic Rock Sample in Simulated Seawater
#Using Example 5 from PHREEQC Version 3 as a base/template
#https://phreeqcusers.org/index.php/topic,2316.0.html -- this is a good post
SOLUTION 1 Seawater
    units   mol/kgw
    temp    25.0
    pH      7.0
KINETICS 1 dissolution over ~40 day
Forsterite
-m0 0.01
-parms 6.3 0.3
Lizardite
-m0 0.1
-parms 6.3 0.3
-steps 3456000 in 5
RATES
Forsterite
#from _________________________________________________
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters #k1 rate for olivine dissolution under acidic condition
11 a1=8.38E+04 #a1 temperature independan pre-exponential factor
12 E1=67206 #E1 is the molar activation energy
13 n1=0.470 #n1 is the reaction order for acidic conditions
20 rem neutral solution parameters #k2 rate constant
21 a2=1.58E+03 #a2 temperature independant pre-exponential factor
22 E2=79000 #E2 is the molar activation energy
30 rem base solution parameters #k3 rate constant
31 a3=1.00E-07 #a3 temperature independant pre-exponential factor
32 E3=56637 #E3 is the molar activation energy
33 n2=-0.600 #n2 is the reaction order for basic conditions
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Forsterite") #saturation ratio (i.e. IAP/K)
41 if (M<0) then goto 200 #if moles of forsterite less than 0 then save moles (prevent negative?)
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 #2/3 or 0.67 is a geometric scaling value assuming spherical grains are decerasing
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("H+")^n2  #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-SR_mineral)*SA*parm(2)
100 moles= Rate*Time
200 save moles
210 PRINT
-end
Lizardite
#from _________________________________________________
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters #k1 rate for olivine dissolution under acidic condition
11 a1=3.36E07 #a1 temperature independan pre-exponential factor
12 E1=75500 #E1 is the molar activation energy
13 n1=0.470 #n1 is the reaction order for acidic conditions
20 rem neutral solution parameters #k2 rate constant
21 a2=3.28E03 #a2 temperature independant pre-exponential factor
22 E2=56600 #E2 is the molar activation energy
30 rem base solution parameters #k3 rate constant *note* lizardite has no reported base solution params
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Lizardite") #saturation ratio (i.e. IAP/K)
41 if (M<0) then goto 200 #if moles of forsterite less than 0 then save moles (prevent negative?)
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 #2/3 or 0.67 is a geometric scaling value assuming spherical grains are decerasing
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
90 Rate=(Rate1+Rate2+Rate3)*(1-SR_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
210 PRINT
-end

SELECTED_OUTPUT
-file beakerdissV4_1.sel
-kinetic_reactants Forsterite Lizardite
-saturation_indices Forsterite Lizardite
-molalities Ca+2 Mg+2
-totals Mg
END


Do you have any idea what could be causing the simulation to not run with the second code?

Best,

[dlparkhurst]

Lizardite equilibrates in 10 s, so you should consider whether the rate is reasonable.

If the equilibration is fast, consider using lizardite in EQUILIBRIUM_PHASES rather than KINETICS.

Your file will run if you use -cvode in KINETICS.

[Bulbasaur]

Hi dlparkhurst,

Thank you for pointing that out. I am finding your feedback very helpful and I think I am getting a handle on the software. I agree that it seems to fast so I defined the rate using information from another publication specifically about lizardite dissolution.

I noticed that I was actually getting precipitation of lizardite rather than dissolution. Because of this I decided to add a gas phase to the model.

I noticed with this I get very negative SIs for forsterite and lizardite but also a huge decrease in the pH of the system. Is there something wrong with the way I have added a gas phase to the model?

Thank you,

Code: [Select]

TITLE Beaker Dissolution of Ultramafic Rock Sample in Simulated Seawater
#Using Example 5 from PHREEQC Version 3 as a base/template
#https://phreeqcusers.org/index.php/topic,2316.0.html -- this is a good post
SOLUTION 1 Ocean_Water
    temp      25
    pH        8.1
    pe        4
    redox     pe
    units     mol/kgw
    density   1.025
    Na        0.48555
    K         0.0102
    Ca        0.01028
    Mg        0.0528
    Cl        0.564
    S(6)      0.02824 as SO4
    C(4)      0.00206 as HCO3
    B         0.000416 as B(OH)3
    Si        0.00006 as H4SiO4
GAS_PHASE 1 --  Air
-fixed_pressure
  -pressure       1.0
  -volume         1.0
  -temperature    25.0
CH4(g)          0.0
CO2(g)          0.000316
O2(g)           0.2
  N2(g)           0.78
 KINETICS 1 dissolution over ~40 day
#Forsterite
# -m0 0.01
# -parms 6.3 1
Lizardite
-m0 0.1
-parms 6.3 1
-steps 3456000 in 30
-cvode true
RATES
Forsterite
#from _________________________________________________
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters #k1 rate for olivine dissolution under acidic condition
11 a1=8.38E+04 #a1 temperature independan pre-exponential factor
12 E1=67206 #E1 is the molar activation energy
13 n1=0.470 #n1 is the reaction order for acidic conditions
20 rem neutral solution parameters #k2 rate constant
21 a2=1.58E+03 #a2 temperature independant pre-exponential factor
22 E2=79000 #E2 is the molar activation energy
30 rem base solution parameters #k3 rate constant
31 a3=1.00E-07 #a3 temperature independant pre-exponential factor
32 E3=56637 #E3 is the molar activation energy
33 n2=-0.600 #n2 is the reaction order for basic conditions
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Forsterite") #saturation ratio (i.e. IAP/K)
41 if (M<0) then goto 200 #if moles of forsterite less than 0 then save moles (prevent negative?)
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 #2/3 or 0.67 is a geometric scaling value assuming spherical grains are decerasing
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("H+")^n2  #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-SR_mineral)*SA*parm(2)
100 moles= Rate*Time
200 save moles
210 PRINT
-end
Lizardite
#from Daval et al. 2013
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters #k1 rate for olivine dissolution under acidic condition
11 a1=5.37E-03 #a1 temperature independan pre-exponential factor
12 E1=42000 #E1 is the molar activation energy
13 n1=0.53 #n1 is the reaction order for acidic conditions
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Lizardite") #saturation ratio (i.e. IAP/K)
41 if (M<0) then goto 200 #if moles of forsterite less than 0 then save moles (prevent negative?)
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 #2/3 or 0.67 is a geometric scaling value assuming spherical grains are decerasing
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #single rate expression
90 Rate=(Rate1)*(1-SR_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
210 PRINT
-end



SELECTED_OUTPUT
-file beakerdissV4_gas1.sel
-kinetic_reactants Forsterite Lizardite
-saturation_indices Forsterite Lizardite
-molalities Ca+2 Mg+2
-totals Mg Si
USER_GRAPH 1
-chart_title "ultramafic dissolution"
-headings Lizardite pH
-axis_titles Days "mol"
-initial_solutions true
-start
10 PLOT_XY total_time/86400, KIN_DELTA("Lizardite"), color = RED, symbol = Square, symbol_size = 6, y-axis = 1
20 PLOT_XY  total_time/86400, -la("H+"), color = Blue, symbol = Square, symbol_size = 6, y-axis = 2
-end
USER_GRAPH 2
-headings Mg
-chart_title "ion concentrations"
-axis_titles days "mmol"
-initial_solution true
-start
10 graph_x total_time/86400
20 graph_y TOTMOLE("Mg")*1000
-end


[dlparkhurst]

Yes, there is something wrong with the gas phase.

By default, PHREEQC will calculate redox equilibrium among redox elements. So, when you add O2 and N2, the oxygen will react with dissolved nitrogen to form nitrate, and the pH will be very low. First question is whether you really need the gas phase; often you simply want to fix the partial pressure of gas phases, in which case you can use EQUILIBRIUM_PHASES. You still must not add both N2 and O2. Phreeqc.dat, pitzer.dat, and phreeqc_rates.dat have non reactive versions of N2--Ntg, which will not react with O2. Ntg has only one redox state and cannot react in redox reactions. It would partition between the gas phase and the solution. So, if you really need a gas phase, you probably want to use Ntg. If you don't need a gas phase, don't define N2 or Ntg.