Redox Kinetics of Aqueous Pu+3 oxidation to Pu+4

[mcomins]

Greetings PHREEQC+ Community,

I am trying to model oxidation of aqueous Pu+3 to Pu+4 in the presence of citrate. However, I can't get the software to "use" the rate expression to calculate the molality of the total of Pu3 or Pu4 species. I am wondering if part of the problem is the presence of aqueous Pu3-citrate and Pu4-citrate species, because I've only defined redox reactions with respect to Pu+3 and Pu+4 aquo ions. Please see my code below. (Note that I can provide additional information regarding aqueous Pu-citrate species, if it might be helpful.)

Code: [Select]

SOLUTION_SPECIES

1.000Pu+4     + 1.000e-     = Pu+3
     log_k    17.69   


RATES

Pu3_ox_to_Pu4_cit
-start
10 Pu3 = TOT("Pu3")
20 if (Pu3 <= 0) then goto 200
30 Pu4 = TOT("Pu4")
40 k1 = (1.7e-7)
50 rf = -k1*Pu3*1e-3
60 logmoles = (-k1 * TIME) + logM0
70 moles = EXP(logmoles)
200 SAVE moles
-end

SOLUTION_SPECIES

PHASES

CALCULATE_VALUES
    pHm
-start
  10 pHm = -LM('H+')
  20 SAVE pHm
-end
    log_Pu
-start
  10 log_Pu = log10(TOT"Pu")
  20 save log_Pu
-end
    log_Pu(III)
-start
  10 log_Pu(III) = log10(TOT"Pu(3)")
  20 save log_Pu(III)
-end
    log_Pu(IV)
-start
  10 log_Pu(IV) = log10(TOT"Pu(4)")
  20 save log_Pu(IV)
-end
   kinPu3
-start
  10 kinPu3 = KIN("Pu3_ox_to_Pu4_cit")
  20 SAVE kinPu3
-end
    d_kinPu3
-start
  10 d_kinPu3 = KIN_DELTA("Pu3_ox_to_Pu4_cit")
  20 SAVE d_kinPu3
-end

SELECTED_OUTPUT 1
    -file                 kinetics_Pu(III)to(IV)_cit.txt
    -high_precision       true
    -reset                false
    -ionic_strength       true
    -charge_balance       true
    -totals               Pu  Cit  Pu(3)  Pu(4)
    -molalities           Pu+3  Pu(Cit)-  Pu(CitH)  Pu(HCitH)+
                          PuH(CitH)2-2  Pu(HCitH)2-  Pu(CitH)2-3  Pu(OH)+2
                          Pu(OH)2+  Pu(OH)3  Pu2(OH)2(CitH)2-2  Pu3(OH)4(CitH)4-7
    -saturation_indices   Pu(OH)3(am)
    -kinetic_reactants    Pu3_ox_to_Pu4_cit
    -calculate_values     pHm  equi_Pu3solid  kinPu3  d_kinPu3

SOLUTION 1
    temp      25
    pH        2.5
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Cl        0.1
    Na        0.1
    Pu(3)     7.67521E-05 Pu(3)/Pu(4)
    Pu(4)     1e-015 Pu(3)/Pu(4)
    Cit       0.001
    -water    1 # kg

KINETICS 1
Pu3_ox_to_Pu4_cit
    -formula  Pu3  -1 Pu4  1
    -m        7.67521e-05
    -m0       7.67521e-05
    -tol      1e-08
-steps       15379200 in 20 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
USE SOLUTION 1
INCREMENTAL_REACTIONS true
Any guidance would be immensely appreciated.

[dlparkhurst]

Here is a file that oxidizes Pu(3) to Pu(4). It is done by adding oxygen to a Pu(3) solution.

Pu(5) and Pu(6) are suppressed, so that oxidation does not proceed past Pu(4). If you do not want to suppress these two redox states, remove the definitions in SOLUTION_SPECIES in the script.

I have hacked your kinetic reaction to be a simple exponential conversion. This may not be what you want, but you can redefine it in a way that you want. I have changed your calculate values to be consistent with the definitions in SELECTED_OUTPUT.

