Processes > Dissolution and precipitation
Correct way to model porosity change
Nourah:
Dear David
I am working on the effects of CO2 on various rock types.
I am especially interested in the evolution of porosity as a function of time and transport (advection and diffusion).
I have come across two different ways to determine porosity, one assumes the pore volume remains at 1000 cm3 (the phreeqc default water volume), the other allows porosity to change (increase or decrease) as minerals dissolve or precipitate.
To start with, I have produced some very simple code that simulates calcite dissolution at high CO2 partial pressure.
I have written code for two user graphs, based on the two different assumptions.
The two different user graphs result in subtly different initial porosity values and quite different degrees of porosity change.
The one that accounts for changed (in this case, decreased) calcite volume results in a four-fold relative increase in porosity compared to the case that assumes fluid volume remains constant at 1000 cm3.
Before I go on to simulate more complex systems with clay and other minerals, as well as diffusion and advection, I thought I would reach out to get your views on the correct way to determine porosity change during simple kinetic batch reactions.
Many thanks,
Nourah AlNajdi
--- Code: ---[RATES
Calcite
-start
1 rem M = current number of moles of calcite
2 rem M0 = number of moles of calcite initially present
3 rem PARM(1) = A/V, cm^2/L
4 rem PARM(2) = exponent for M/M0
10 si_cc = SI("Calcite")
20 if (M <= 0 and si_cc < 0) then goto 200
30 k1 = 10^(0.198 - 444.0 / TK )
40 k2 = 10^(2.84 - 2177.0 / TK)
50 if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
60 if TC > 25 then k3 = 10^(-1.1 - 1737.0 / TK )
70 t = 1
80 if M0 > 0 then t = M/M0
90 if t = 0 then t = 1
100 area = PARM(1) * (t)^PARM(2)
110 rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120 rem 1e-3 converts mmol to mol
130 rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
140 moles = rate * TIME
200 SAVE moles
-end
KINETICS 1
Calcite
-formula CaCO3 1
-m 81
-m0 81
-parms 5 0.3
-tol 1.e-8
-steps 172800 in 48 steps
SOLUTION 1
temp 100
pressure 100
pH 8
pe 4
redox pe
units mol/kgw
density 1
Na 1
Cl 1
-water 1
GAS_PHASE 1
-fixed_volume
-volume 0.1
-pressure 300
-temperature 100
CO2(g) 300
H2O(g) 0.0
SELECTED_OUTPUT 1
-file Calcite_water_kinetics.txt
-reset false
-step true
-pH true
-reaction true
-totals Ca C(4)
-saturation_indices Calcite
-kinetic_reactants Calcite
USER_GRAPH 1 #starts at higher porosity than USER_GRAPH 2 and has a bigger porosity difference - which is correct? 1 or 2
-axis_titles "Time (hrs)" "Porosity" ""
-chart_title "Limestone dissolution: pore volume changed by mineral dissol/pptn"
-axis_scale y_axis 0.25055 0.25074 auto auto
-initial_solutions false
-connect_simulations true
-plot_concentration_vs t
-start
10 new_por = 1000 + 36.93 * (81 - KIN("Calcite"))
20 x = TOTAL_TIME/3600
30 PLOT_XY x, new_por /(new_por + 36.93 * KIN("Calcite"))
-end
-active true
USER_GRAPH 2 #starts at lower porosity than USER_GRAPH 1 and has a smaller porosity difference - which is correct? 1 or 2
-axis_titles "Time (hrs)" "Porosity" ""
-chart_title "Limestone dissolution: porosity kept at 1000 cm3 as per phreeqc volume"
-axis_scale y_axis 0.25055 0.25074 auto auto
-initial_solutions false
-connect_simulations true
-plot_concentration_vs t
-start
20 x = TOTAL_TIME/3600
30 PLOT_XY x, 1000 /(1000 + 36.93 * KIN("Calcite"))
-end
-active true
END]
--- End code ---
dlparkhurst:
--- Code: ---SOLUTION 1
-units mol/kgw
pH 3 charge
Cl 1
-water 0.25
END
EQUILIBRIUM_PHASES
# porosity 0.25 (L water/L aquifer), 0.25 L water, 0.75 L calcite
# calcite density = 2.7 kg/L, 0.75 * 2.7 = 2.0 kg calcite per L aquifer
# gram formula weight of calcite 100 g/mol, 2000 g/100 = 20 mol
# molar volume = 100 g/mol / 2.7 g/cm3 = 37 cm^3/mol
Calcite 0 20
END
USE solution 1
USE EQUILIBRIUM_PHASES 1
USER_PRINT
-start
10 delta_rock_vol = EQUI_DELTA("Calcite")*37/1000
20 new_rock_vol = EQUI("Calcite")*37/1000
30 new_porosity = 1-new_rock_vol
40 PRINT "Delta rock volume: ", STR_F$(delta_rock_vol, 10, 4)
50 PRINT "New rock volume: ", STR_F$(new_rock_vol, 10, 4)
60 PRINT "New porosity: ", STR_F$(new_porosity, 10, 4)
70 END
-end
--- End code ---
Nourah:
Dear Dr Parkhurst
I have been going in loops with porosity calculation in the reactive transport model. I still can't figure out why the original porosity calculation, calculated in USER_GRAPH 4, differs from User_Graph 8. The porosity linked to GET_POR seems to be double accounting for the original porosity, making it unrealistic.
I have been trying a simple basic equilibrium diffusion model of quartz and calcite
with CO2 being introduced in GAS_PHASE
Please help me understand why the two porosity calculations don't match and what GET_POR is doing. I would appreciate your reply.
--- Code: ---# diffusion equilibrium only
PHASES
CO2(g)
CO2 = CO2
-log_k -1.468
-delta_h -4.776 kcal
-analytic 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5
-T_c 304.2 # critical T, K
-P_c 72.86 # critical P, atm
-Omega 0.225 # acentric factor
SiO2(a)
SiO2 + 2 H2O = H4SiO4
-log_k -2.71
-delta_h 3.340 kcal
-analytic -0.26 0.0 -731.0
SOLUTION 0
temp 65
pH 7
pe 4
redox pe
units mg/l
density 1
Al 0.5
C 230 as HCO3
Ca 1300
Cl 32330 charge
K 110
Mg 310
Na 19350
Si 13
Fe 5
-water 1
GAS_PHASE 0
-fixed_pressure
-pressure 300
-volume 1000
-temperature 65
CO2(g) 300
SAVE solution 0
END
SOLUTION 1-50 # Seal
temp 65
pH 7
pressure 300
units mg/l
Al 0.5
C 230 as HCO3
Ca 1300
Cl 32330 charge
K 110
Mg 310
Na 19350
Si 13
Fe 5
-water 1
SAVE SOLUTION 1-50
EQUILIBRIUM_PHASES 1-50 #seal
Calcite 0 10
Quartz 0 95
CO2(g) 2.25 10
USE solution 1
USE equilibrium_phases 1
END
TRANSPORT
-cells 5
-shifts 10
-time_step 10.00 yr
-flow_direction diffusion_only
-boundary_conditions constant closed
-diffusion_coefficient 5e-11
-lengths 5*0.5
-punch_frequency 1
-multi_D true 1e-9 0.27 0 1
-porosities 5*0.27
-punch_cells 1
-print_cells 1
SELECTED_OUTPUT 1
-file diffusion_Rodcell.txt
-reset false
-distance true
-time true
-pH true
-saturation_indices Quartz Calcite
- equilibrium_phases Quartz Calcite
-totals Na Cl Ca K Mg
USER_PUNCH
-headings water_density Pressure mole_fraction_CO2 partial_pressure_CO2 TIME Porosity pH
-start
10 PUNCH RHO
20 PUNCH PRESSURE
30 PUNCH MOL("CO2")
40 PUNCH PR_P("CO2(g)")
50 PUNCH TOTAL_TIME/3.15e7
80 PUNCH CELL_NO
60 PUNCH GET_POR(CELL_NO)
70 PUNCH -La("H+")
100 REM calculate the porosity
120 dv = 23.69 * (EQUI_DELTA("Quartz")) + 36.93 * (EQUI_DELTA("Calcite"))
130 Vtot = EQUI("Calcite") * 36.93 + EQUI("Quartz") * 23.69 + 1000
140 dPor = dv/Vtot
150 new_por = GET_POR(CELL_NO) - dPor
160 CHANGE_POR(new_por,cell_no)
-end
USER_GRAPH 4 # works correctly
-headings Time Por
-axis_titles "Time years" "fractional porosity" ""
-chart_title "CO2 rich water diffuses through mudstone"
-axis_scale x_axis auto auto auto auto log
-initial_solutions false
-connect_simulations false
-plot_concentration_vs t
-start
10 vtot = (EQUI("Calcite") * 36.93) + (EQUI("Quartz") * 23.69)/1000
20 dv = 23.69 * (EQUI_DELTA("Quartz")) + 36.93 * (EQUI_DELTA("Calcite"))
30 new_por = 1000 + dv
40 graph_x TOTAL_TIME/3.15e7
50 graph_y new_por / (new_por + (36.93 * EQUI("Calcite"))+ (23.69 * EQUI("Quartz")))
-active true
USER_GRAPH 8
-headings Time Por
-axis_titles "Time (years)" "Porosity"
-chart_title "Porosity Evolution for All Cells"
-axis_scale x_axis auto auto auto auto log
-initial_solutions false
-connect_simulations false
-plot_concentration_vs t
-start
10 GRAPH_X TOTAL_TIME /3.15e7
20 GRAPH_Y GET_POR(CELL_NO)
-end
-active true
END
TRANSPORT
-shifts 1000
-time_step 10 yr
END
--- End code ---
dlparkhurst:
I think you problem is related to your choice of -punch_cells in TRANSPORT. The USER_PUNCH 1 definition is only calculated for cell number 1. You will have to set -punch_cells to 1-5 to make sure the porosity is redefined in each cell.
Then you need to add some definitions in the USER_GRAPH definitions so you don't get more data plotting than you want. Here is an example that you can modify.
--- Code: ---# diffusion equilibrium only
PHASES
CO2(g)
CO2 = CO2
-log_k -1.468
-delta_h -4.776 kcal
-analytic 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5
-T_c 304.2 # critical T, K
-P_c 72.86 # critical P, atm
-Omega 0.225 # acentric factor
SiO2(a)
SiO2 + 2 H2O = H4SiO4
-log_k -2.71
-delta_h 3.340 kcal
-analytic -0.26 0.0 -731.0
SOLUTION 0
temp 65
pH 7
pe 4
redox pe
units mg/l
density 1
Al 0.5
C 230 as HCO3
Ca 1300
Cl 32330 charge
K 110
Mg 310
Na 19350
Si 13
Fe 5
-water 1
GAS_PHASE 0
-fixed_pressure
-pressure 300
-volume 1000
-temperature 65
CO2(g) 300
SAVE solution 0
END
SOLUTION 1-50 # Seal
temp 65
pH 7
pressure 300
units mg/l
Al 0.5
C 230 as HCO3
Ca 1300
Cl 32330 charge
K 110
Mg 310
Na 19350
Si 13
Fe 5
-water 1
SAVE SOLUTION 1-50
EQUILIBRIUM_PHASES 1-50 #seal
Calcite 0 10
Quartz 0 95
CO2(g) 2.25 10
USE solution 1
USE equilibrium_phases 1
END
TRANSPORT
-cells 5
-shifts 10
-time_step 10.00 yr
-flow_direction diffusion_only
-boundary_conditions constant closed
-diffusion_coefficient 5e-11
-lengths 5*0.5
-punch_frequency 1
-multi_D true 1e-9 0.27 0 1
-porosities 5*0.27
-punch_cells 1-5
-print_cells 1
SELECTED_OUTPUT 1
-file diffusion_Rodcell.txt
-reset false
-distance true
-time true
-pH true
-saturation_indices Quartz Calcite
- equilibrium_phases Quartz Calcite
-totals Na Cl Ca K Mg
USER_PUNCH 1
-headings water_density Pressure mole_fraction_CO2 partial_pressure_CO2 TIME Porosity pH
-start
10 PUNCH RHO
20 PUNCH PRESSURE
30 PUNCH MOL("CO2")
40 PUNCH PR_P("CO2(g)")
50 PUNCH TOTAL_TIME/3.15e7
80 PUNCH CELL_NO
60 PUNCH GET_POR(CELL_NO)
70 PUNCH -La("H+")
100 REM calculate the porosity
120 dv = 23.69 * (EQUI_DELTA("Quartz")) + 36.93 * (EQUI_DELTA("Calcite"))
130 Vtot = EQUI("Calcite") * 36.93 + EQUI("Quartz") * 23.69 + 1000
140 dPor = dv/Vtot
150 new_por = GET_POR(CELL_NO) - dPor
160 CHANGE_POR(new_por,cell_no)
-end
USER_GRAPH 4 # works correctly
-headings Time Por
-axis_titles "Time years" "fractional porosity" ""
-chart_title "CO2 rich water diffuses through mudstone"
-axis_scale x_axis auto auto auto auto log
-initial_solutions false
-connect_simulations false
-plot_concentration_vs t
-start
5 IF CELL_NO <> 1 THEN GOTO 100
10 vtot = (EQUI("Calcite") * 36.93) + (EQUI("Quartz") * 23.69)/1000
20 dv = 23.69 * (EQUI_DELTA("Quartz")) + 36.93 * (EQUI_DELTA("Calcite"))
30 new_por = 1000 + dv
40 graph_x TOTAL_TIME/3.15e7
50 graph_y new_por / (new_por + (36.93 * EQUI("Calcite"))+ (23.69 * EQUI("Quartz")))
100 END
-active true
END
USER_GRAPH 8
-headings Time Por
-axis_titles "Time (years)" "Porosity"
-chart_title "Porosity Evolution for All Cells"
-axis_scale x_axis auto auto auto auto log
-initial_solutions false
-connect_simulations false
-plot_concentration_vs t
-start
5 IF STEP_NO <> 500 THEN GOTO 100
10 GRAPH_X DIST
20 GRAPH_Y GET_POR(CELL_NO)
100 END
-end
-active true
TRANSPORT
-shifts 500 #1000
-time_step 10 yr
END
--- End code ---
Nourah:
Dear Dr Parkhurst,
Thank you very much for your reply. I have tried running the code, but the software keeps crashing and closing every time I run it.
I am using phreeqc the latest version.
Maybe Line number 5 is what causing the code to crash ?
Navigation
[0] Message Index
[#] Next page
Go to full version