Reactive Transport Modelling - Saturation Indices

[AntonioG]

Hello,

I trying to simulate the impact of a geothermal fluid heated to 95C and reinjected back into the save reservoir at 50C.

The model was first developed as a batch reaction following these steps:
1- Recalculate the reservoir water composition (at 50C) with pH, Eh, gas pressure and GLR measured at the production bore head.
2- Equilibrate the composition with the mineral phases detected in the reservoir.
3- Raise the temperature to 95C and identify key mineral phases likely to dissolve/PP.
4- Estimate degree of potential clogging/void generated under equilibrium and rate limited diss/PP conditions.

Now I'd like to estimate the potential impact on the reservoir within 55 m around the injection bore (assuming initial reservoir composition at 50C etc...) for a 6-month period.

The temperature of the infill solution raised to 95C initially returns the same saturation indices as for the batch reactions (calcite SI = 0.61, ...).

In the transport simulation I would have expected, for example SI of calcite to decrease from 0.61 (95C) to 0 (50C) with the increasing distance.
However, SI at the infill boundary equal to zero from the first shift and SI in the other cells are negligible and/or negative as the warm fluid infiltrates.
It doesn't seem to make sense as no reactions are included yet. I'm a bit lost, any suggestion/explanation? Thanks.'

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

#DATABASE database\pitzer.dat
# DATABASE database\minteq.v4.dat
#DATABASE database\core10.dat
# DATABASE database\phreeqc.dat
# DATABASE database\LLNL.dat
DATABASE PHREEQC_ThermoddemV1.10_15Dec2020.dat

PHASES
Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
     log_k 1.54
     delta_h 0
     -analytical_expression -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
     
pe_fix
    e- = e-
    log_k     0
ph_fix
    H+ = H+
    log_k     0
   
SOLUTION_MASTER_SPECIES
Ac       Ac-           0.0   60.00      60.00   

SOLUTION_SPECIES
#Acetic acid
Ac- = Ac-
        log_k   0.0;          -gamma  1e7   0.0
Ac- + H2O = HAc + OH-
        log_k   4.75;         -gamma  1e7   0.0
 

SELECTED_OUTPUT 1
    -file                 outputs/selected_output.xls
    -reset                true
    -simulation           false
    -state                false
    -solution             false
    -reaction             false
    # -equilibrium_phases   Calcite Dolomite Magnesite(Natur) Ankerite Chalcedony Quartz(alpha) Phillipsite(K) Albite(low) Kaolinite Illite(Mg) Montmorillonite(MgK) Montmorillonite(MgNa) Sudoite Boehmite Diaspore Gibbsite Actinolite Arsenopyrite Marcassite Pyrite Anhydrite Gypsum
    -saturation_indices  Calcite Dolomite Magnesite(Natur) Ankerite Chalcedony Quartz(alpha) Phillipsite(K) Albite(low) Kaolinite Illite(Mg) Montmorillonite(MgK) Montmorillonite(MgNa) Sudoite Boehmite Diaspore Gibbsite Actinolite Arsenopyrite Marcassite Pyrite Anhydrite Gypsum

