changing moles of a reactant

[Cristina Domenech]

Dear
We are considering a system with several processes controlled by kinetic reactions. When the run starts, current moles of reactants change, in our case, diminishe. Thats ok. But at a given time we would like to increase the moles of one of the kinetic reactants while keeping the moles of the other at that time (simulating, for example an addition of that reactant in a lab experiment).
Can we change "M" value within a rate expression, for example    If total_time> 10000 M=0.5? We have tried it and it doesn't work.
We would appreciate your contributions
thanks

[dlparkhurst]

No, you cannot change M directly with M = x. You could SAVE 0.5 - M. RATES is called multiple times over a complete time step, so you would have to do this for all of the sub-time steps. I think this could work, but there may be complications with the sub time steps (the TOTAL_TIME is different at each sub-time step), or it is possible there would be some numerical issues.

Probably better would be to use KINETICS_MODIFY to change the number of moles of one reactant, before or after RunCells.

[Cristina Domenech]

Thank you very much. We will try it!

[Barbara]

Hi!
I think we are having a similar problem set. In our model scenario we also consider a system with several processes that are controlled by kinetic reactions.
The model involves the kinetic dissolution of a basaltic rock, thus the moles of the kinetic reactants of interest (e.g., S(-2), S(6), C(-4), C(4) ) in the solution are continuously changing over time. If we understood correctly from the above problem, we cannot directly change M in the KINETICS data block into a code that would allow PHREEQC to read values of M of the kinetic reactants as a function of time, right?
I just would like to double check: Would this mean that we need to select time steps (or cell runs) in our basalt dissolution model after which we “manually” change M for the kinetic reactants by using for example, KINETICS_MODIFY? Or is there another more “automated” way for this scenario to be set up in PHREEQC that we are unaware of?
Any insights would be helpful.
Thank you!

[dlparkhurst]

The moles of a KINETIC reactant (M) is adjusted automatically during kinetic integration. So, if the rate integrates to 1 mol over a time step; 1.0 moles of reactant has been added to the solution (actual moles of elements added to solution depends on the stoichiometry of -formula), but M is decreased by 1.0 moles.

So, in general, M does not need to be adjusted externally. If you are dissolving basalt, the moles of reactant will decrease as basalt dissolves.

The previous post had a less common situation to "increase the moles of one of the kinetic reactants while keeping the moles of the other at that time (simulating, for example an addition of that reactant in a lab experiment). " If you really need to change the number of moles externally, you would have to either define an entirely new KINETICS block or use KINETICS_MODIFY to adjust the number of moles of a reactant.

[Barbara]

Thank you very much for your quick response! This cleared things up.

[padhi]

Hello

I have a similar problem. I want to model CO2 dissolution kinetics, i.e., pure water is injected with CO2 at a particular pressure (10 atm for example).

Here is the script.


SOLUTION 1
    temp      25
    pressure 10
    pH        7 
 END

GAS_PHASE 1
    -fixed_pressure
    -pressure 10
    -volume   1
    -temperature 25.0
    CO2(g) 10
END

RATES
    CO2_kinetic
-start
 10   k = 3.0*(10^(-2))
 #20  conc_co2=TOT("CO2(g)")
 #20  if (conc_co2<= 0) then goto 200
 20  if (SI("CO2(g)")<= -8) then goto 200
 30  rateco2 = k*ACT("CO2")*ACT("H2O")
 40  moles = rateco2 * TIME
 200 SAVE moles
-end

KINETICS 1
CO2_kinetic
      -formula CO2 1.0
      -m0       0.1
      -tol      1e-08
-steps 25 50 100 200 300 400 500 600 700 
-step_divide 1e-4
-cvode true

INCREMENTAL_REACTIONS True

REACTION_TEMPERATURE 1
   25

USE SOLUTION 1
USE GAS_PHASE 1
USE KINETICS 1
USE REACTION_TEMPERATURE 1

USER_GRAPH 1 CO2 dissolution kinetics at 25C
    -headings               time TDIC CO2(available)  pH
    -axis_titles            "Time in hours" "CO2" "pH"
    -chart_title            "CO2 dissolution at 25C"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y ACT("CO2"), TOT("C(4)"), KIN("CO2_kinetic")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

USER_GRAPH 2 CO2 dissolution kinetics, CO2(g)
    -headings               time CO2(g) in atm
    -axis_titles            "Time in hours" "CO2(g)"
    -chart_title            "CO2 dissolution kinetics-CO2(g)"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y PR_P("CO2(g)")
  -end

