Dolomite and calcite precipitation in geothermal system

[AntonioG]

Hello,

I have been trying to estimate the amount/volume of mineral phases which can potentially dissolve and precipitate when a geothermal water from a closed system is re-injected in the same a calcareous aquifer (deep Dogger Limestone) at lower temperatures. I have the following data:
- groundwater composition measured in the lab at ~25C and Eh measured at the bore head.
- xrd data of the aquifer mineral composition (calcite 71%, dolomite 11%, Albite 6%, Qz 6%, K-Feld 6%)
- gas phase pressure (175 bar) composition (Ar 0.45%, CH4 15%, CO2 16.9%, H2S 1.14%) and gas liquid ratio measured (0.0495) measured at the bore head at 57.5C.

First I equilibrate the groundwater composition with measured Eh (pe), mineral phases, pressure and partial pressures of the gas at the initial reservoir temperature (57.5C). Next I incrementally decrease the temperature to 20C.

As expected silicates precipitate at lower temperatures, but some carbonate minerals like calcite tend to also precipitate with the decreasing temperature, while dolomite dissolves. ?????? I have tried phreeqc, minteqv4, thermoden and llnl database but I still get the same results. I can't identify what I have done wrong, any ideas? Thanks.


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

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

SOLUTION 1 #INITIAL Reservoir Composition !!!!!!! No iron !!!!!!
units mg/L
density 1 calculate
temp 57.5
pe -3.81
redox pe
Alkalinity 415 as HCO3
pH   6.0
Al   0.01
B   12.6
Ba   0.4
Ca   907
Cl   14117 #charge
F   3.6
# C(4)   415
Li   2.5
Mg   343
Mn   0.03
Na   8126
Si   26.8
S(6)   686
Sr   80
water 1


EQUILIBRIUM_PHASES 1
pe_fix    3.81  O2    1000 #Eh = -250 mV
Calcite 0 117
Dolomite 0 9.8
Magnesite 0 0
Quartz 0 16.5
Albite_low 0 3.8
# Illite 0.0 2.5
K-Feldspar 0 3.6


GAS_PHASE 1 #Fixed volume gas phase
   -fixed_volume
   -pressure 298 #INITIAL PRESSURE ADJUSTED TO REACH EQUILIBRIUM PRESSURE of 172.71 atm
   #########-pressure 172.71 #PRESSURE AT EQUILIBRIUM (measured)
   -temperature 57.5
        -volume   0.0498 #equal to GLR for a solution of 1L at 60C
   ########## measured Equilibrium pressures
   ########## Ar(g)   1.28
        ########## CH4(g)   14.28
   ########## CO2(g)   16.37
        ########## H2S(g)   1.04
   ########## He(g)      1.08
        ########## N2(g)      143.16
   ########## Initial partial pressures to reach measured Equilibrium pressures
   Ar(g)   2.180
   CH4(g)   24.282
   CO2(g)   27.824
   H2S(g)   1.767
   He(g)   1.833
   N2(g)   243.380
   #Ntg(g)   243.380
   
SAVE SOLUTION 1   
END

USER_GRAPH 4
-headings Dolomite Calcite Strontianite Magnesite Huntite
-axis_scale y_axis auto #-40000 40000
-axis_scale x_axis auto #22 52

10 plot_xy tc, EQUI_DELTA("Dolomite")*184.4*1000,  line_width = 2,symbol_size = 4, color = Red
15 plot_xy tc, (EQUI_DELTA("Calcite")*100*1000),  line_width = 2,symbol_size = 4, color = magenta
25 plot_xy tc, (EQUI_DELTA("Strontianite")*147.73*1000),  line_width = 2,symbol_size = 4, color = green
35 plot_xy tc, (EQUI_DELTA("Magnesite")*84.1*1000),  line_width = 2,symbol_size = 4, color = blue
45 plot_xy tc, (EQUI_DELTA("Huntite")*353*1000),  line_width = 2,symbol_size = 4, color = orange

USER_GRAPH 5
-headings Anhydrite Barite Celestite
-axis_scale y_axis auto #-4 4
-axis_scale x_axis auto #22 52

10 plot_xy tc, EQUI_DELTA("Anhydrite")*102*1000,  line_width = 2,symbol_size = 4, color = green
15 plot_xy tc, (EQUI_DELTA("Barite")*233.4*1000),  line_width = 2,symbol_size = 4, color = blue
25 plot_xy tc, (EQUI_DELTA("Celestite")*183.8*1000),  line_width = 2,symbol_size = 4, color = red