USER_PUNCH 1
   -headings    Time_(Days) pH _Temperature Alkalinity(mg/L_CaCO3) Calcite Dolomite Magnesite Ankerite Chalcedony Quartz K-Feldspar Albite Kaolinite Illite Montmorillonite Sudoite Boehmite Diaspore Gibbsite Actinolite Al As Br Ca Cl F K Fe Mg Mn Na Si SO4 CO2
   
   
   Ba Ca Cl Fe K Mg Mn Na O SO4 Si_ Sr_  CO2
   -start
           1 PUNCH TOTAL_TIME/86400
      2 PUNCH -LA("H+")
      3 PUNCH tc
      4 PUNCH alk*100*1000
      # 110 PUNCH EQUI_DELTA("Calcite")*100*1000/2.71
      110 PUNCH KIN_DELTA("Calcite")*100*1000/2.71
      # 115 PUNCH EQUI_DELTA("Dolomite")*184.4*1000/2.86
      115 PUNCH KIN_DELTA("Dolomite")*184.4*1000/2.86
      # 120 PUNCH EQUI_DELTA("Magnesite(Natur)")*84.1*1000/3.05
      120 PUNCH KIN_DELTA("Magnesite(Natur)")*84.1*1000/3.05      
      130 PUNCH EQUI_DELTA("Ankerite")*206.39*1000/3.05
      140 PUNCH EQUI_DELTA("Chalcedony")*60*1000/2.61
      # 150 PUNCH EQUI_DELTA("Quartz(alpha)")*60*1000/2.65
      150 PUNCH KIN_DELTA("Quartz(alpha)")*60*1000/2.65
      #155 PUNCH EQUI_DELTA("Phillipsite(K)")*661.10*1000/2.2
      155 PUNCH KIN_DELTA("Phillipsite(K)")*661.10*1000/2.2
      160 PUNCH EQUI_DELTA("Albite(low)")*263*1000/2.62
      # 165 PUNCH EQUI_DELTA("Kaolinite")*258.1*1000/2.65
      165 PUNCH KIN_DELTA("Kaolinite")*258.1*1000/2.65
      #175 PUNCH EQUI_DELTA("Illite(Mg)")*389.34*1000/2.75
      175 PUNCH KIN_DELTA("Illite(Mg)")*389.34*1000/2.75
      180 PUNCH EQUI_DELTA("Montmorillonite(MgK)")*549*1000/2.5
      185 PUNCH EQUI_DELTA("Sudoite")*546.77*1000/2.65
      190 PUNCH EQUI_DELTA("Bohemite")*60*1000/3.03
      195 PUNCH EQUI_DELTA("Diaspore")*60*1000/3.4
      200 PUNCH EQUI_DELTA("Gibbsite")*78*1000/2.41
      205 PUNCH EQUI_DELTA("Actinolite")*853.16*1000/3.02
      210 PUNCH TOT("Al")*137.2*1000
      215 PUNCH TOT("As")*137.2*1000
      220 PUNCH TOT("Br")*137.2*1000
      225 PUNCH TOT("Ca")*40*1000
      230 PUNCH TOT("Cl")*35.5*1000
      235 PUNCH TOT("F")*137.2*1000
      240 PUNCH TOT("K")*39*1000
      245 PUNCH TOT("Fe")*56*1000
      250 PUNCH TOT("Mg")*24*1000
      255 PUNCH TOT("Mn")*55*1000    
      260 PUNCH TOT("Na")*23*1000      
      265 PUNCH TOT("Si")*28*1000
      270 PUNCH TOT("S(6)")*96*1000
      275 PUNCH MOL("CO2")*44*1000   

SOLUTION 0 #INITIAL Reservoir Composition
-units mg/L
-density 1 calculate
-temp 50
pe -3.61
-redox pe
Alkalinity    459 as HCO3
pH   7.0
Al   0.006
As   0.018
Br      0.71
Ca   19.0
Cl   199 charge
F   4.10
K    9.0
Fe   0.4
Mg   4.2
Mn   0.046
Na   341
N(-3)   0.28   # as NH3
Si   26.9
S(6)   138
-water 1

EQUILIBRIUM_PHASES 0
pe_fix    3.61  O2        100
Calcite 0 22.6
Kaolinite 0 7.26
Illite(Mg) 0 2.16
Quartz(alpha) 0 80.6
Phillipsite(K) 0 0.86
Actinolite 0 0
Albite(low) 0 0
Arsenopyrite 0 0
Marcassite 0 0
Pyrite 0 0
Diaspore 0 0
Gibbsite 0 0
Dolomite 0 0
Chalcedony 0 0
Boehmite 0 0
Montmorillonite(MgK) 0 0
Montmorillonite(MgNa) 0 0
Magnesite(Natur) 0 0
Ankerite 0 0
Sudoite 0 0
Anhydrite 0 0