END


My concerns are:

1. Is it possible to consider the amount of CO2 available to react kinetically to be defined in the gas phase (using the above code considers -m0 as the amount of CO2 available). Or else, is it that initial CO2 concentration needs to be explicitly calculated and defined in -m0?

2. In example 9 of PHREEQC manual for Fe oxidation, concentration of Fe_di was defined in the solution to be 1e-4 moles. By default, 1 mole of reactant was added during kinetic simulation (-m or -m0 was not defined). Why the initial concentration of Fe_di was 1e-4 moles, though 1 mole of reactant is being added in the kinetic simulation?
P.S. if I change -m0 drastically and simulate for longer time for this example, it does affect Fe_di concentration though.

Thank you

[dlparkhurst]

You may or may not use M directly in the RATES expression. Often reactions are proportional to M or there is a term like M/M0 to account for changing surface area. In your case, the rate depends explicitly on ACT("CO2") and ACT("H2O"). Although not explicitly used in the rate expression, M is important if it goes to zero. When M is zero, the reaction stops.

So, in the case of Fe_di, M is not explicitly used in the rates expression and M is a relatively large number, so it does not go to zero and has no real effect on the calculation.

You have found a bug in the Peng-Robinson calculation when using INCREMENTAL_REACTIONS. I'll have to fix this bug in the next release, but for now, do not use INCREMENTAL_REACTIONS true with Peng Robinson, at least don't use it when M goes to zero; it was correct up until that point.

I'm not sure what you are trying to model. You are adding CO2 to your system with your kinetic reaction, in addition to the amount in solution and the gas phase. The only effect is to increase the volume and moles of gas in the gas phase. The solution composition is constant. Because of your choice of M0 (really, the initial M), you only allowed the kinetic reaction to proceed until M decreased from 0.1 to zero. After that the kinetic reaction stopped and the system continued unchanged for the remainder of the time intervals.

Here is a version of your simulation without the error with INCREMENTAL_REACTIONS. The amount of gas in the gas phase is shown in the second plot.

Code: [Select]

SOLUTION 1
    temp      25
    pressure 10
    pH        7
 END

GAS_PHASE 1
    -fixed_pressure
    -pressure 10
    -volume   1
    -temperature 25.0
    CO2(g) 10
END

RATES
    CO2_kinetic
-start
 10   k = 3.0*(10^(-2))
 #20  conc_co2=TOT("CO2(g)")
 #20  if (conc_co2<= 0) then goto 200
 20  if (SI("CO2(g)")<= -8) then goto 200
 30  rateco2 = k*ACT("CO2")*ACT("H2O")
 40  moles = rateco2 * TIME
 200 SAVE moles
-end

KINETICS 1
CO2_kinetic
      -formula CO2 1.0
      -m0       0.1
      -tol      1e-08
-steps 25 in 25 # 25 50 100 200 300 400 500 600 700
-step_divide 1e-4
-cvode true

#INCREMENTAL_REACTIONS True

REACTION_TEMPERATURE 1
   25

USE SOLUTION 1
USE GAS_PHASE 1
USE KINETICS 1
USE REACTION_TEMPERATURE 1

USER_GRAPH 1 CO2 dissolution kinetics at 25C
    -headings               time ACT(CO2) TOT(CO2) KIN(CO2)  pH
    -axis_titles            "Time in hours" "CO2" "pH"
    -chart_title            "CO2 dissolution at 25C"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y ACT("CO2"), TOT("C(4)"), KIN("CO2_kinetic")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

USER_GRAPH 2 CO2 dissolution kinetics, CO2(g)
    -headings               time P(CO2(g)) GAS(CO2)
    -axis_titles            "Time in hours" "P[CO2(g)]" "CO2(g), Moles"
    -chart_title            "CO2 dissolution kinetics-CO2(g)"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y PR_P("CO2(g)")
30 GRAPH_SY GAS("CO2(g)")
  -end

END


[padhi]

Thank you very much Dr. Parkhurst.

I want to kinetically dissolve CO2, where the partial pressure or concentration of CO2 is defined. Having both gas phase and kinetics makes the solution to react with the gas phase too changing the pH of the solution according to CO2 in the gas phase.
 
I tried the following ways:

1. Defining small amount of C(4) in the solution, removing the gas phase and defining initial moles of CO2 as -m0 in the kinetics data block.

