Co2 mineralization

[Amr khan]

Hello, I am new to PHREEQC and still developing my understanding of its concepts. I am working on the effect of storage efficiency due to changes in porosity and permeability caused by carbon mineralization in a basalt reservoir. However, when I run my input file, I encounter an error related to porosity calculation (not giving any change in output also).

Code: [Select]

PHASES
Magnesite
    MgCO3 = Mg+2 + CO3-2
    log_k     -8.029
    delta_h   -6.169 kcal

# Simulation at 25?C
SOLUTION 1  # Initial formation water at 25?C
    temp      25
    pH        7.5
    pe        4
    units     mol/kgw
    Ca        0.0132
    Mg        0.0165
    Fe        0.0003
    Al        1e-6
    Si        0.0006
    Na        0.0054
    Cl        0.0148
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    CO2(g)    1.398  10.0    # 25 bar PCO2
    Calcite   0.0    0.0
    Magnesite 0.0    0.0
    Dolomite  0.0    0.0
   
KINETICS 1
    -steps    0 1 2 4 8 24 48 72 days  # Time steps in days
   
    Plagioclase
        -formula  CaAl2Si2O8
        -m0       0.45        # Initial moles
        -m        0.45
        -parms    65.0  1e-12 # Ea (kJ/mol), k25 (mol/m2/s)
        -tol      1e-8
   
    Pyroxene
        -formula  MgSiO3
        -m0       0.35
        -m        0.35
        -parms    40.0  3e-11
        -tol      1e-8
   
    Olivine
        -formula  Mg2SiO4
        -m0       0.18
        -m        0.18
        -parms    67.2  6.9e-10
        -tol      1e-8