GAS_PHASE 0 #Fixed volume gas phase
   fixed_pressure
   pressure 2.5
   temperature 50
   volume   0.0344
   Ar(g)   0.0304
   CH4(g)   0.0957
   CO2(g)   1.4085
   He(g)   0.0027
   N2(g)   0.9625
   O2(g)   0.0183

# GAS_PHASE 1 #Fixed volume gas phase
   # fixed_pressure
   # pressure 158
   # temperature 50
   # volume  0.0344
   # Ar(g)   1.94
   # CH4(g) 6.087
   # CO2(g) 89.58
   # He(g)   0.172
   # N2(g)   61.211
   # O2(g)   1.164

SAVE SOLUTION 0
COPY SOLUTION 0 1-182
END
            
# SOLUTION_MODIFY 0
# -temp 95
# SAVE SOLUTION 0
# END
# USE SOLUTION 0

USE SOLUTION 0
REACTION_TEMPERATURE 0
95

# END

TRANSPORT 1-182
    -cells                 182
    -shifts                9
    -time_step             21 day
    -lengths               0.3
    -flow_direction     forward
    -boundary_conditions  flux flux #constant flux
    -dispersivities        0.01
    -correct_disp          true
    -diffusion_coefficient 1.25e-09
    -thermal_diffusion     1#2
    -punch_cells           0-5
    -punch_frequency       1
END

[dlparkhurst]

I think by defining EQUILIBRIUM_PHASES 0 and GAS_PHASE 0, initial solution 0 is always equilibrated with these reactants before entering the column. You should either define them as say user number 1000 so they do not further affect the TRANSPORT calculation, or DELETE them before the TRANSPORT calculation.

Code: [Select]

PHASES
Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
     log_k 1.54
     delta_h 0
     -analytical_expression -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
     
pe_fix
    e- = e-
    log_k     0
ph_fix
    H+ = H+
    log_k     0
   
SOLUTION_MASTER_SPECIES
Ac       Ac-           0.0   60.00      60.00   

SOLUTION_SPECIES
#Acetic acid
Ac- = Ac-
        log_k   0.0;          -gamma  1e7   0.0
Ac- + H2O = HAc + OH-
        log_k   4.75;         -gamma  1e7   0.0
 
END

SOLUTION 0 #INITIAL Reservoir Composition
-units mg/L
-density 1 calculate
-temp 50
pe -3.61
-redox pe
Alkalinity    459 as HCO3
pH   7.0
Al   0.006
As   0.018
Br      0.71
Ca   19.0
Cl   199 charge
F   4.10
K    9.0
Fe   0.4
Mg   4.2
Mn   0.046
Na   341
N(-3)   0.28   # as NH3
Si   26.9
S(6)   138
-water 1
END
EQUILIBRIUM_PHASES 0
pe_fix    3.61  O2        100
Calcite 0 22.6
Kaolinite 0 7.26
Illite(Mg) 0 2.16
Quartz(alpha) 0 80.6
Phillipsite(K) 0 0.86
Actinolite 0 0
Albite(low) 0 0
Arsenopyrite 0 0
Marcassite 0 0
Pyrite 0 0
Diaspore 0 0
Gibbsite 0 0
Dolomite 0 0
Chalcedony 0 0
Boehmite 0 0
Montmorillonite(MgK) 0 0
Montmorillonite(MgNa) 0 0
Magnesite(Natur) 0 0
Ankerite 0 0
Sudoite 0 0
Anhydrite 0 0
END
GAS_PHASE 0 #Fixed volume gas phase
   fixed_pressure
   pressure 2.5
   temperature 50
   volume   0.0344
   Ar(g)   0.0304
   CH4(g)   0.0957
   CO2(g)   1.4085
   He(g)   0.0027
   N2(g)   0.9625
   O2(g)   0.0183
END
USE solution 0
USE equilibrium_phases 0
USE gas_phase 0
SAVE SOLUTION 0-182
END
USE SOLUTION 0
REACTION_TEMPERATURE 0
95
SAVE solution 0
END
DELETE
    -equilibrium_phases 0
    -gas_phase 0
    -reaction_temperature 0
END

