Reservoir Modeling.

[serhatttt]

Hi, Dear Prof. Parkhurst

I want to simulate a reservoir based on different years using phreeqc. I prepared the code, but the calcite, anorthite, and alunite minerals go to zero after 100 years. And also, Kaolinite goes to zero after 10000 years. Is the code correct? I have one more question. Also, in this geothermal well, we observe calcite scaling at about 140 m. The depth of the geothermal well is 580 m. The reservoir temperature is 180 degrees, and the reservoir pressure is 50 atm at 580 m. In this well, I want to simulate the thickness of the calcite deposit depending on temperature. Is it possible to simulate this situation? Which parameters do I need to change in the code?

Thank you for helping.
Best regards,

Here is the code;

DATABASE c:\phreeqc\database\llnl.dat

SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
     log_k 16.681
     delta_h -5.738 kcal
    -analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9
    -dw 1.92e-9
    -Vm 7.29 0.92 2.07 -1.23 -1.60 # ref. 1 + McBride et al. 2015, JCED 60, 171
   
2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T
     log_k -1.8
    -analytical_expression 8.68 -0.0103 -2190
    -Vm 14.58 1.84 4.14 -2.46 -3.20
END

PHASES
CO2(g)
CO2 = CO2
     log_k -1.468
     delta_h -4.776 kcal
     -analytical_expression 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
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
END

SOLUTION 1 Geothermal Well
      temp 180
      pH 5.6
      pe 4
      redox pe
      units mg/l
      density 1
      C(4) 135
      Ca 2920
      Cl 34560
      Mg 122
      Na 17250
      S(6) 206.6
      K 1815
      F 4.22
      B 25.2
      Si 195
      Li 21.6
      Fe 5.49
      water 1 # kg
      pressure 50 atm #Reservoir pressure
EQUILIBRIUM_PHASES 1
      CO2(g) 1.6989 0.002 #Partial pressure of carbondioxide at 50 atm, and injected mol
RATES
Quartz
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_qtz = SI("Quartz")
30 if (M <= 0 and si_qtz < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_qtz > 0) then SA = 1e-05 #nucleation
60 k_acid = 0
70 k_neut = 10^(-13.99)*EXP(-87.60e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_qtz))
190 moles = r * TIME
200 SAVE moles
   -end
Anorthite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_anort = SI("Anorthite")
30 if (M <= 0 and si_anort < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_anort > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.5)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.411)
70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_anort))
190 moles = r * TIME
200 SAVE moles
    -end
Albite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alb = SI("Albite")
30 if (M <= 0 and si_alb < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_alb > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))
70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA *(1-(10^si_alb))
190 moles = r * TIME
200 SAVE moles
    -end
K-Feldspar
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kfeld = SI("K-Feldspar")
30 if (M <= 0 and si_kfeld < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_kfeld > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_kfeld))
190 moles = r * TIME
200 SAVE moles
    -end
Corundum
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_cornd = SI("Corundum")
30 if (M <= 0 and si_cornd < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_cornd > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-12.53)*EXP(-47.61e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.150)
70 k_neut = 10^(-15.70)*EXP(-47.61e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-11.75)*EXP(-47.61e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.550)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_cornd))
190 moles = r * TIME
200 SAVE moles
    -end
Enstatite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_enst = SI("Enstatite")
30 if (M <= 0 and si_enst < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_enst > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-9.02)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.600)
70 k_neut = 10^(-12.72)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_enst))
190 moles = r * TIME
200 SAVE moles
    -end
Ilmenite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_ilmen = SI("Ilmenite")
30 if (M <= 0 and si_ilmen < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_ilmen > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-8.35)*EXP(-37.90e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.421)
70 k_neut = 10^(-11.16)*EXP(-37.90e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_ilmen))
190 moles = r * TIME
200 SAVE moles
   -end
Fluorapatite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_fluorap = SI("Fluorapatite")
30 if (M <= 0 and si_fluorap < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_fluorap > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.73)*EXP(-25.00e+04/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.613)
70 k_neut = 10^(-8.00)*EXP(-25.00e+04/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_fluorap))
190 moles = r * TIME
200 SAVE moles
   -end
Kaolinite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kaol = SI("Kaolinite")
30 if (M <= 0 and si_kaol < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_kaol > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-11.31)*EXP(-65.90e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.777)
70 k_neut = 10^(-13.18)*EXP(-22.20e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-17.05)*EXP(-17.90e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.472)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_kaol))
190 moles = r * TIME
200 SAVE moles
   -end