RATES
    Plagioclase
        -start
        10  if TK <= 0 then TK = 273.15 + 25  # Ensure TK is initialized
        20  R = 1e-12 * exp(-65000/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.5
        30  moles = R * TIME
        40  save moles
        -end
   
    Pyroxene
        -start
        10  if TK <= 0 then TK = 273.15 + 25  # Ensure TK is initialized
        20  R = 3e-11 * exp(-40000/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.7
        30  moles = R * TIME
        40  save moles
        -end
   
    Olivine
        -start
        10  if TK <= 0 then TK = 273.15 + 25  # Ensure TK is initialized
        20  R = 6.9e-10 * exp(-67200/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.5
        30  moles = R * TIME
        40  save moles
        -end

USER_PUNCH
    -headings Time pH SI_Cal SI_Mag SI_Dol Por Perm Storage CO2_min
    -start
    10  REM Initial parameters
    20  phi_0 = 0.15          # Initial porosity
    30  k_0 = 1e-13          # Initial permeability (m^2)
    40  rho_CO2 = 65.7       # CO2 density at 25?C, 25 bar (kg/m^3)
   
    50  REM Calculate mineral volumes
    60  vol_min = TOT("Calcite")*36.9 + TOT("Magnesite")*28.0 + TOT("Dolomite")*64.3
   
    70  REM Update porosity
    80  new_phi = phi_0 - vol_min/1000
   
    90  REM Kozeny-Carman permeability
    100 new_k = k_0 * (new_phi/phi_0)^3 * ((1-phi_0)/(1-new_phi))^2
   
    110 REM CO2 storage capacity
    120 storage = new_phi * rho_CO2 + (1-new_phi) * 0.3 * rho_CO2
   
    130 REM CO2 mineralized
    140 co2_min = (TOT("Calcite") + TOT("Magnesite") + 2*TOT("Dolomite")) * 44.01
   
    150 REM Output results
    160 punch TOTAL_TIME
    170 punch -LA("H+")
    180 punch SI("Calcite")
    190 punch SI("Magnesite")
    200 punch SI("Dolomite")
    210 punch new_phi
    220 punch new_k
    230 punch storage
    240 punch co2_min
    -end

SELECTED_OUTPUT
    -file                 results_25C.xls
    -reset               false
    -time               true
    -pH                 true
    -alkalinity         true
    -ionic_strength     true
    -totals             Ca Mg Fe Al Si C(4)
    -saturation_indices Calcite Magnesite Dolomite
    -kinetic_reactants  Plagioclase Pyroxene Olivine




[dlparkhurst]

There are several posts in the forum related to porosity changes.

The Basic function TOT applies to the total molality of an element in solution. EQUI applies to equilibrium phases and KIN applies to kinetic reactants. Look at the section The Basic Interpreter in the manual for more details and more Basic functions.

[Amr khan]

Okay... Thank you.

[Amr khan]

I amended a few steps but am still not getting the right result. If you can help me with this i will be grateful
Thank you 

Code: [Select]

PHASES
Magnesite
    MgCO3 = Mg+2 + CO3-2
    log_k     -8.029
    delta_h   -6.169 kcal

# Simulation at 25?C
SOLUTION 1  # Initial formation water at 25?C
    temp      25
    pH        7.5
    pe        4
    units     mol/kgw
    Ca        0.0132
    Mg        0.0165
    Fe        0.0003
    Al        1e-6
    Si        0.0006
    Na        0.0054
    Cl        0.0148
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    CO2(g)    1.398  10.0    # 25 bar PCO2
    Calcite   0.0    0.0
    Magnesite 0.0    0.0
    Dolomite  0.0    0.0
   
KINETICS 1
    -steps    0 1 2 4 8 24 48 72 days  # Time steps in days
   
    Plagioclase
        -formula  CaAl2Si2O8
        -m0       0.45        # Initial moles
        -m        0.45
        -parms    65.0  1e-12 # Ea (kJ/mol), k25 (mol/m2/s)
        -tol      1e-8
   
    Pyroxene
        -formula  MgSiO3
        -m0       0.35
        -m        0.35
        -parms    40.0  3e-11
        -tol      1e-8
   
    Olivine
        -formula  Mg2SiO4
        -m0       0.18
        -m        0.18
        -parms    67.2  6.9e-10
        -tol      1e-8

RATES
    Plagioclase
        -start
        10  if TK <= 0 then TK = 273.15 + 25  # Ensure TK is initialized
        20  R = 1e-12 * exp(-65000/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.5
        30  moles = R * TIME
        40  save moles
        -end
   
    Pyroxene
        -start
        10  if TK <= 0 then TK = 273.15 + 25  # Ensure TK is initialized
        20  R = 3e-11 * exp(-40000/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.7
        30  moles = R * TIME
        40  save moles
        -end
   
    Olivine
        -start
        10  if TK <= 0 then TK = 273.15 + 25  # Ensure TK is initialized
        20  R = 6.9e-10 * exp(-67200/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.5
        30  moles = R * TIME
        40  save moles
        -end

USER_PUNCH
    -headings Time pH SI_Cal SI_Mag SI_Dol Por Perm Storage CO2_min
    -start
    10  REM Initial parameters
    20  phi_0 = 0.15          # Initial porosity
    30  k_0 = 1e-13          # Initial permeability (m^2)
    40  rho_CO2 = 65.7       # CO2 density at 25?C, 25 bar (kg/m^3)
   
    50  REM Calculate mineral volumes
    60  vol_min = EQUI("Calcite")*36.9 + EQUI("Magnesite")*28.0 + EQUI("Dolomite")*64.3
   
    70  REM Update porosity
    80  new_phi = phi_0 - vol_min/1000
   
    90  REM Kozeny-Carman permeability
    100 new_k = k_0 * (new_phi/phi_0)^3 * ((1-phi_0)/(1-new_phi))^2
   
    110 REM CO2 storage capacity
    120 storage = new_phi * rho_CO2 + (1-new_phi) * 0.3 * rho_CO2
   
    130 REM CO2 mineralized
    140 co2_min = (EQUI("Calcite") + EQUI("Magnesite") + 2*EQUI("Dolomite")) * 44.01
   
    150 REM Output results
    160 punch TOTAL_TIME
    170 punch -LA("H+")
    180 punch SI("Calcite")
    190 punch SI("Magnesite")
    200 punch SI("Dolomite")
    210 punch new_phi
    220 punch new_k
    230 punch storage
    240 punch co2_min
    -end

SELECTED_OUTPUT
    -file                 results_123.xls
    -reset               false
    -time               true
    -pH                 true
    -alkalinity         true
    -ionic_strength     true
    -totals             Ca Mg Fe Al Si C(4)
    -saturation_indices Calcite Magnesite Dolomite
    -kinetic_reactants  Plagioclase Pyroxene Olivine



[dlparkhurst]

First, you need to decide on the volume of your cell. Two possibilities are (1) you have 1 L of solution, or (2) 1 L of porous media. If you have porosity of 0.15, then in case (1) you have 1/0.15 = 6.67 L of porous media; 5.57 L of rock, and 1 L of water. In case (2) you have 1 L of porous media, 0.85 L of rock, and 0.15 L of water.

Your input file as 1 kg of water, so about 1 L of water, which would translate to a representative volume (RV) of 6.67 L.

However you define RV, you can use the following definitions to calculate the change in porosity from before reaction (p) to after reaction (p').

Let dn(i) represent the change in moles of mineral i due to reaction,
vm(i) the molar volume of mineral i, cm^3/mol,

then the change in porosity is

Code: [Select]

p' = p - [sum(dn(i)*vm(i))/1000)]/RV

If you consider only reactions with equilibrium phases, each dn(i) is given by the function EQUI_DELTA. If you are also considering kinetic reactions, then dn for kinetic reactions is given by KIN_DELTA, and you need to account for all of the minerals, both equilibrium and kinetic, that are reacting.

[Amr khan]

Hi, Thankyou for your guidance.
I added this porosity calculation but still in my result there is no change in porosity, permeability, and storage, although mineralization occurs. The porosity (Por) column shows all zeros (0.00E+00)
The permeability (Perm) column shows all zeros (0.00E+00)
The storage capacity (Storage) stays constant at 19.7 kg/m?

Code: [Select]

PHASES
Magnesite
    MgCO3 = Mg+2 + CO3-2
    log_k     -8.029
    delta_h   -6.169 kcal

# Simulation at 50?C with improved conditions
SOLUTION 1
    temp      50
    pH        8.5
    pe        4
    units     mol/kgw
    Ca        0.015      # Slightly increased for carbonate formation
    Mg        0.018
    Fe        0.0003
    Al        1e-6
    Si        0.0006
    Na        0.0054
    Cl        0.0148
    C(4)      2e-3      # Increased bicarbonate to buffer pH
    -water    1

EQUILIBRIUM_PHASES 1
    CO2(g)    1.398  10.0  # 25 bar PCO2
    Calcite   0.0    0.1   # Ensure presence for precipitation
    Magnesite 0.0    0.1
    Dolomite  0.0    0.1
   
KINETICS 1
    -steps    0 86400 172800 345600 691200 2073600 4147200 6220800 8208000 10368000 13392000  # Converted days to seconds

    Plagioclase
        -formula  CaAl2Si2O8
        -m0       0.45 
        -m        0.45
        -parms    65.0  5e-13  # Reduced dissolution rate
        -tol      1e-8

    Pyroxene
        -formula  MgSiO3
        -m0       0.35
        -m        0.35
        -parms    40.0  1e-11  # Adjusted reaction rate
        -tol      1e-8

    Olivine
        -formula  Mg2SiO4
        -m0       0.18
        -m        0.18
        -parms    67.2  3.5e-10 # Adjusted to promote precipitation
        -tol      1e-8

RATES
    Plagioclase
        -start
        10 if TK <= 0 then TK = 273.15 + 25
        20 R = 5e-13 * exp(-65000/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.5
        30 moles = R * TIME
        40 save moles
        -end

    Pyroxene
        -start
        10 if TK <= 0 then TK = 273.15 + 25
        20 R = 1e-11 * exp(-40000/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.7
        30 moles = R * TIME
        40 save moles
        -end

    Olivine
        -start
        10 if TK <= 0 then TK = 273.15 + 25
        20 R = 3.5e-10 * exp(-67200/8.314*(1/TK - 1/298.15)) * (1 + 10^(-pH))^0.5
        30 moles = R * TIME
        40 save moles
        -end

USER_PUNCH
    -headings Time pH SI_Cal SI_Mag SI_Dol Por Perm Storage CO2_min
    -start
    10  REM Initial conditions
    20  phi_0 = 0.15          # Initial porosity
    30  k_0 = 1e-13           # Initial permeability (m^2)
    40  rho_CO2 = 65.7        # CO2 density at 50?C, 25 bar (kg/m^3)
    50  RV = 1000             # Representative volume in cm3 (1 L)

    60  REM Molar volumes (cm3/mol)
    70  vm_olivine    = 43.0
    80  vm_magnesite  = 28.0
    90  vm_calcite    = 36.9
    100 vm_dolomite   = 64.3
    110 vm_plagio     = 100.0
    120 vm_pyroxene   = 31.0

    130 REM Calculate mineral volumes - use TOT() for absolute amounts
    140 vol_calc = EQUI("Calcite")*vm_calcite
    150 vol_mag = EQUI("Magnesite")*vm_magnesite
    160 vol_dol = EQUI("Dolomite")*vm_dolomite
    170 vol_olv = KIN("Olivine")*vm_olivine
    180 vol_pyr = KIN("Pyroxene")*vm_pyroxene
    190 vol_pla = KIN("Plagioclase")*vm_plagio
   
    200 REM Sum all mineral volumes
    210 vol_total = vol_calc + vol_mag + vol_dol + vol_olv + vol_pyr + vol_pla
   
    220 REM Update porosity - scale properly with representative volume
    230 new_phi = max(0.01, phi_0 - vol_total/RV)
   
    240 REM Kozeny-Carman permeability
    250 new_k = k_0 * (new_phi/phi_0)^3 * ((1-phi_0)/(1-new_phi))^2
   
    260 REM CO2 storage capacity
    270 storage = new_phi * rho_CO2 + (1-new_phi) * 0.3 * rho_CO2
   
    280 REM CO2 mineralized
    290 co2_min = (TOT("Calcite") + TOT("Magnesite") + 2*TOT("Dolomite")) * 44.01
   
    300 REM Output results
    310 punch TOTAL_TIME
    320 punch -LA("H+")
    330 punch SI("Calcite")
    340 punch SI("Magnesite")
    350 punch SI("Dolomite")
    360 punch new_phi
    370 punch new_k
    380 punch storage
    390 punch co2_min
    -end

SELECTED_OUTPUT
    -file               results_94D.xls
    -reset              false
    -time               true
    -pH                 true
    -alkalinity         true
    -ionic_strength     true
    -totals             Ca Mg Fe Al Si C(4)
    -saturation_indices Calcite Magnesite Dolomite
    -kinetic_reactants  Plagioclase Pyroxene Olivine

PRINT
    -saturation_indices Calcite Magnesite Dolomite

[dlparkhurst]

There is no Basic function MAX; I think your usage is probably interpreted as a 2-dimensional array with a value of 0.

You can use PRINT statements to debug your RATES and USER_PUNCH calculations.

[Amr khan]

Thank you. Do I need to use different CO₂ densities and pressure when running simulations at different temperatures?

[dlparkhurst]

You are specifying the CO2 partial pressure.

The temperature and total pressure of the system is the value defined in solution (pressure default 1 atm), or you can set the temperature and pressure with REACTION_TEMPERATURE and REACTION_PRESSURE. The aqueous density, log Ks, and activity coefficients of the solution will adjust based on these temperature and pressure values (for phreeqc.dat, Amm.dat or pitzer.dat).

There is also a GAS_PHASE capability with either a fixed_volume or fixed_pressure gas phase. Pressure is adjusted for fixed volume and volume is adjusted for fixed pressure as the gas phase equilibrates with the solution.

[Amr khan]

Thank you for your time.