High pCO2 values usage in CO2 capture script

[Erdinc]

Hello PHREEQC gurus,

I have been working on script to calculate the dissolution of minerals and CO2 capture through weathering of slag mineralization. I have run my simulations from 0.1% to up to 99% CO2. from 0.1 to around 6 % makes some sense but anything around 50% and 99% shows erratic readings which I believe could be PHREEQC sensitivity. Would you please let me know your thoughts? Many thanks.


Best Regards,
Erdinc

Code: [Select]

TITLE DATABASE C:\Users\ErdincCosgun\Documents\Projects\Databases\carbfix_2022.dat
TITLE TESTING KINETICS OF SPECIFIC SOLIDS
TITLE 100 CELLS, 1CM HIGH
TITLE FLOW-THROUGH TEST, 3 LITRE
TITLE === FULL ANTIMAGNETIC SUITE ===
TITLE KARDEMIR PROCESS WATER MEAN VALUES
TITLE REDOX CONTROLLED BY IRON SPECIES
TITLE CHEMICAL SPECIATION
TITLE PHASES FROM Carbfix2022
TITLE CALCITE DISSOLVES AS KINETIC PHASE
TITLE CALCITE PRECIPITATION AS EQUILIBRIUM PHASE
TITLE CO2 INJECTED FROM DAC 50% TO SOLUTION, log(0.5):-0.301
TITLE Li ADDED IN TRACE TO KARDEMIR SUPPLIED WATER

PHASES # Minerals, Directly taken from Carbfix2022.dat PHASES.

Calcite_diss # Calcite that dissolves only. Alphanumeric name of phase; no spaces are allowed.
CaCO3 + H+ = Ca+2 + HCO3- # Dissolution reaction
log_k 1.8487 # Log K at 25 ?C for the reaction
-delta_H -25.7149 kJ/mol # Enthalpy of reaction at 25 ?C for the reaction
# deltafH -288.552 kcal/mol
-analytic -1.4978e2 -4.8370e-2 4.8974e3 6.0458e1 7.6464e1 # Identifier for coefficients for an analytical expression for the temperature dependence of log K.
# Range 0-350
-Vm 36.934 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=100.09, Density (g/ml or g/cm3)= 2.71 therfore, Vm= 100.09 / 2.71 = 36.934
# Extrapol supcrt92
# Ref HDN+78

Calcite_ppt # Calcite that percipitates only
CaCO3 + H+ = Ca+2 + HCO3-
log_k 1.8487
-delta_H -25.7149 kJ/mol
# deltafH -288.552 kcal/mol
-analytic -1.4978e2 -4.8370e-2 4.8974e3 6.0458e1 7.6464e1
# Range 0-350
-Vm 36.934 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=100.09, Density (g/ml or g/cm3)= 2.71 therfore, Vm= 100.09 / 2.71 = 36.934
# Extrapol supcrt92
# Ref HDN+78

Brucite
Mg(OH)2 + 2 H+ = Mg+2 + 2 H2O
log_k 16.2980
-delta_H -111.34 kJ/mol
# deltafH -221.39 kcal/mol
-analytic -1.0280e2 -1.9759e-2 9.0180e3 3.8282e1 1.4075e2
# Range 0-350
-Vm 24.63 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=58.32, Density (g/ml or g/cm3)= 2.39 therfore, Vm= 58.32 / 2.39 = 24.40
# Extrapol supcrt92
# Ref HDN+78

Larnite
Ca2SiO4 + 4 H+ = SiO2 + 2 Ca+2 + 2 H2O
log_k 38.4665
-delta_H -227.061 kJ/mol
# deltafH -551.74 kcal/mol
-analytic 2.6900e1 -2.1833e-3 1.0900e4 -9.5257 -7.2537e4
# Range 0-300
-Vm 51.6 # HDN+78 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=172.24, Density (g/ml or g/cm3)= 3.28 therfore, Vm= 172.24 / 3.28 = 52.512
# Extrapol Cp integration
# Ref 82sar/bar

Portlandite
Ca(OH)2 + 2 H+ = Ca+2 + 2 H2O
log_k 22.5552
-delta_H -128.686 kJ/mol
#     deltafH -986.074 kJ/mol
-analytic -8.3848e1 -1.8373e-2 9.3154e3 3.2584e1 1.4538e2
# Range 0-300
-Vm 33.056 # thermo.com.V8.R6+.tdat # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=74.09, Density (g/ml or g/cm3)= 2.23 therfore, Vm= 74.09 / 2.23 = 33.224
# Extrapol Cp integration
# Ref RHF79
# All other solids are defined in the chosen database