Alunite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alun = SI("Alunite")
30 if (M <= 0 and si_alun < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_alun > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.65)*EXP(-32.30e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.133)
70 k_neut = 10^(-12.53)*EXP(-39.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-12.53)*EXP(-39.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.194)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_alun))
190 moles = r * TIME
200 SAVE moles
   -end
Calcite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_calc = SI("Calcite")
30 if (M <= 0 and si_calc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_calc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_calc))
190 moles = r * TIME
200 SAVE moles
   -end
Anhydrite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_anhyd = SI("Anhydrite")
30 if (M <= 0 and si_anhyd < 0) then goto 200
40 SA = PARM(1)*M
50 if (M = 0 and si_anhyd > 0) then SA = 1e-05 #nucleation
60 k_acid = 0
70 k_neut = 10^(-3.19)*EXP(-14.30e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_anhyd))
190 moles = r * TIME
200 SAVE moles
   -end
Gypsum
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_gyps = SI("Gypsum")
30 if (M <= 0 and si_gyps < 0) then goto 200
40 SA = PARM(1)*M
50 if (M = 0 and si_gyps > 0) then SA = 1e-05 #nucleation
60 k_acid = 0
70 k_neut = 10^(-2.79)*EXP(-0.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_gyps))
190 moles = r * TIME
200 SAVE moles
   -end
Hematite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_hema = SI("Hematite")
30 if (M <= 0 and si_hema < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_hema > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-9.39)*EXP(-66.20e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-14.60)*EXP(-66.20e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_hema))
190 moles = r * TIME
200 SAVE moles
    -end
Dolomite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_dolo = SI("Dolomite")
30 if (M <= 0 and si_dolo < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_dolo > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_dolo))
190 moles = r * TIME
200 SAVE moles
    -end
   
KINETICS 1
Quartz
    -formula SiO2 1
    -m 41.28
    -m0 41.28
    -parms 1.36
    -tol 1e-08
   
Anorthite
    -formula CaAl2(SiO4)2 1
    -m 11.35
    -m0 11.35
    -parms 6.03
    -tol 1e-08

Albite
    -formula NaAlSi3O8 1
    -m 8.55
    -m0 8.55
    -parms 6.02
    -tol 1e-08

K-Feldspar
    -formula KAlSi3O8 1
    -m 45.65
    -m0 45.65
    -parms 6.523
    -tol 1e-08

Corundum
    -formula Al2O3 1
    -m 16.52
    -m0 16.52
    -parms 1.537
    -tol 1e-08
   
Enstatite
    -formula MgSiO3 1
    -m 9.621
    -m0 9.621
    -parms 1.6235
    -tol 1e-08

Ilmenite
    -formula FeTiO3 1
    -m 0.6551
    -m0 0.6551
    -parms 1.916
    -tol 1e-08
   
Hematite
     -formula Fe2O3 1
     -m 9.60
     -m0 9.60
     -parms 1.825
     -tol 1e-08   
   
Kaolinite
    -formula Al2Si2O5(OH)4 1
    -m 10
    -m0 10
    -parms 5.98
    -tol 1e-08

Alunite
    -formula KAl3(OH)6(SO4)2 1
    -m 10
    -m0 10
    -parms 9.07
    -tol 1e-08

Calcite
    -formula CaCO3 1
    -m 10
    -m0 10
    -parms 2.2661
    -tol 1e-08
   
Dolomite
     -formula CaMg(CO3)2 1
     -m 10
     -m0 10
     -parms 3.8685
     -tol 1e-08
     
Anhydrite
    -formula CaSO4 1
    -m 10
    -m0 10
    -parms 2.80
    -tol 1e-08
    -steps 3 30 300 3000 30000 300000 3000000 30000000 300000000 3000000000 30000000000 300000000000
    -cvode true
   
KNOBS
    -convergence_tolerance 1e-10

INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings Time Quartz Anorthite Albite K-Feldspar Corundum Kaolinite Alunite Calcite Anhydrite Hematite Dolomite
    -axis_titles "Log10 Time" "Mineral (mol)"
    -initial_solutions false
    -connect_simulations true
    -plot_concentration_vs x
    -start
10 GRAPH_X log10(total_time)
20 GRAPH_Y kin("Quartz"), kin("Anorthite"), kin("Albite"), kin("K-Feldspar"), kin("Corundum"), kin("Kaolinite"), kin("Alunite"), kin("Calcite"), kin("Anhydrite"), kin("Hematite"), kin("Dolomite")
    -end
   
-active true

