CO2 Capture by Steel Slag Heap

[Erdinc]

Dear PHREEQC users,

I am currently working on 1D Reactive Transport Model. Intention is to weather slag heap and find potential carbon capture through weathering. I have calculated Total Cumulative Carbon capture by using HCO3 determined from the dissolution of slag minerals which are subtracted from changes in CaCO3 precipitation. Looks like I have Cumulative Carbon removal of 29.62 kg/tonneslag or 0.02962 Ton/Ton. Since this is simulation, I wanted to verify the answer with amount of mineral dissolved (e.g. k_dissolved mineral). From output I can see some changes on some minerals (e.g. calcite_dissolved) but cannot verify clearly as the end of simulation value came back to initial mols available. Can anyone please help me to understand or guide me if I am doing the right approach on comparing? My script is below.


[code]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/GRID 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 TO SOLUTION
TITLE LI ADDED IN TRACE TO KARDEMIR SUPPLIED GRID WATER


# Information following a # is not read

########################## SWITCH OFF OUTPUT TO file.pqo ##########################
###################################################################################
# Need to increase speed of program and to cut down on the megabytes of output created

# PRINT
# -reset FALSE # suspend all print output

# This can always be changed if specific runs need to be investigated.Using PRINT, I can modify the print to the output tab by setting output block identifiers true or false. These identifiers can be clicked and copied to the input tab after expanding keyword PRINT in the right window. The option -reset false is rather essential when doing transport simulations.This keyword data block is used to select which results are written to the output file. In addition, this data block enables or disables writing results to the selected-output file and writing a status line to the screen, which monitors the type of calculation being performed.
# If used, this identifier should be the first identifier of the data block. If true, all data blocks described on Lines are printed to the output file; if false, these data blocks are excluded from the output file. If neither true nor false is entered on the line, true is assumed.
# For long TRANSPORT and ADVECTION calculations with KINETICS, printing the status line [-status true (default)] may cause a significant increase in run time. In this case, setting -status false will speed up run times. Below status set to false.

######################################################################################################################### DEFINE SOLID PHASES IN SLAG SYSTEM #################################################################################################################################################
########################################################################################################################################################################################################################################################################################################

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

# PHASES keyword data block is used to define a name, chemical reaction, log K, and temperature dependence of log K for each gas component and mineral that can be used for speciation, batch-reaction and transport calculations. In addition, molar volumes can be defined for solids, and the critical temperature and pressure and the acentric factor can be defined for gases. Normally, this data block is included in the database file and only additions and modifications are included in the input file.
# There are 4 mineral phases on Antimagnetic Slag that I need to write below so that model can identify the phases. These are Brucite, Calcite_diss (Slag calcite dissolves only), Calcite_ppt (Slag calcite precipitate only), Larnite and Portlandite.

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- # 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

Brucite
   Mg(OH)2 + 2 H+ = Mg+2 + 2 H2O # Dissolution reaction
   log_k      16.2980 # Log K at 25 ?C for the reaction
   -delta_H   -111.34      kJ/mol # Enthalpy of reaction at 25 ?C for the reaction
#   deltafH      -221.39      kcal/mol
   -analytic   -1.0280e2 -1.9759e-2 9.0180e3 3.8282e1 1.4075e2 # Identifier for coefficients for an analytical expression for the temperature dependence of log K.
#   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 # Dissolution reaction
   log_k      38.4665 # Log K at 25 ?C for the reaction
   -delta_H   -227.061   kJ/mol # Enthalpy of reaction at 25 ?C for the reaction
#   deltafH      -551.74      kcal/mol
   -analytic   2.6900e1 -2.1833e-3 1.0900e4 -9.5257 -7.2537e4 # Identifier for coefficients for an analytical expression for the temperature dependence of log K.
#   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 # Dissolution reaction
   log_k      22.5552 # Log K at 25 ?C for the reaction
   -delta_H   -128.686   kJ/mol # Enthalpy of reaction at 25 ?C for the reaction
#     deltafH      -986.074   kJ/mol
   -analytic   -8.3848e1 -1.8373e-2 9.3154e3 3.2584e1 1.4538e2 # Identifier for coefficients for an analytical expression for the temperature dependence of log K.
#   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



########################################################################################################### GEOCHAMICAL NODE 1 - MAKE UP WATER/GRID WATER OR WASTEWATER ###############################################################################################################################
########################################################################################################################################################################################################################################################################################################

# Process water/grid water / makeup water is pumped from the Kardemir steel plant. This conceptual model is the Kardemir process water in the hypothetical scenario of no contact with atmosphere, still in the pumping/piping system. Mean Values are used.
# Process SOLUTION Use of Kardemir water25062022 - MEAN CONCENTRATIONS (AVG) + Tracer = Li, Chloride charge balance.


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