SOLUTION 1 Kardemir Process/grid water # Directly taken from KARDEMIR Steel Manufacturer Plant. Define the composition of an aqueous solution. Speciate an analysis, calculate saturation indices.
    temp      25 # temperature in degrees Celcius, default = 25 C, Kardemir 25C
    pH        8.25 # +-0.1,default pH = 7, Kardemir 8.5
    pe        4 # pe = -log(electron activity), default = 4.
    redox     pe # reduction-oxidation reactions
    units     ppm # same as mg/l, default units ppm, mmol/kg water
    density   1 # density, default = 1 kg/L, Kardemir 1
    Ca        160 # +-0.1,total dissolved Ca, Kardemir 140-180
    #Mg        3.4 # +-0.1,total dissolved Mg # Speciate that can be in grid water. No presence of Mg in Kardemir water
    #N(5)      0.84 # 0.84% ppm or mg/l of Nitrogen # Speciate that can be in grid water. No presence of N(5) in Kardemir water
    N(-3)      0.35 # Kardemir NH3 0.2-0.5
    #Na        22 # 11-33 mg/l,total dissolved Na # Speciate that can be in grid water. No presence of Na in Kardemir water
    #K         0.82 # total dissolved K # Speciate that can be in grid water. No presence of K in Kardemir water
    P         0.22 # 0.3% ppm or mg/l of Phosforous, Kardemir P2O5 0.15-0.30
    Fe(2)        0.25 # total dissolved Fe,generally <0.01, Kardemir 0.2-0.3, assumed all to be as ferrous
    Fe(3)        1e-9 # ferric iron assumed as present in negligible quantities
    Alkalinity 240 # 30-400ppm reasonable range alkalinity for domestic drinking water,low alkalinity 0.4- 0.8, Kardemir 210-270 CaCo3
    C(4)          4.620e-03 mol/kgs # equivalent of Alkalinity 240 and pH 8.25
    Cl        85 charge # total dissolved Cl, Kardemir 50-120. Charged balanced.
    S(6)      70 # ppm or mg/l of Sulfate, total dissolved Sulfate, Kardemir 60-80
    #Li        6.941E-03 # Lithium (1umol/L) is a tracer for Process Water
    #Si        2.80855 # mg/L = 0.1 mmol/L
    -water    1 # kg water, default = 1 kg

SAVE SOLUTION 1

END


SOLUTION 2  Median Black Sea rainfall concentrations from Alagha and Tunkel (2003). CO2(atm), O2(atm) added. Chloride charge balance
    temp      25
    pH        5.15
    pe        4
    redox     pe
    units     ppm
    density   1
    Al        3.37
    Ca        1.8
    Cl        1.01 charge # adapt Cl to obtain charge balance
    C(4)      1   CO2(g)   -3.373 # equilibrium with atmospheric CO2
    Fe        2.82
    K         1.72
    Mg        0.88
    Mn        0.07
    N(-3)     0.56
    N(5)      1.5
    Na        2
    S(6)      2.52
    O(0)      0.01 O2(g) -0.678 # equilibrium with atmospheric O2
    -water    1 # kg

END


EQUILIBRIUM_PHASES 1-100

#    Equilibrate SOLUTION 1  water with calcite precipitation

#    Brucite   0 4.733 dissolve_only # 4.733 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (4.733 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 4.733/100=0.04733 cells will give mols /cell.See supporting document
#    Calcite_diss   0 9.49 dissolve_only # Slag calcite dissolves, 9.49 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (9.49 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 9.49/100=0.0949 cells will give mols /cell.See supporting document
     Calcite_ppt   0 0 # Switch on calcite as a phase that can precipitate within days
#    Larnite   0 9.87 dissolve_only # 9.87 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (9.87 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 9.87/100=0.0987 cells will give mols /cell.See supporting document
#    Portlandite 0 14.43 dissolve_only # 14.43 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (14.43 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 14.43/100=0.1443 cells will give mols /cell.See supporting document

END


MIX 1 Combine rainwater with (process/grid water + DAC)
    1    0.88 # 1-(Rainfall_rate/(Rainfall_rate+Procwtr_rate), use fraction 0.88 of solution 1
    2    0.12 # 1-(Procwtr_rate/(Rainfall_rate+Procwtr_rate), use fraction 0.12 of solution 2


USE equilibrium_phases 1

SAVE SOLUTION 0

END


SOLUTION 1-101 Simulated mixture of (process water, DAC) and rainwater.
    temp      25
    pH        6.963
    pe        13.64
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
C(4) 5.06E-03
Ca 3.52E-03
Cl 2.12E-03
Fe 3.94E-06
K 1.07E-06
Mg 1.93E-06
N 2.79E-05
Na 1.07E-05
P 6.26E-06
S(6) 6.45E-04

END


REACTION 0

# Remove / Evaporate water in steps. All solutions defined by SOLUTION input are scaled to have exactly 1 kg (approximately 55.5 mol) of water, unless -water identifier is used.
    H2O        -1
# Evaporate 11.1012 moles(= 0.2 * 55.506), 20% evaporation from pond in 4 step. # 11.1012 mole in 4 step. # The first step removes(-) 2.775e+00 moles of reaction (units are ?moles?) to the initial solution; the second step removes 5.551e+00 moles of reaction to the initial solution; the third 8.326e+00 moles; and the fourth 1.110e+01 moles; each reaction step begins with the same initial solution and removes only the amount of reaction specified.
    11.1012 moles in 4 steps


EQUILIBRIUM_PHASES 101 DAC CO2 capture to 50% #  BOTTOM CELL IS THE EFFLUENT FROM THE COLUMN, AS IT EXPERIENCES CONTACT WITH DAC-AMENDED AIR
    CO2(g)    -0.301 10 # Calculated CO2 SI= -0.301, Amount of 10 mol available
    O2(g)     -0.679 10 # Calculated  O2 SI= -0.679, Amount of 10 mol available
#   N2(g)     -0.107 10
    Calcite_ppt   0 0