SOLUTION 1
    temp      25
    pressure 10
    pH        7
    C(4)      0.0001
 END

GAS_PHASE 1
    -fixed_pressure
    -pressure 10
    -volume   1
    -temperature 25.0
    CO2(g) 10
END

RATES
    CO2_kinetic
-start
 10   k = 3.0*(10^(-2))
 #20  conc_co2=TOT("CO2(g)")
 #20  if (conc_co2<= 0) then goto 200
 20  if (SI("CO2(g)")<= -8) then goto 200
 30  rateco2 = k*ACT("CO2")*ACT("H2O")
 40  moles = rateco2 * TIME
 200 SAVE moles
-end

KINETICS 1
CO2_kinetic
      -formula CO2 1.0
      -m0       0.1
      -tol      1e-08
-steps 1000 in 100 # 25 50 100 200 300 400 500 600 700
-step_divide 1e-4
-cvode true

#INCREMENTAL_REACTIONS True

REACTION_TEMPERATURE 1
   25

USE SOLUTION 1
#USE GAS_PHASE 1
USE KINETICS 1
USE REACTION_TEMPERATURE 1

USER_GRAPH 1 CO2 dissolution kinetics at 25C
    -headings               time ACT(CO2) TOT(CO2) KIN(CO2)  pH
    -axis_titles            "Time in hours" "CO2" "pH"
    -chart_title            "CO2 dissolution at 25C"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y ACT("CO2"), TOT("C(4)"), KIN("CO2_kinetic")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

USER_GRAPH 2 CO2 dissolution kinetics, CO2(g)
    -headings               time P(CO2(g)) GAS(CO2)
    -axis_titles            "Time in hours" "P[CO2(g)]" "CO2(g), Moles"
    -chart_title            "CO2 dissolution kinetics-CO2(g)"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y PR_P("CO2(g)")
30 GRAPH_SY GAS("CO2(g)")
  -end
END
 

With this option, CO2 concentration has to be explicitly defined (as -m0) or in the solution. As I want to quantify CO2 in solution too at high pressure in gas phase, without defining GAS PHASE, it seems not feasible. In the above code, GAS(CO2) is zero.


2. Defining another gas phase, which reacts kinetically to produce CO2. By defining a faster rate, I assume CO2 formation, which will dissolve kinetically in water.

Code:
SOLUTION_MASTER_SPECIES
    Cog       Cog         0     Cog         44.1
SOLUTION_SPECIES
Cog = Cog
   -log_k   0.0   
PHASES
Cog(g)
   Cog = Cog
   -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
END

SOLUTION 1
    temp      25
    pressure  10
    pH        7 
    units     mmol/kgw
    C(4)    0.0001 
    Cog      0.0001   
END

GAS_PHASE 1
    -fixed_pressure
    -pressure 10
    -volume   1
    -temperature 25.0
    Cog(g)   10
    CO2(g)    0
END

RATES
 CO2_formation
 -start
 20 rate = 1 #atm/sec
 30 moles = rate* PR_P("Cog(g)")
 40 SAVE moles
  -end   

 CO2_kinetic
 -start
 10   k = 3.0*(10^(-2))
 #20  conc_co2=TOT("CO2(g)")
 #20  if (conc_co2<= 0) then goto 200
 20  if (SI("CO2(g)")<= -8) then goto 200
 30  rateco2 = k*ACT("CO2")*ACT("H2O")
 40  moles = rateco2 * TIME
 200 SAVE moles
-end

KINETICS 1
  CO2_formation
      -formula CO2(g) 1 Cog(g) -1
      -m0      1
      -tol    1e-08

  CO2_kinetic
      -formula CO2 1.0
      -m0      1
     -tol      1e-08
-steps 100 in 100 steps
-step_divide 1e-4
-cvode true

#INCREMENTAL_REACTIONS True

REACTION_TEMPERATURE 1
   25
END

USE SOLUTION 1
USE KINETICS 1
USE GAS_PHASE 1
USE REACTION_TEMPERATURE 1

USER_GRAPH 1 CO2 dissolution kinetics at 25C
    -headings               time ACT(CO2) TOT(CO2) KIN(CO2) pH
    -axis_titles            "Time in hours" "CO2" "pH"
    -chart_title            "CO2 dissolution at 25C"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y ACT("CO2"), TOT("C(4)"), KIN("CO2_kinetic")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