# The solution is charge balanced. In Ionic solution, the sum of negative electrical charges equals to sum of positive electrical charges (C+ = A-).(-4.15 meaning more negative ions) and balance adjustment may be required for future calculations, hence charge balance. Percentage error simply describes Cations and Anions balance. It should be noted that Percent error, 100*(Cat-|An|)/(Cat+|An|) is unitless and it is important to adjust the system to become electrically neutral/ balanced because PHREEQC will look for demand of equilibrium and will force charge balance with one of the components (that can be designated), as they solve the matrix of non-linear equations.
# Only analyses of water samples with an ion-balance error < 5% are discussed to be tolerable (Appelo  and Postma, 2005) and , (Wisotzky et al., 2018). If errors exceed 5% either not all ionic compounds were reflected in the analysis, or individual compound concentrations may be over- or underestimated or individual analysis results are incorrect (Appelo  and Postma, 2005).
# When the simulation run without charge balance, Electrical balance (eq)  =  -6.707e-16 is pretty close to 0 so it seems ok to use, no further adjustment needed
# Values given in solution now converted in molality/moles. For example given Alkalinity (In PHREEQC the alkalinity is always stated as CaCO 3 equivalents) 240ppm is also equilvalent to 4.801e-03 moles.
# All carbon in the form of carbonic acid species are given as C (4) which specifies element valence with brackets C (4)= 4.624e-03.Carbon speciation contrubution is HCO3- (4.246E-03/4.624e-03= 92%), CaCO3 (1.628E-04/4.624e-03= 3.5%), CaHCO3+ (1.18E-04/4.624e-03= 2.55%), CO2   (4.99E-05/4.624e-03= 1%) and CO3-2   (4.66E-05/4.624e-03= 1%).
# The Ca concentrations is 3.996e-03 molar. It is listed in the order of Magnitude, most abandon is Ca+2 with 3.581e-03 molar.
# Phase saturations, Calcite   1.37   : Supersaturated(expected to form), Fe(OH)3  -7.39 : Undersaturated(expected to dissolve) , CO2(g)  -2.84 : Undersaturated(expected to dissolve), Portlandite  -8.68 : Undersaturated(expected to dissolve) and O2(g)  -34.12  : Undersaturated(expected to dissolve).     
# Iron is present in solution, such that ferrihydrite (Fe(OH)3) is understaurated. This is probably due to the pe of the system being arbitrarily set at 4.0.Iron redox behaviour is important because iron forms solid precipitates such as Fe(OH)3, often called ferrihydrite or amorphous ferric hydroxide, when a reduced solution is oxidized (usually by contact with the atmosphere). Ferrihydrite is a good adsorption surface and will remove many metals from solution. Ferrihydrite is gradually converted by de-hydration to limonite, goethite and finally, hematite.Ferrihydrite readily dissolves in the Kardemir process water system. It is likely that the redox status of the system is such that ferrous ion is stabilised in solution. Ferrous ion is much less likely to precipitate as an oxide.
# It is also possible that calcite is remaining supersaturated during the lifetime of the pumping of process water. If this is true, there will be no scaling potential.Calcite is allowed to precipitate, though none of the other phases are allowed to do so due to kinetic reaction.



SAVE SOLUTION 1

# SAVE keyword data block is used to save the composition of a solution, exchange assemblage, gas phase, equilibrium-phase assemblage, solid-solution assemblage, or surface assemblage following a batch-reaction calculation. The composition is stored internally in computer memory and can be retrieved subsequently with the USE keyword during the remainder of the computer run. SAVE affects only the internal storage of chemical-composition information during the current run; it does not save information between PHREEQC runs. The SAVE data block applies only at the end of batch-reaction calculations and has no effect following initial solution or transport calculations.

END

# Ends the program. The input file is then read and processed simulation by simulation (as defined by END keywords) until the end of the file.


########################################################################################################### GEOCHAMICAL NODE 2 - RAINFALL AND HEAP LEACH IRRIGATION ####################################################################################################################################
########################################################################################################################################################################################################################################################################################################

# The chemistry of the rainfall entering on the Kardemir slag heap is taken as the median concentrations from (Alagha and Tunkel, 2003) (DOI: 10.1023/A:1026049027446).
# The alkalinity values for these chemical analyses are not reported. Therefore, C(4) is introduced at atmospheric concentration. Chloride is here chosen as the charge-balancing anion.
# The pH drop of Kardemir make up water due to CO2 introduction promotes the dissolution of alkali earth metals such as calcium and magnesium. These metals react with the dissolved carbon bearing species such as bicarbonate or carbonate which result in the precipitation of carbonate minerals such as calcite. As an example, the reaction between carbonate and calcium leads to the precipitation of calcite. Interaction among CO2- water-slag heap produces a complex series of reactions that varies depending on the chemical composition of the fluid and the mineralogy of the slag. (Rosenbauer et al., 2005) defined similar complex reactions between CO2-brine-rock. These complex reactions may promote the change in porosity of Slag heap. Excess dissolution of slag heap mineral will possibly lead to porosity increase which is desirable since it increases the CO2 storage capacity of a slag. When precipitation of new or existing mineral is dominant, porosity decrease will be observed.

# USE SOLUTION 1

# USE keyword data block is used to specify explicitly which solution is to be used in the batch-reaction calculation of a simulation. USE also can specify kinetically controlled reactions (KINETICS data block) and mixing parameters (MIX data block) to be used in a batch-reaction calculation.

# SOLUTION 1 is the process/grid water that is used to irrigate the slag heap

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

# Fugacity of CO2(g) in the atmosphere is approximately 0.000424, or 0.04%. Saturation Index of the CO2(g) in the atmosphere is therefore Saturation Index CO2 (-pCO2)= SI(CO2(g) = log10(fugacity) = log10(0.000424) = -3.373. The gas pressure can be varied, setting SI. The SI is equal to log(gas pressure / 1 atm).This is the number used in the PHREEQC calculations.
# The saturation state of O2 is more complex. Oxygen readily dissolves in aqueous solutions exposed to atmosphere. Equilibrate water with O2(g) 209460ppm = 0.20946= SI(CO2(g) = log10(fugacity=0.20946) =-0.678. (https://en.wikipedia.org/wiki/Atmosphere_of_Earth). Thus, the anticipated SI of oxygen is log10(0.209) = -0.678. This is the number used in the PHREEQC calculations.