SELECTED_OUTPUT 0
    -file                 230828_1125_Antimag_1cm_101cell_DAC_FlowThru_12pw-flush_dt-dx_Cal-KinEqu_CV5_Si_Init-solns20March2025V0_50pc500Days.tsv
    -ionic_strength       true
    -distance             true
    -totals               Ca C(4) Si
    -molalities           Ca+2 CaCO3  CO3-2  HCO3-  CO2
    -saturation_indices   Calcite  Larnite  Brucite  Portlandite CO2(g)
    -equilibrium_phases   Calcite Calcite_ppt CO2(g)  O2(g)
    -kinetic_reactants    Calcite_diss Larnite  Brucite  Portlandite

END


USER_PRINT

10 p = SYS("equi", count , name$ , type$ , moles )
20 FOR i = 1 to count
30   dv = PHASE_VM(name$(i))*EQUI_DELTA(name$(i))
40   delta_v = delta_v + dv
50   PRINT name$(i), dv, " cm^3"
60 NEXT i
70 PRINT "Total change in mineral volume: ", delta_v, " cm^3"
80 END


RATES
    Calcite_diss # calcite dissolving only
-start # Identifier marks the beginning of a Basic program by which the moles of reaction for a time subinterval are calculated.
 5  REM Formula = CaCO3
 10 REM Rate law modified from Zhang et al. 2019
 20 REM Rate law from Marty et al 2015
 30 REM unit should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Calcite") # Calculate calcite saturation ratio
 95 si_mineral = SI("Calcite") # Calculate calcite saturation index
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 320
130 a0 = PARM(1)
132 REM Volume is 1 litre solution
135 v = 1
140 temp = 1/TK
145 REM acid solution parameters
150 k_acid = 0
160 eapp_acid = 0
170 n_acid = 0
175 REM neutral solution parameters
180 k_neut = 6.59E+04
190 eapp_neut = 66000
195 REM base solution parameters
200 k_base = 1.04E+09
210 eapp_base = 67000
220 n_base = 1.6
225 REM Rate Law in TST form
230 hplus = ACT("H+")
240 r = 8.314
250 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
260 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
270 REM Base in terms of bicarbonate
280 r_base = k_base*EXP((-eapp_base/r)*temp)*(ACT("HCO3-")^n_base)
290 r_all = r_acid+r_neut+r_base
300 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2) # Calculate overall rate, factor
310 moles = rate*TIME # Calculate moles of reaction over time interval given by TIME. Note that the multiplication of the rate by TIME must be present in one of the Basic lines.
320 SAVE moles # Return moles of reaction for time subinterval with ?SAVE?. A SAVE statement must always be present in a rate program.
-end # Identifier marks the end of a Basic program by which the number of moles of a reaction for a time subinterval is calculated.

    Brucite
-start
 5  REM Formula = Mg(OH)2
 10 REM Rate parameters from Zhang et al. 2019 database
 20 REM Rate law from Palandri and Kharaka (2004)
 25 REM experimental condition range T=25-75C, pH=1-5
 30 REM unit should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Brucite")
 95 si_mineral = SI("Brucite")
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 380
130 a0 = PARM(1)
140 REM Volume is 1 litre solution
150 v = 1
160 temp = 1/TK
170 REM acid solution parameters
180 k_acid = 4.00E+05
190 eapp_acid = 59000
200 n_acid = 0.500
210 REM neutral solution parameters
220 k_neut = 1.30E-01
230 eapp_neut = 42000
240 REM base solution parameters
250 k_base = 0
260 eapp_base = 0
270 n_base = 0
280 REM Rate Law in TST form
290 hplus = ACT("H+")
300 r = 8.314
310 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
320 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
340 r_base = k_base*EXP((-eapp_base/r)*temp)*(hplus^n_base)
350 r_all = r_acid+r_neut+r_base
360 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2)
370 moles = rate*TIME
380 SAVE moles
-end

    Larnite
-start
 5  REM Formula = Ca2SiO4
 10 REM Rate parameters from Carbfix_2022 database
 20 REM Rate law from Palandri and Kharaka (2004)
 30 REM unit should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Larnite")
 95 si_mineral = SI("Larnite")
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 380
130 a0 = PARM(1)
140 REM Volume is 1 litre solution
150 v = 1
160 temp = 1/TK
170 REM acid solution parameters
180 k_acid = 5.25e8
190 eapp_acid = 70400
200 n_acid = 0.44
210 REM neutral solution parameters
220 k_neut = 0
230 eapp_neut = 0
240 REM base solution parameters
250 k_base = 8.25e5
260 eapp_base = 60900
270 n_base = 0.22
280 REM Rate Law in TST form
290 hplus = ACT("H+")
300 r = 8.314
310 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
320 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
330 IF(k_base = 0) THEN GOTO 350
340 r_base = k_base*EXP((-eapp_base/r)*temp)*(hplus^n_base)
350 r_all = r_acid+r_neut+r_base
360 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2)
370 moles = rate*TIME
375 REM PRINT "SI(Larnite) = " si_mineral ", Molar adjustment = " moles
380 SAVE moles
-end

    Portlandite
