Negative mols of Si error

[Mgill9460]

Hi All,

I get an error when running this simulation saying negative moles of Si in solution. I have tried adjusting the Si in initial solution but still get the same error.

Any input is appreciated.

Code: [Select]

TITLE Beaker Dissolution Using carbfix_kin.dat
#the carbfix_kin database has predefined rates blocks within the database
SOLUTION 1 Ocean_Water
    temp      70
    pH        8.19
    pe        10
    redox     pe
    units     mol/l
    density   1
    Mg        0.058
    Si        4.9e-06
    Ca        0.008
    Fe        1.2e-07
    Cl        0.54 CHARGE
    Na        0.28
    C         0.0034
    S(6)      0.027
    K         0.0097
    -water    1 # kg


KINETICS 1 dissolution over ~40 day
Forsterite
-m0 0.32
-parms 0 40 0
#Parm(1) = specifies if specific surface area is m2/g of rock (0) or m2 per kg of water (1)
#Parm(2) = specifies the specific surface area either m2 or m2/kgw based on Parm(1)
#Parm(3) = specifies how surface area changes during dissolution (only avaliable when Parm(1) = 0
#      (0) surface area changes linearly with moles of mineral present
#      (1) surface area changes according to the geometry  of dissolving cubes or spheres

Lizardite
-m0 0.15
-parms 0 40 0
-steps 3197000 in 7
-cvode true
SAVE Solution 1
RATES
#note - in carbfix_kin.dat rate for forsterite is already defined
#note - there is not rate block for lizardite so it is defined here
#Dissolution parameters for Lizardite in Palandri and Kharaka 2004
Lizardite
-start
1 name$ = "Lizardite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =3.36e7# mol.m-2.s-1
1001 Ab  =3.28e-03# mol.m-2.s-1
1002 Ea  =75500# J/mol
1003 Eb  =56600# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.8
1008 nb = 0.
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk))) * S
2009 rplus = rplusa + rplusb
4000 rate = rplus * (1 - SR ("Lizardite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end

SELECTED_OUTPUT
-file CC1-FT-144hrv3.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 concentration Mg"
-axis_titles days "mmol"
-initial_solution true
-start
10 graph_x total_time/86400
20 graph_y TOTMOLE("Mg")*1000
-end
USER_GRAPH 3
-headings Si
-chart_title "ion concentration Si"
-axis_titles days "mmol"
-initial_solution true
-start
10 graph_x total_time/86400
20 graph_y TOTMOLE("Si")*1000
-end

[dlparkhurst]

Do you get an ERROR or a WARNING? If you get only WARNINGs, then the calculation should run to completion and the results should be correct. An ERROR will cause the run to abort and not finish.

I tried downloading CARBFIX_KIN, but there was not RATE definition for Forsterite. If you need more help, add the Forsterite RATES definition to your input file and repost it.

[Mgill9460]

Hi David, just a warning actually.

Strange that your Carb_Kin.dat doesn't have forsterite. I added it here.

I also added brucite to the kinetic reactants. I notice that when I run the model now precipitates until reaching and SI of 0.

Is there a way to block the precipitation of this mineral if I want but still allow for dissolution?

Code: [Select]

TITLE Beaker Dissolution Using carbfix_kin.dat
#the carbfix_kin database has predefined rates blocks within the database
SOLUTION 1 Ocean_Water
    temp      70
    pH        8.19
    pe        10
    redox     pe
    units     mol/l
    density   1
    Mg        0.058
    Si        4.9e-06
    Ca        0.008
    Fe        1.2e-07
    Cl        0.54 CHARGE
    Na        0.28
    C         0.0034
    S(6)      0.027
    K         0.0097
    -water    1 # kg


KINETICS 1 dissolution over ~40 day
#Parm(1) = specifies if specific surface area is m2/g of rock (0) or m2 per kg of water (1)
#Parm(2) = specifies the specific surface area either m2 or m2/kgw based on Parm(1)
#Parm(3) = specifies how surface area changes during dissolution (only avaliable when Parm(1) = 0
#      (0) surface area changes linearly with moles of mineral present
#      (1) surface area changes according to the geometry  of dissolving cubes or spheres
Forsterite
-m0 0.32
-parms 0 10 0

Lizardite
-m0 0.15
-parms 0 10 0
Brucite
-m0 0.036
-parms 0 10 0
-steps 3197000 in 7
-cvode true
SAVE Solution 1
RATES
#note - in carbfix_kin.dat rate for forsterite is already defined
#note - there is not rate block for lizardite so it is defined here
#Dissolution parameters for Lizardite in Palandri and Kharaka 2004

Forsterite   #Mg2SiO4, M 140.692 g/mol   
-start
1 name$ = "Forsterite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =14.8e4# mol.m-2.s-1
1001 Ac  =220# mol.m-2.s-1
1002 Ea  =70400# J/mol
1003 Ec  =60900# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.44
1008 nc = 0.22
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusc = Ac* (exp(-Ec/ (R * Tk)))* (act("H+")^nc) * S
2009 rplus = rplusa + rplusc
4000 rate = rplus * (1 - SR ("Forsterite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end

Lizardite
-start
1 name$ = "Lizardite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =3.36e7# mol.m-2.s-1
1001 Ab  =3.28e-03# mol.m-2.s-1
1002 Ea  =75500# J/mol
1003 Eb  =56600# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.8
1008 nb = 0
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk))) * S
2009 rplus = rplusa + rplusb
4000 rate = rplus * (1 - SR ("Lizardite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end

Brucite
-start
1 name$ = "Brucite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =4.04e05# mol.m-2.s-1
1001 Ab  =1.31e-01# mol.m-2.s-1
1002 Ea  =59000# J/mol
1003 Eb  =42000# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.5
1008 nb = 0
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk))) * S
2009 rplus = rplusa + rplusb
4000 rate = rplus * (1 - SR ("Brucite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end
SELECTED_OUTPUT
-file 2_CC1-FT-T70_SA10.sel
-kinetic_reactants Forsterite Lizardite Brucite
-saturation_indices Forsterite Lizardite Brucite
-molalities Ca+2 Mg+2
-totals Mg Si Ca
USER_GRAPH 1
-chart_title "ultramafic dissolution"
-headings Forsterite pH
-axis_titles Days "mol"
-initial_solutions true
-start
10 PLOT_XY total_time/86400, KIN_DELTA("Forsterite"), 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 concentration Mg"
-axis_titles days "mmol"
-initial_solution true
-start
10 graph_x total_time/86400
20 graph_y TOTMOLE("Mg")*1000
-end


[Mgill9460]

Hi David,

Another follow up question actually.

Is there a way to buffer against large pH changes in this type of simulation?

The pH jumps from 8 to 12 in the first step of the simulation. Is there some sort of phase that can be added as a buffer?

[dlparkhurst]

You can increase the log K for Brucite or change the kinetic RATES definition so that the reaction stops at a different point.

To end up at a lower pH, you must either decrease the extent of reaction, or perhaps add more CO2.

[Mgill9460]

Hi David, is there a way to fix pH or another way to buffer without adding CO2 to equilibrium phases? When I add the CO2 phase the Mg concentration increases too quickly.

Code: [Select]

TITLE Beaker Dissolution Using carbfix_kin.dat
#the carbfix_kin database has predefined rates blocks within the database
SOLUTION 1 Ocean_Water
    temp      55
    pH        8.19
    pe        8.451
    redox     pe
    units     mol/l
    density   1.02
    Mg        0.058
    Si        4.9e-06
    Ca        0.008
    Fe        1.2e-07
    Cl        0.54
    Na        0.28
    C         0.0034
    S(6)      0.027
    K         0.0097
    -water    1 # kg

EQUILIBRIUM_PHASES 1
CO2(g) -3.5 1

KINETICS 1 dissolution over ~40 day
#Parm(1) = specifies if specific surface area is m2/g of rock (0) or m2 per kg of water (1)
#Parm(2) = specifies the specific surface area either m2 or m2/kgw based on Parm(1)
#Parm(3) = specifies how surface area changes during dissolution (only avaliable when Parm(1) = 0
#      (0) surface area changes linearly with moles of mineral present
#      (1) surface area changes according to the geometry  of dissolving cubes or spheres
Forsterite
-m0 0.32
-parms 0 25 0

Lizardite
-m0 0.15
-parms 0 25 0


-steps 3197000 in 7
-cvode true
SAVE Solution 1
RATES
#note - in carbfix_kin.dat rate for forsterite is already defined
#note - there is not rate block for lizardite so it is defined here
#Dissolution parameters for Lizardite in Palandri and Kharaka 2004

Forsterite   #Mg2SiO4, M 140.692 g/mol   
-start
1 name$ = "Forsterite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =14.8e4# mol.m-2.s-1
1001 Ac  =220# mol.m-2.s-1
1002 Ea  =70400# J/mol
1003 Ec  =60900# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.44
1008 nc = 0.22
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusc = Ac* (exp(-Ec/ (R * Tk)))* (act("H+")^nc) * S
2009 rplus = rplusa + rplusc
4000 rate = rplus * (1 - SR ("Forsterite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end

Lizardite
-start
1 name$ = "Lizardite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =3.36e7# mol.m-2.s-1
1001 Ab  =3.28e-03# mol.m-2.s-1
1002 Ea  =75500# J/mol
1003 Eb  =56600# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.8
1008 nb = 0
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk))) * S
2009 rplus = rplusa + rplusb
4000 rate = rplus * (1 - SR ("Lizardite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end

Brucite
-start
1 name$ = "Brucite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =4.04e05# mol.m-2.s-1
1001 Ab  =1.31e-01# mol.m-2.s-1
1002 Ea  =59000# J/mol
1003 Eb  =42000# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.5
1008 nb = 0
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk))) * S
2009 rplus = rplusa + rplusb
4000 rate = rplus * (1 - SR ("Brucite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end
SELECTED_OUTPUT
-file 11_CC1-FT-T55_SA25.sel
-kinetic_reactants Forsterite Lizardite
-saturation_indices Forsterite Lizardite
-molalities Ca+2 Mg+2
-totals Mg Si Ca
USER_GRAPH 1
-chart_title "ultramafic dissolution"
-headings Forsterite pH
-axis_titles Days "mol"
-initial_solutions true
-start
10 PLOT_XY total_time/86400, KIN_DELTA("Forsterite"), 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 concentration Mg"
-axis_titles days "mmol"
-initial_solution true
-start
10 graph_x total_time/86400
20 graph_y TOTMOLE("Mg")*1000
-end

Cheers,

[dlparkhurst]

I think you need a process. Just forcing the pH to a value doesn't give me much confidence that you are really modeling your system.

That said, you could account for buffering with SURFACE reactions.

[Mgill9460]

Hi David,

Thank you.

I don't think that surface_master_species are defined in the CARBFIX_KIN.dat.

I've copied over he SURFACE_MASTER_SPECIES code from the generic PHREEQC.dat

Code: [Select]

SURFACE_MASTER_SPECIES
Hfo_s Hfo_sOH
Hfo_w Hfo_wOH
SURFACE_SPECIES
# All surface data from
# Dzombak and Morel, 1990
#
#
# Acid-base data from table 5.7
#
# strong binding site--Hfo_s,

Hfo_sOH = Hfo_sOH
log_k 0.0

Hfo_sOH + H+ = Hfo_sOH2+
log_k 7.29 # = pKa1,int

Hfo_sOH = Hfo_sO- + H+
log_k -8.93 # = -pKa2,int

# weak binding site--Hfo_w

Hfo_wOH = Hfo_wOH
log_k 0.0

Hfo_wOH + H+ = Hfo_wOH2+
log_k 7.29 # = pKa1,int

Hfo_wOH = Hfo_wO- + H+
log_k -8.93 # = -pKa2,int
###############################################
# CATIONS #
###############################################
#
# Cations from table 10.1 or 10.5
#
# Calcium
Hfo_sOH + Ca+2 = Hfo_sOHCa+2
log_k 4.97

Hfo_wOH + Ca+2 = Hfo_wOCa+ + H+
log_k -5.85
# Strontium
Hfo_sOH + Sr+2 = Hfo_sOHSr+2
log_k 5.01

Hfo_wOH + Sr+2 = Hfo_wOSr+ + H+
log_k -6.58

Hfo_wOH + Sr+2 + H2O = Hfo_wOSrOH + 2H+
log_k -17.6
# Barium
Hfo_sOH + Ba+2 = Hfo_sOHBa+2
log_k 5.46

Hfo_wOH + Ba+2 = Hfo_wOBa+ + H+
log_k -7.2 # table 10.5
#
# Cations from table 10.2
#
# Cadmium
Hfo_sOH + Cd+2 = Hfo_sOCd+ + H+
log_k 0.47

Hfo_wOH + Cd+2 = Hfo_wOCd+ + H+
log_k -2.91
# Zinc
Hfo_sOH + Zn+2 = Hfo_sOZn+ + H+
log_k 0.99

Hfo_wOH + Zn+2 = Hfo_wOZn+ + H+
log_k -1.99
# Copper
Hfo_sOH + Cu+2 = Hfo_sOCu+ + H+
log_k 2.89

Hfo_wOH + Cu+2 = Hfo_wOCu+ + H+
log_k 0.6 # table 10.5
# Lead
Hfo_sOH + Pb+2 = Hfo_sOPb+ + H+
log_k 4.65

Hfo_wOH + Pb+2 = Hfo_wOPb+ + H+
log_k 0.3 # table 10.5
#
# Derived constants table 10.5
#
# Magnesium
Hfo_wOH + Mg+2 = Hfo_wOMg+ + H+
log_k -4.6
# Manganese
Hfo_sOH + Mn+2 = Hfo_sOMn+ + H+
log_k -0.4 # table 10.5

Hfo_wOH + Mn+2 = Hfo_wOMn+ + H+
log_k -3.5 # table 10.5
# Iron, strong site: Appelo, Van der Weiden, Tournassat & Charlet, subm.
Hfo_sOH + Fe+2 = Hfo_sOFe+ + H+
log_k -0.95
# Iron, weak site: Liger et al., GCA 63, 2939, re-optimized for D&M
Hfo_wOH + Fe+2 = Hfo_wOFe+ + H+
log_k -2.98

Hfo_wOH + Fe+2 + H2O = Hfo_wOFeOH + 2H+
log_k -11.55
###############################################
# ANIONS #
###############################################
#
# Anions from table 10.6
#
# Phosphate
Hfo_wOH + PO4-3 + 3H+ = Hfo_wH2PO4 + H2O
log_k 31.29

Hfo_wOH + PO4-3 + 2H+ = Hfo_wHPO4- + H2O
log_k 25.39

Hfo_wOH + PO4-3 + H+ = Hfo_wPO4-2 + H2O
log_k 17.72
#
# Anions from table 10.7
#
# Borate
Hfo_wOH + H3BO3 = Hfo_wH2BO3 + H2O
log_k 0.62
#
# Anions from table 10.8
#
# Sulfate
Hfo_wOH + SO4-2 + H+ = Hfo_wSO4- + H2O
log_k 7.78

Hfo_wOH + SO4-2 = Hfo_wOHSO4-2
log_k 0.79
#
# Derived constants table 10.10
#
Hfo_wOH + F- + H+ = Hfo_wF + H2O
log_k 8.7

Hfo_wOH + F- = Hfo_wOHF-
log_k 1.6
#
# Carbonate: Van Geen et al., 1994 reoptimized for D&M model
#
Hfo_wOH + CO3-2 + H+ = Hfo_wCO3- + H2O
log_k 12.56
 
Hfo_wOH + CO3-2 + 2H+= Hfo_wHCO3 + H2O
log_k 20.62
I used code from example 8 to include HFo as a SURFACE species but get the error:

Code: [Select]

ERROR: Elements in species have not been tabulated, Ba+2.
ERROR: Reaction for species has not been defined, Ba+2.
ERROR: Elements in species have not been tabulated, Cd+2.
ERROR: Reaction for species has not been defined, Cd+2.
ERROR: Elements in species have not been tabulated, H3BO3.
ERROR: Reaction for species has not been defined, H3BO3.
ERROR: Elements in species have not been tabulated, Pb+2.
ERROR: Reaction for species has not been defined, Pb+2.
ERROR: Elements in species have not been tabulated, Sr+2.
ERROR: Reaction for species has not been defined, Sr+2.
ERROR: Calculations terminating due to input errors.

Also is it possible that defining alkalinity in the solution could help with large pH values? I have tried to include it but noticed no change.

Code: [Select]

TITLE Beaker Dissolution Using carbfix_kin.dat
#the carbfix_kin database has predefined rates blocks within the database
SOLUTION 1 Ocean_Water
    temp      55
    pH        8.19
    pe        8.451
    redox     pe
    units     mol/l
    density   1.02
    Alkalinity 180 mg/L       as HCO3
    C         0.0034
    Ca        0.008
    Cl        0.54
    Fe        1.2e-07
    K         0.0097
    Mg        0.058
    Na        0.28
    S(6)      0.027
    Si        4.9e-06
    -water    1 # kg

KINETICS 1 dissolution over ~40 day
#Parm(1) = specifies if specific surface area is m2/g of rock (0) or m2 per kg of water (1)
#Parm(2) = specifies the specific surface area either m2 or m2/kgw based on Parm(1)
#Parm(3) = specifies how surface area changes during dissolution (only avaliable when Parm(1) = 0
#      (0) surface area changes linearly with moles of mineral present
#      (1) surface area changes according to the geometry  of dissolving cubes or spheres
Forsterite
-m0 0.32
-parms 0 25 0

Lizardite
-m0 0.15
-parms 0 25 0


-steps 3197000 in 7
-cvode true
#----------------------------SURFACE-----------------------------#
SURFACE_SPECIES
        Hfo_sOH  + H+ = Hfo_sOH2+
        log_k  7.18
SURFACE 1
        Hfo_sOH        1e-3    25    0.1

SAVE Solution 1
RATES
#note - in carbfix_kin.dat rate for forsterite is already defined
#note - there is not rate block for lizardite so it is defined here
#Dissolution parameters for Lizardite in Palandri and Kharaka 2004

Forsterite   #Mg2SiO4, M 140.692 g/mol   
-start
1 name$ = "Forsterite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =14.8e4# mol.m-2.s-1
1001 Ac  =220# mol.m-2.s-1
1002 Ea  =70400# J/mol
1003 Ec  =60900# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.44
1008 nc = 0.22
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusc = Ac* (exp(-Ec/ (R * Tk)))* (act("H+")^nc) * S
2009 rplus = rplusa + rplusc
4000 rate = rplus * (1 - SR ("Forsterite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end

Lizardite
-start
1 name$ = "Lizardite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =3.36e7# mol.m-2.s-1
1001 Ab  =3.28e-03# mol.m-2.s-1
1002 Ea  =75500# J/mol
1003 Eb  =56600# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.8
1008 nb = 0
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk))) * S
2009 rplus = rplusa + rplusb
4000 rate = rplus * (1 - SR ("Lizardite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end

Brucite
-start
1 name$ = "Brucite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa  =4.04e05# mol.m-2.s-1
1001 Ab  =1.31e-01# mol.m-2.s-1
1002 Ea  =59000# J/mol
1003 Eb  =42000# J/mol
1004 R =  8.314 #J.deg-1.mol-1
1006 Sig = 1   
1007 na = 0.5
1008 nb = 0
#Rate Equation
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk))) * S
2009 rplus = rplusa + rplusb
4000 rate = rplus * (1 - SR ("Brucite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end
SELECTED_OUTPUT
-file 11_CC1-FT-T55_SA25.sel
-kinetic_reactants Forsterite Lizardite
-saturation_indices Forsterite Lizardite
-molalities Ca+2 Mg+2
-totals Mg Si Ca
USER_GRAPH 1
-chart_title "ultramafic dissolution"
-headings Forsterite pH
-axis_titles Days "mol"
-initial_solutions true
-start
10 PLOT_XY total_time/86400, KIN_DELTA("Forsterite"), 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 concentration Mg"
-axis_titles days "mmol"
-initial_solution true
-start
10 graph_x total_time/86400
20 graph_y TOTMOLE("Mg")*1000
-end


Best,

[dlparkhurst]

I think this is just causing you more trouble without getting to your essential question. So, I give you this with reservations.

As you change the number of SurfOH sites from zero to 10 or more moles, you will be adding more and more buffering capacity. I think there are some quirky things with activity coefficients and charge balance, but it will maintain a pH closer and closer to 8.19 as you add more sites.

I would much rather see an inverse model that accounts for the changes, than invoking this mechanism without more evidence.

Code: [Select]

SURFACE_MASTER_SPECIES
Surf SurfOH
SURFACE_SPECIES
SurfOH = SurfOH
-log_k 0

SurfOH + H+ = SurfOH2+
-log_k 7.29 # = pKa1,int

SurfOH = SurfO- + H+
-log_k -8.93 # = -pKa2,int
END
SURFACE
SurfOH 10 1 600
-no_edl
-eq 1