EQUILIBRIUM_PHASES 1-100

#    Equilibrate SOLUTION 1  water with calcite precipitation

#    Brucite   0 4.733 dissolve_only # 4.733 (mol_mineral/L_pores)calculated, see supporting document
#    Calcite_diss   0 9.49 dissolve_only # Slag calcite dissolves, 9.49 (mol_mineral/L_pores)calculated, 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, see supporting document
#    Portlandite 0 14.43 dissolve_only # 14.43 (mol_mineral/L_pores)calculated, see supporting document

END

############################################################################################### CALCULATE FLOW RATES THROUGH CELLS, AS RETENTION TIMES AND DEPTHS FOR SSANTIMAGNETIC #####################################################################################################################
##########################################################################################################################################################################################################################################################################################################


# Depth or height of cell (-lengths)   0.01   metres = 1cm
# Number of horizontal slices   100
# No ponding on the surface of the slag heap is assumed, therefore the advective ingress of rainfall may be estimated as follows. Total annual millimetres of rainfall can be extracted to be 522.1 (mm/a)(https://en.wikipedia.org/wiki/Karab%C3%BCk). A broad average of rainfall, in terms of metres per second is thus calculated:
# Total annual millimetres of rainfall 522.1   mm/a
# Convert to m/a= 522.1/1000= 0.522   m/a
# Year = 365.25   d/a
# Convert to m/d= 0.522/ 365.25 = 1.429E-03   m/d
# Convert to s/d = 60x60x24 = 86400   s/d
# Convert to m/s = 1.429E-03 / 86400 = 1.654E-08   m/s. Annual average precipitation is thus 1.654E-08 m/s.
# Slag porosity for Antimagnetic is calculated at: 46% pore volume/slag bulk volume and 56% pore volume on Magnetic.see supporting document

# Advective velocity of rainfall through the Antimagnetic slag heap is thus 1.654E-08/0.46 m/s. ########## Advective velocity of rainfall through the slag heap is 3.596E-08 m/s. and 2.953E-08 (Magnetic) ##########
# The slag heap is 1 metre high and is constantly irrigated with Process-water/grid water and rainwater. The irrigation rate is 10 millimetres per day. The velocity of the irrigation water before it enters the heap is 10   mm_irrigation/day.
# Convert to milimetres per second   10mm/86400sec = 1.160E-04 mm_irrigation/second.
# Convert to metres per second irrigation   1.160E-04/ 1000 = 1.160E-07   m_irrgation/second.
# Porosity of Anti Magnetic slag 0.46 L_pores/L_slag. ########## Thus the velocity of pore waters is   1.160E-07/0.46 = 2.520E-07 m/s. and 2.07E-07 (Magnetic)##########
# Combined velocity of pore waters = 3.596E-08 m/s + 2.520E-07 m/s = ########## 2.88E-07 m/s. and 2.366E-07 (Magnetic) ##########
# Proportion  of total for Rain water = 3.60E-08 / 2.88E-07 = ########## 0.125. and 0.124. (Magnetic) ##########
# Proportion  of total for Process water = 2.52E-07 / 2.88E-07 = ########## 0.875. and 0.875.(Magnetic) ##########
# Retention time for rainwater leach only = 0.01m / 3.596E-08 = ########## 278086.7. and 338638.6 (Magnetic) ##########
# Retention time for process water leach only = 0.01m / 2.520E-07 = ########## 39682.5. and 48309.2.(Magnetic)##########
# Retention time for combined water leach = 0.01m / 2.88E-07 = ########## 34722.2. and 42265.4.(Magnetic) ##########


#SS Antimagnetic   Velocity in pores (m_pw/s)   Proportion of total   Retention time in cells (-time_step/s)
#Rain water      3.60E-08            0.12            278086.7   # This is retention time for rainwater leach only
#Process water   2.52E-07            0.88            39682.5   # This is fictitious, for calculation only
#Combined      2.88E-07            1            34722.2   # This is retention time for Rainwater + DAC




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

# MIXing two or more aqueous solutions in specified proportions. This procedure effectively increases the total mass (or volume) of the aqueous phase but maintains the same concentrations. Often the question of the compatibility of the supplied water with used water needs to be addressed. Supplied water can be mixed it with used water. Mixing results in a batch-reaction calculation, which produces aqueous equilibrium, including redox equilibrium. Mixing may be used alone, in combination with additional reactions, or during advection or transport calculations. All applications of MIX result in a batch-reaction calculation that produces aqueous equilibrium, including redox equilibrium. Mixing solution represents 1 L for each solution meaning 50% of each. In mixing, each solution is multiplied by its mixing fraction and a new solution is calculated by summing over all the fractional solutions. The mass of water is multiplied by the same fractions.

USE equilibrium_phases 1

SAVE SOLUTION 0

END

## USE SOLUTION 0

# SOLUTION 0, the leaching solution (https://speciation.net/Glossary/column-leaching-;938), is combination: rainwater with (process/grid water + DAC)

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

# Evaporation of mixture of (process water, DAC) and rainwater

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

END

############################################################################################################################# GEOCHAMICAL NODE 3 - DAC INTRODUCTION AND AMBIENT AIR #########################################################################################################################
#############################################################################################################################################################################################################################################################################################################