-start
 5  REM Formula = Ca(OH)2
 10 REM Rate law modified from Zhang et al. 2019
 20 REM Rate law from marty et al 2015
 25 REM experimental condition range T=25-80C, pH=5-7
 30 REM units should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Portlandite")
 95 si_mineral = SI("Portlandite")
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 320
130 a0 = PARM(1)
132 REM Volume is 1 litre solution
135 v = 1
140 temp = 1/TK
145 REM acid solution parameters
150 k_acid = 1.10E+10
160 eapp_acid = 75000
170 n_acid = 0.600
175 REM neutral solution parameters
180 k_neut = 3.04E+05
190 eapp_neut = 75000
195 REM base solution parameters
200 k_base = 0
210 eapp_base = 0
220 n_base = 0
225 REM Rate Law in TST form
230 hplus = ACT("H+")
240 r = 8.314
250 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
260 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
270 REM Base in terms of bicarbonate
280 r_base = k_base*EXP((-eapp_base/r)*temp)*(ACT("OH-")^n_base)
290 r_all = r_acid+r_neut+r_base
300 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2)
310 moles = rate*TIME
320 SAVE moles
-end



kin_HCO3
-start
10 REM
20 kin_ca=1*KIN_DELTA("Portlandite")+2*KIN_DELTA("Larnite")+0*KIN_DELTA("Brucite")+(0.5)*KIN_DELTA("Calcite_diss")
30 kin_mg=0*KIN_DELTA("Portlandite")+0*KIN_DELTA("Larnite")+1*KIN_DELTA("Brucite")+0*KIN_DELTA("Calcite_diss")
60 moles = -2*(kin_ca+kin_mg)
70 print kin_ca, kin_mg, moles
320 SAVE moles
-end
KINETICS 1-100
Portlandite
    -formula  Ca(OH)2  1
    -m        0.1443 # 14.43/100 cells = 0.1443 cells will give mols /cell. Initial moles of reactant. Mass conc port vs pores/ Molecular-mass_port
    -m0       0.1443 # Default is equal to initial moles
    -parms    2.303 1 # Bulk area * Area% Port
    -tol      1e-08 # Default is 1?10-8 mol
Larnite
    -formula  Ca2SiO4  1
    -m        0.0987 # 9.87/100 cells = 0.0987 cells will give mols /cell
    -m0       0.0987
    -parms    2.499 1
    -tol      1e-08
Calcite_diss
    -formula  CaCO3  1
    -m        0.0949 # 9.49/100 cells = 0.0949 cells will give mols /cell
    -m0       0.0949
    -parms    1.682 1
    -tol      1e-08
Brucite
    -formula  Mg(OH)2  1
    -m        0.04733 # 4.733/100 cells = 0.04733 cells will give mols /cell
    -m0       0.04733
    -parms    0.556 1
    -tol      1e-08

-steps       1
-step_divide 0.1
-runge_kutta 6
-bad_step_max 500
-cvode true
-cvode_steps 200
-cvode_order 5

END

USE SOLUTION 0
USE REACTION 0
USE equilibrium_phases 101
USE KINETICS 1
SAVE equilibrium_phases 1-10
SAVE solution 1-10



TRANSPORT

    -cells                 101  # 100 cells slag heap, last cell closed pond (1 litre)
  # -shifts                303 # 3 passes of pore water through the slag heap
    -shifts                1245 # 12.326 passes of pore water through the slag heap
  #  -shifts               101 # 1 passes of pore water through the slag heap
    -lengths               0.01 # 1cm cells
    -time_step             34722.2 # seconds retention time in cells - Antimagnetic = 0.401 days. Shift (101) x Timestep(34722.2) = 3506942.2 sec = 3506942.2 /60/ 60/ 24 = The total time of the simulation (40.5896 days).Shift (303) x Timestep(34722.2) = 10520826.6 sec = 10520826.6/60/ 60/ 24 = The total time of the simulation (121.768 days).
  #  Shift (1245) x Timestep(34722.2) = 43200000sec = 43200000 /60/ 60/ 24 = The total time of the simulation (500 days)
    -print_cells           1, 100-101
  # -print_cells           10-11
    -print_frequency       1245 # Print only at end
  # -print_frequency       33 # Print only at end



USER_PUNCH 1 # This keyword data block is used to define Basic programs that print user-defined quantities to the selected-output file. Any Basic ?PUNCH? statement will write to the selected-output file.
    -headings kin_hco3 Pore_vol # Cum-C-heading
    -start
 7 REM IF(CELL_NO < 101) THEN GOTO 100
10 IF (TOTAL_TIME = 0) THEN PUT(0,1)
20 kin_ca=1*KIN_DELTA("Portlandite")+2*KIN_DELTA("Larnite")+0*KIN_DELTA("Brucite")+(0.5)*KIN_DELTA("Calcite_diss")
30 kin_mg=0*KIN_DELTA("Portlandite")+0*KIN_DELTA("Larnite")+1*KIN_DELTA("Brucite")+0*KIN_DELTA("Calcite_diss")
60 kin_hco3 = -2*(kin_ca+kin_mg)
65 PUT(GET(1)+kin_hco3,1)
67 cum_c = GET(1)
68 PUNCH kin_hco3
#70 PUNCH cum_c
75 REM punch pore volumes, no_cells = 100
77 no_cells = 100
80 PUNCH (STEP_NO + .5) / no_cells
100 REM END
    -end