USER_GRAPH 6
-headings Chalcedony Illite KFeldspar Albite Muscovite Kaolinite Quartz
-axis_scale y_axis auto #-150000 150000
-axis_scale x_axis auto #22 52

10 plot_xy tc, EQUI_DELTA("Chalcedony")*60*1000,  line_width = 2,symbol_size = 4, color = black
15 plot_xy tc, (EQUI_DELTA("Illite")*389.3*1000),  line_width = 2,symbol_size = 4, color = red
35 plot_xy tc, (EQUI_DELTA("K-Feldspar")*278*1000),  line_width = 2,symbol_size = 4, color = blue
45 plot_xy tc, (EQUI_DELTA("Albite_low")*263*1000),  line_width = 2,symbol_size = 4, color = green
55 plot_xy tc, (EQUI_DELTA("Muscovite")*398.7*1000),  line_width = 2,symbol_size = 4, color = orange
65 plot_xy tc, (EQUI_DELTA("Kaolinite")*258.1*1000),  line_width = 2,symbol_size = 4, color = magenta
75 plot_xy tc, (EQUI_DELTA("Quartz")*60*1000),  line_width = 2,symbol_size = 4, color = yellow
85 plot_xy tc, (EQUI_DELTA("Montmor-Na")*549*1000),  line_width = 2,symbol_size = 4, color = red


SELECTED_OUTPUT 1
    -file                 selected_output_mass_and_volumes_v2.xls
    # -file                 selected_output_SI_v2.xls
    -reset                true
    -simulation           false
    -state                false
    -solution             false
    -reaction             false
    -equilibrium_phases   Dolomite Calcite Strontianite Magnesite Witherite Barite Celestite Illite Chalcedony Quartz K-Feldspar Albite_low Muscovite Kaolinite Montmor-Na Gibbsite Siderite Anhydrite Gypsum
    -saturation_indices   Dolomite Calcite Strontianite Magnesite Witherite Barite Celestite Illite Chalcedony Quartz K-Feldspar Albite_low Muscovite Kaolinite Montmor-Na Gibbsite Siderite Anhydrite Gypsum

USER_PUNCH 1
   -headings    pH _Temperature Alkalinity(mg/L_CaCO3) Dolomite Calcite Witherite Strontianite Magnesite Barite Celestite Illite Chalcedony Quartz K-Feldspar Albite Kaolinite Montmorillonite-Na Gibbsite Siderite Dolomite Calcite Witherite Strontianite Magnesite Barite Celestite Illite Chalcedony Quartz K-Feldspar Albite Kaolinite Montmorillonite-Na Gibbsite Siderite Ba Ca Cl Fe K Mg Mn Na O SO4 Si_ Sr_  CO2
   -start
           1 PUNCH -LA("H+")
      2 PUNCH tc
      3 PUNCH alk*100*1000
      ## masses in mg/kgw
      10 PUNCH EQUI_DELTA("Dolomite")*184.4*1000
      15 PUNCH EQUI_DELTA("Calcite")*100*1000
      17 PUNCH EQUI_DELTA("Witherite")*197.37*1000   
      20 PUNCH EQUI_DELTA("Strontianite")*147.74*1000
      24 PUNCH EQUI_DELTA("Magnesite")*84.1*1000
      35 PUNCH EQUI_DELTA("Barite")*233.4*1000
      40 PUNCH EQUI_DELTA("Celestite")*183.8*1000
      45 PUNCH EQUI_DELTA("Illite")*389.3*1000
      50 PUNCH EQUI_DELTA("Chalcedony")*60*1000
      52 PUNCH EQUI_DELTA("Quartz")*60*1000
      55 PUNCH EQUI_DELTA("K-Feldspar")*278*1000
      57 PUNCH EQUI_DELTA("Albite_low")*263*1000
      60 PUNCH EQUI_DELTA("Kaolinite")*258.1*1000
      61 PUNCH EQUI_DELTA("Montmor-Na")*549*1000
      62 PUNCH EQUI_DELTA("Gibbsite")*78*1000
      63 PUNCH EQUI_DELTA("Siderite")*115.8*100
      ## volumes in mL/kgw
      65 PUNCH EQUI_DELTA("Dolomite")*184.4*1000/2.86
      70 PUNCH EQUI_DELTA("Calcite")*100*1000/2.71
      72 PUNCH EQUI_DELTA("Witherite")*197.37*1000/3.75
      75 PUNCH EQUI_DELTA("Strontianite")*147.74*1000/3.78
      76 PUNCH EQUI_DELTA("Magnesite")*84.1*1000/3.05
      85 PUNCH EQUI_DELTA("Barite")*233.4*1000/4.48
      90 PUNCH EQUI_DELTA("Celestite")*183.8*1000/3.95
      95 PUNCH EQUI_DELTA("Illite")*389.3*1000/2.75
      100 PUNCH EQUI_DELTA("Chalcedony")*60*1000/2.61
      101 PUNCH EQUI_DELTA("Quartz")*60*1000/2.65
      105 PUNCH EQUI_DELTA("K-Feldspar")*278*1000/2.60
      110 PUNCH EQUI_DELTA("Albite_low")*263*1000/2.62
      120 PUNCH EQUI_DELTA("Kaolinite")*258.1*1000/2.65
      122 PUNCH EQUI_DELTA("Montmor-Na")*549*1000/2.5
      125 PUNCH EQUI_DELTA("Gibbsite")*78*1000/2.41
      127 PUNCH EQUI_DELTA("Siderite")*115.8*1000/3.96
      150 PUNCH TOT("Ba")*137.2*1000
      210 PUNCH TOT("Ca")*40*1000
      215 PUNCH TOT("Cl")*35.5*1000
      220 PUNCH TOT("Fe")*56*1000
      240 PUNCH TOT("K")*39*1000
      255 PUNCH TOT("Mg")*24*1000
      260 PUNCH TOT("Mn")*55*1000    
      265 PUNCH TOT("Na")*23*1000   
      270 PUNCH TOT("O(0)")*16*1000
      275 PUNCH TOT("S(6)")*96*1000
      277 PUNCH TOT("Si")*28*1000
      280 PUNCH TOT("Sr")*87.62*1000   
      300 PUNCH MOL("CO2")*44*1000   