# DAC is expected to raise the concentration in the gas from 0.04% atmospheric CO2(atm) to >3% CO2(g) .Fugacity of CO2(g) in this gas is at minimum 0.03. Gas phase equilibrium incorporated into calculation using DAC avg 3% CO2, ppm_CO2/Air log10(fugacity) = log10(0.03)= -1.5230. These are the numbers used in the PHREEQC calculations.
# The saturation state of oxygen is more problematic. Oxygen readily dissolves in aqueous solutions exposed to atmosphere. Equilibrate water with O2(g) 209460ppm = 0.20946= SI(CO2(g) = log10(fugacity=0.20946) =-0.679. (https://en.wikipedia.org/wiki/Atmosphere_of_Earth). Thus, the anticipated SI of oxygen is log10(0.209) = -0.679.
# Addition of DAC to process and rain water

EQUILIBRIUM_PHASES 101 DAC CO2 capture to 3%
    CO2(g)    -1.523 10 # Calculated CO2 SI= -1.523, 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 # Calculated  N2 , log 10 (0.7809), SI= -0.107, Amount of 10 mol available.
    Calcite_ppt   0 0

# REACTION 1 # ferric iron assumed as present in negligible quantities
#    Fe+3      1
#    1e-08 moles

###################################################################################################### SETUP OUTPUT0 PUNCH TSV FILE ANTIMAGNETIC ###################################################################################################
###################################################################################################################################################################################################################################################



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

END

#REACTION 1 # ferric iron assumed as present in negligible quantities
#    Fe+3      1
#    1e-08 moles


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

# -Vm, Identifier for the molar volume of the solid phase. It identifies for parameters used to calculate the specific volume (cm3/mol) of aqueous species with a Redlich-type equation (Redlich and Meyer, 1964). The molecular weight divided by the density of the solid at 25 ?C. Mineral volume is directly disproportional with Porosity which means if slag property volume decreases (-) , porosity increases.
# Porosity is defined as the ratio of the volume of voids (the spaces between the grains) in a rock to the bulk volume (the overall volume of the voids plus the grains) of that rock (Amyx et al., 1960). In this case rock is replaced by Slag heap.  Porocity (%) = (1-bulk density slag/particle density slag) ? 100
# The mineral volume increased by Calcite_ppt   1.5205e-03  cm^3. Total change in mineral volume:    1.5205e-03  cm^3    = 0.0015 cm^3 which means the porosity decreases. Increase in volume is  0.0015 / 1000= 0.0000015 L is very small.     

############################################################################################################################ DEFINE RATES IN SLAG SYSTEM #######################################################################################################################################################
################################################################################################################################################################################################################################################################################################################

# RATES
# The rate equation for Larnite in carbfix is reproduced here for study purposes.
# Since the rate is embedded in the carbfix2022.dat it does not ordinary have to be reproduced.

# Larnite   #Ca2SiO4;M 172.237 g/mol, from chemical composition database.

##------------------Kinetic calculation---------------------##
      # parameters
#1000 Aa  = k_acid = 5.25e8 # mol.m-2.s-1
#1001 Ac  = k_base = 8.25e5 # mol.m-2.s-1
#1002 An  = k_neut = 0      # mol.m-2.s-1
#1003 Ea  = eapp_acid = 70400 # J/mol
#1004 Ec  = eapp_base = 60900 # J/mol
#1005 En  = eapp_neut = 0     # J/mol
#1006 R =  r = 8.314 #J.deg-1.mol-1
#1007 na = n_acid = 0.44
#1008 nc = n_base = 0.22
#1009 temp = 1/TK
#1010 hplus = ACT("H+")
#1011 a0 = PARM(1)# parm(1) is surface area in the unit of m2/kgw
#1012 parm(2) is a correction factor# 1e-3 converts mmol to mol
#1013 v = 1
#1014 M is current moles of minerals
#1015 M0 is the initial moles of minerals
#1016 sr_mineral = SR("Larnite")
#1017 si_mineral = SI("Larnite")

      #rate equations
#2000 r_acid = k_acid*EXP((-eapp_acid/r)*temp)*(hplus^n_acid)# acid rate expression
#2001 r_neut = k_neut*EXP((-eapp_neut/r)*temp) # neutral rate expression
#2002 r_base = k_base*EXP((-eapp_base/r)*temp)*(hplus^n_base)# base rate expression
#2003 r_all = r_acid+r_neut+r_base
#2004 rate = (a0/v)*(M/M0)^0.67*r_all*(1-sr_mineral)*PARM(2) # 2/3*si_mineral)=0.67
#2005 moles = rate*TIME
#2006 SAVE moles
#-end

# KINETICS 1-100
#-------------------------Mineral---------------------------------
# Larnite   #
             #-formula  Ca2SiO4 1
#----------------------------Mass---------------------------------
             #-m0       9.883 #moles Larnite per litre pore-water
#-------------------------Surface area----------------------------
# The mineral surface area is specified using 3 parameters.
# The first parameter specifies if the specific surface area is calculated per g of rock (0) or kg of water (1)
# The second parameter specifies the specific surface area of the mineral in m2/g
# The third parameters define how the surface area changes during dissolution and has three possible values:
#   0: The surface area is constant
#   1: The surface area changes linearly with the moles of the mineral present
#   2: The surface area changes according to the geometry of dissolving cubes
#   3: The surface area changes according to the geometry of dissolving spheres

# Here, parameter #1 is 1, so that parameter #2 is the surface area
# in m2 per kg of water. See readme file for details.

          #-parms    1 2.5 2# 2.500 m2_Larnite/L_pores - calcs from Erdinc's data
             #-tol      1e-08
