Equilibrium problem in the kinetic modeling

[Hamed]

Dear all,
I am working on a model considering the reactive transport of CO2 in porous media. Specifically, I want to check the porosity change. To obtain the brine composition, I first equilibrated it with 1.7 M NaCl and minerals present in the medium. I have both anorthite and laumontite in the mineralogy, and this causes the undersaturation of anorthite in the initial equilibrium. Thus, in the reactive transport model, the porosity of the cells that are not exposed to CO2 yet does not reach the initial porosity value of 0.1 over time. I don't know how I can manage this problem. I would appreciate any guidance related to this.

Code: [Select]

PRINT; -reset false

PHASES


Dawsonite
        NaAlCO3(OH)2 +3.0000 H+  =  + 1.0000 Al+++ + 1.0000 HCO3- + 1.0000 Na+ + 2.0000 H2O
        log_k           4.3464
-delta_H -76.3549 kJ/mol # Calculated enthalpy of reaction Dawsonite
        -analytic -1.1393e+002 -2.3487e-002 7.1758e+003 4.0900e+001 1.2189e+002
        -Vm  59

Smectite
       Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12 +8.0000 H+ + 2.0000 H2O =  + 0.0250 Ca++ + 0.1000 Na+ + 0.2000 Fe+++ + 0.2000 K+ + 0.5000 Fe++ + 1.1500 Mg++ + 1.2500 Al+++ + 3.5000H4SiO4
        log_k           17.4200
-delta_H -199.841 kJ/mol # Calculated enthalpy of reaction
        -analytic -9.6102e+000 1.2551e-003 1.8157e+004 -7.9862e+000 -1.3005e+006
         -Vm  140

Laumontite
        CaAl2Si4O12:4H2O +8.0000 H+  =  + 1.0000 Ca++ + 2.0000 Al+++ + 4.0000 H4SiO4 #+ 8.0000 H2O
        log_k           13.6667
-delta_H -184.657 kJ/mol # Calculated enthalpy of reaction Laumontite
        -analytic 1.1904e+000 8.1763e-003 1.9005e+004 -1.4561e+001 -1.5851e+006
        -Vm  207.53

CO2(g)
CO2 = CO2
-log_k -1.468
-delta_h -4.776 kcal
-analytic   10.5624  -2.3547e-2  -3972.8  0  5.8746e5  1.9194e-5
-T_c  304.2  # critical T, K
-P_c   72.86 # critical P, atm
-Omega 0.225 # acentric factor

END

RATES

    Laumontite