USE SOLUTION 1
EQUILIBRIUM_PHASES 2
Calcite 0 117
Dolomite 0 9.8
Strontianite 0 0.0
Magnesite 0 0
######### Quartz 0 16.5
#######  Huntite 0.0 0.0
Witherite 0 0
#######  Anhydrite 0 0.0
Barite 0 0
Celestite 0.0 0.0
Albite_low 0.0 3.8
Gibbsite 0 0
Diaspore 0 0
Chalcedony 0 0.0
Illite 0  0 #2.5
K-Feldspar 0 3.6

Kaolinite 0.0 0.0
Montmor-Na 0 0
Siderite 0 0
pe_fix    3.81  O2     1000
REACTION_TEMPERATURE 1
57.5 20 in 20 steps
END

 

[dlparkhurst]

I would try, as you have, to recreate the formation water. The analysis you have is nearly in equilibrium with several of the minerals--calcite, dolomite, others, which seems reasonable. However, if you add the gas back into that solution, you will be undersaturated with the carbonates. I don't know for sure how degassing and precipitation has affected your analyses. You will have to consider how to best recreate the formation water.

As for the gas phase, I think you should define a fixed-pressure gas phase with the percentages at the temperature, pressure, and volume of measurement. PHREEQC will calculate a number of moles of each component under those conditions. Then you can react at the temperature and pressure of the formation. The GAS_PHASE will then adjust in volume to equilibrate at the formation temperature and pressure.

Once you have the formation water, my first approach would be to look at precipitation of carbonates, sulfates, chalcedony, and kaolinite as you lower the temperature. Define these minerals in EQUILIBRIUM_PHASES with zero moles. That would give you an idea of which scales might form during the re-injection process. Longer term, you can consider interaction of the water with the formation minerals.

[AntonioG]

Hi dlparkhurst
Thanks for your help, the main issue was that I got confused between fixed volume and fixed pressure in the definition of the gas phase.
I might be wrong but I realised that a kinetic approach should be more appropriate because PP/diss kinetics of the silicates is actually slower than the carbonates and the flow velocity in the pore space around the injection bore is high. Do you think it would be relevant to include the following rates with the llnl database? Thanks,

RATES
    Calcite