#-------------------------Mineral---------------------------------
 #Brucite
             #-formula  Mg(OH)2 1
#----------------------------Mass---------------------------------
             #-m0       4.733 #moles Brucite per litre pore-water

#-------------------------Surface area----------------------------
          #-parms    1 0.556 2# 0.556 m2_Brucite/L_pores - calcs from Erdinc's data
             #-tol      1e-08
#-------------------------Mineral---------------------------------
 #Calcite_diss
             #-formula  CaCO3 1
#----------------------------Mass---------------------------------
             #-m0       9.46 #moles Calcite per litre pore-water

#-------------------------Surface area----------------------------
          #-parms    1 1.683 2# 1.683 m2_Calcite/L_pores - calcs from Erdinc's data
             #-tol      1e-08
#-------------------------Mineral---------------------------------
 #Portlandite
             #-formula  Ca(OH)2 1

#----------------------------Mass---------------------------------
             #-m0       14.389 #moles Calcite per litre pore-water

#-------------------------Surface area----------------------------
          #-parms    1 2.303 2# 2.303 m2_Calcite/L_pores - calcs from Erdinc's data
             #-tol      1e-08


              #-steps       1 # Time steps over which to integrate the rate expressions(s). The -steps identifier is used only during batch-reaction calculations; it is not needed for transport calculations. By default, the list of time steps are considered to be independent times, all starting from zero.
                #-step_divide 0.1 # This parameter affects integration by the Runge-Kutta solver. If step_divide is greater than 1.0, the first time interval of each integration is set to time step / step_divide; at least two time intervals must be integrated to reach the total time of time step?0 to time step / step_divide and time step / step_divide to time step. If step_divide is less than 1.0, then step_divide is the maximum moles of reaction that can be added during a kinetic integration subinterval. Frequently reaction rates are fast initially, thus requiring small time intervals to produce an accurate integration of the rate expressions.
                #-runge_kutta 6 # This parameter affects integration by the Runge-Kutta solver. It designates the preferred number of time subintervals to use when integrating rates and is related to the order of the integration method. A value of 6 specifies that a 5th order embedded Runge-Kutta method, which requires six intermediate rate evaluations, will be used for all integrations. For values of 1, 2, or 3, the program will try to limit the rate evaluations to this number.
                #-bad_step_max 500 # Defines the maximum number of attempts at integrating a set of kinetic reactions over a time step. For the Runge-Kutta method, it is the maximum number of times the integration fails in integrating over a time interval. For the CVODE method, it is the number of times that CVODE is invoked in integrating over a time interval. Default is 500.
                #-cvode true # Specifies whether to use the explicit Runge-Kutta method or the implicit CVODE method (Cohen and Hindmarsh, 1996) to integrate the kinetic rate equations. Default is false. A value of true indicates the CVODE implicit integration method is used.
                #-cvode_steps 100 # Specifies the maximum number of steps that will be taken during one invocation of CVODE. Default is 100 if -cvode_steps is not included.
                #-cvode_order 5 # Specifies the number of terms to use in the extrapolation of rates when using the CVODE method. Default is 5.


 
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


# Moles of a kinetic reactant(HCO3).Kinetic_reactants?Identifier allows definition of a list of kinetically controlled reactants for which two values are written to the selected-output file: (1) the current moles of the reactant, and (2) the moles transferred of the reactant. Optionally, kin, -k[inetics. Carbon capture by charge balance (as generation of bicarbonate ion = HCO3.  Bicarbonate (HCO3?) is most commonly employed as the electrolyte anion because it is known to facilitate CO2 electroreduction.In particular,  HCO3-(Seifritz, 1990; Lackner et al., 1995; Schuiling and Krijgsman, 2006; K?hler et al., 2010; Kelemen et al., 2011; Hartmann et al., 2013). This process of increasing
# HCO3- export is known as enhanced weathering (Schuiling and Krijgsman, 2006). The majority of studies concerning enhanced weathering have evaluated the potential of mafic silicate rocks (e.g., basalt) and minerals (e.g., forsterite, Mg2SiO4 and wollastonite, CaSiO3) as reactants (Schuiling and Tickell, 2010; Moosdorf et al., 2014; Renforth et al., 2015; Taylor et al., 2016; Griffioen, 2017; Montserrat et al., 2017; Rigopoulos et al., 2018; Andrews and Taylor, 2019; Haque et al., 2019)

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


############################################################################################################################################################### DEFINE KINETIC PHASES IN ANTIMAGNETIC SLAG SYSTEM ###############################################################################################
#################################################################################################################################################################################################################################################################################################################

# Here, the initial amounts of minerals in the slag system are defined
# The KINETICS keyword specifies which of the cells in the 1D flow path contain the kinetic minerals
# Kinetics blocks refer to the RATES expressions.
# -parms: parm(1) is surface area in the unit of m2/kgw, parm(2)is an adjustment factor on the rate expression
# CLOSED POND


KINETICS 1-100
Portlandite
    -formula  Ca(OH)2  1
    -m        14.43
    -m0       14.43
    -parms    2.303 1
    -tol      1e-08
Larnite
    -formula  Ca2SiO4  1
    -m        9.87
    -m0       9.87
    -parms    2.499 1
    -tol      1e-08
Calcite_diss
    -formula  CaCO3  1
    -m        9.49
    -m0       9.49
    -parms    1.682 1
    -tol      1e-08
