Time dependent dissolution of CO2 gas phase

[smc]

Hi,
I am modelling a system of solid minerals, formation water, and CO2. I'd like to track the change in pH with time. Currently the dissolution of the solid minerals is set up through KINETICS, however I am not sure how to make CO2 dissolution time dependent. Currently i am using GAS_PHASE, which obviously equilibrates with the formation water at step one and thus changes it instantaneously. Is there a way to set up the dissolution of CO2 into water in a time dependent manner?

[dlparkhurst]

It's just another kinetic reaction. You'll have to decide whether kinetics or equilibrium is more appropriate.

Here's a kinetics example, where the target log partial pressure CO2 and rate constant are parameters in the KINETICS definition.

Code: [Select]

RATES
CO2_dissolution
20 k = parm(2) / 86400
30 eq_CO2 = 10^PARM(1)*10^LK_PHASE("CO2(g)")
40 act_CO2 = ACT("CO2")
50 moles = k * (eq_CO2 - act_CO2) * TIME
60 SAVE moles

END
SOLUTION
KINETICS
CO2_dissolution
-formula CO2 1
-parm 0 1 # target log PCO2(g), rate constant 1/day
-time 500000 in 10
USER_GRAPH 1
    -axis_titles            "Time, in days" "LOG P(CO2(g))" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/86400
20 GRAPH_Y SI("CO2(g)")
  -end


[MichaelZ20]

Hi David!
When I run your code, CO2(g) dissolves, but its pressure increases.
What may be a reason?

[dlparkhurst]

The original post wanted to be able to add CO2 to a solution kinetically. The script I wrote adds CO2, which will necessarily increase the partial pressure, but it needed some limit. The script is designed to add CO2 until the partial pressure reaches a target log partial pressure given by the first parameter in kinetics. So, in this case, the kinetic reaction stops when the partial pressure reaches 1 atmosphere. How fast it reaches 1 atm is determined by the second parameter--the rate constant, with units of 1/d.

[MichaelZ20]

David, thank you!
Actually, you set in KINETICS the target SI("CO2(g))").

[dlparkhurst]

The Henry's law relationship is

Code: [Select]

logK(CO2(g)) = a(CO2(aq)) / f(CO2(g))
or
a(CO2(aq)) = logK(CO2(g))*f(CO2(g))
and
f(CO2(g)) is approximately P(CO2(g))

LogK(CO2(g)) is given by the Basic function LK_PHASE("CO2(g)").

So, really, we are setting a target f(CO2(g)) of 10^0 (which is close to P(CO2(g)) at low pressure), the activity of CO2(aq) would be logk(CO2(g))*1. The RATES definition calculates the difference between the current activity and the target activity to determine a rate. If the activity is equal to the target activity, the rate goes to zero. When the activity is the target activity, the fugacity is the target fugacity.

[smc]

Thank you for the code. Just to clarify, this code is for the conversion of CO2(g) to CO2 (aq) correct, not the hydration of CO2 to H2CO3. If I wanted to represent that reaction kinetically, would I just set it up in the same way mineral kinetics are set up in the database?

[dlparkhurst]

Yes, the kinetic reaction dissolves CO2(g).

In phreeqc.dat, Amm.dat, wateq4f.dat, and others, CO2(aq) is actually H2CO3*, which is the sum of CO2(aq) and H2CO3(aq). My recollection is that CO2(aq) is predominant over H2CO3(aq). If you want to treat the species separately (using phreeqc.dat as an example), you would need to define H2CO3(aq) in SOUTION_SPECIES as CO2 + H2O -> H2CO3, with an appropriate log K. However, you would also need to adjust the log Ks for other carbonate reactions for HCO3- and, possibly CO3-2. Logically, I think you may also need to adjust the Henry's law constant for CO2(g).

[smc]

I am trying to implement reactions 1 and 2 from this paper: https://www.researchgate.net/publication/256178256_The_kinetics_of_the_reaction_CO2_H2O_H_HCO3-_as_one_of_the_rate_limiting_steps_for_the_dissolution_of_calcite_in_the_system_H2OCO2CaCO3

reaction 1:CO2 + H20 -> H2C03 -> H+ + HCO3- (acidic conditions)
reaction 2: CO2 + H20 -> CO2 + H+ + OH- -> H+ + HCO3- (basic conditions)
**these reactions dont need to both exist in the same simulation, i would just implement one or the other based on initial conditions
so i think i can bypass adding h2co3 explicitly. The paper above has the rate constants for the attached reactions, do you think there are any adjustments i'd need to make to implement the reactions kinetically? From my understanding, right now im dissolving CO2, and then speciating it according to equilibrium, and I'd like for the speciation of CO2 -> HCO3- to occur in accordance with the rates, and then HCO3- <-> CO3(2-) can happen according to equilibrium but perhaps I am misunderstanding the operation of the software
Thanks!

[dlparkhurst]

I think this may take a couple of exchanges to figure out, but here is the first try.

Here is reaction 1:

Code: [Select]

CO2 + H2O -> H2CO3 -> H+ + HCO3-

