Issues with Simulating Dolomite Kinetics

[davidblevy]

Hello: I am trying to simulate the kinetic dissolution of quartz and dolomite under elevated PCO2(g) and total P. The file was running fine with quartz dissolution kinetics but not after adding the dolomite kinetics. I have tried decreasing time step and increasing bad_step_max but with no success. I don't have as much experience with modeling kinetics, so hopefully I am not missing something basic here. Any information would be appreciated, thank you. David

#Kinetic Modeling for RGN.

RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
 -start
  1  A0 = parm(1)
  2  V =  parm(2)
  3  dif_temp = 1/TK - 1/298
  4  pk_w = 13.7 + 4700.4 * dif_temp
  10 rate = (A0 / V) * (m/m0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
  20 save rate * time
 -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
Quartz
 -formula SiO2
 -m0  138                       
 -parms 23.13  0.16
 
Dolomite
 -formula  CaMg(CO3)2 
 -m0       3.97
 -parms    33.192
   
-steps   5 year in 50

INCREMENTAL_REACTIONS true

USER_GRAPH
  -chart_title "Quartz & Dolomite Dissolution"
 -axis_titles "Years " "[Si] mmol / L" "[Ca] mmol/L"
 -axis_scale x_axis 0 5
 -axis_scale y_axis 0.20 0.25
 -axis_scale sy_axis 0 auto
 -initial_solutions  true

10 plot_xy total_time / 3.1536e7,tot("Si") * 1e3, y-axis = 1
20 plot_xy total_time / 3.1536e7, tot("Ca")*1e3, y-axis = 2

SOLUTION 1 RGN
units      mg/L
pressure   97.4
temp      55.1 #131.1 F
pe      -0.78
pH      7.2
Al      0.4
Ba      0.43
Br      9.0
Sr      13.5
Alkalinity   248 as HCO3
Cl      10700
S(6)      2680 charge
F      1.0
Ca      606
Fe(2)      78.5
Mg      110
K      130
Si      11.2
Na      6410

GAS_PHASE 1
-fixed_pressure
-pressure      97.4
CO2(g)      0.999955
N2(g)         0.00000133
O2(g)         0.00001

END

[dlparkhurst]

I would use Ntg(g) instead of N2(g), which would make dissolved nitrogen inert and avoid potential redox reactions.

The script ran when I used -cvode in the KINETICS definition, but I did not look at the results carefully.

[davidblevy]

Thank you very much, the file is running now. I just removed the N2(g) since I am treating it as being inert. I did not see much change in the dissolved Ca concentration vs. time, so I reduced the time to 5 seconds in 500 steps and modified the file to just plot Ca concentration. I think what this is showing is that dolomite dissolution to reach equilibrium occurs within less than 1 second. Again, I have not worked much with kinetics, and I wanted to ask if this seems reasonable to you? Thanks again, David.

Kinetic Modeling for RGN.

RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
 -start
  1  A0 = parm(1)
  2  V =  parm(2)
  3  dif_temp = 1/TK - 1/298
  4  pk_w = 13.7 + 4700.4 * dif_temp
  10 rate = (A0 / V) * (m/m0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
  20 save rate * time
 -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
Quartz
 -formula SiO2
 -m0  138                       
 -parms 23.13  0.16
 
Dolomite
 -formula  CaMg(CO3)2 
 -m0       3.97
 -parms    33.192
   
-steps   5 in 500
-cvode   true

INCREMENTAL_REACTIONS true

USER_GRAPH
  -chart_title "Quartz & Dolomite Dissolution"
 -axis_titles "Seconds " "[Ca] mmol / L"
 -axis_scale x_axis 0 auto
 -axis_scale y_axis 15 auto
 
 -initial_solutions  true

10 plot_xy total_time,tot("Ca") * 1e3, y-axis = 1

SOLUTION 1 RGN
units      mg/L
pressure   97.4
temp      55.1 #131.1 F
pe      -0.78
pH      7.2
Al      0.4
Ba      0.43
Br      9.0
Sr      13.5
Alkalinity   248 as HCO3
Cl      10700
S(6)      2680 charge
F      1.0
Ca      606
Fe(2)      78.5
Mg      110
K      130
Si      11.2
Na      6410

GAS_PHASE 1
-fixed_pressure
-pressure      97.4
CO2(g)      0.999955
O2(g)         0.00001

END

[dlparkhurst]

Let's first recap the simulation. A SOLUTION is defined at temperature 55.1 C and pressure 97.4 atm. A GAS_PHASE is defined at 97.4 atm and 25 C; however, in the initial gas-phase calculation, the total pressure is not used. The amount of CO2 in the gas phase is (PV = nRT) about 1 * 1 / (0.082 * 298.15) = 0.04 mol.

 The initial solution is supersaturated with dolomite with a saturation index of 1.16. When the gas phase is reacted with this solution, all of the gas dissolves at a pressure of 97 atm and 55 C. Effectively, 0.04 mol of CO2 is added to the solution, the pH decreases to about 4, and dolomite is undersaturated with a saturation index of about -3. At this point the kinetic reaction kicks in and dissolves dolomite. The rate is indeed calculated to be very fast, on the order of a fraction of a second to reach equilibrium.

Just a couple comments. Make sure that the GAS_PHASE has the composition you expect and is reacting as you expected in producing the low pH. As for the rate, assuming you defined the rates correctly (perhaps from Palandri and Kharaka) you probably do not want to change any of the rate constants. Note that PR_P("CO2(g)") is zero because the gas phase doe not exist. That leaves the surface area. You will have to decide whether 33 m^2/mol is reasonable (total of about 130 m^2 given 3.97 mol of dolomite). The rate will be directly proportional to the surface area.

[davidblevy]

Thank you again and I do have a few follow-up questions:

You mentioned total pressure is not used in the initial GAS_PHASE calculation and it appears the default of 1 atm is being used. I defined the pressure as 97.4 atm under GAS_PHASE, so how would I modify the file to have the total pressure used in the initial GAS_PHASE calculation?

And with a different approach - Rather than having all of the gas dissolve leaving no gas phase, I would like the solution to be equilibrated with CO2(g) at 97.4 atm. So, I removed the GAS_PHASE and added CO2(g) under EQUILIBRIUM_PHASES (log PCO2(g) = 1.99 and 1000 moles to provide excess CO2(g). In that case, it appears that dolomite precipitates rather than dissolves, and now the Ca concentrations decrease with time. This seems counter-intuitive at high pressure but I assume the excess of carbon going into solution keeps the solution saturated with respect to dolomite and it precipitates. If I am trying to simulate water-mineral interactions over time under these constant conditions of pressure and elevated CO2(g), is the EQUILIBRIUM_PHASES approach more appropriate than the GAS_PHASE approach?

(There was a question regarding output for moles of minerals reacted for a kinetic simulation- I overlooked KIN_DELTA in the manual, so that should solve the third question I had).

Thank you again, David

[dlparkhurst]

You defined the partial pressure of CO2 to be 0.999955 atm in GAS_PHASE. The number of moles of CO2 is then calculated from this number, volume 1 L (default), and 25 C (default). When the GAS_PHASE reacts with a solution, the pressure of the system will be set to 97.4 (unless overridden by REACTION_PRESSURE). You could define the partial pressure of CO2(g) to be 97.4 in GAS_PHASE. Then, it would use that pressure in PV =nRT to determine the number of moles of CO2(g) in the initial gas phase. The number of moles would change when the GAS_PHASE reacted with the solution as would the volume of gas. If CO2(g) is the only component in the gas phase, then its partial pressure would remain 97.4 unless the gas dissolved completely.

Dolomite should still dissolve. I think perhaps there is an END before your definition of EQUILIBRIUM_PHASES so that dolomite precipitates from the initial solution that is supersaturated with dolomite.

EQUILIBRIUM_PHASES is more straightforward in that it maintains a fixed partial pressure. In GAS_PHASE, the bubble may dissolve and disappear, or the partial pressures may vary because of reactions involving the other components of the GAS_PHASE. The ratio of gas to liquid might be hard to determine.

There are functions KIN, which gives the total number of moles of kinetic reaction (M) at any given point, and KIN_DELTA, which, when using INCREMENTAL_REACTIONS, will give the change in M over the last time step. If INCREMENTAL_REACTIONS is false, KIN_DELTA gives the change in moles from the beginning of the kinetic calculation.

[davidblevy]

Thank you again, this is very helpful. I went back to my initial equilibrium calculations using the approach of defining CO2(g) under EQUILIBRIUM_PHASES and this seems to be the best approach. Revisiting the kinetic calculations, things seems to be working fine (so far!). David