-start
10 SR_mineral=SR("Laumontite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Laumontite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


    Anorthite
-start
10 SR_mineral=SR("anorthite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Anorthite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Albite
-start
10 SR_mineral=SR("Albite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*(act("H+")^(0.457))  + k_neut + k_base*act("H+")^(-0.572)
100 r = (k_rateconst) * SA * (1-SR("Albite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Kaolinite
-start
10 SR_mineral=SR("Kaolinite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 4.9e-12*EXP((-65900/8.3145)*(1/Tk-1/298.15))
70 k_neut = 6.92e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
80 k_base = 8.91e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
100 r = k_rateconst * SA * (1-SR("Kaolinite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

   Illite
-start
10 SR_mineral=SR("Illite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.05e-11*EXP(-23600/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+")^0.34
100 r = k_rateconst * SA * (1-SR("Illite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


   Smectite
-start
10 SR_mineral=SR("Smectite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
100 r = 5.62e-14* SA * (1-SR("Smectite"))*PARM(2)*act("H+")^0.5
150 moles = r * TIME
200 SAVE moles
-end

  Dawsonite
-start
10 SR_mineral=SR("Dawsonite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.457e-4*EXP(-36100/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 1.26e-9*EXP(-62760/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut 
100 r = k_rateconst * SA * (1-SR("Dawsonite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


END

SOLUTION 0

-temperature 60
-pH 7
-units   mol/kgw

GAS_PHASE 0
      -fixed_pressure
      -pressure      250
      -volume        2.34
      -temperature   60           
     CO2(g)          250

SAVE solution 0
END

SOLUTION 1-50     

-temperature 60
-pH 7
-units   mol/kgw

Al  6.370e-04
Ca  7.924e-04
Cl  1.700e+00  charge
Fe  5.638e-09
K   7.000e-01
Mg  9.263e-09
Na  9.998e-01
Si  1.845e-04

KINETICS 1-50

Albite
    -m        1.6
    -m0       1.6
    -parms    29.23 0.001
    -tol      1e-06

Anorthite
    -m        1.75
    -m0       1.75
    -parms    35 0.001
    -tol      1e-06

Laumontite
    -m        3.35
    -m0       3.35
    -parms    110.45 0.001
    -tol      1e-06


Smectite
    -m        21.81
    -m0       21.81
    -parms    8724  0.001
    -tol      1e-06

Kaolinite
    -m        6.71
    -m0       6.71
    -parms    2008 0.001
    -tol      1e-06

Dawsonite
    -m       1e-9
    -m0      1e-9
    -parms   1 1
    -tol      1e-06

-cvode true
-bad_step_max 5000
INCREMENTAL_REACTIONS true
END

   
TRANSPORT
-cells 5
-shifts 5
-time_step  10 yr
-flow_direction diffusion_only
-boundary_conditions constant constant
-diffusion_coefficient  5e-10
-lengths  5*1
-punch_frequency   1
-fix_current 1e-8
-multi_D true 1e-9 0.10 0.01 1 true
-porosities    5*0.10
-implicit true 1
-punch 1-5

SELECTED_OUTPUT 1
         -reset           false

USER_PUNCH 1
-headings TIME Cell Porosity pH
-start
10 PUNCH TOTAL_TIME
20 PUNCH CELL_NO
30 PUNCH GET_POR(CELL_NO) 
40 PUNCH -La("H+")

100 REM calculate the porosity
110 dv = KIN_DELTA("Anorthite")*PHASE_VM("Anorthite") + KIN_DELTA("Albite")*PHASE_VM("Albite") + KIN_DELTA("Laumontite")*PHASE_VM("Laumontite")+ KIN_DELTA("Kaolinite")*PHASE_VM("Kaolinite") + KIN_DELTA("Smectite")*PHASE_VM("Smectite")+ KIN_DELTA("Dawsonite")*PHASE_VM("Dawsonite")
1110 Vtot = 9000
1120 dPor =  dv/Vtot           
1130 new_por = GET_POR(CELL_NO) - dPor                 
1140 CHANGE_POR(new_por,cell_no)                                                       
 -end
 
USER_GRAPH 1
    -headings  Dis Por
    -axis_titles            "Distance" "Porosity"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 GRAPH_X  DIST
20 GRAPH_Y  GET_POR(CELL_NO)
  -end
    -active                 true
END

[dlparkhurst]

First, your initial solutions 1-50 are at 1 atm pressure. Shouldn't the pressure be 250 atm?

Second, the initial solution does not seem to be in equilibrium with the minerals of the formation, even though you state "I first equilibrated it with 1.7 M NaCl and minerals present in the medium".

If the minerals are not in equilibrium with the initial solution, then your kinetic reactions will generate a porosity change whether CO2 is introduced or not. The following is the way I interpret your description, but you can tell me what is wrong with it.

Code: [Select]

#PRINT; -reset false

PHASES


Dawsonite
        NaAlCO3(OH)2 +3.0000 H+  =  + 1.0000 Al+++ + 1.0000 HCO3- + 1.0000 Na+ + 2.0000 H2O
        log_k           4.3464
-delta_H -76.3549 kJ/mol # Calculated enthalpy of reaction Dawsonite
        -analytic -1.1393e+002 -2.3487e-002 7.1758e+003 4.0900e+001 1.2189e+002
        -Vm  59

Smectite
       Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12 +8.0000 H+ + 2.0000 H2O =  + 0.0250 Ca++ + 0.1000 Na+ + 0.2000 Fe+++ + 0.2000 K+ + 0.5000 Fe++ + 1.1500 Mg++ + 1.2500 Al+++ + 3.5000H4SiO4
        log_k           17.4200
-delta_H -199.841 kJ/mol # Calculated enthalpy of reaction
        -analytic -9.6102e+000 1.2551e-003 1.8157e+004 -7.9862e+000 -1.3005e+006
         -Vm  140

Laumontite
        CaAl2Si4O12:4H2O +8.0000 H+  =  + 1.0000 Ca++ + 2.0000 Al+++ + 4.0000 H4SiO4 #+ 8.0000 H2O
        log_k           13.6667
-delta_H -184.657 kJ/mol # Calculated enthalpy of reaction Laumontite
        -analytic 1.1904e+000 8.1763e-003 1.9005e+004 -1.4561e+001 -1.5851e+006
        -Vm  207.53

CO2(g)
CO2 = CO2
-log_k -1.468
-delta_h -4.776 kcal
-analytic   10.5624  -2.3547e-2  -3972.8  0  5.8746e5  1.9194e-5
-T_c  304.2  # critical T, K
-P_c   72.86 # critical P, atm
-Omega 0.225 # acentric factor

END

RATES

    Laumontite
-start
10 SR_mineral=SR("Laumontite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Laumontite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


    Anorthite
-start
10 SR_mineral=SR("anorthite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Anorthite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Albite
-start
10 SR_mineral=SR("Albite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*(act("H+")^(0.457))  + k_neut + k_base*act("H+")^(-0.572)
100 r = (k_rateconst) * SA * (1-SR("Albite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Kaolinite
-start
10 SR_mineral=SR("Kaolinite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 4.9e-12*EXP((-65900/8.3145)*(1/Tk-1/298.15))
70 k_neut = 6.92e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
80 k_base = 8.91e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
100 r = k_rateconst * SA * (1-SR("Kaolinite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

   Illite
-start
10 SR_mineral=SR("Illite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.05e-11*EXP(-23600/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+")^0.34
100 r = k_rateconst * SA * (1-SR("Illite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


   Smectite
-start
10 SR_mineral=SR("Smectite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
100 r = 5.62e-14* SA * (1-SR("Smectite"))*PARM(2)*act("H+")^0.5
150 moles = r * TIME
200 SAVE moles
-end

  Dawsonite
-start
10 SR_mineral=SR("Dawsonite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.457e-4*EXP(-36100/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 1.26e-9*EXP(-62760/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut
100 r = k_rateconst * SA * (1-SR("Dawsonite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


END

SOLUTION 0

-temperature 60
-pH 7
-units   mol/kgw

GAS_PHASE 0
      -fixed_pressure
      -pressure      250
      -volume        2.34
      -temperature   60           
     CO2(g)          250

SAVE solution 0
END

SOLUTION 1     

-temperature 60
-pressure 250
-pH 7
-units   mol/kgw

Al  6.370e-04
Ca  7.924e-04
Cl  1.700e+00  charge
Fe  5.638e-09
K   7.000e-01
Mg  9.263e-09
Na  9.998e-01
Si  1.845e-04
END
USE solution 1
EQUILIBRIUM_PHASES 101
Albite
Anorthite
Laumontite
Smectite
Kaolinite
SAVE solution 1-50
END
KINETICS 1-50

Albite
    -m        1.6
    -m0       1.6
    -parms    29.23 0.001
    -tol      1e-06

Anorthite
    -m        1.75
    -m0       1.75
    -parms    35 0.001
    -tol      1e-06

Laumontite
    -m        3.35
    -m0       3.35
    -parms    110.45 0.001
    -tol      1e-06


Smectite
    -m        21.81
    -m0       21.81
    -parms    8724  0.001
    -tol      1e-06

Kaolinite
    -m        6.71
    -m0       6.71
    -parms    2008 0.001
    -tol      1e-06

Dawsonite
    -m       1e-9
    -m0      1e-9
    -parms   1 1
    -tol      1e-06
-cvode true
-bad_step_max 5000
INCREMENTAL_REACTIONS true
END

   
TRANSPORT
-cells 5
-shifts 2
-time_step  10 yr
-flow_direction diffusion_only
-boundary_conditions constant constant
-diffusion_coefficient  5e-10
-lengths  5*1
-punch_frequency   1
-fix_current 1e-8
-multi_D true 1e-9 0.10 0.01 1 true
-porosities    5*0.10
-implicit true 1
-punch_cells 1-5
-print_cells 1-5

SELECTED_OUTPUT 1
         -reset           false

USER_PUNCH 1
-headings TIME Cell Porosity pH
-start
10 PUNCH TOTAL_TIME
20 PUNCH CELL_NO
30 PUNCH GET_POR(CELL_NO)
40 PUNCH -La("H+")

100 REM calculate the porosity
110 dv = KIN_DELTA("Anorthite")*PHASE_VM("Anorthite") + KIN_DELTA("Albite")*PHASE_VM("Albite") + KIN_DELTA("Laumontite")*PHASE_VM("Laumontite")+ KIN_DELTA("Kaolinite")*PHASE_VM("Kaolinite") + KIN_DELTA("Smectite")*PHASE_VM("Smectite")+ KIN_DELTA("Dawsonite")*PHASE_VM("Dawsonite")
1110 Vtot = 9000
1120 dPor =  dv/Vtot           
1130 new_por = GET_POR(CELL_NO) - dPor                 
1140 CHANGE_POR(new_por,cell_no)                                                       
 -end
 
USER_GRAPH 1
    -headings  Dis Por
    -axis_titles            "Distance" "Porosity"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 GRAPH_X  DIST
20 GRAPH_Y  GET_POR(CELL_NO)
  -end
    -active                 true
END

Note that you should put an END statement between your SOLUTION and KINETICS definitions. Otherwise SOLUTION 1 will react with KINETICS 1, and, although SOLUTION 1 will not change, KINETICS 1 will change.

[Hamed]

Thank you, Dr. Parkhust. In the previous code, I did not include all the minerals. In the attached code, I have included all the minerals, and if you run that, you can see that porosity changes even without minerals being introduced to CO2. In the equilibrium phases, all anorthite and illite dissolve. I think the problem is related to these minerals. So, I commented them. Any guidance to manage this problem is appreciated.
Thank you in advance for your time,

Code: [Select]

#PRINT; -reset false

PHASES

Dawsonite
        NaAlCO3(OH)2 +3.0000 H+  =  + 1.0000 Al+++ + 1.0000 HCO3- + 1.0000 Na+ + 2.0000 H2O
        log_k           4.3464
-delta_H -76.3549 kJ/mol # Calculated enthalpy of reaction Dawsonite
        -analytic -1.1393e+002 -2.3487e-002 7.1758e+003 4.0900e+001 1.2189e+002
        -Vm  59

Smectite
       Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12 +8.0000 H+ + 2.0000 H2O =  + 0.0250 Ca++ + 0.1000 Na+ + 0.2000 Fe+++ + 0.2000 K+ + 0.5000 Fe++ + 1.1500 Mg++ + 1.2500 Al+++ + 3.5000H4SiO4
        log_k           17.4200
-delta_H -199.841 kJ/mol # Calculated enthalpy of reaction
        -analytic -9.6102e+000 1.2551e-003 1.8157e+004 -7.9862e+000 -1.3005e+006
         -Vm  140

Laumontite
        CaAl2Si4O12:4H2O +8.0000 H+  =  + 1.0000 Ca++ + 2.0000 Al+++ + 4.0000 H4SiO4 #+ 8.0000 H2O
        log_k           13.6667
-delta_H -184.657 kJ/mol # Calculated enthalpy of reaction Laumontite
        -analytic 1.1904e+000 8.1763e-003 1.9005e+004 -1.4561e+001 -1.5851e+006
        -Vm  207.53

Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
 log_k -17.4
 delta_h 6.98 kJ
-analytic -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
 -Vm  69

CO2(g)
CO2 = CO2
-log_k -1.468
-delta_h -4.776 kcal
-analytic   10.5624  -2.3547e-2  -3972.8  0  5.8746e5  1.9194e-5
-T_c  304.2  # critical T, K
-P_c   72.86 # critical P, atm
-Omega 0.225 # acentric factor

END

RATES

    Laumontite
-start
10 SR_mineral=SR("Laumontite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Laumontite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


    Anorthite
-start
10 SR_mineral=SR("anorthite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Anorthite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

  Ankerite
-start
10 SR_mineral=SR("Ankerite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.59e-4*EXP((-45000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 1.26e-9*EXP((-62800/8.3145)*(1/Tk-1/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut
100 r = (k_rateconst) * SA * (1-SR("Ankerite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

    K-feldspar
-start
10 SR_mineral=SR("K-feldspar")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))
70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.823)
100 r = k_rateconst * SA * (1-SR("K-Feldspar"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Albite
-start
10 SR_mineral=SR("Albite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*(act("H+")^(0.457))  + k_neut + k_base*act("H+")^(-0.572)
100 r = (k_rateconst) * SA * (1-SR("Albite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Kaolinite
-start
10 SR_mineral=SR("Kaolinite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 4.9e-12*EXP((-65900/8.3145)*(1/Tk-1/298.15))
70 k_neut = 6.92e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
80 k_base = 8.91e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
100 r = k_rateconst * SA * (1-SR("Kaolinite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

   Illite
-start
10 SR_mineral=SR("Illite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.05e-11*EXP(-23600/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+")^0.34
100 r = k_rateconst * SA * (1-SR("Illite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


   Smectite
-start
10 SR_mineral=SR("Smectite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
100 r = 5.62e-14* SA * (1-SR("Smectite"))*PARM(2)*act("H+")^0.5
150 moles = r * TIME
200 SAVE moles
-end


Siderite
-start
10 SR_mineral=SR("Siderite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.59e-4*EXP((-45000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 1.26e-9*EXP((-62800/8.3145)*(1/Tk-1/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut
100 r = (k_rateconst) * SA * (1-SR("Siderite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

    Chalcedony
-start
10 SR_mineral=SR("Chalcedony")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
70 k_neut = 1.26e-14*EXP((-87500/8.3145)*(1/TK-1/298.15))
90 k_rateconst = k_neut
100 r = k_rateconst * SA * (1-SR("Chalcedony"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


Calcite
-start
10 SR_mineral=SR("Calcite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 5e-1*EXP(-14400/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 1.55e-6*EXP(-23500/8.314*(1.0/TK-1.0/298.15))
80 k_base = 3.31e-4*EXP(-35400/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+") + k_neut + k_base*act("H+")
100 r = k_rateconst * SA * (1-SR("Calcite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

  Pyrite
-start
10 SR_mineral=SR("Pyrite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
70 k_acid = 4e-11*EXP(-62760/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+")^0.5
100 r = k_rateconst * SA * (1-SR("Pyrite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

  Dawsonite
-start
10 SR_mineral=SR("Dawsonite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.457e-4*EXP(-36100/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 1.26e-9*EXP(-62760/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut
100 r = k_rateconst * SA * (1-SR("Dawsonite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


END

SOLUTION 0

-temperature 60
-pH 7
-units   mol/kgw

GAS_PHASE 0
      -fixed_pressure
      -pressure      250
      -volume        2.34
      -temperature   60           
     CO2(g)          250

SAVE solution 0
END

SOLUTION 1     

-temperature 60
-pressure 250
-pH 7
-units   mol/kgw

Na 1.7
Cl 1.7 charge

END
USE solution 1
EQUILIBRIUM_PHASES 101
Albite
#Anorthite
Laumontite
Chalcedony
Pyrite
#Illite
Smectite
Kaolinite
Calcite
Siderite
K-feldspar
SAVE solution 1-50
END
KINETICS 1-50

Albite
    -m        1.6
    -m0       1.6
    -parms    29.23 0.001
    -tol      1e-06

Anorthite
    -m        1.75
    -m0       1.75
    -parms    35 0.001
    -tol      1e-06

Laumontite
    -m        3.35
    -m0       3.35
    -parms    110.45 0.001
    -tol      1e-06

Chalcedony
    -m       62.77
    -m0      62.77
    -parms    258.34 0.001
    -tol      1e-06

Illite
    -m        14.44
    -m0       14.44
    -parms    2628 0.001
    -tol      1e-06

Pyrite
    -m        6.09
    -m0       6.09
    -parms    26.52 0.001
    -tol      1e-06

Smectite
    -m        21.81
    -m0       21.81
    -parms    8724  0.001
    -tol      1e-06

Calcite
    -m        1.5
    -m0       1.5
    -parms    10.5 0.001
    -tol      1e-06

Siderite
    -m        9.96
    -m0       9.96
    -parms    53.24 0.001
    -tol      1e-06

Kaolinite
    -m        6.71
    -m0       6.71
    -parms    2008 0.001
    -tol      1e-06

K-feldspar
    -m        0.87
    -m0       0.87
    -parms    17.19 0.001
    -tol      1e-06

Dawsonite
    -m       1e-9
    -m0      1e-9
    -parms   1 1
    -tol      1e-06

Ankerite
    -m       1e-9
    -m0      1e-9
    -parms   1 1
    -tol      1e-06

-cvode true
-bad_step_max 5000
INCREMENTAL_REACTIONS true
END

   
TRANSPORT
-cells 5
-shifts 3
-time_step  5 yr
-flow_direction diffusion_only
-boundary_conditions constant constant
-diffusion_coefficient  5e-10
-lengths  5*1
-punch_frequency   1
-fix_current 1e-8
-multi_D true 1e-9 0.10 0.01 1 true
-porosities    5*0.10
-implicit true 1
-punch_cells 1-5
-print_cells 1-5

SELECTED_OUTPUT 1
         -reset           false

USER_PUNCH 1
-headings TIME Cell Porosity pH
-start
10 PUNCH TOTAL_TIME
20 PUNCH CELL_NO
30 PUNCH GET_POR(CELL_NO)
40 PUNCH -La("H+")

100 REM calculate the porosity
110 dv = KIN_DELTA("Ankerite")*PHASE_VM("Ankerite")+ KIN_DELTA("Anorthite")*PHASE_VM("Anorthite") + KIN_DELTA("Siderite")*PHASE_VM("Siderite") + KIN_DELTA("K-feldspar")*PHASE_VM("K-feldspar") + KIN_DELTA("Albite")*PHASE_VM("Albite") + KIN_DELTA("Laumontite")*PHASE_VM("Laumontite")+ KIN_DELTA("Kaolinite")*PHASE_VM("Kaolinite") + KIN_DELTA("Illite")*PHASE_VM("Illite") + KIN_DELTA("Calcite")*PHASE_VM("Calcite") + KIN_DELTA("Pyrite")*PHASE_VM("Pyrite") + KIN_DELTA("Smectite")*PHASE_VM("Smectite")+ KIN_DELTA("Chalcedony")*PHASE_VM("Chalcedony") + KIN_DELTA("Dawsonite")*PHASE_VM("Dawsonite")
1110 Vtot = 8000
1120 dPor =  dv/Vtot           
1130 new_por = GET_POR(CELL_NO) - dPor                 
1140 CHANGE_POR(new_por,cell_no)                                                       
 -end
 
USER_GRAPH 1
    -headings  Dis Por
    -axis_titles            "Distance" "Porosity"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 GRAPH_X  DIST
20 GRAPH_Y  GET_POR(CELL_NO)
  -end
    -active                 true
END

[dlparkhurst]

If your initial water composition for transport is not in equilibrium with anorthite and laumontite, then they will react whether CO2 is added or not. With your rate definitions, the anorthite reacts quickly in the absence of CO2, which causes the change in porosity.

If I include all of the minerals in equilibrium phases,  I find that chalcedony and illite are removed completely (using phreeqc.dat). If you equilibrate with all but those two, then the porosity changes much more slowly in the absence of CO2.

I'm not saying there is a complete answer to your problem. You have more minerals than can coexist in thermodynamic equilibrium according to the Gibbs' Phase Rule. If you use the whole set, two minerals will dissolve completely. So, regardless of the set that you initially equilibrate that is less than the whole set, the other minerals will be out of equilibrium and will react kinetically causing other minerals to react kinetically. Regardless of your choices, the porosity will change. I have just switched faster reacting minerals for slower reacting minerals.

Code: [Select]

#PRINT; -reset false

PHASES

Dawsonite
        NaAlCO3(OH)2 +3.0000 H+  =  + 1.0000 Al+++ + 1.0000 HCO3- + 1.0000 Na+ + 2.0000 H2O
        log_k           4.3464
-delta_H -76.3549 kJ/mol # Calculated enthalpy of reaction Dawsonite
        -analytic -1.1393e+002 -2.3487e-002 7.1758e+003 4.0900e+001 1.2189e+002
        -Vm  59

Smectite
       Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12 +8.0000 H+ + 2.0000 H2O =  + 0.0250 Ca++ + 0.1000 Na+ + 0.2000 Fe+++ + 0.2000 K+ + 0.5000 Fe++ + 1.1500 Mg++ + 1.2500 Al+++ + 3.5000H4SiO4
        log_k           17.4200
-delta_H -199.841 kJ/mol # Calculated enthalpy of reaction
        -analytic -9.6102e+000 1.2551e-003 1.8157e+004 -7.9862e+000 -1.3005e+006
         -Vm  140

Laumontite
        CaAl2Si4O12:4H2O +8.0000 H+  =  + 1.0000 Ca++ + 2.0000 Al+++ + 4.0000 H4SiO4 #+ 8.0000 H2O
        log_k           13.6667
-delta_H -184.657 kJ/mol # Calculated enthalpy of reaction Laumontite
        -analytic 1.1904e+000 8.1763e-003 1.9005e+004 -1.4561e+001 -1.5851e+006
        -Vm  207.53

Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
 log_k -17.4
 delta_h 6.98 kJ
-analytic -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
 -Vm  69

CO2(g)
CO2 = CO2
-log_k -1.468
-delta_h -4.776 kcal
-analytic   10.5624  -2.3547e-2  -3972.8  0  5.8746e5  1.9194e-5
-T_c  304.2  # critical T, K
-P_c   72.86 # critical P, atm
-Omega 0.225 # acentric factor

END

RATES

    Laumontite
-start
10 SR_mineral=SR("Laumontite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Laumontite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


    Anorthite
-start
10 SR_mineral=SR("anorthite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 3.16e-4*EXP(-16600/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 7.59e-10*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^1.411 + k_neut
100 r = k_rateconst * SA * (1-SR("Anorthite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

  Ankerite
-start
10 SR_mineral=SR("Ankerite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.59e-4*EXP((-45000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 1.26e-9*EXP((-62800/8.3145)*(1/Tk-1/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut
100 r = (k_rateconst) * SA * (1-SR("Ankerite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

    K-feldspar
-start
10 SR_mineral=SR("K-feldspar")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))
70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.823)
100 r = k_rateconst * SA * (1-SR("K-Feldspar"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Albite
-start
10 SR_mineral=SR("Albite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*(act("H+")^(0.457))  + k_neut + k_base*act("H+")^(-0.572)
100 r = (k_rateconst) * SA * (1-SR("Albite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

Kaolinite
-start
10 SR_mineral=SR("Kaolinite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 4.9e-12*EXP((-65900/8.3145)*(1/Tk-1/298.15))
70 k_neut = 6.92e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
80 k_base = 8.91e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
100 r = k_rateconst * SA * (1-SR("Kaolinite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

   Illite
-start
10 SR_mineral=SR("Illite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.05e-11*EXP(-23600/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+")^0.34
100 r = k_rateconst * SA * (1-SR("Illite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


   Smectite
-start
10 SR_mineral=SR("Smectite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
100 r = 5.62e-14* SA * (1-SR("Smectite"))*PARM(2)*act("H+")^0.5
150 moles = r * TIME
200 SAVE moles
-end


Siderite
-start
10 SR_mineral=SR("Siderite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 1.59e-4*EXP((-45000/8.3145)*(1/Tk-1/298.15))
70 k_neut = 1.26e-9*EXP((-62800/8.3145)*(1/Tk-1/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut
100 r = (k_rateconst) * SA * (1-SR("Siderite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

    Chalcedony
-start
10 SR_mineral=SR("Chalcedony")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
70 k_neut = 1.26e-14*EXP((-87500/8.3145)*(1/TK-1/298.15))
90 k_rateconst = k_neut
100 r = k_rateconst * SA * (1-SR("Chalcedony"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


Calcite
-start
10 SR_mineral=SR("Calcite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 5e-1*EXP(-14400/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 1.55e-6*EXP(-23500/8.314*(1.0/TK-1.0/298.15))
80 k_base = 3.31e-4*EXP(-35400/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+") + k_neut + k_base*act("H+")
100 r = k_rateconst * SA * (1-SR("Calcite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

  Pyrite
-start
10 SR_mineral=SR("Pyrite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
70 k_acid = 4e-11*EXP(-62760/8.314*(1.0/TK-1.0/298.15))
90 k_rateconst = k_acid*act("H+")^0.5
100 r = k_rateconst * SA * (1-SR("Pyrite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end

  Dawsonite
-start
10 SR_mineral=SR("Dawsonite")
20 if (M<0) then goto 200
30 if (M=0 and SR_mineral<1) then goto 200
40 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 k_acid = 6.457e-4*EXP(-36100/8.314*(1.0/TK-1.0/298.15))
70 k_neut = 1.26e-9*EXP(-62760/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid*act("H+")^0.5 + k_neut
100 r = k_rateconst * SA * (1-SR("Dawsonite"))*PARM(2)
150 moles = r * TIME
200 SAVE moles
-end


END

SOLUTION 0

-temperature 60
-pH 7
-units   mol/kgw

GAS_PHASE 0
      -fixed_pressure
      -pressure      250
      -volume        2.34
      -temperature   60           
     CO2(g)          250

SAVE solution 0
END

SOLUTION 1     

-temperature 60
-pressure 250
-pH 7
-units   mol/kgw

Na 1.7
Cl 1.7 charge

END
USE solution 1
EQUILIBRIUM_PHASES 101
Albite
Anorthite
Laumontite
#Chalcedony
Pyrite
#Illite
Smectite
Kaolinite
Calcite
Siderite
K-feldspar
SAVE solution 1-50
END
KINETICS 1-50

Albite
    -m        1.6
    -m0       1.6
    -parms    29.23 0.001
    -tol      1e-06

Anorthite
    -m        1.75
    -m0       1.75
    -parms    35 0.001
    -tol      1e-06

C
    -m        3.35
    -m0       3.35
    -parms    110.45 0.001
    -tol      1e-06

Chalcedony
    -m       62.77
    -m0      62.77
    -parms    258.34 0.001
    -tol      1e-06

Illite
    -m        14.44
    -m0       14.44
    -parms    2628 0.001
    -tol      1e-06

Pyrite
    -m        6.09
    -m0       6.09
    -parms    26.52 0.001
    -tol      1e-06

Smectite
    -m        21.81
    -m0       21.81
    -parms    8724  0.001
    -tol      1e-06

Calcite
    -m        1.5
    -m0       1.5
    -parms    10.5 0.001
    -tol      1e-06

Siderite
    -m        9.96
    -m0       9.96
    -parms    53.24 0.001
    -tol      1e-06

Kaolinite
    -m        6.71
    -m0       6.71
    -parms    2008 0.001
    -tol      1e-06

K-feldspar
    -m        0.87
    -m0       0.87
    -parms    17.19 0.001
    -tol      1e-06

Dawsonite
    -m       1e-9
    -m0      1e-9
    -parms   1 1
    -tol      1e-06

Ankerite
    -m       1e-9
    -m0      1e-9
    -parms   1 1
    -tol      1e-06

-cvode true
-bad_step_max 5000
INCREMENTAL_REACTIONS true
END

   
TRANSPORT
-cells 5
-shifts 5
-time_step  5 yr
-flow_direction diffusion_only
-boundary_conditions constant constant
-diffusion_coefficient  5e-10
-lengths  5*1
-punch_frequency   1
-fix_current 1e-8
-multi_D true 1e-9 0.10 0.01 1 true
-porosities    5*0.10
-implicit true 1
-punch_cells 1-5
-print_cells 1-5

SELECTED_OUTPUT 1
         -reset           false

USER_PUNCH 1
-headings TIME Cell Porosity pH
-start
10 PUNCH TOTAL_TIME
20 PUNCH CELL_NO
30 PUNCH GET_POR(CELL_NO)
40 PUNCH -La("H+")

100 REM calculate the porosity
110 dv = KIN_DELTA("Ankerite")*PHASE_VM("Ankerite")+ KIN_DELTA("Anorthite")*PHASE_VM("Anorthite") + KIN_DELTA("Siderite")*PHASE_VM("Siderite") + KIN_DELTA("K-feldspar")*PHASE_VM("K-feldspar") + KIN_DELTA("Albite")*PHASE_VM("Albite") + KIN_DELTA("Laumontite")*PHASE_VM("Laumontite")+ KIN_DELTA("Kaolinite")*PHASE_VM("Kaolinite") + KIN_DELTA("Illite")*PHASE_VM("Illite") + KIN_DELTA("Calcite")*PHASE_VM("Calcite") + KIN_DELTA("Pyrite")*PHASE_VM("Pyrite") + KIN_DELTA("Smectite")*PHASE_VM("Smectite")+ KIN_DELTA("Chalcedony")*PHASE_VM("Chalcedony") + KIN_DELTA("Dawsonite")*PHASE_VM("Dawsonite")
1110 Vtot = 8000
1120 dPor =  dv/Vtot           
1130 new_por = GET_POR(CELL_NO) - dPor                 
1140 CHANGE_POR(new_por,cell_no)                                                       
 -end
 
USER_GRAPH 1
    -headings  Dis Por
    -axis_titles            "Distance" "Porosity"
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 GRAPH_X  DIST
20 GRAPH_Y  GET_POR(CELL_NO)
  -end
    -active                 true
END