Processes > Dissolution and precipitation
Geothermal Bore Scaling
(1/1)
AntonioG:
I'm trying estimate the impact of barite, quartz and dolomite precipitation on scaling in a geothermal bore as a function of temperature (20-66C) and pressure (2-220 bar). I did a first estimation under equilibrium conditions and now I'm trying to include the kinetics assuming none of these mineral phases are initially present, but I don't understand why barite and quartz don't precipitate despite high surface areas, positive SI values and SR > 1. Any suggestions? Thanks.
--- Code: ---#DATABASE database\phreeqc.dat
DATABASE database\PITZER.dat
PHASES
pe_fix
e- = e-
log_k 0
ph_fix
H+ = H+
log_k 0
RATES
Barite #BaSO4 ; M 233.404 g/mol
-start
1 name$ = "Barite"
2 if (PARM(1) = 0) then goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 100
5 S = PARM(2)*TOT("water")
##-----------------Dissolution and precipitation options---------------------##
100 if (PARM(4) = 0) then goto 1000 else goto 110
110 if (PARM(4) = 1) Then GoTo 150 else goto 200 #
150 If (SR (name$) < 1) Then GoTo 5000 else GoTO 1000 # warning no dissolution reaction
200 If (SR (name$) > 1) Then GoTo 5000 else GoTO 1000 # warning no precipitation reaction
##------------------Kinetic calculation---------------------##
# parameters
1000 Aa = 2.5e-3 #mol.m-2.s-1
1003 Ea = 26000 #J.mol-1
1006 R = 8.314 #J.deg-1.mol-1
1007 ACTI = ACT ("H+")
1008 na = 0.11
1010 Sig = 1
#rate equation
2000 rplusa = Aa * ACTI^na * exp(-Ea/ (R * Tk)) * S
2003 rplus = rplusa
3000 rate = rplus * (1 - (SR ("Barite")^(1/Sig)))
4000 moles = rate * time
5000 save moles
-end
# Barite
# -start
# 20 si_ = SI("Barite")
# 30 IF (M <= 0 and si_ < 0) THEN GOTO 200
# 40 e1 = 30.8
# 50 k1 = 10^-6.9 * EXP(-e1*1000/8.314*(1/(TC+273.15)-1/298.15))
# 60 n1 = 0.22
# 70 e2 = 30.8
# 80 k2 = 10^-7.9* EXP(-e2*1000/8.314*(1/(TC+273.15)-1/298.15)) #in mol.m-2.s-1
# 90 e3 = 0
# 100 k3 = 10^-30* EXP(-e3*1000/8.314*(1/(TC+273.15)-1/298.15))
# 110 n3 = 0
# 120 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M # m2 = m2.mol-1 x mol
# 130 rate = area * ( k2 ) # mol.s-1 = m2 x mol.m-2.s-1
# 140 rate = rate * (1 - SR("Barite"))
# # 200 moles = rate * 0.001 * TIME # convert from mmol to mol
# 200 moles = rate * TIME
# 210 SAVE moles
# -end
Quartz#SiO2; M 60.08 g/mol
-start
1 name$ = "Quartz"
2 if (PARM(1) = 0) then goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 100
5 S = PARM(2)*TOT("water")
##----------------- Dissolution only or precipitation only option---------------------##
100 if (PARM(4) = 0) then goto 1000 else goto 110
110 if (PARM(4) = 1) Then GoTo 150 else goto 200 #
150 if (SR (name$) < 1) Then GoTo 5000 else goto 1000# warning no dissolution reaction
200 if (SR (name$) > 1) Then GoTo 5000 else goto 1000 # warning no precipitation reaction
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa = 4.03e-4#mol/m2/s
1001 Ab = 0.105#mol/m2/s
1002 na = 0.309
1003 nb = -0.41
1004 Ea = 45600
1005 Eb = 80000
1006 R = 8.314 #J.deg-1.mol-1
1007 ACTI = act("H+")
1009 Sig = 1
#rate equations
2002 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(ACTI^na)* S
2003 rplusb = Ab* (exp(-Eb/ (R * Tk)))*(ACTI^nb)* S
2009 rplus = rplusa + rplusb
3000 rate = rplus * (1 - (SR("Quartz")^(1/Sig)))
4000 moles = rate * time
5000 save moles
-end
Dolomite # CaMg(CO3)2: M 184.40 g/mol !!!
-start
1 name$ = "Dolomite"
2 if (PARM(1) = 0) then goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 100
5 S = PARM(2)*TOT("water")
##-----------------Dissolution and precipitation options---------------------##
100 if (PARM(4) = 0) then goto 1000 else goto 110
110 if (PARM(4) = 1) Then GoTo 150 else goto 200 #
150 If (SR (name$) < 1) Then GoTo 5000 else GoTO 1000 # warning no dissolution reaction
200 If (SR (name$) > 1) Then GoTo 5000 else GoTO 1000 # warning no precipitation reaction#
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa =1.2e-3# mol.m-2.s-1
1001 Ac = 650 # mol.m-2.s-1
1002 Ea =10000# J/mol
1003 Eac =65000# J/mol
1004 R = 8.314 #J.deg-1.mol-1
1006 Sig = 1.9
1007 na =0.5
1008 kc =160
1009 act_c = act("HCO3-")+act("CO3-2")
1010 carb_term = 1-(kc*act_c)/(1+kc*act_c)
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2001 rplusc = Ac* (exp(-Eac/ (R * Tk)))*carb_term
2002 rplus = rplusa + rplusc
4000 rate = rplus * (1 - SR("Dolomite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end
##################### WINNIPEG
SOLUTION 1
-units mg/L
-density 1.0 calculate
-temp 52
pe -0.310
-redox pe
pH 7.87
Alkalinity 378
Al 1.800
B 4.445
Ba 0.447
Br 23
Ca 1474
Cl 71333 charge
Fe 1.83
Hg 7.23
K 1272
Li 7.600
Mg 330
Mn 0.913
Na 44806
Pb 1.225
Se 0.700
Si 15.53
S(6) 4199
Sr 38.0
U 0.003
# Zn 1.625
#-water 1
GAS_PHASE 1
-pressure 215.23
-temperature 52 #winnipeg
-volume 0.016
-fixed_pressure
CO2(g) 0.0494
Ntg(g) 207.13
Oxg(g) 0.803
H2Sg(g) 0.637
CO2(g) 8.810
EQUILIBRIUM_PHASES 1
# pe_fix 0.310 O2 100 #H2Sg 100
# ph_fix -5.2 CO2 100
Quartz 0 17.97
# # #K-feldspar 0 0.059
Dolomite 0 0.031
Halite 0 0.394
Anhydrite 0 0.021
# # #Illite 0 0.039
# # #Pyrite 0 0.192
# # #Chlorite(14A) 0 0.085
Chalcedony 0 0
Barite 0 0
SAVE EQUILIBRIUM_PHASES 1
SAVE SOLUTION 1
END
##################### DEADWOOD
SOLUTION 2
-units mg/L
-density 1.0 calculate
-temp 66.67
pe -1.4
pH 7.55
-redox pe
Alkalinity 354.5
Al 4.53
B 6.6425
Ba 0.23
Br 29
Ca 1648.5
Cl 56386 charge
Fe 0.73
Hg 7.55
K 2046
Li 9.8
Mg 368.75
Mn 0.98
Na 49306
Se 4.53
Si 10.05
S(6) 4236
Sr 47.4
Rb 0.9
Zn 3.53
GAS_PHASE 2
-pressure 221.46 #deadwood
-temperature 66.67#deadwood
-volume 0.016 #GLR
-fixed_pressure
Ntg(g) 210
Oxg(g) 0.815
H2Sg(g) 0.646
CO2(g) 8.938
EQUILIBRIUM_PHASES 2
# pe_fix 1.409 H2Sg 100
# ph_fix -5.3 CO2 100
Dolomite 0 0.056
Halite 0 0.532
Anhydrite 0 0.726
# # # Illite 0 0.080
# # # Ca-montmorillonite 0 0.029
# # # siderite 0 0.089 ##in reservoir
Quartz 0 13.93
# # #K-feldspar 0 0.158
Chalcedony 0 0
Barite 0 0
SAVE EQUILIBRIUM_PHASES 2
SAVE SOLUTION 2
END
####### GET AVERAGE PRODUCTION COMPOSITION - mixing ratios based on relative net thicknesses of formations
MIX SOLUTION 3
1 0.279 #fraction from winnipeg formation
2 0.721 #fraction from deadwood formation
SAVE SOLUTION 3
END
REACTION_TEMPERATURE 3
66 20 in 27 steps
END
EQUILIBRIUM_PHASES 3
# pe_fix 1.409 H2Sg 100
# ph_fix -5.3 CO2 100
# Dolomite 0 0
# Halite 0 0
# Anhydrite 0 0
# Quartz 0 0
# Calcite 0 0
# Chalcedony 0 0
# Gypsum 0 0
# Barite 0 0
END
KINETICS 3
Barite
-formula BaSO4 1
-m0 0
-m 0
-parms 0 1.07 0 0#0 0.000107 0 0 #Castillo 2015, Kazmierczak 2022
-tol 1e-08
Dolomite
-formula CaMg(CO3)2 #1
-m0 0
-m 0
-parms 0 0.0012 0 0
-tol 1e-08
Quartz
-formula SiO2 #1
-m0 0
-m 0
-parms 0 6.8 0 0 #0 0.00068 0 0
-tol 1e-08
-steps 10 day
# -bad_step_max 1000
-cvode true
-cvode_steps 100
-cvode_order 5
END
SELECTED_OUTPUT 1
-file outputs/Regina_output.xls
-reset true
-simulation false
-state false
-solution false
-reaction false
-saturation_indices Barite Quartz Anhydrite Dolomite
USER_PUNCH 1
-headings pH _Temperature Pressure Alkalinity(mg/L_CaCO3) Calcite Halite Anhydrite Gypsum KIN_barite KIN_quartz KIN_dolomite SR_Barite SR_QUARTZ SR_DOLOMITE
-start
1 PUNCH -LA("H+")
2 PUNCH tc
3 PUNCH pressure
4 PUNCH alk*50*1000
60 PUNCH EQUI_DELTA("Calcite")*100*1000/2.71
100 PUNCH EQUI_DELTA("Halite")*58.4*1000/2.16
110 PUNCH EQUI_DELTA("Anhydrite")*136.14*1000/2.96
160 PUNCH EQUI_DELTA("Gypsum")*172.2*1000/2.96
######################
202 PUNCH KIN_DELTA("Barite")*233.39*1000/4.5
206 PUNCH KIN_DELTA("Quartz")*60*1000/2.65
210 PUNCH KIN_DELTA("Dolomite")*184.4*1000/2.84
550 PUNCH SR("Barite")
650 PUNCH SR("Quartz")
750 PUNCH SR("Dolomite")
#################################################
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
2
USE KINETICS 3
USE EQUILIBRIUM_PHASES 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
10
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
20
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
40
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
70
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
100
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
130
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
160
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
190
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
USE SOLUTION 3
USE REACTION_TEMPERATURE 3
REACTION_PRESSURE
220
USE EQUILIBRIUM_PHASES 3
USE KINETICS 3
END
--- End code ---
dlparkhurst:
You can use PRINT statements in RATES definitions to try to debug rate issues. The PRINT statements will be encountered many times because RATES is evaluated multiple times per time step, but it allows you to evaluate how the rate is calculated.
Hint: Your surface area is zero.
AntonioG:
I got it, thanks for your reply dlparkhurst (sorry for the disturbance).
Navigation
[0] Message Index
Go to full version