TRANSPORT
    -cells                 18 #182
    -shifts                9
    -time_step             21 day
    -lengths               0.3
    -flow_direction     forward
    -boundary_conditions  flux flux #constant flux
    -dispersivities        0.01
    -correct_disp          true
    -diffusion_coefficient 1.25e-09
    -thermal_diffusion     1#2
    #-punch_cells           0-5
    -punch_frequency       9 #1
USER_GRAPH 1
    -headings               dist TC SI(Calcite)
    -axis_titles            "Distance, meters" "Celsius" "SI(Calcite)"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X DIST
20 GRAPH_Y TC
30 GRAPH_SY SI("Calcite")
  -end
    -active                 true
END
As you have defined your gas phase, you probably get unwanted redox reactions. The O2, N2, CH4, and CO2 will react to redox equilibrium; all of the O2 will be instantly consumed. So, look at the output for the solution formed when SOLUTION 0 is reacted with equilibrium phases and the gas phase that is saved as SOLUTION 0. Phreeqc.dat and Amm.dat have unreactive versions Oxg(g), Ntg(g), and Mtg(g) that you may want to use phreeqc.dat/Amm.dat or add the defiitions with SOLUTION_MASTER_SPECIES, SOLUTION_SPECIES and PHASES data blocks. You probably want Ntg to remain inert, but you might want to react Oxg with Mtg or other dissolved constituents with KINETICS. Initially, and for simplicity, you may want to define the GAS_PHASE with the unreactive versions; these gases will partition between the gas phase and solution, but will not participate in any other reactions. 

[AntonioG_v2]

Thanks for your reply dlparkhurst,

I just have an additional question, EQUILIBRIUM_PHASES was included in cells 1-20 at equilibrium with the initial solution at the reservoir temperature (50C).
When I increase the temperature from 50 to 95C quartz precipitates and calcite dissolves, but in the batch reaction approach calcite precipitate and quartz dissolves (as I would expect). I don't understand why the transport approach gives different results as no reaction is involved.


# DATABASE phreeqc.dat
# DATABASE Amm.dat
DATABASE PHREEQC_ThermoddemV1.10_15Dec2020.dat

PHASES   
pe_fix
    e- = e-
    log_k     0
ph_fix
    H+ = H+
    log_k     0
END

SOLUTION 0 #INITIAL Reservoir Composition
-units mg/L
-density 1 calculate
-temp 50
pe -3.61
-redox pe
Alkalinity    459 as HCO3
pH   7.0
Al   0.006
As   0.018
Br      0.71
Ca   19.0
Cl   199 charge
F   4.10
K    9.0
Fe   0.4
Mg   4.2
Mn   0.046
Na   341
N(-3)   0.28   # as NH3
Si   26.9
S(6)   138
-water 1
# COPY SOLUTION 0 1-100
END

EQUILIBRIUM_PHASES 0
pe_fix    3.61  O2        100
Calcite 0 22.6
# Kaolinite 0 7.26
# Illite(Mg) 0 2.16
Quartz(alpha) 0 80.6
# Phillipsite(K) 0 0.86
# Actinolite 0 0
# Albite(low) 0 0
# Arsenopyrite 0 0
# Marcassite 0 0
# Pyrite 0 0
# Diaspore 0 0
# Gibbsite 0 0
# Dolomite 0 0
# Chalcedony 0 0
# Boehmite 0 0
# Montmorillonite(MgK) 0 0
# Montmorillonite(MgNa) 0 0
# Magnesite(Natur) 0 0
# Ankerite 0 0
# Sudoite 0 0
# Anhydrite 0 0
# Goethite 0 0
# Hematite 0 0
END

GAS_PHASE 0 #Fixed volume gas phase
   fixed_pressure
   pressure 2.5
   temperature 50
   volume   0.0344
   Ar(g)   0.0304
   CO2(g)   1.4085
   He(g)   0.0027
   O2(g)   0.0183
   CH4(g)   0.0957
   N2(g)   0.9625     
   # Oxg(g)  0.0183 #O2
   # Mtg(g)  0.0957 #CH4
   # Ntg(g)  0.9625 #N2
