Forsterite rates for kinetic modelling

[maria]

Dear David,

I would very much appreciate your help on my question. I am trying to set up a kinetic modelling for basalts samples to simulate geothermal experiments; I have calculated the minerals moles for my basalt sample, consisting of augite 2.52 moles, albite of 1.27 moles, forsterite of 1.44 moles and hematite of 0.53 moles.

I also model calcite as a secondary mineral with 0 moles. I was able to find diopside parameters so I used those for the augite; I was also able to find the ones for albite, hematite and calcite at the same structure for Rates. I am proving those below. Unfortunately, I couldn't find the ones for forsterite similar to the form/ structure as others are.

Code: [Select]

RATES
Augite
start
10 REM PARM(1) = MSA (Molar surface area)
I would very much thank you if you could  [m^2/mol]
20 si_augi = SI("Augite")
30 if (M <= 0 and si_augi < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_augi > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-6.36)*EXP(-96.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.710)
70 k_neut = 10^(-11.11)*EXP(-40.60e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_augi))
190 moles = r * TIME
200 SAVE moles
end

Albite
start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alb = SI("Albite")
30 if (M <= 0 and si_alb < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_alb > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))
70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA *(1-(10^si_alb))
190 moles = r * TIME
200 SAVE moles
end

Hematite
start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_hema = SI("Hematite")
30 if (M <= 0 and si_hema < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_hema > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-9.39)*EXP(-66.20e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-14.60)*EXP(-66.20e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_hema))
190 moles = r * TIME
200 SAVE moles
end

Calcite
start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_calc = SI("Calcite")
30 if (M <= 0 and si_calc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_calc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_calc))
190 moles = r * TIME
200 SAVE moles
end

KINETICS 1
Augite
formula FeMgSi2O6 1
m 2.52
m0 2.52
parms 3.974
tol 1e-09

Albite
formula NaAlSi3O8 1
m 1.27
m0 1.27
parms 5.26
tol 1e-08

Hematite
formula Fe2O3 1
m 0.53
m0 0.53
parms 3672.87
tol 1e-08

Calcite
formula CaCO3 1
m 0
m0 0
parms 15.0135
tol 1e-08

I was able to find from another post the following Rates code for forsterite but it is not the same and I don't know how to convert it to my code.

Code: [Select]

RATES
Forsterite
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
11 a1=8.38E+04
12 E1=67206
13 n1=0.470
20 rem neutral solution parameters
21 a2=1.58E+03
22 E2=79000
30 rem base solution parameters
31 a3=1.00E-07
32 E3=56637
33 n2=-0.600
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("forsterite")
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("H+")^n2         #base rate expression
90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
I would very much appreciate your help on that so I can set up the code in the proper way.

Thank you very much in advance.

Kind regards,
Maria

[dlparkhurst]

These are all questions that you need to answer. However, because new capabilities for rate calculations were introduced in PHREEQC version 3.8, I will expand a bit.

First, the database KInec_v3.dat was released by Oelkers and others. This database has RATES definitions for many primary and secondary minerals. For this database, you do not have to write your own RATES definitions; you just need to write the KINETICS data block.

Second, the database phreeqc_rates.dat includes the new keywords RATE_PARAMETERS_PK and RATE_PARAMETERS_HERMANSKA (references in the database comments). The rate parameter tabulations for may minerals are included in these data blocks. However, it is necessary to write your own RATES definitions. The rate in mol/m^2/sec can be defined by the Basic functions RATE_PK and RATE_HERMANSKA, but it is necessary for you to add the definitions to account for surface area and affinity (often 1 - SR("Forsterite"), for example). If you are working at pressure, you should include the molar volume for the minerals of interest so that the equilibrium constant is adjusted correctly for pressure effects.

Here is an example for Forsterite dissolution at 100 atm of CO2(g) pressure. You can adjust the script as necessary for your needs. The RATE_PARAMETERS_PK data block is included for clarity, but is not necessary if you want the parameters defined in the database. See the primary reference for the meaning of the parameters. Rates calculated by the Hermanska parameters may be defined analogously using the Basic function RATE_HERMANSKA. More examples are included in the comments in the phreeqc_rates.dat database.

Code: [Select]

PHASES
Forsterite
Mg2SiO4 + 4 H+ = 2 Mg+2 + H4SiO4
log_k 28.609
delta_h -217.115 #kJ/mol #Internal calculation
-analytic -1.0983766E3 -1.5385695E-1 7.321503E4 3.91599E2 -3.7061609E6
      [b]-Vm 42.6 # mw/rho = 140.7 / 3.3 = 42.6[/b]
RATE_PARAMETERS_PK
#               Acid                    Neutral         Base
#               log K   E       n(H+)   log K   E       log K   E       n(OH-)
#               ================================================================
Forsterite      -6.85   67.2    0.47    -10.64  79      -30     0       0       # Table 23
RATES
Forsterite_PK # Palandri and Kharaka, 2004
-start
5  REM PARMS: 1 m^2/mol, 2 exponent
10 affinity = 1 - SR("Forsterite")
20 rate = RATE_PK("Albite")
30 IF M > 0 THEN area = M * parm(1) * (M/M0)^parm(2) ELSE area = 0
40 SAVE area * rate * affinity * TIME
 -end
END
SOLUTION 1
C(4) 1 CO2(g) 2.0 # 100 atm
END
REACTION_PRESSURE 1
100 # atm
END
KINETICS 1
Forsterite_PK
  -formula Mg2SiO4
  -parms 1 0.67
  -m0 1
  -time 15*3.15e10
END
INCREMENTAL_REACTIONS true
END
RUN_CELLS
-cell 1
USER_GRAPH 1
    -headings               time Mg SI(Forsterite)
    -axis_titles            "Time, years" "Molality" "SI"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3.15e7
20 GRAPH_Y TOT("Mg")
30 GRAPH_SY SI("Forsterite")
40 END
  -end
    -active                 true
END