Geothermal Bore Scaling

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

#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


[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).