END

USE solution 0
USE equilibrium_phases 0
USE gas_phase 0
SAVE SOLUTION 0
SAVE EQUILIBRIUM_PHASES 0
SAVE GAS_PHASE 0
COPY EQUILIBRIUM_PHASES 0 1-20
COPY SOLUTION 0 1-20
COPY GAS_PHASE 0 1-20
END

USE SOLUTION 0
REACTION_TEMPERATURE 0
95
SAVE solution 0
END

DELETE
    -equilibrium_phases 0
    -gas_phase 0
    -reaction_temperature 0
END
   
SELECTED_OUTPUT 1      
    -file                 outputs/transport.xls      
    -reset                true      
    -simulation           false      
    -state                false      
    -solution             false      
    -reaction             false       
    -saturation_indices  Calcite Quartz(alpha)
    -equilibrium_phases  Calcite Quartz(alpha)
USER_PUNCH 1   
     -headings   Tempp Calcite Quartz
     -start     
   2  PUNCH tc
   30 PUNCH EQUI_DELTA("Calcite")*100/2.71/0.00435
   40 PUNCH EQUI_DELTA("Quartz(alpha)")*60/2.65/0.00435
   # 50 PUNCH EQUI_DELTA("Kaolinite")*258.1/2.65/0.00435
   # 60 PUNCH EQUI_DELTA("Illite(Mg)")*389.34/2.75/0.00435
   # 70 PUNCH EQUI_DELTA("Phillipsite(K)")*661.10/2.2/0.00435
   # 80 PUNCH EQUI_DELTA("Hematite")*159.7/5.3/0.00435
   # 90 PUNCH EQUI_DELTA("Goethite")*88.85/3.8/0.00435   
      
TRANSPORT
    cells                 20#100      
    shifts                10#60      
    time_step             0.5 day#3 day      
    lengths               0.2      
    flow_direction     forward      
    boundary_conditions  flux flux #constant flux      
    dispersivities        0.1      
    correct_disp          true      
    diffusion_coefficient 1.25e-09      
    thermal_diffusion    1#2      
    punch_frequency      1 #1      
    porosities         0.24

    # USER_GRAPH 1      
    # -headings               dist TC SI(Quartz)      
    # -axis_titles            "Distance, meters" "Celsius" "SI(Quartz)"      
    # -initial_solutions      false      
    # -connect_simulations    false      
    # -plot_concentration_vs  x      
  # -start      
# 10 GRAPH_X DIST      
# 20 GRAPH_Y TC      
# 40 GRAPH_SY SI("Quartz")      
  # -end      
    # -active                 true
END      


######################## BATCH
SELECTED_OUTPUT 1      
    -file                 outputs/batch.xls      
    -reset                true      
    -simulation           false      
    -state                false      
    -solution             false      
    -reaction             false       
    -saturation_indices  Calcite Quartz(alpha)
    -equilibrium_phases  Calcite Quartz(alpha)
USER_PUNCH 1   
     -headings   Tempp Calcite Quartz
     -start     
   2  PUNCH tc
   30 PUNCH EQUI_DELTA("Calcite")*100/2.71/0.00435
   40 PUNCH EQUI_DELTA("Quartz(alpha)")*60/2.65/0.00435


USE SOLUTION 0
USE EQUILIBRIUM_PHASES 1
REACTION_TEMPERATURE 0
50 95 in 10 steps
SAVE solution 0
END

# USE SOLUTION 0
# USE EQUILIBRIUM_PHASES 1
# REACTION_TEMPERATURE 0
# 50 10 in 10 steps
# SAVE solution 0
# END

[dlparkhurst]

Both transport and batch give the same results. Most of the reaction in the transport calculation occurs in the first cell, when a solution at equilibrium with the equilibrium phases and gas phase at 50 C (but assigned a temperature of 95) enters the first cell and equilibrates with equilibrium phases and gas phase at 95 C. Calcite precipitates and quartz dissolves. The batch reaction results are plotted at distance 0 meters and show the same mole transfer on heating to 95C.