SELECTED_OUTPUT 1
    -file                 230828_1125_Antimag_1cm_101cell_DAC_FlowThru_12pw-flush_dt-dx_Cal-KinEqu_CV5_Si_adv-react20March2025V0_50pc500Days.tsv
    -ionic_strength       true
    -distance             true
    -totals               Li Ca C(4) Si   # Lithium is a tracer for Process Water
    -molalities           Ca+2 CaCO3  CO3-2  HCO3-  CO2
    -saturation_indices   Calcite  Larnite  Brucite  Portlandite CO2(g)
    -equilibrium_phases   Calcite Calcite_ppt CO2(g)  O2(g)
    -kinetic_reactants    Calcite_diss Larnite  Brucite  Portlandite kin_HCO3


END



[dlparkhurst]

I don't understand what you are simulating. I have simplified your script and only simulated 20 cells.

You have evaporated 20% of the water in the first 10 cells. However, that evaporated water advects through the column and will be gone after 20 shifts (for my simulation), after which, there is no effect of evaporated water. Calcite is always in equilibrium.

You will have explain what you mean by "erratic readings". If possible, make a batch reaction where a change in CO2 concentration demonstrates the erratic readings. Otherwise, please use a column with fewer cells that runs quickly and explain what results are unexpected.

The following shows the "bubble" of evaporated water moving down the column after 5 shifts and the equilibrium with calcite.

Code: [Select]

TITLE DATABASE C:\Users\ErdincCosgun\Documents\Projects\Databases\carbfix_2022.dat
TITLE TESTING KINETICS OF SPECIFIC SOLIDS
TITLE 100 CELLS, 1CM HIGH
TITLE FLOW-THROUGH TEST, 3 LITRE
TITLE === FULL ANTIMAGNETIC SUITE ===
TITLE KARDEMIR PROCESS WATER MEAN VALUES
TITLE REDOX CONTROLLED BY IRON SPECIES
TITLE CHEMICAL SPECIATION
TITLE PHASES FROM Carbfix2022
TITLE CALCITE DISSOLVES AS KINETIC PHASE
TITLE CALCITE PRECIPITATION AS EQUILIBRIUM PHASE
TITLE CO2 INJECTED FROM DAC 50% TO SOLUTION, log(0.5):-0.301
TITLE Li ADDED IN TRACE TO KARDEMIR SUPPLIED WATER
RATES
    Calcite_diss # calcite dissolving only
-start # Identifier marks the beginning of a Basic program by which the moles of reaction for a time subinterval are calculated.
 5  REM Formula = CaCO3
 10 REM Rate law modified from Zhang et al. 2019
 20 REM Rate law from Marty et al 2015
 30 REM unit should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Calcite") # Calculate calcite saturation ratio
 95 si_mineral = SI("Calcite") # Calculate calcite saturation index
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 320
130 a0 = PARM(1)
132 REM Volume is 1 litre solution
135 v = 1
140 temp = 1/TK
145 REM acid solution parameters
150 k_acid = 0
160 eapp_acid = 0
170 n_acid = 0
175 REM neutral solution parameters
180 k_neut = 6.59E+04
190 eapp_neut = 66000
195 REM base solution parameters
200 k_base = 1.04E+09
210 eapp_base = 67000
220 n_base = 1.6
225 REM Rate Law in TST form
230 hplus = ACT("H+")
240 r = 8.314
250 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
260 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
270 REM Base in terms of bicarbonate
280 r_base = k_base*EXP((-eapp_base/r)*temp)*(ACT("HCO3-")^n_base)
290 r_all = r_acid+r_neut+r_base
300 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2) # Calculate overall rate, factor
310 moles = rate*TIME # Calculate moles of reaction over time interval given by TIME. Note that the multiplication of the rate by TIME must be present in one of the Basic lines.
320 SAVE moles # Return moles of reaction for time subinterval with ?SAVE?. A SAVE statement must always be present in a rate program.
-end # Identifier marks the end of a Basic program by which the number of moles of a reaction for a time subinterval is calculated.

    Brucite
-start
 5  REM Formula = Mg(OH)2
 10 REM Rate parameters from Zhang et al. 2019 database
 20 REM Rate law from Palandri and Kharaka (2004)
 25 REM experimental condition range T=25-75C, pH=1-5
 30 REM unit should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Brucite")
 95 si_mineral = SI("Brucite")
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 380
130 a0 = PARM(1)
140 REM Volume is 1 litre solution
150 v = 1
160 temp = 1/TK
170 REM acid solution parameters
180 k_acid = 4.00E+05
190 eapp_acid = 59000
200 n_acid = 0.500
210 REM neutral solution parameters
220 k_neut = 1.30E-01
230 eapp_neut = 42000
240 REM base solution parameters
250 k_base = 0
260 eapp_base = 0
270 n_base = 0
280 REM Rate Law in TST form
290 hplus = ACT("H+")
300 r = 8.314
310 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
320 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
340 r_base = k_base*EXP((-eapp_base/r)*temp)*(hplus^n_base)
350 r_all = r_acid+r_neut+r_base
360 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2)
370 moles = rate*TIME
380 SAVE moles
-end

    Larnite