The study cites a forward and a back rate constant, but H2CO3, H+, and HCO3- are assumed to be at equilibrium. So, it seems to me, you must consider CO2 and H2CO3 as two separate species. So, there must be an equilibrium constant for CO2 + H2O -> H2CO3 that represents equilibrium, and it is implicit in the forward and backward rate constants.

As for the second reaction, the following reaction makes more sense to me. Again the forward and backward rate constants define the equilibrium constant.

Code: [Select]

CO2 + OH- -> HCO3-

Here is a script that you can consider. I've assumed the rate constants apply to the activities on the left- and right-hand-sides, which may not be correct. You could combine the two RATES definitions because the -formula definitions are the same. I'm not sure that the equilibrium constants implied by the rate constants are consistent with typical carbonate equilibrium constants, but that's on you. I had to redefine some of the carbonate system to account for separate CO2 and H2CO3 species.

Code: [Select]

SOLUTION_MASTER_SPECIES
    [CO2]         [CO2]            0     CO2             44
SOLUTION_SPECIES
[CO2] = [CO2]
log_k 0

CO3-2 + 2 H+ = H2CO3
-log_k 16.681
-delta_h -5.738 kcal
-analytic 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
2H2CO3 = (CO2)2 + 2H2O # 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
PHASES
CO2(g)
CO2 + H2O = H2CO3
-log_k -1.468
-delta_h -4.776 kcal
-analytic   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
RATES
CO2_hydration
-start
# CO2 + H2O -> H2CO3
10 kf = 3e-2 # 1/s
20 kb = 12   # 1/s
30 rf = kf*ACT("[CO2]")*ACT("H2O")^2
40 rb = kb*ACT("H2CO3")
50 rate = (rf - rb)
60 moles = rate * TIME
70 SAVE moles
-end
CO2_hydroxylation
-start
# CO2 + OH- = HCO3-
10 kf = 8.5e-3 # mol-1 s-1
20 kb = 2e-4  # 1/s
30 rf = kf*ACT("[CO2]")*ACT("OH-")
40 rb = kb*ACT("HCO3-")
50 rate = (rf - rb)
60 moles = rate * TIME
70 SAVE moles
-end
END
SOLUTION
-pH 9
[CO2] 1
KINETICS
CO2_hydration
-formula [CO2] -1 CO2 +1
CO2_hydroxylation
-formula [CO2] -1 CO2 +1
-step 2 in 10
USER_GRAPH 1
    -headings               time [CO2] H2CO3 HCO3- pH
    -axis_titles            "Seconds" "log molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y LOG10[TOT("[CO2]")], LM("H2CO3"), LM("HCO3-")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
END


[dlparkhurst]

I had some errors in the forward and backward rate expressions in the previous post.

Code: [Select]

SOLUTION_MASTER_SPECIES
    [CO2]         [CO2]            0     CO2             44
SOLUTION_SPECIES
[CO2] = [CO2]
log_k 0

CO3-2 + 2 H+ = H2CO3
-log_k 16.681
-delta_h -5.738 kcal
-analytic 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
2H2CO3 = (CO2)2 + 2H2O # 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
PHASES
CO2(g)
CO2 + H2O = H2CO3
-log_k -1.468
-delta_h -4.776 kcal
-analytic   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
RATES
CO2_hydration
-start
# CO2 + H2O -> H2CO3
10 kf = 3e-2 # 1/s
20 kb = 12   # 1/s
30 rf = kf*ACT("[CO2]")*ACT("H2O")
40 rb = kb*ACT("H2CO3")
50 rate = (rf - rb)
60 moles = rate * TIME
70 SAVE moles
-end
CO2_hydroxylation
-start
# CO2 + OH- = HCO3-
10 kf = 8.5e-3 # mol-1 s-1
20 kb = 2e-4  # 1/s
30 rf = kf*ACT("[CO2]")*ACT("H2O")
40 rb = kb*ACT("HCO3-")*ACT("H+")
50 rate = (rf - rb)
60 moles = rate * TIME
70 SAVE moles
-end
END
SOLUTION
-pH 9
[CO2] 1
KINETICS
CO2_hydration
-formula [CO2] -1 CO2 +1
CO2_hydroxylation
-formula [CO2] -1 CO2 +1
-step 2 in 10
USER_GRAPH 1
    -headings               time [CO2] H2CO3 HCO3- pH
    -axis_titles            "Seconds" "log molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y LOG10[TOT("[CO2]")], LM("H2CO3"), LM("HCO3-")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
END

[smc]

Hi,
Thanks for the code. When I set an equilibrium gas phase of CO2, it still automatically speciates and changes the pH of the solution. Essentially what I want is for CO2 to enter the aqueous phase in an equilibrium manner, but the speciation to occur based on the kinetic hydration reactions. i.e. we assume that the diffusion of CO2 into CO2(aq) is much faster than the hydration process. What would I need to change to get this result?
These reactions also seem to be having a lot of trouble with the step time, they seem to need a step time that is much smaller than what I can give. what would be a good way to troubleshoot this?

[dlparkhurst]

