Serpentinization

[Olivine]

Hi,

I am trying to model the Serpentinization of Olivine (at 200-300C, 500bar) to assess the role of temperature and grain size, and estimate the precipitation of products to compare to experimental data.

Since i am new to this, i have first attempted to model the process at standard conditions so i could take advantage of the Phases and Rates blocks provided by (Zhang et al. 2019). My attempt is below, however i end up with "Warning negative moles in solution" and it will not finish Kinetic step 1.

Any advice is much appreciated!

Code: [Select]

#Phases and kinetics from (Zhang et al. 2019)

PHASES 1
Forsterite
    Mg2Si1O4 + 4H+ = 2H2O + 2Mg+2 + SiO2
    -analytical_expression -4970.994 -1.223246 243993.9 1880.021 -12003840 0.0003271594
    -Vm       43.66 cm3/mol
Fayalite
    Fe2Si1O4 + 4H+ = 2Fe+2 + 2H2O + SiO2
    -analytical_expression -3539.104 -0.8316186 178241 1330.28 -8977947 0.000207672
    -Vm       46.31 cm3/mol
Lizardite
    Mg3Si2O9H4 + 6H+ = 5H2O + 3Mg+2 + 2SiO2
    -analytical_expression -6396.55 -1.496386 319194.2 2405.399 -16184800 0.0003771849
    -Vm       106.45 cm3/mol
Magnetite
    Fe3O4 + 8H+ = Fe+2 + 2Fe+3 + 4H2O
    -analytical_expression -17446.4 -4.73879 781485.5 6686.528 -37407500 0.00140278
    -Vm       44.52 cm3/mol
Brucite
    Mg1O2H2 + 2H+ = 2H2O + Mg+2
    -analytical_expression -3594.082 -0.962404 165332.5 1376.033 -7736330 0.0002818074

RATES 1
   Forsterite
-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
 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 
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

   Fayalite
-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
 11 a1=5.48E+11
 12 E1=94400
 13 n1=0
 20 rem neutral solution parameters
 21 a2=5.48E+02
 22 E2=94400
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Fayalite")
 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 
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

    Lizardite
-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
 11 a1=3.39E+07
 12 E1=75516
 13 n1=0.8
 20 rem neutral solution parameters
 21 a2=3.33E-03
 22 E2=56637
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Lizardite")
 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 
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

    Brucite
-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
 11 a1=4.00E+05 
 12 E1=59000
 13 n1=0.500
 20 rem neutral solution parameters
 21 a2=1.30E-01
 22 E2=42000
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("brucite")
 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 
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

    Magnetite
-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
 11 a1=4.66E-06
 12 E1=18600
 13 n1=0.279
 20 rem neutral solution parameters
 21 a2=3.01E-08
 22 E2=18600
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Magnetite")
 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 
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

SOLUTION 1
 temp 70
 units mol/kgw
 ph 6
 Mg 0
 Fe 0

KINETICS 1
Forsterite
      -M0 0.1
      -parms 1e3 0.1
Fayalite
      -M0 0.9
-parms 1e3 0.1
Brucite
      -M0 0
-parms 1e3 0.1
Magnetite
      -M0 0
-parms 1e3 0.1
Lizardite
      -M0 0
      -parms 1e3 0.1
-steps 2000 year in 25 steps
  -step_divide 100


INCREMENTAL_REACTIONS TRUE

USER_GRAPH 1
 -chart_title "F080 dissolution"
 -axis_titles Time "mol", "Year"
 -initial_solutions false
 -start
 10 graph_x total_time
 20 graph_y KIN("Fayalite")
 20 graph_y KIN("Forsterite")
 30 graph_y KIN("Brucite")
 40 graph_y KIN("Lizardite")
 50 graph_y KIN("Magnetite")
-end
END



[dlparkhurst]

You need to use the -cvode option in KINETICS.

Your reaction reaches equilibrium in about 20 days with the RATES and KINETICS definitions that you have. Note that the use of both units and steps in KINETICS does not seem to work; best to simply define the time steps in seconds.

Code: [Select]

#Phases and kinetics from (Zhang et al. 2019)

PHASES 1
Forsterite
    Mg2Si1O4 + 4H+ = 2H2O + 2Mg+2 + SiO2
    -analytical_expression -4970.994 -1.223246 243993.9 1880.021 -12003840 0.0003271594
    -Vm       43.66 cm3/mol
Fayalite
    Fe2Si1O4 + 4H+ = 2Fe+2 + 2H2O + SiO2
    -analytical_expression -3539.104 -0.8316186 178241 1330.28 -8977947 0.000207672
    -Vm       46.31 cm3/mol
Lizardite
    Mg3Si2O9H4 + 6H+ = 5H2O + 3Mg+2 + 2SiO2
    -analytical_expression -6396.55 -1.496386 319194.2 2405.399 -16184800 0.0003771849
    -Vm       106.45 cm3/mol
Magnetite
    Fe3O4 + 8H+ = Fe+2 + 2Fe+3 + 4H2O
    -analytical_expression -17446.4 -4.73879 781485.5 6686.528 -37407500 0.00140278
    -Vm       44.52 cm3/mol
Brucite
    Mg1O2H2 + 2H+ = 2H2O + Mg+2
    -analytical_expression -3594.082 -0.962404 165332.5 1376.033 -7736330 0.0002818074

RATES 1
   Forsterite
-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
 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
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

   Fayalite
-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
 11 a1=5.48E+11
 12 E1=94400
 13 n1=0
 20 rem neutral solution parameters
 21 a2=5.48E+02
 22 E2=94400
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Fayalite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

    Lizardite
-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
 11 a1=3.39E+07
 12 E1=75516
 13 n1=0.8
 20 rem neutral solution parameters
 21 a2=3.33E-03
 22 E2=56637
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Lizardite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

    Brucite
-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
 11 a1=4.00E+05
 12 E1=59000
 13 n1=0.500
 20 rem neutral solution parameters
 21 a2=1.30E-01
 22 E2=42000
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("brucite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

    Magnetite
-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
 11 a1=4.66E-06
 12 E1=18600
 13 n1=0.279
 20 rem neutral solution parameters
 21 a2=3.01E-08
 22 E2=18600
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Magnetite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)               
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2   
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

SOLUTION 1
 temp 70
 units mol/kgw
 ph 6
 Mg 0
 Fe 0

KINETICS 1
Forsterite
      -M0 0.1
      -parms 1e3 0.1
Fayalite
      -M0 0.9
-parms 1e3 0.1
Brucite
      -M0 0
-parms 1e3 0.1
Magnetite
      -M0 0
-parms 1e3 0.1
Lizardite
      -M0 0
      -parms 1e3 0.1
-steps 20*86400
-cvode true


INCREMENTAL_REACTIONS TRUE

USER_GRAPH 1
 -chart_title "F080 dissolution"
 -axis_titles Day "mol"
 -initial_solutions false
 -start
 10 graph_x total_time/86400
 20 graph_y KIN_DELTA("Fayalite")
 20 graph_y KIN_DELTA("Forsterite")
 30 graph_y KIN_DELTA("Brucite")
 40 graph_y KIN_DELTA("Lizardite")
 50 graph_y KIN_DELTA("Magnetite")
-end
END