Brucite
    -formula  Mg(OH)2  1
    -m        4.8
    -m0       4.8
    -parms    0.556 1
    -tol      1e-08
kin_HCO3
    -formula  HCO3-1  1
    -m        0
    -m0       0
    -tol      1e-08

# Define the Differential Equation solver, with parameters

-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



############## GRAPHICAL OUTPUT GENERATED WHILST PHREEQC RUNS FOR ANTIMAGNETIC ##############
#############################################################################################

### THE FOLLOWING GRAPHS ARE PRESENTED WHILST PHREEQC RUNS
### GRAPH 1 AND GRAPH 2 ARE THE QUANTITY OF SLAG MATERIALS DISSOLVED IN EACH STEP, AND THE SI IN EACH STEP.
### GRAPH 3 IS THE pH OVER THE ENTIRE COLUMN. THESE GRAPHS ARE EXPECTED TO HAVE A "SAWTOOTH" PROFILE
### GRAPH 4 AND 5 DEFINING MOLES CO2 CAPTURED VS DAYS AND TIME
### KIN_DELTA is the change in the amount of each kinetic reactant. So, the function KIN("Portlandite") is the total amount of Portlandite available at any time; its initial value is set with -M in KINETICS (1 mole by default). If portlandite decreases, KIN("Portlandite") decreases and KIN_DELTA is negative. The reaction will stop when there is no kinetic reactant goes to zero.

USER_GRAPH 1
    -headings               d_Portlandite d_Larnite d_Brucite d_Calcite
    -axis_titles            "Days" "Dissolved mnrl (mol/L)" ""
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - DISSOLVED MINERALS"
    -axis_scale y_axis      auto 0.005 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
  1 REM do not Save as emf
  5 REM IF(TOTAL_TIME = 0) THEN GOTO 100
  7 REM (CELL_NO < 11) THEN GOTO 100
 10 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Portlandite"), color = Green, line_width = 2, y-axis = 1
 20 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Larnite"), color = Blue, line_width = 2, y-axis = 1
 30 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Brucite"), color = Orange, line_width = 2, y-axis = 1
 40 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Calcite_diss"), color = Magenta, line_width = 2, y-axis = 1
100 REM END
  -end
    -active                 true


###################################################################################################
###################################################################################################


USER_GRAPH 2
    -headings               SI_Portlandite SI_Larnite SI_Brucite SI_Calcite
    -axis_titles            "Days" "Saturation Index"
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - SATURATION INDICES"
    -axis_scale y_axis      -20 20 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 1 REM do not Save as emf
 5 REM IF(TOTAL_TIME = 0) THEN GOTO 100
 7 REM IF(CELL_NO < 11) THEN GOTO 100
10 PLOT_XY TOTAL_TIME/86400, SI("Portlandite"), color = Green, line_width = 2, y-axis = 1
20 PLOT_XY TOTAL_TIME/86400, SI("Larnite"), color = Blue, line_width = 2, y-axis = 1
30 PLOT_XY TOTAL_TIME/86400, SI("Brucite"), color = Orange, line_width = 2, y-axis = 1
40 PLOT_XY TOTAL_TIME/86400, SI("Calcite"), color = Magenta, line_width = 2, y-axis = 1
100 REM END
  -end
    -active                 true


###################################################################################################
###################################################################################################


USER_GRAPH 3
    -headings               pH
    -axis_titles            "Days" "pH"
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - pH EVOLUTION"
    -axis_scale y_axis     auto auto auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 1 REM do not Save as emf
 5 REM IF(TOTAL_TIME = 0) THEN GOTO 100
90 PLOT_XY TOTAL_TIME/86400, -LA("H+"), color = Red, line_width = 2, y-axis = 1
100 REM END
  -end
    -active                 true


###################################################################################################
###################################################################################################


USER_GRAPH 4
    -headings               time HCO3
    -axis_titles            "Time, seconds" "Moles CO2 captured" ""
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - MOLES CO2 CAPTURED"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
  1 REM do not Save as emf
  7 REM IF(CELL_NO < 11) 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")+1*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)
100 GRAPH_X TOTAL_TIME # Time interval for which moles of reaction are calculated in rate programs, automatically set in the time-step algorithm of the numerical integration method, in seconds.
110 GRAPH_Y kin_hco3
  -end
    -active                 true


###################################################################################################
###################################################################################################


USER_GRAPH 5
    -headings               Depth HCO3
    -axis_titles            "Distance, in meters" "Moles CO2 captured" ""
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - MOLES CO2 CAPTURED"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
  1 REM do not Save as emf
  7 REM IF(CELL_NO < 11) THEN GOTO 100
 10 IF (DIST = 0) THEN PUT(0,1)
 20 kin_ca=1*KIN_DELTA("Portlandite")+2*KIN_DELTA("Larnite")+0*KIN_DELTA("Brucite")+1*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)
100 GRAPH_X DIST
110 GRAPH_Y kin_hco3
  -end
    -active                 true