-start
 5  REM Formula = Ca2SiO4
 10 REM Rate parameters from Carbfix_2022 database
 20 REM Rate law from Palandri and Kharaka (2004)
 30 REM unit should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Larnite")
 95 si_mineral = SI("Larnite")
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 380
130 a0 = PARM(1)
140 REM Volume is 1 litre solution
150 v = 1
160 temp = 1/TK
170 REM acid solution parameters
180 k_acid = 5.25e8
190 eapp_acid = 70400
200 n_acid = 0.44
210 REM neutral solution parameters
220 k_neut = 0
230 eapp_neut = 0
240 REM base solution parameters
250 k_base = 8.25e5
260 eapp_base = 60900
270 n_base = 0.22
280 REM Rate Law in TST form
290 hplus = ACT("H+")
300 r = 8.314
310 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
320 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
330 IF(k_base = 0) THEN GOTO 350
340 r_base = k_base*EXP((-eapp_base/r)*temp)*(hplus^n_base)
350 r_all = r_acid+r_neut+r_base
360 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2)
370 moles = rate*TIME
375 REM PRINT "SI(Larnite) = " si_mineral ", Molar adjustment = " moles
380 SAVE moles
-end

    Portlandite
-start
 5  REM Formula = Ca(OH)2
 10 REM Rate law modified from Zhang et al. 2019
 20 REM Rate law from marty et al 2015
 25 REM experimental condition range T=25-80C, pH=5-7
 30 REM units should be J, m2, mol,kgw-1 and second-1
 40 REM parm(1) is surface area in the unit of m2/kgw
 50 REM calculation of surface area can be found in the note
 60 REM M is current moles of minerals. M0 is the initial moles of minerals
 70 REM parm(2) is a correction factor
 80 REM Affinity ON
 90 sr_mineral = SR("Portlandite")
 95 si_mineral = SI("Portlandite")
100 moles = 0
110 REM rate=0 if mineral supersaturated in solution
120 IF(M <= 0 OR sr_mineral > 1) THEN GOTO 320
130 a0 = PARM(1)
132 REM Volume is 1 litre solution
135 v = 1
140 temp = 1/TK
145 REM acid solution parameters
150 k_acid = 1.10E+10
160 eapp_acid = 75000
170 n_acid = 0.600
175 REM neutral solution parameters
180 k_neut = 3.04E+05
190 eapp_neut = 75000
195 REM base solution parameters
200 k_base = 0
210 eapp_base = 0
220 n_base = 0
225 REM Rate Law in TST form
230 hplus = ACT("H+")
240 r = 8.314
250 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)
260 r_neut = k_neut*EXP((-eapp_neut/r)*temp)
270 REM Base in terms of bicarbonate
280 r_base = k_base*EXP((-eapp_base/r)*temp)*(ACT("OH-")^n_base)
290 r_all = r_acid+r_neut+r_base
300 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2)
310 moles = rate*TIME
320 SAVE moles
-end



kin_HCO3
-start
10 REM
20 kin_ca=1*KIN_DELTA("Portlandite")+2*KIN_DELTA("Larnite")+0*KIN_DELTA("Brucite")+(0.5)*KIN_DELTA("Calcite_diss")
30 kin_mg=0*KIN_DELTA("Portlandite")+0*KIN_DELTA("Larnite")+1*KIN_DELTA("Brucite")+0*KIN_DELTA("Calcite_diss")
60 moles = -2*(kin_ca+kin_mg)
70 print kin_ca, kin_mg, moles
320 SAVE moles
-end
PHASES # Minerals, Directly taken from Carbfix2022.dat PHASES.

Calcite_diss # Calcite that dissolves only. Alphanumeric name of phase; no spaces are allowed.
CaCO3 + H+ = Ca+2 + HCO3- # Dissolution reaction
log_k 1.8487 # Log K at 25 ?C for the reaction
-delta_H -25.7149 kJ/mol # Enthalpy of reaction at 25 ?C for the reaction
# deltafH -288.552 kcal/mol
-analytic -1.4978e2 -4.8370e-2 4.8974e3 6.0458e1 7.6464e1 # Identifier for coefficients for an analytical expression for the temperature dependence of log K.
# Range 0-350
-Vm 36.934 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=100.09, Density (g/ml or g/cm3)= 2.71 therfore, Vm= 100.09 / 2.71 = 36.934
# Extrapol supcrt92
# Ref HDN+78

Calcite_ppt # Calcite that percipitates only
CaCO3 + H+ = Ca+2 + HCO3-
log_k 1.8487
-delta_H -25.7149 kJ/mol
# deltafH -288.552 kcal/mol
-analytic -1.4978e2 -4.8370e-2 4.8974e3 6.0458e1 7.6464e1
# Range 0-350
-Vm 36.934 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=100.09, Density (g/ml or g/cm3)= 2.71 therfore, Vm= 100.09 / 2.71 = 36.934
# Extrapol supcrt92
# Ref HDN+78

Brucite
Mg(OH)2 + 2 H+ = Mg+2 + 2 H2O
log_k 16.2980
-delta_H -111.34 kJ/mol
# deltafH -221.39 kcal/mol
-analytic -1.0280e2 -1.9759e-2 9.0180e3 3.8282e1 1.4075e2
# Range 0-350
-Vm 24.63 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=58.32, Density (g/ml or g/cm3)= 2.39 therfore, Vm= 58.32 / 2.39 = 24.40
# Extrapol supcrt92
# Ref HDN+78

Larnite
Ca2SiO4 + 4 H+ = SiO2 + 2 Ca+2 + 2 H2O
log_k 38.4665
-delta_H -227.061 kJ/mol
# deltafH -551.74 kcal/mol
-analytic 2.6900e1 -2.1833e-3 1.0900e4 -9.5257 -7.2537e4
# Range 0-300
-Vm 51.6 # HDN+78 # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=172.24, Density (g/ml or g/cm3)= 3.28 therfore, Vm= 172.24 / 3.28 = 52.512
# Extrapol Cp integration
# Ref 82sar/bar