SELECTED_OUTPUT 1
    -file TwellKineticyear.sel
    -reset false
    -time true
    -step true
    -ph true
    -pe true
    -totals C(4)
    -molalities K+ Na+ Ca+2 Mg+2 HCO3- SO4-2 Cl- CO3-2 CO2
    -kinetic_reactants Quartz Anorthite Albite K-feldspar Corundum Kaolinite Alunite Calcite Anhydrite Hematite Dolomite
END

[dlparkhurst]

Is it correct? I would say no. The rates seem too fast, and you are precipitating K-spar and other minerals that I think should be dissolving.

Surface area is usually the fitting factor for kinetic reactions, and it is difficult to assess. I don't have much hope that you will be able to predict the kinetic evolution of the formation.

Another problem is determining the water composition at depth if reactions are occurring in the well bore. If the water has been in the formation for geologic time, I would assume that calcite, dolomite, anhydrite, and perhaps other minerals are in equilibrium with the formation water. Your water composition is supersaturated with these minerals, so I think there may be a systematic error in your estimation of the formation water composition.

Assuming your water composition is correct, it is supersaturated with calcite at depth. As the water moves up the well with decreasing temperature and pressure, the water becomes undersaturated with calcite. With the calculation below, CO2 alone does not produce a gas phase, but methane has not been taken into account. I used phreeqc.dat because it has the Peng-Robinson data for gases.

Code: [Select]

SOLUTION 1 Geothermal Well
      temp 180
      pH 5.6
      pe 4
      redox pe
      units mg/l
      density 1
      C(4) 135
      Ca 2920
      Cl 34560
      Mg 122
      Na 17250
      S(6) 206.6
      K 1815
      F 4.22
      B 25.2
      Si 195
      Li 21.6
      Fe 5.49
      water 1 # kg
      pressure 50 atm #Reservoir pressure
END
GAS_PHASE
-fixed_pressure
CO2(g) 0
END
USE solution 1
USE gas_phase 1
EQUILIBRIUM_PHASES
Calcite 0 0
REACTION_TEMPERATURE
180 25 in 10
REACTION_PRESSURE
50 2 in 10
USER_GRAPH 1
    -headings               depth Calcite Temp Pressure
    -axis_titles            "Depth" "Calcite moles" "T and P"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 5 depth =  - (TC - 25)/(180 - 25) * 580
10 GRAPH_X depth
20 GRAPH_Y EQUI("Calcite")
30 GRAPH_SY TC, PRESSURE
  -end
    -active                 true
END


[serhatttt]

Dear Prof. Parkhurst,

Thank you so much for your kind help. Yes, this field is a little bit complex. Electric conductivity values of the geothermal water can reach up to 80.000-90.000 uS/cm in this geothermal field.
This simulation you provided is very useful for me. Many thanks for this.

How can I simulate the thickness (in mm) of the calcite deposit with respect to depth in this well? Is it possible to simulate the thickness of the deposit in a geothermal well using phreeqc? The depth of the well is 580 m. The total flow rate of geothermal water in this well is 360 cubic meters per hour, and the diameter of the well is 50 cm. I couldn't find any example that includes this data in phreeqc.

[dlparkhurst]

Given your formation solution composition, the following calculates the amount of calcite that could precipitate. I leave it to you to decide how to distribute the volume precipitated and account for the flow rate and the surface area of the pipe.

Code: [Select]

SOLUTION 1 Geothermal Well
      temp 180
      pH 5.6
      pe 4
      redox pe
      units mg/l
      density 1
      C(4) 135
      Ca 2920
      Cl 34560
      Mg 122
      Na 17250
      S(6) 206.6
      K 1815
      F 4.22
      B 25.2
      Si 195
      Li 21.6
      Fe 5.49
      water 1 # kg
      pressure 50 atm #Reservoir pressure
END
USE solution 1
EQUILIBRIUM_PHASES
Calcite 0 0
USER_PRINT
10 PRINT "Moles of calcite per liter fm water:        ", EQUI("Calcite")
20 PRINT "Molar volume calcite, cm^3/mol:             ", PHASE_VM("Calcite")
30 PRINT "Volume of calcite per liter fm water, cm^3: ", EQUI("Calcite")*PHASE_VM("Calcite")
END

That said, I still do not think that the formation water is supersaturated with calcite, so I am skeptical of the above calculation. If the well produces gas, it could be that calcite precipitates because of loss of CO2. There is a slight decrease in solubility from a pressure decrease alone, but an increase in solubility with a temperature decrease. I think you need to consider why it is that calcite is precipitating in your well.

[serhatttt]

Dear Prof. Parkhurst,

Many thanks for your kind effort and time.