################## TRANSPORT - THROUGH SLAG HEAP PHREEQC RUNS FOR ANTIMAGNETIC ##########################################################
#########################################################################################################################################
# TRANSPORT keyword data block is used to simulate one-dimensional (1D) transport of solutes, water, colloids, and heat due to the processes of advection and dispersion, diffusion, and diffusion into stagnant zones adjacent to the 1D flow system. Purely advective transport plus reactions?without diffusion, dispersion, or stagnant zones?can be simulated with the ADVECTION data block.
# The model can be envisioned as a series of cells each of which has the same pore volume, filled with an initial fluid, which then is shifted (transported) from the first cell into the adjacent cell and so forth down the line of cells. In transport simulations the three most important variables are the time step, shifts and cells. ?The time step is the amount of time it takes for water to travel from one cell to the next; the ?shifts? define how many times water travels (shifts) from one cell to the next; and the length(s) define the length of the cells. For each cell, the velocity is the cell length divided by the time step. It takes the same length of time to traverse any of the cells, but the cells may have different lengths; therefore, the velocities may change as you move through the column.
# Need to run at least as many shifts as possible to have cells for the front to move through the entire column? (Parkhurst, 2002).
# A column with 10 cells is filled with a Kardemir Make Up and Rainwater solution that equilibrates with the CO2 and O2. Then, a solution flushes three pore volumes through the column, which requires 33 shifts.
# PHREEQC calculation first does advection, which is simply moving the water from cell 0 to cell 1, cell 1 to cell 2, ... cell n-1 to cell n. Following advection and aqueous equilibration is performed. The surface water will propagate through the system in time steps


# Depth or height of cell (-lengths)   0.01   metres = 1cm
# Number of horizontal slices   100
# No ponding on the surface of the slag heap is assumed, therefore the advective ingress of rainfall may be estimated as follows. Total annual millimetres of rainfall can be extracted to be 522.1 (mm/a)(https://en.wikipedia.org/wiki/Karab%C3%BCk). A broad average of rainfall, in terms of metres per second is thus calculated:
# Total annual millimetres of rainfall 522.1   mm/a
# Convert to m/a= 522.1/1000= 0.522   m/a
# Year = 365.25   d/a
# Convert to m/d= 0.522/ 365.25 = 1.429E-03   m/d
# Convert to s/d = 60x60x24 = 86400   s/d
# Convert to m/s = 1.429E-03 / 86400 = 1.654E-08   m/s. Annual average precipitation is thus 1.654E-08 m/s.
# Slag porosity for Antimagnetic is calculated at: 46% pore volume/slag bulk volume and 56% pore volume on Magnetic.see supporting document

# Advective velocity of rainfall through the Antimagnetic slag heap is thus 1.654E-08/0.46 m/s. ########## Advective velocity of rainfall through the slag heap is 3.596E-08 m/s. and 2.953E-08 (Magnetic) ##########
# The slag heap is 1 metre high and is constantly irrigated with Process-water/grid water and rainwater. The irrigation rate is 10 millimetres per day. The velocity of the irrigation water before it enters the heap is 10   mm_irrigation/day.
# Convert to milimetres per second   10mm/86400sec = 1.160E-04 mm_irrigation/second.
# Convert to metres per second irrigation   1.160E-04/ 1000 = 1.160E-07   m_irrgation/second.
# Porosity of Anti Magnetic slag 0.46 L_pores/L_slag. ########## Thus the velocity of pore waters is   1.160E-07/0.46 = 2.520E-07 m/s. and 2.07E-07 (Magnetic)##########
# Combined velocity of pore waters = 3.596E-08 m/s + 2.520E-07 m/s = ########## 2.88E-07 m/s. and 2.366E-07 (Magnetic) ##########
# Proportion  of total for Rain water = 3.60E-08 / 2.88E-07 = ########## 0.125. and 0.124. (Magnetic) ##########
# Proportion  of total for Process water = 2.52E-07 / 2.88E-07 = ########## 0.875. and 0.875.(Magnetic) ##########
# Retention t

[dlparkhurst]

Sorry, there is no way I will work through your lengthy script. I will answer specific questions about simple scripts.

Also, please use the # button between the forum text box and the subject to define scripts.

[Erdinc]

Hello David,

Thanks for your prompt response. The reason I attached the whole script because you may want to see how I built.

I am trying to create CO2 capture through weathering steel slag. wanted to verify the answer with amount of mineral dissolved (e.g. k_dissolved mineral). From output I can see some changes on some minerals (e.g. calcite_dissolved) but cannot verify clearly as the end of simulation value came back to initial mols available. Can you please help me understand or guide me if I am doing the right approach on comparing? I am trying to attach xls file but keeps giving me error.

[dlparkhurst]

I will look at a script that is less than 200 lines and that is entered with the # button above the text box. Eliminate all SELECTED_OUTPUT and USER_PUNCH data blocks. Do not define minerals that are in the database. Eliminate all the comments.

PHREEQC maintains mole balance; I am 99.9% sure that there is no loss of mass in the calculations.

Looking at the first few lines, I don't think you need to define Calcite_diss and Calcite_ppt. There is a calcite definition in the database. Use either EQUILIBRIUM_PHASES or KINETICS for a mineral, but not both.

[Erdinc]

Thanks David,

Here is the shortened script with #. The reason I split Calcite is that I wanted to define calcite in precipitation and how much in dissolution. Thanks in advance.

Code: [Select]


SOLUTION 1 Kardemir Process/grid water
    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 #
    #N(5)      0.84 # 0.84% ppm or mg/l of Nitrogen
    N(-3)      0.35 # Kardemir NH3 0.2-0.5
    #Na        22 # 11-33 mg/l,total dissolved Na
    #K         0.82 # total dissolved K
    P         0.22 # 0.3% ppm or mg/l of Phosforous, Kardemir P2O5 0.15-0.30
    Fe(2)        0.25
    Fe(3)        1e-9 # ferric iron assumed as present in negligible quantities
    Alkalinity 240
    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

     Calcite_ppt   0 0 # Switch on calcite as a phase that can precipitate within days

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

    H2O        -1
    11.1012 moles in 4 steps