Do not forget the previous comments about redox reactions among N2, CH4, O2, etc.

Code: [Select]

# DATABASE phreeqc.dat
# DATABASE Amm.dat
#DATABASE PHREEQC_ThermoddemV1.10_15Dec2020.dat

PHASES   
pe_fix
    e- = e-
    log_k     0
ph_fix
    H+ = H+
    log_k     0
END

SOLUTION 0 #INITIAL Reservoir Composition
-units mg/L
-density 1 calculate
-temp 50
pe -3.61
-redox pe
Alkalinity    459 as HCO3
pH   7.0
Al   0.006
As   0.018
Br      0.71
Ca   19.0
Cl   199 charge
F   4.10
K    9.0
Fe   0.4
Mg   4.2
Mn   0.046
Na   341
N(-3)   0.28   # as NH3
Si   26.9
S(6)   138
-water 1
# COPY SOLUTION 0 1-100
END

EQUILIBRIUM_PHASES 0
pe_fix    3.61  O2        100
Calcite 0 22.6
# Kaolinite 0 7.26
# Illite(Mg) 0 2.16
Quartz(alpha) 0 80.6
# Phillipsite(K) 0 0.86
# Actinolite 0 0
# Albite(low) 0 0
# Arsenopyrite 0 0
# Marcassite 0 0
# Pyrite 0 0
# Diaspore 0 0
# Gibbsite 0 0
# Dolomite 0 0
# Chalcedony 0 0
# Boehmite 0 0
# Montmorillonite(MgK) 0 0
# Montmorillonite(MgNa) 0 0
# Magnesite(Natur) 0 0
# Ankerite 0 0
# Sudoite 0 0
# Anhydrite 0 0
# Goethite 0 0
# Hematite 0 0
END

GAS_PHASE 0 #Fixed volume gas phase
   fixed_pressure
   pressure 2.5
   temperature 50
   volume   0.0344
   Ar(g)   0.0304
   CO2(g)   1.4085
   He(g)   0.0027
   O2(g)   0.0183
   CH4(g)   0.0957
   N2(g)   0.9625     
   # Oxg(g)  0.0183 #O2
   # Mtg(g)  0.0957 #CH4
   # Ntg(g)  0.9625 #N2
END

USE solution 0
USE equilibrium_phases 0
USE gas_phase 0
SAVE SOLUTION 0
SAVE EQUILIBRIUM_PHASES 0
SAVE GAS_PHASE 0
COPY EQUILIBRIUM_PHASES 0 1-20
COPY SOLUTION 0 1-20
COPY GAS_PHASE 0 1-20
END

USE SOLUTION 0
REACTION_TEMPERATURE 0
95
SAVE solution 0
END

DELETE
    -equilibrium_phases 0
    -gas_phase 0
    -reaction_temperature 0
END

     
TRANSPORT
    cells                 20#100     
    shifts                10#60     
    time_step             0.5 day#3 day     
    lengths               0.2     
    flow_direction     forward     
    boundary_conditions  flux flux #constant flux     
    dispersivities        0.1     
    correct_disp          true     
    diffusion_coefficient 1.25e-09     
    thermal_diffusion    1#2     
    punch_frequency      10 #1     
    punch_cells          1-20
    porosities         0.24

USER_GRAPH 1
    -headings               dist TC Quartz_delta Calcite_delta
    -axis_titles            "Distance, meters" "Celsius" "Moles"
    -axis_scale x_axis      -1 auto auto auto
    -initial_solutions      false
    -connect_simulations    false
    -plot_concentration_vs  x
  -start
10 GRAPH_X DIST     
20 GRAPH_Y TC     
40 GRAPH_SY EQUI_DELTA("Quartz(alpha)"), EQUI_DELTA("Calcite")
  -end
    -active                 true

END     


######################## BATCH


USE SOLUTION 0
USE EQUILIBRIUM_PHASES 1
USE gas_phase 1
REACTION_TEMPERATURE 0
95 #50 95 in 10 steps
SAVE solution 0
END


[AntonioG_v2]

Thanks for your reply.