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: [Select]

[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]

[dlparkhurst]

Code: [Select]

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



[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: [Select]

# 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

[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: [Select]

# 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


[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 ?

[dlparkhurst]

Runs for me with PhreeqcI version 3.8.5.

Are you running PhreeqcI, batch, or Tony's version from Hydrochemistry.edu?

Latest versions of codes are now being released at https://github.com/usgs-coupled/. For each program, the Release is located on the right-hand-side of the page.

[Nourah]

Dear Dr Parkhurst,

I have downloaded the latest version of Phreeqc which is 3.8.5
The code runs with the new version. But I seem to have the same issue with get por.

Am I doing something wrong with my code or calculation ?
The original calculation calculated based on number of moles is correct, the original porosity is almost 27% with the 10 moles of calcite and 95 quartz, the porosity increases to 31 , as calcite dissolved and part if quartz. But the porosity in get pore goes up to almoat 38 which is unreal.

[dlparkhurst]

I did not check your logic, but it is probably right.

I think these are approximately the initial conditions.

95*23.7/1000 = 2.25 L quartz
10*36.9/1000 = 0.37 L calcite
                         1.    L water

Initially, calculated porosity is 1 / 3.62 = 0.28.

With enough CO2, all the calcite dissolves, and the porosity is 1.37/3.62 = 0.38.

Here is a PHREEQC calculation that shows the same result.

Code: [Select]

# 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 3                                 
-shifts 300                             
-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-3
-punch_frequency 300
-print_cells 1


USER_GRAPH 4                # works correctly
   -headings                Time Porosity Calcite_volume Water_volume
    -axis_titles            "Time years" "Porosity" "Calcite volume, L"
    -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  x
  -start

110 rock_v_original = (10 * 36.93 + 95 * 23.69)/1000   
120 rock_v_current = (EQUI("Calcite") * 36.93 + EQUI("Quartz") * 23.69)/1000   
130 cell_v = 1.0 + rock_v_original
140 water_v_current = 1 + (rock_v_original - rock_v_current)
150 porosity_current = water_v_current / cell_v
155 calcite_v_current = EQUI("Calcite") * 36.93 / 1000
160 graph_x DIST
170 graph_y porosity_current
180 GRAPH_SY calcite_v_current, water_v_current
1000 END
    -active                 true
END