Portlandite
Ca(OH)2 + 2 H+ = Ca+2 + 2 H2O
log_k 22.5552
-delta_H -128.686 kJ/mol
#     deltafH -986.074 kJ/mol
-analytic -8.3848e1 -1.8373e-2 9.3154e3 3.2584e1 1.4538e2
# Range 0-300
-Vm 33.056 # thermo.com.V8.R6+.tdat # molar_volume, the molecular weight divided by the density of the solid at 25 ?C. Ex: Molecular Weight/Mass (g/mole)=74.09, Density (g/ml or g/cm3)= 2.23 therfore, Vm= 74.09 / 2.23 = 33.224
# Extrapol Cp integration
# Ref RHF79
# All other solids are defined in the chosen database


SOLUTION 1 Kardemir Process/grid water # Directly taken from KARDEMIR Steel Manufacturer Plant. Define the composition of an aqueous solution. Speciate an analysis, calculate saturation indices.
    temp      25 # temperature in degrees Celcius, default = 25 C, Kardemir 25C
    pH        8.25 # +-0.1,default pH = 7, Kardemir 8.5
    pe        4 # pe = -log(electron activity), default = 4.
    redox     pe # reduction-oxidation reactions
    units     ppm # same as mg/l, default units ppm, mmol/kg water
    density   1 # density, default = 1 kg/L, Kardemir 1
    Ca        160 # +-0.1,total dissolved Ca, Kardemir 140-180
    #Mg        3.4 # +-0.1,total dissolved Mg # Speciate that can be in grid water. No presence of Mg in Kardemir water
    #N(5)      0.84 # 0.84% ppm or mg/l of Nitrogen # Speciate that can be in grid water. No presence of N(5) in Kardemir water
    N(-3)      0.35 # Kardemir NH3 0.2-0.5
    #Na        22 # 11-33 mg/l,total dissolved Na # Speciate that can be in grid water. No presence of Na in Kardemir water
    #K         0.82 # total dissolved K # Speciate that can be in grid water. No presence of K in Kardemir water
    P         0.22 # 0.3% ppm or mg/l of Phosforous, Kardemir P2O5 0.15-0.30
    Fe(2)        0.25 # total dissolved Fe,generally <0.01, Kardemir 0.2-0.3, assumed all to be as ferrous
    Fe(3)        1e-9 # ferric iron assumed as present in negligible quantities
    Alkalinity 240 # 30-400ppm reasonable range alkalinity for domestic drinking water,low alkalinity 0.4- 0.8, Kardemir 210-270 CaCo3
    C(4)          4.620e-03 mol/kgs # equivalent of Alkalinity 240 and pH 8.25
    Cl        85 charge # total dissolved Cl, Kardemir 50-120. Charged balanced.
    S(6)      70 # ppm or mg/l of Sulfate, total dissolved Sulfate, Kardemir 60-80
    #Li        6.941E-03 # Lithium (1umol/L) is a tracer for Process Water
    #Si        2.80855 # mg/L = 0.1 mmol/L
    -water    1 # kg water, default = 1 kg

END


SOLUTION 2  Median Black Sea rainfall concentrations from Alagha and Tunkel (2003). CO2(atm), O2(atm) added. Chloride charge balance
    temp      25
    pH        5.15
    pe        4
    redox     pe
    units     ppm
    density   1
    Al        3.37
    Ca        1.8
    Cl        1.01 charge # adapt Cl to obtain charge balance
    C(4)      1   CO2(g)   -3.373 # equilibrium with atmospheric CO2
    Fe        2.82
    K         1.72
    Mg        0.88
    Mn        0.07
    N(-3)     0.56
    N(5)      1.5
    Na        2
    S(6)      2.52
    O(0)      0.01 O2(g) -0.678 # equilibrium with atmospheric O2
    -water    1 # kg

END


EQUILIBRIUM_PHASES 1-100

#    Equilibrate SOLUTION 1  water with calcite precipitation

#    Brucite   0 4.733 dissolve_only # 4.733 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (4.733 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 4.733/100=0.04733 cells will give mols /cell.See supporting document
#    Calcite_diss   0 9.49 dissolve_only # Slag calcite dissolves, 9.49 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (9.49 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 9.49/100=0.0949 cells will give mols /cell.See supporting document
     Calcite_ppt   0 0 # Switch on calcite as a phase that can precipitate within days
#    Larnite   0 9.87 dissolve_only # 9.87 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (9.87 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 9.87/100=0.0987 cells will give mols /cell.See supporting document
#    Portlandite 0 14.43 dissolve_only # 14.43 (mol_mineral/L_pores)calculated. PHREEQC output assigns this value (14.43 )to each cell rather than the total depth of the core therfore dividing the moles of mineral by the number of cells is required, 14.43/100=0.1443 cells will give mols /cell.See supporting document

END


MIX 1 Combine rainwater with (process/grid water + DAC)
    1    0.88 # 1-(Rainfall_rate/(Rainfall_rate+Procwtr_rate), use fraction 0.88 of solution 1
    2    0.12 # 1-(Procwtr_rate/(Rainfall_rate+Procwtr_rate), use fraction 0.12 of solution 2


USE equilibrium_phases 1

SAVE SOLUTION 0

END