I did not have your Cit complexing definitions, so they have been ignored.

Code: [Select]

SOLUTION_SPECIES

1.000Pu+4     + 1.000e-     = Pu+3
     log_k    17.69   
# Suppress Pu(5)
1.000PuO2+2     + 0.500H2O     - 1.000H+     - 0.250O2     = PuO2+
     log_k    -20
# Suppress Pu(6)
1.000PuO2+     + 3.000H+     - 1.500H2O     - 0.250O2     = Pu+4
     log_k    20
RATES
Pu3_ox_to_Pu4_cit
-start
10 Pu3 = TOT("Pu(3)")
40 k1 = (1.7e-7)
60  moles = (k1 * Pu3 * TIME)
200 SAVE moles
-end

CALCULATE_VALUES
    pHm
-start
  10 pHm = -LM('H+')
  20 SAVE pHm
-end
    log_Pu
-start
  10 log_Pu = log10(TOT"Pu")
  20 save log_Pu
-end
    log_Pu(III)
-start
  10 log_Pu3 = log10(TOT"Pu(3)")
  20 save log_Pu3
-end
    log_Pu(IV)
-start
  10 log_Pu4 = log10(TOT"Pu(4)")
  20 save log_Pu4
-end
   kinPu3
-start
  10 kinPu3 = KIN("Pu3_ox_to_Pu4_cit")
  20 SAVE kinPu3
-end
    d_kinPu3
-start
  10 d_kinPu3 = KIN_DELTA("Pu3_ox_to_Pu4_cit")
  20 SAVE d_kinPu3
-end

SELECTED_OUTPUT 1
    -file                 kinetics_Pu(III)to(IV)_cit.txt
    -high_precision       true
    -reset                false
    -ionic_strength       true
    -charge_balance       true
    -totals               Pu  Cit  Pu(3)  Pu(4)
    -molalities           Pu+3  Pu(Cit)-  Pu(CitH)  Pu(HCitH)+
                          PuH(CitH)2-2  Pu(HCitH)2-  Pu(CitH)2-3  Pu(OH)+2
                          Pu(OH)2+  Pu(OH)3  Pu2(OH)2(CitH)2-2  Pu3(OH)4(CitH)4-7
    -saturation_indices   Pu(OH)3(am)
    -kinetic_reactants    Pu3_ox_to_Pu4_cit
    -calculate_values     pHm  log_Pu log_Pu(III) log_Pu(IV) kinPu3  d_kinPu3

SOLUTION 1
    temp      25
    pH        2.5
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Cl        0.1
    Na        0.1
    Pu(3)     7.67521E-05 Pu(3)/Pu(4)
    Pu(4)     1e-015 Pu(3)/Pu(4)
    Cit       0.001
    -water    1 # kg

KINETICS 1
Pu3_ox_to_Pu4_cit
    -formula  O2 0.25 #Pu3  -1 Pu4  1
    -m        7.67521e-05
    -m0       7.67521e-05
    -tol      1e-08
-steps       15379200 in 20 steps # seconds
-cvode
USE SOLUTION 1
INCREMENTAL_REACTIONS true
USER_GRAPH 1
    -headings               time Pu Pu(3) Pu(4)
    -axis_titles            "Time, seconds" "Molalit" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3.15e7
20 GRAPH_Y TOT("Pu"), TOT("Pu(3)"), TOT("Pu(4)")
  -end
    -active                 true


[mcomins]

This looks great! Many thanks for your insight.

[mcomins]

I have a follow-up question with a similar problem. I'm trying to model reduction of Pu(VI) to Pu(V) and Pu(V) to Pu(IV), so that two separate reduction reactions are proceeding simultaneously (the Pu(VI) reduction being slightly faster than Pu(V) reduction) (please see code below).

Code: [Select]

SOLUTION_SPECIES
 
1.000PuO2+     - 1.000e-     = PuO2+2
     log_k     -15.82     #18TAS/GAO