USER_GRAPH 2 CO2 dissolution kinetics, CO2(g)
    -headings               time P(CO2(g)) GAS(CO2)
    -axis_titles            "Time in hours" "P[CO2(g)]" "CO2(g), Moles"
    -chart_title            "CO2 dissolution kinetics-CO2(g)"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y PR_P("CO2(g)")
30 GRAPH_SY GAS("CO2(g)")
  -end

END

With this code, I wonder with the following outputs:
1. when Cog(g) partial pressure is 10 atm, CO2 gas formed is larger than CO2 formed when Cog(g) partial pressure is 100 atm.
2. increasing kinetic reaction time makes the program run for very long time.
3. is there any other way to define CO2 concentration or partial pressure directly, that will be available for the kinetic reaction?

Thank you.

[dlparkhurst]

I still don't think you are posing your question correctly. Equilibrium between gas and water is pretty quick. If you include a fixed-pressure GAS_PHASE with sufficient gas, you are fixing the solution composition; pH, total C, activity of CO2(aq) are all fixed. There is little in the way of kinetics about it. Pure water equilibrium with CO2 gas at a given temperature and pressure gives you a unique solution composition.

Perhaps this script will help you understand. It adds CO2 to water by a REACTION, but you could make it KINETICS if you want. A gas phase is allowed to form when the partial pressure reaches 200 atm. pH and dissolved C increase as the reaction proceeds. Once the gas phase forms, the solution composition is constant. By choosing different pressures or temperatures you will have graphs that look similar, but will become constant at different points.

You could also consider the simulation to be reaction of a gas phase with successively larger volumes of CO2(g) at 200 atm. If the volume of gas is insufficient, it will dissolve completely, producing the increasing parts of the curves. With a sufficient amount of gas, the system produces gas-water equilibrium at a fixed gas pressure and solution composition.

Code: [Select]

SOLUTION
REACTION 1
CO2(g) 1
2 in 20 steps
GAS_PHASE
-pressure 200
CO2(g) 0

USER_GRAPH 1
    -headings               moles TOT(C) pH f(CO2(g))
    -axis_titles            "CO2 added, moles" "Total C and pH" "Fugacity of CO2(g)"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN
20 GRAPH_Y TOT("C")
30 GRAPH_Y -LA("H+")
40 GRAPH_SY 10^SI("CO2(g)")
  -end
    -active                 true
END


[padhi]

Dear Dr. Parkhurst,

Thank you very much for the explanation. It cleared my doubts.
I modified your previous code by adding kinetics.
SOLUTION
C(4)    0.0001 
#REACTION 1
#CO2(g) 1
#2 in 20 steps
GAS_PHASE
-pressure 200
CO2(g) 0

RATES
    CO2_kinetic
-start
 10   k = 3.0*(10^(-2))
 20  if (SI("CO2(g)")<= -8) then goto 200
 30  rateco2 = k*ACT("CO2")*ACT("H2O")
 40  moles = rateco2 * TIME
 200 SAVE moles
-end

KINETICS 1
CO2_kinetic
      -formula CO2 1.0
      -m0       2
      -tol      1e-08
-steps 1000 in 100 steps # 50 100 200 300 400 500 600 700 
#-steps 25 50 100 200 300 400 500 600 700 
#-steps 25 50 100 200 300 400 500 700 1000 1500 2000
-step_divide 1e-4
-cvode true
REACTION_TEMPERATURE 1
   25

USER_GRAPH 1
    -headings               moles TOT(C) pH f(CO2(g))
    -axis_titles            "time (hours)" "Total C and pH" "Fugacity of CO2(g)"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
#10 GRAPH_X RXN
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y TOT("C")
30 GRAPH_Y -LA("H+")
40 GRAPH_SY 10^SI("CO2(g)")
  -end
    -active                 true
END


Here, 2 moles of CO2 is added (-m0 n kinetics). If in the experiment, the volume of CO2 injected is v at a partial pressure of p, should I calculate the moles of CO2 and input as -m0?

Thank you so much.

[dlparkhurst]

I think you should try to model your experiment. If your experiment is to take a volume of water, a volume of gas at a specified pressure, and then bring the two together, you should simply define a GAS_PHASE with P and V and equilibrate your solution with the GAS_PHASE.

You seem to want to use KINETICS, so maybe your experiment is different.

[padhi]

Dear Dr. Parkhurst,

Thank you for your suggestion and your time.