SOLUTION 1-101 Simulated mixture of (process water, DAC) and rainwater.
    temp      25
    pH        6.963
    pe        13.64
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
C(4) 5.06E-03
Ca 3.52E-03
Cl 2.12E-03
Fe 3.94E-06
K 1.07E-06
Mg 1.93E-06
N 2.79E-05
Na 1.07E-05
P 6.26E-06
S(6) 6.45E-04

END


REACTION 0

# Remove / Evaporate water in steps. All solutions defined by SOLUTION input are scaled to have exactly 1 kg (approximately 55.5 mol) of water, unless -water identifier is used.
    H2O        -1
# Evaporate 11.1012 moles(= 0.2 * 55.506), 20% evaporation from pond in 4 step. # 11.1012 mole in 4 step. # The first step removes(-) 2.775e+00 moles of reaction (units are ?moles?) to the initial solution; the second step removes 5.551e+00 moles of reaction to the initial solution; the third 8.326e+00 moles; and the fourth 1.110e+01 moles; each reaction step begins with the same initial solution and removes only the amount of reaction specified.
    11.1012 moles #in 4 steps


EQUILIBRIUM_PHASES 101 DAC CO2 capture to 50% #  BOTTOM CELL IS THE EFFLUENT FROM THE COLUMN, AS IT EXPERIENCES CONTACT WITH DAC-AMENDED AIR
    CO2(g)    -0.301 10 # Calculated CO2 SI= -0.301, Amount of 10 mol available
    O2(g)     -0.679 10 # Calculated  O2 SI= -0.679, Amount of 10 mol available
#   N2(g)     -0.107 10
    Calcite_ppt   0 0

END

KINETICS 1-100
Portlandite
    -formula  Ca(OH)2  1
    -m        0.1443 # 14.43/100 cells = 0.1443 cells will give mols /cell. Initial moles of reactant. Mass conc port vs pores/ Molecular-mass_port
    -m0       0.1443 # Default is equal to initial moles
    -parms    2.303 1 # Bulk area * Area% Port
    -tol      1e-08 # Default is 1?10-8 mol
Larnite
    -formula  Ca2SiO4  1
    -m        0.0987 # 9.87/100 cells = 0.0987 cells will give mols /cell
    -m0       0.0987
    -parms    2.499 1
    -tol      1e-08
Calcite_diss
    -formula  CaCO3  1
    -m        0.0949 # 9.49/100 cells = 0.0949 cells will give mols /cell
    -m0       0.0949
    -parms    1.682 1
    -tol      1e-08
Brucite
    -formula  Mg(OH)2  1
    -m        0.04733 # 4.733/100 cells = 0.04733 cells will give mols /cell
    -m0       0.04733
    -parms    0.556 1
    -tol      1e-08

-steps       1
-step_divide 0.1
-runge_kutta 6
-bad_step_max 500
-cvode true
-cvode_steps 200
-cvode_order 5

END

USE SOLUTION 0
USE REACTION 0
USE equilibrium_phases 101
USE KINETICS 1
SAVE equilibrium_phases 1-10
SAVE solution 1-10

END

TRANSPORT
    -cells                 20 # 101  # 100 cells slag heap, last cell closed pond (1 litre)
    -shifts                5 # 1245 # 12.326 passes of pore water through the slag heap
    -lengths               0.01 # 1cm cells
    -time_step             34722.2 # seconds retention time in cells - Antimagnetic = 0.401 days. Shift (101) x Timestep(34722.2) = 3506942.2 sec = 3506942.2 /60/ 60/ 24 = The total time of the simulation (40.5896 days).Shift (303) x Timestep(34722.2) = 10520826.6 sec = 10520826.6/60/ 60/ 24 = The total time of the simulation (121.768 days).
    -print_cells           0-20
    -punch_cells           1-20
    -print_frequency       5 # 1245 # Print only at end
    -punch_frequency       5

USER_GRAPH 1
    -headings               dist Portlandite Larnite Calcite Brucite H2O_kg
    -axis_titles            "Distance" "SI" "Water, kg"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start

 10 graph_x DIST
 20 graph_y SI("Portlandite"), SI("Larnite"), SI("Calcite_diss"), SI("Brucite")
 30 graph_sy TOT("water")
100 END
  -end
END

[Erdinc]

Thanks David. I have calculated the shift by number of days I want to simulate (500 days in my case) to understand the dissolution of the mineral kinetics. I have run the simulations from different percentages of pCO2 values and from 0.1% up to 6% there is good indication on how much of the minerals dissolve during the life of the 500 days. When I switch the pCO2 value 7% and above up to 99%, I see the simulations stops completing all the way to 500 days and therefore same for minerals dissolutions and CO2 capture potentials are not known. I was thinking the effect of accurate thermodynamic solubility data, the non-ideal behavior of CO2 fugacity, and potential limitations in the way PHREEQC handles gas-water interactions would be a reason for the high pCO2 values. Is this something you come across before? Unfortunately, I do not have options to attach what I am simulating in regard to graphs but this is what I can explain to the best of my knowledge.

Many Thanks for the support

[dlparkhurst]

Okay, I ran your (very) long simulation with llnl.dat and a version of carbfix.dat. Both ran to completion.

What error am I supposed to see?

Will a 20 cell model with no mixruns demonstrate the problem? The current run time is too long to debug usefully.

I'll send you a message and you can send me a script or graph that shows me what you are talking about. I don't see anything wrong.