1.000Pu+4     + 2.000H2O     - 1.000e-     - 4.000H+     = PuO2+
     log_k     -17.45     #18TAS/GAO

RATES

Pu6_red_to_Pu5_cit
-start
10 Pu6 = TOT("Pu(6)")
20 k1 = (2.5e-7)
30 moles = (k1 * Pu6 * TIME)
200 SAVE moles
-end

Pu5_red_to_Pu4_cit
-start
10 Pu5 = TOT("Pu(5)")
20 k1 = (9.0e-9)
30 moles = (k1 * Pu5 * TIME)
200 SAVE moles
-end

CALCULATE_VALUES
    pHm
-start
  10 pHm = -LM('H+')
  20 SAVE pHm
-end
    log_Pu
-start
  10 log_Pu = log10(TOT"Pu")
  20 save log_Pu
-end
    log_Pu(VI)
-start
  10 log_Pu6 = log10(TOT"Pu(6)")
  20 save log_Pu6
-end
log_Pu(V)
-start
  10 log_Pu5 = log10(TOT"Pu(5)")
  20 save log_Pu5
-end
    log_Pu(IV)
-start
  10 log_Pu4 = log10(TOT"Pu(4)")
  20 save log_Pu4
-end
    kinPu6
-start
  10 kinPu6 = KIN("Pu6_red_to_Pu5_cit")
  20 SAVE kinPu6
-end
    d_kinPu6
-start
  10 d_kinPu6 = KIN_DELTA("Pu6_red_to_Pu5_cit")
  20 SAVE d_kinPu6
-end
    rate_Pu6_red
-start
  10 rate_Pu6_red = KIN_DELTA("Pu6_red_to_Pu5_cit")/KIN_TIME
  20 SAVE rate_Pu6_red
-end
    kinPu5
-start
  10 kinPu5 = KIN("Pu5_red_to_Pu4_cit")
  20 SAVE kinPu5
-end
    d_kinPu5
-start
  10 d_kinPu5 = KIN_DELTA("Pu5_red_to_Pu4_cit")
  20 SAVE d_kinPu5
-end
    rate_Pu5_red
-start
  10 rate_Pu5_red = KIN_DELTA("Pu5_red_to_Pu4_cit")/KIN_TIME
  20 SAVE rate_Pu5_red
-end

SELECTED_OUTPUT 1
    -file                 kinetics_Pu(VI) reduction_cit.txt
    -high_precision       true
    -reset                false
    -ionic_strength       true
    -charge_balance       true
    -totals               Pu  Cit  Pu(6)  Pu(5)  Pu(4)
    -saturation_indices   PuO2(OH)2:H2O(cr)  PuO2(OH)(am)  PuO2:2H2O(am)
    -kinetic_reactants    Pu6_red_to_Pu5_cit  Pu5_red_to_Pu4_cit
    -calculate_values     pHm  log_Pu  log_Pu(VI)  log_Pu(V)
                          log_Pu(IV)  kinPu6  d_kinPu6  rate_Pu6_red
                          kinPu5  d_kinPu5  rate_Pu5_red

SOLUTION 1
    temp      25
    pH        2.6
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Cit       0.001
    Cl        0.1
    Na        0.1
    Pu(4)     1.3e-06
    Pu(5)     3.4e-05 Pu(4)/Pu(5)
    Pu(6)     0.0001 Pu(5)/Pu(6)
    -water    1 # kg

KINETICS 1
Pu6_red_to_Pu5_cit
    -formula  O2  -0.25 Pu5  1
    -m        0.0001
    -m0       0.0001
    -tol      1e-08
Pu5_red_to_Pu4_cit
    -formula  O2  -0.25 Pu5  -1 Pu4  1
    -m        3.4e-05
    -m0       3.4e-05
    -tol      1e-08
-steps       5184000 in 20 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
USE SOLUTION 1

INCREMENTAL_REACTIONS true
USER_GRAPH 1
    -headings               time Pu Pu(6) Pu(5) Pu(4)
    -axis_titles            "Time, days" "Molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/86400