I think the script I gave you is close to what you want. If you remove the KINETICS block, then the solution contains only [CO2], which is unhydrated carbon dioxide. It is only by the kinetic reaction that [CO2] is transferred to an equilibrium carbon "C" pool. The kinetic reaction should stop when the forward and backward reactions are equal.

Considering further, I have been using phreeqc.dat, and I should have removed CO2(aq). As it stands, CO2(aq) is double accounted for with [CO2](aq) and CO2(aq). I also disabled the (CO2)2 species, which is only important at high pressures.

Code: [Select]

SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
log_k -100
2H2CO3 = (CO2)2 + 2H2O # activity correction for CO2 solubility at high P, T
log_k -100

However, to use this approach correctly, you need to change the log_k for the reaction CO3-2 + 2 H+ = H2CO3 because the log K in the database applies to H2CO3*, not H2CO3(aq).

You should be able to get the equilibrium constant from the forward and backward rate constants.

Code: [Select]

At equilibrium,
kf[CO2][H2O] = kb[H2CO3]

kf/kb = [H2CO3]/([CO2][H2O]) = Keq = 2.5e-3

CO2 + H2O = H2CO3
log_k -2.6

So then, assuming that H2CO3(aq) is negligible in H2CO3*, you can consider H2CO3* ~ CO2(aq), and you get

Code: [Select]

CO2 + H2O = H2CO3
log_k -2.6

CO3-2 + 2 H+ = H2CO3* ~ CO2(aq) (from phreeqc.dat)
-log_k 16.681

CO3-2 + 2H+ = H2CO3
log_k = 16.681 + (-2.6) = 14.081

Fine and dandy.

I have trouble with the second reaction. Using the same reasoning

Code: [Select]

CO2 + H2O = HCO3- + H+
kf[CO2][H2O] = kr[HCO3-][H+]

kf/kr = 8.5e-3/2e-4 = [HCO3-][H+]/([CO2][H2O]) = Keq = 42.5

CO2 + H2O = HCO3- + H+
log_k   1.63

Compare that to the accepted log K for that reaction:

Code: [Select]

CO3-2 + H+ = HCO3-
-log_k 10.329
CO3-2 + 2 H+ = CO2 + H2O
-log_k 16.681

CO2 + H2O = HCO3- + H+
log_k = 10.329 - 16.681 = -6.352

I don't see how the rate constants for the second reaction can be consistent with accepted equilibrium constants.

If you ignore the second reaction, this script should allow [CO2](aq) (unhydrated dissolved CO2) to react kinetically to H2CO3(aq). In this example, the partial pressure of [CO2](g) is fixed, so the activity of [CO2](aq) is fixed. I reduced some tolerances to smooth the curves. As for pH, you cannot expect the pH to remain constant. You can think of the pH as the ratio of HCO3-:H2CO3(aq) (in your system), and you are constantly decreasing this ratio as you hydrate CO2(aq) to H2CO3(aq), so pH decreases in my example script until equilibrium is reached.

Code: [Select]

KNOBS
-conv 1e-13
SOLUTION_MASTER_SPECIES
    [CO2]         [CO2]            0     CO2             44

SOLUTION_SPECIES
[CO2] = [CO2]
log_k 0

CO3-2 + 2 H+ = CO2 + H2O
log_k -100

2H2CO3 = (CO2)2 + 2H2O # activity correction for CO2 solubility at high P, T
log_k -100

CO3-2 + 2 H+ = H2CO3
-log_k      14.081
# -log_k 16.681
# -delta_h -5.738 kcal
# -analytic 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

PHASES
[CO2](g)
[CO2] = [CO2]
-log_k -1.468
-delta_h -4.776 kcal
-analytic   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
RATES
CO2_hydration
-start
# CO2 + H2O -> H2CO3
10 kf = 3e-2 # 1/s
20 kb = 12   # 1/s
30 rf = kf*ACT("[CO2]")*ACT("H2O")
40 rb = kb*ACT("H2CO3")
50 rate = (rf - rb)
60 moles = rate * TIME
70 SAVE moles
-end
CO2_hydroxylation
-start
# CO2 + H2O = HCO3- + H+
10 kf = 8.5e-3 # mol-1 s-1
20 kb = 2e-4  # 1/s
30 rf = kf*ACT("[CO2]")*ACT("H2O")
40 rb = kb*ACT("HCO3-")*ACT("H+")
50 rate = (rf - rb)
60 moles = rate * TIME
70 SAVE moles
-end
END
SOLUTION
-pH 9
Na 1 charge
#[CO2] 1
KINETICS
CO2_hydration
-formula [CO2] -1 CO2 +1
-tol 1e-12
#CO2_hydroxylation
#-formula [CO2] -1 CO2 +1
-step 40 in 10
EQUILIBRIUM_PHASES
[CO2](g) -3.4 10
USER_GRAPH 1
    -headings               time [CO2] H2CO3 HCO3- pH
    -axis_titles            "Seconds" "log molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y LOG10[TOT("[CO2]")], LM("H2CO3"), LM("HCO3-")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
END