-start
 10 si_ = SI("Calcite")
 20 IF (M <= 0  and si_ < 0) THEN GOTO 200
 30 e1 = 14.4
 40 k1 = 10^-0.3 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
 50 n1 = 1
 60 e2 = 23.5
 70 k2 = 10^-5.81* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15))
 80 e3 = 35.4
 90 k3 = 10^-3.48* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
100 n3 = 1
110 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
120 rate = area * ( k2  )
130 rate = rate * (1 - SR("Calcite"))
140 moles = rate * 0.001 * TIME
200 SAVE moles
-end
    Dolomite-dis
-start
 10 si_ = SI("Dolomite-dis")
 20 IF (M <= 0  and si_ < 0) THEN GOTO 200
 30 e1 = 14.4
 40 k1 = 10^-0.3 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
 50 n1 = 1
 60 e2 = 36.1
 70 k2 = 10^-3.19* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15))
 80 e3 = 35.4
 90 k3 = 10^-3.48* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
100 n3 = 1
110 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
120 rate = area * ( k2  )
130 rate = rate * (1 - SR("Dolomite-dis"))
140 moles = rate * 0.001 * TIME
200 SAVE moles
-end
    Albite_low
-start
 20 si_ = SI("Albite_low")
 30 IF (M <= 0  and si_ < 0) THEN GOTO 200
 40 e1 = 14.4
 50 k1 = 10^-0.3 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
 60 n1 = 1
 70 e2 = 14.4
 80 k2 = 10^-10.56* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15))
 90 e3 = 69.8
100 k3 = 10^-3.48* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
110 n3 = 1
120 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
130 rate = area * ( k2  )
140 rate = rate * (1 - SR("Albite_low"))
200 moles = rate * 0.001 * TIME
210 SAVE moles
-end
    Illite
-start
 10 si_ = SI("Illite")
 20 IF (M <= 0  and si_ < 0) THEN GOTO 200
 30 e1 = 14.4
 40 k1 = 10^-0.3 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
 50 n1 = 1
 60 e2 = 62.76
 70 k2 = 1e-13* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15))
 80 e3 = 35.4
 90 k3 = 10^-3.48* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
100 n3 = 1
110 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
120 rate = area * ( k2  )
130 rate = rate * (1 - SR("Illite"))
140 moles = rate * 0.001 * TIME
200 SAVE moles
-end
    Anhydrite
-start
 10 si_ = SI("Anhydrite")
 20 IF (M <= 0  and si_ < 0) THEN GOTO 200
 30 e1 = 14.4
 40 k1 = 10^-0.3 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
 50 n1 = 1
 60 e2 = 14.3
 70 k2 = 10^-3.79* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15))
 80 e3 = 35.4
200 k3 = 10^-3.48* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
210 n3 = 1
220 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
230 rate = area * ( k2  )
240 rate = rate * (1 - SR("Anhydrite"))
250 moles = rate * 0.001 * TIME
260 SAVE moles
-end
    Chalcedony
-start
 10 si_ = SI("Chalcedony")
 20 IF (M <= 0  and si_ < 0) THEN GOTO 150
 30 e1 = 14.4
 40 k1 = 10^-0.3 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
 50 n1 = 1
 60 e2 = 87.6
 70 k2 = 10^-13.99* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15))
 80 e3 = 35.4
 90 k3 = 10^-3.48* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
100 n3 = 1
110 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
120 rate = area * ( k2  )
130 rate = rate * (1 - SR("Chalcedony"))
140 moles = rate * 0.001 * TIME
150 SAVE moles
-end
    K-Feldspar
-start
 10 si_ = SI("K-Feldspar")
 20 IF (M <= 0  and si_ < 0) THEN GOTO 150
 30 e1 = 14.4
 40 k1 = 10^-0.3 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
 50 n1 = 1
 60 e2 = 38
 70 k2 = 10^-12.41* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15))
 80 e3 = 35.4
 90 k3 = 10^-3.48* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
100 n3 = 1
110 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
120 rate = area * ( k2  )
130 rate = rate * (1 - SR("K-Feldspar"))
140 moles = rate * 0.001 * TIME
150 SAVE moles
-end

END

[dlparkhurst]

You haven't said whether you are most concerned with scaling on injection or long-term changes in the aquifer, or both.

You have some rate expressions, so you are in a position to simulate either equilibrium or kinetics and compare results. I don't know whether dolomite will be calculated to precipitate, but you may want to restrict dolomite to dissolution because it does not usually form under low-temperature conditions.

[AntonioG]

Hi dlparkhurst thanks for your help.