20 GRAPH_Y TOT("Pu"), TOT("Pu(6)"), TOT("Pu(5)"), TOT("Pu(4)")
  -end
    -active                 true

Two strange things are happening: The first is that total [Pu] is increasing as the reaction progresses. I can't find where in the script I've caused that to happen. The second is that all Pu is fully reduced to Pu(IV) by the first step, so the rates function isn't being utilized. Please let me know what I might be doing wrong. Thanks!

[dlparkhurst]

You are actually using equilibrium thermodynamics to control the Pu system. As you remove O2, the Pu species are all in equilibrium. You should not include Pu in the kinetics reaction when you are simulating reduction in this way. If you remove the Pu terms from -formula, Pu will be conserved.

One consequence of this approach is that Pu(3) will also be in equilibrium. The revised script shows that Pu(4) is never a significant redox state; Pu(3) is formed instead.

I'm fine with the equilibrium approach in most cases. If you really want to force the states to evolve kinetically, you need to split each Pu redox state into a new "element". You could define Pu_six, Pu_five, Pu_four, and Pu_three in SOLUTION_MASTER_SPECIES and add all of the aqueous reactions for each redox state in SOLUTION_SPECIES. See example 9 in the user's manual for an example with Fe(2) and Fe(3). Be forewarned that you will have to ensure that the kinetic transformations occur in a thermodyanamically reasonable way. The RATES definition will transfer Pu from one state to another regardless of the solution redox environment unless you do something like inhibit Pu(4) formation until after Pu(6) is consumed.

Code: [Select]

SOLUTION_SPECIES
 
1.000PuO2+     - 1.000e-     = PuO2+2
     log_k     -15.82     #18TAS/GAO

1.000Pu+4     + 2.000H2O     - 1.000e-     - 4.000H+     = PuO2+
     log_k     -17.45     #18TAS/GAO

RATES

Pu6_red_to_Pu5_cit
-start
10 Pu6 = TOT("Pu(6)")
20 k1 = (2.5e-7)
30 moles = (k1 * Pu6 * TIME)
200 SAVE moles
-end

Pu5_red_to_Pu4_cit
-start
10 Pu5 = TOT("Pu(5)")
20 k1 = (9.0e-9)
30 moles = (k1 * Pu5 * TIME)
200 SAVE moles
-end

SOLUTION 1
    temp      25
    pH        2.6
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Cit       0.001
    Cl        0.1
    Na        0.1
    Pu(4)     1.3e-06
    Pu(5)     3.4e-05 Pu(4)/Pu(5)
    Pu(6)     0.0001 Pu(5)/Pu(6)
    -water    1 # kg

KINETICS 1
Pu6_red_to_Pu5_cit
    -formula  O2  -0.25 #Pu5  1
    -m        0.0001
    -m0       0.0001
    -tol      1e-08
Pu5_red_to_Pu4_cit
    -formula  O2  -0.25 #Pu5  -1 Pu4  1
    -m        3.4e-05
    -m0       3.4e-05
    -tol      1e-08
-steps       10*51840. 10*519400 10*5194000 # seconds
-cvode true

USE SOLUTION 1

INCREMENTAL_REACTIONS true
USER_GRAPH 1
    -headings               time Pu Pu(6) Pu(5) Pu(4) Pu(3)
    -axis_titles            "Time, days" "Molality" ""
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/86400
20 GRAPH_Y TOT("Pu"), TOT("Pu(6)"), TOT("Pu(5)"), TOT("Pu(4)"), TOT("Pu(3)")
  -end
    -active                 true



[mcomins]

Thank you for your quick reply. I may have been making a mistake in my previous script because I hadn't seen any Pu(3) formation. I was only getting [Pu(4)] that was close to TOT("Pu") for each step.

I think it may be worthwhile to use kinetics instead of thermodynamics so that I can model Pu(4) and Pu(5) formation simultaneously, but I am encountering issues when trying to define these new Pu species. The SOLUTION_MASTER_SPECIES and related identity species below seem fine. It's when I try to define redox reactions with them that I get errors saying that each new Pu species is not balanced.