END

EQUILIBRIUM_PHASES 101 DAC CO2 capture to 3%
    CO2(g)    -1.523 10 # Calculated CO2 SI= -1.523, Amount of 10 mol available
    O2(g)     -0.679 10 # Calculated  O2 SI= -0.679, Amount of 10 mol available
    Calcite_ppt   0 0



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        14.43
    -m0       14.43
    -parms    2.303 1
    -tol      1e-08
Larnite
    -formula  Ca2SiO4  1
    -m        9.87
    -m0       9.87
    -parms    2.499 1
    -tol      1e-08
Calcite_diss
    -formula  CaCO3  1
    -m        9.49
    -m0       9.49
    -parms    1.682 1
    -tol      1e-08
Brucite
    -formula  Mg(OH)2  1
    -m        4.8
    -m0       4.8
    -parms    0.556 1
    -tol      1e-08
kin_HCO3
    -formula  HCO3-1  1
    -m        0
    -m0       0
    -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


USER_GRAPH 1
    -headings               d_Portlandite d_Larnite d_Brucite d_Calcite
    -axis_titles            "Days" "Dissolved mnrl (mol/L)" ""
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - DISSOLVED MINERALS"
    -axis_scale y_axis      auto 0.005 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
  1 REM do not Save as emf
  5 REM IF(TOTAL_TIME = 0) THEN GOTO 100
  7 REM (CELL_NO < 11) THEN GOTO 100
 10 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Portlandite"), color = Green, line_width = 2, y-axis = 1
 20 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Larnite"), color = Blue, line_width = 2, y-axis = 1
 30 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Brucite"), color = Orange, line_width = 2, y-axis = 1
 40 PLOT_XY TOTAL_TIME/86400, KIN_DELTA("Calcite_diss"), color = Magenta, line_width = 2, y-axis = 1
100 REM END
  -end
    -active                 true


USER_GRAPH 2
    -headings               SI_Portlandite SI_Larnite SI_Brucite SI_Calcite
    -axis_titles            "Days" "Saturation Index"
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - SATURATION INDICES"
    -axis_scale y_axis      -20 20 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 1 REM do not Save as emf
 5 REM IF(TOTAL_TIME = 0) THEN GOTO 100
 7 REM IF(CELL_NO < 11) THEN GOTO 100
10 PLOT_XY TOTAL_TIME/86400, SI("Portlandite"), color = Green, line_width = 2, y-axis = 1
20 PLOT_XY TOTAL_TIME/86400, SI("Larnite"), color = Blue, line_width = 2, y-axis = 1
30 PLOT_XY TOTAL_TIME/86400, SI("Brucite"), color = Orange, line_width = 2, y-axis = 1
40 PLOT_XY TOTAL_TIME/86400, SI("Calcite"), color = Magenta, line_width = 2, y-axis = 1
100 REM END
  -end
    -active                 true

USER_GRAPH 3
    -headings               pH
    -axis_titles            "Days" "pH"
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - pH EVOLUTION"
    -axis_scale y_axis     auto auto auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
 1 REM do not Save as emf
 5 REM IF(TOTAL_TIME = 0) THEN GOTO 100
90 PLOT_XY TOTAL_TIME/86400, -LA("H+"), color = Red, line_width = 2, y-axis = 1
100 REM END
  -end
    -active                 true


USER_GRAPH 4
    -headings               time HCO3
    -axis_titles            "Time, seconds" "Moles CO2 captured" ""
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - MOLES CO2 CAPTURED"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
  1 REM do not Save as emf
  7 REM IF(CELL_NO < 11) 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)
100 GRAPH_X TOTAL_TIME # Time interval for which moles of reaction are calculated in rate programs, automatically set in the time-step algorithm of the numerical integration method, in seconds.
110 GRAPH_Y kin_hco3
  -end
    -active                 true


USER_GRAPH 5
    -headings               Depth HCO3
    -axis_titles            "Distance, in meters" "Moles CO2 captured" ""
    -chart_title            "KARDEMIR ANTIMAGNETIC SLAG FLOW THROUGH TESTS - MOLES CO2 CAPTURED"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
  1 REM do not Save as emf
  7 REM IF(CELL_NO < 11) THEN GOTO 100
 10 IF (DIST = 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)
100 GRAPH_X DIST
110 GRAPH_Y kin_hco3
  -end
    -active                 true

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                101 # 1 passes of pore water through the slag heap
    -lengths               0.01 # 1cm cells
    -time_step             34722.2 The total time of the simulation (40.5896 days).
    -print_cells           1, 100-101
  # -print_cells           10-11
    -print_frequency       303 # Print only at end
  # -print_frequency       33 # Print only at end


[dlparkhurst]

Your file does not run because you have not defined Calcite_ppt in PHASES, but before I go on, this kinetic reaction is incorrect.

Code: [Select]

kin_HCO3
    -formula  HCO3-1  1
    -m        0
    -m0       0
    -tol      1e-08

Charge is ignored in a -formula for kinetics, so the effect of this reaction is to add, for example
CO2, 0.5H2O, and 0.5O2.
So, it will be an oxidation reaction.

The -formula should be charge balanced, and it represents elements that are added or removed from the system. I don't know what you have in mind with this reaction. 

Once you have fixed the HCO3 kinetic reaction and have a file that runs, please use just one or two graphs to show me what your specific question is.