Code: [Select]

SOLUTION_MASTER_SPECIES

Pu_three      Pu_three+3     0            Pu_three      242
Pu_four       Pu_four+4      0            Pu_four       242
Pu_five       Pu_fiveO2+     0            Pu_five       242
Pu_six        Pu_sixO2+2     0            Pu_six        242

1.000Pu_three+3     = Pu_three+3
     log_k     0.0

1.000Pu_four+4     = Pu_four+4
     log_k     0.0

1.000Pu_fiveO2+     = Pu_fiveO2+
     log_k     0.0

1.000Pu_sixO2+2     = Pu_sixO2+2
     log_k     0.0

1.000Pu_three+3     - 1.000e-     = Pu_four+4
     log_k     -17.69     
     delta_h   -51.9     #kJ/mol 

1.000Pu_four+4      + 2.000H2O      - 4.000H+     - 1.000e-     = Pu_fiveO2+
     log_k     -17.45     

1.000Pu_fiveO2+     - 1.000e-     = Pu_sixO2+2
     log_k     -15.82

RATES

Pu6_red_to_Pu5_cit
-start
10 Pu_six = TOT("Pu_six")
20 k1 = (2.5e-7)
30 moles = (k1 * Pu_six * TIME)
200 SAVE moles
-end

Pu5_red_to_Pu4_cit
-start
10 Pu_five = TOT("Pu_five")
20 k1 = (9.0e-9)
30 moles = (k1 * Pu_five * TIME)
200 SAVE moles
-end

KINETICS 1
Pu6_red_to_Pu5_cit
    -formula  O2  -0.25
    -m        0.0001
    -m0       0.0001
    -tol      1e-08
Pu5_red_to_Pu4_cit
    -formula  O2  -0.25 Pu_four  1
    -m        1e-05
    -m0       1e-05
    -tol      1e-08
-steps       5184000 in 20 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5


[dlparkhurst]

You are mixing the two approaches. If you want to use KINETICS-only for the transformations between redox states, do not include any Pu reactions with e-. Remove your last three SOLUTION_SPECIES equations. Each redox state is completely independent from the other redox states except for the transformations you impose with KINETICS, for example:

-formula Pu_six -1 Pu_five +1

Adding O2 will have no effect unless you include it as part of your RATES equations, for example, a factor

rate = rate * MOL("O2")/(1e-5 + MOL("O2")),

which would inhibit the rate after O2 is consumed.

Or

rate = rate*(1e-5/(1e-5 + MOL("O2")),

which would inhibit a rate while O2 is present. You may want to do the same thing to turn off further reduction until, say, Pu_six is small.

[mcomins]

Ah, I see, that makes much more sense. I hope I can ask you one more troubleshooting question. When I try to run the kinetic version, it's giving me:

WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-13. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-14. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-15. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-16. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-17. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-18. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-19. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-20. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-21. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-22. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-23. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-24. Recovering...
WARNING: Negative moles in solution 1 for Pu_six, -1.000000e-25. Recovering...

These are my solution conditions:

Code: [Select]

SOLUTION 1
    temp      25
    pH        2.6
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Cit       0.001
    Cl        0.1
    Na        0.1
    Pu(4)     1.3e-06
    Pu(5)     3.4e-05
    Pu(6)     0.0001
    -water    1 # kg

[dlparkhurst]

If you do not get an "ERROR", you should be OK. The integration should recover from bad steps that lead to negative concentrations.

If the integration fails you can try -cvode if you have not already. Also you may want to add a factor like

TOT("Pu_six")/(1e-10 + TOT("Pu_six")) so that the concentrations don't get quite so small.

[mcomins]

Thank you for your feedback and for keeping this site running. I realized part of the issue came from my SOLUTION_MASTER_SPECIES where I used the primary master species "Pu" in both thermodynamic reactions (Pu(3), Pu(4), etc.) and kinetic reactions (Pu_three, Pu_four, etc.). I created new primary and secondary master species with Pv instead and got rid of that last error I was getting. Thank you again.