Modeling concentrated aqueous H3PO4 solution with varying temperature

[jungsdao]

Dear developer,
Hello, I'm trying to model concentrated aqueous H3PO4 solution but encountering some error which I don't know how to deal with.
Following is my input code

Code: [Select]

DATABASE c:\phreeqc\database\pitzer.dat

USER_PUNCH
-head Spec_Cond
PUNCH SC
SELECTED_OUTPUT
-ionic_strength true
-pH
-water
-totals water
-charge_balance true
-saturation_indices
-equilibrium_phases

SOLUTION 1
pH 7 charge
temp 25
units mol/kgw
#density 1.1 calculate
water 0.015   # kg
P(5) 57.826 as P

SOLUTION_MASTER_SPECIES
P        HPO4-2         2.0     P               30.9738
P(5)     HPO4-2         2.0     P               30.9738   
P(5)     H2PO4-         1.0     P               
P(5)     H3P2O7-        1.0     P         
P(5)     H3PO4          0.0     P
P(5)     H4P2O7         0.0     P
#F        F-             0.0     F               18.9984
#O(0)     O2             0.0     O


SOLUTION_SPECIES
HPO4-2 =  HPO4-2
-llnl_gamma 4.0000
log_k 0
-delta_H 0 kJ/mol # Calculated enthalpy of reaction HPO4-2

H2PO4- =  H2PO4-
-llnl_gamma 4.0000
log_k 0
-delta_H 0 kJ/mol # Calculated enthalpy of reaction HPO4-2

H3P2O7- = H3P2O7-
-llnl_gamma 4.0000
log_k 0
-delta_H 0 kJ/mol # Calculated enthalpy of reaction HPO4-2

H3PO4 = H3PO4
-llnl_gamma 4.0000
log_k 0
-delta_H 0 kJ/mol # Calculated enthalpy of reaction HPO4-2

#2.0000 H+ + 1.0000 HPO4--  =  H3PO4
H3PO4  =  1.0000 H+ + 1.0000 H2PO4-
        -llnl_gamma           3.0   
        log_k           -2.15
-delta_H -8.29 kJ/mol # Calculated enthalpy of reaction H3PO4
# Enthalpy of formation: -307.92 kcal/mol
        #-analytic 1.8380e+002 6.7320e-002 -3.7792e+003 -7.3463e+001 -5.9025e+001
#       -Range:  0-300

H2PO4-  =  1.0000 H+ + 1.0000 HPO4--
        -llnl_gamma           3.0   
        log_k           -7.21
-delta_H 41.13 kJ/mol # Calculated enthalpy of reaction H3PO4
# Enthalpy of formation: -307.92 kcal/mol
        #-analytic 1.8380e+002 6.7320e-002 -3.7792e+003 -7.3463e+001 -5.9025e+001
#       -Range:  0-300

1.0000 HPO4--  =  PO4--- +1.0000 H+
        -llnl_gamma           4.0   
        log_k           -12.3218
-delta_H 14.7068 kJ/mol # Calculated enthalpy of reaction PO4-3
# Enthalpy of formation: -305.3 kcal/mol
        -analytic -7.6170e+001 -3.3574e-002 1.3405e+002 2.9658e+001 2.1140e+000
#       -Range:  0-300

2.0000 H3PO4 = 1.0000 H4P2O7 + 1.0000 H2O
log_k           -2.80
-delta_H        21.57 kJ/mol

1.0000 H4P2O7 = 1.0000 H+ +1.000 H3P2O7-
log_k           -1.54
-delta_H        -7.9

2.0000 HPO4-- + 1.0000 H+  =  HP2O7--- +1.0000 H2O
        -llnl_gamma           4.0   
        log_k           +5.4498
-delta_H 23.3326 kJ/mol # Calculated enthalpy of reaction HP2O7-3
# Enthalpy of formation: -2274.99 kJ/mol
        -analytic 3.9159e+002 1.5438e-001 -8.7071e+003 -1.6283e+002 -1.3598e+002
#       -Range:  0-300


REACTION_TEMPERATURE 1
25 200 in 50 steps


 USER_GRAPH 1
 -headings "temperature"  "H3PO4" "H2PO4-" "H3P2O7-" "H4P2O7" "pH"
 -axis_titles "Temperature","Molality (m)" ,"pH"
 -chart_title  "Ions in 200degC H3PO4 solution"
      -axis_scale x_axis 0 auto 
      -axis_scale y_axis 1e-2 auto auto auto
      -axis_scale sy_axis -1 auto
 -start
 10 pH=-LA("H+")
 20 GRAPH_X TC
 30 GRAPH_Y MOL("H3PO4"), MOL("H2PO4-"), MOL("H3P2O7-") ,MOL("H4P2O7")
 40 GRAPH_SY pH
 -end
 -active                 true

When I run this, I get following error message.
How can this be circumvented?
Any comment on input file and solution to this problem would be thankful!
Many thanks in advance.

Code: [Select]

ERROR: Could not reduce equation to primary master species, H4P2O7.
ERROR: Could not reduce equation to secondary master species, H2PO4-.
ERROR: Non-master species in secondary reaction, H2PO4-.
ERROR: Could not reduce equation to secondary master species, H3P2O7-.
ERROR: Non-master species in secondary reaction, H3P2O7-.
ERROR: Could not reduce equation to secondary master species, H3PO4.
ERROR: Non-master species in secondary reaction, H3PO4.
ERROR: Every primary master species for a redox element
must also be a secondary master species.
HPO4-2 is the primary master species for element P.
Another entry in SOLUTION_MASTER_SPECIES is needed.
Define species HPO4-2 as a secondary master species for a valence state.
For example:
P(0) HPO4-2 alk gfw
ERROR: Calculations terminating due to input errors.


[dlparkhurst]

First, phosphorus is not included in the pitzer.dat database. The Pitzer approach requires sets of parameters for interactions among species to account for non-ideality (activity coefficients). You cannot simply add aqueous species without including additional parameters in the PITZER keyword data block. The Pitzer approach would add a limited number of aqueous species of phosphorus, and account for additional effects with interaction coefficients.

Ion-association models tend to use aqueous species plus relatively simple activity coefficients to account for non-ideality. You can look at phreeqc.dat, llnl.dat, minteq.v4.dat to see which phosphorus aqueous species have been included in each database. It is a question which if any of these aqueous models are consistent with data for mean-activity coefficients, mineral solubility, and other empirical data.

Finally, if you do add species to a database (which I am not sure you should do), you must add only one line for each redox state of the species in SOLUTION_MASTER_SPECIES. If you define a redox state multiple times, only the last will be used. You must then define reactions for all of the aqueous species of the element in SOLUTION_SPECIES in such a way that all of the aqueous species can be written in terms of the master species defined in SOLUTION_MASTER_SPECIES.

I don't know what you are trying to model, but I would start using phreeqc.dat to determine whether it is satisfactory or not, and why.

[jungsdao]

Dear Mr. Parkhurst

Thank you for your reply, and I fixed code accordingly.
I searched a literature and included Pitzer parameter for phosphate species.
What I want to model is the profile of molality concentration of ionic species at extremely high concentration (50mol%) with increasing temperature.
Following code works fine when initial molality of Phosphoric acid is 1 mol/kgw, but when it is increased, it crashes at higher temperature.
I'm not sure this problem is due to lack of applicability of Pitzer parameter or some numerical noise which can be fixed.
Without inclusion of Pitzer parameter, the model gives prediction for even higher concentration with high temperature up to 200 degC.
Any comment would be appreciated!

Best regards,
Hyunwook

Code: [Select]

DATABASE c:\phreeqc\database\pitzer.dat

USER_PUNCH
-head Spec_Cond
PUNCH SC
SELECTED_OUTPUT
-ionic_strength true
-pH
-water
-totals water
-charge_balance true
-saturation_indices
-equilibrium_phases

LLNL_AQUEOUS_MODEL_PARAMETERS
-temperatures
         0.0100   25.0000   60.0000  100.0000
       150.0000  200.0000  250.0000  300.0000
#debye huckel a (adh)
-dh_a
         0.4939    0.5114    0.5465    0.5995
         0.6855    0.7994    0.9593    1.2180
#debye huckel b (bdh)
-dh_b
         0.3253    0.3288    0.3346    0.3421
         0.3525    0.3639    0.3766    0.3925
-bdot
         0.0374    0.0410    0.0438    0.0460
         0.0470    0.0470    0.0340    0.0000
#cco2   (coefficients for the Drummond (1981) polynomial)
-co2_coefs
        -1.0312              0.0012806
          255.9                 0.4445
      -0.001606

PITZER
# M. Cherif et al. / Fluid Phase Equilibria 175 (2000) 197–212
-B0
 H2PO4- H+ 0.0227
-B1
 H2PO4- H+ -3.4000
-LAMDA
 H3PO4 H+  3.3600
-MU
 H3PO4 H3PO4 H3PO4 0.0135
 
 
SOLUTION 1
pH 7 charge
temp 25
units mol/kgw
#density 1.1 calculate
water 0.015   # kg
P 4#57.826 # as P

SOLUTION_MASTER_SPECIES
#P(5) PO4-3 2.0 P 30.9738
P     HPO4--         2.0     P                  30.9738
#P(5) PO4-3 2.0 P 30.9738
#P        HPO4-2         2.0     P               30.9738
#P PO4-3 2.0 P 30.9738
#F        F-             0.0     F               18.9984
#O(0)     O2             0.0     O


REACTION_TEMPERATURE 1
25 200 in 50 steps

SOLUTION_SPECIES
HPO4-- =  HPO4--
-llnl_gamma 4.0000
log_k 0
-delta_H 0 kJ/mol # Calculated enthalpy of reaction HPO4-2
# Enthalpy of formation: -308.815 kcal/mol
2.0000 H+ + 1.0000 HPO4--  =  H3PO4
        -llnl_gamma           3.0   
        log_k           +9.3751
-delta_H 3.74468 kJ/mol # Calculated enthalpy of reaction H3PO4
# Enthalpy of formation: -307.92 kcal/mol
        -analytic 1.8380e+002 6.7320e-002 -3.7792e+003 -7.3463e+001 -5.9025e+001
#       -Range:  0-300

1.0000 HPO4-- + 1.0000 H+  =  H2PO4-
        -llnl_gamma           4.0   
        log_k           +7.2054
-delta_H -4.20492 kJ/mol # Calculated enthalpy of reaction H2PO4-
# Enthalpy of formation: -309.82 kcal/mol
        -analytic 8.2149e+001 3.4077e-002 -1.0431e+003 -3.2970e+001 -1.6301e+001
#       -Range:  0-300

1.0000 HPO4--  =  PO4--- +1.0000 H+
        -llnl_gamma           4.0   
        log_k           -12.3218
-delta_H 14.7068 kJ/mol # Calculated enthalpy of reaction PO4-3
# Enthalpy of formation: -305.3 kcal/mol
        -analytic -7.6170e+001 -3.3574e-002 1.3405e+002 2.9658e+001 2.1140e+000
#       -Range:  0-300

4.0000 H+ + 2.0000 HPO4--  =  H4P2O7 +1.0000 H2O
        -llnl_gamma           3.0   
        log_k           +15.9263
-delta_H 29.7226 kJ/mol # Calculated enthalpy of reaction H4P2O7
# Enthalpy of formation: -2268.6 kJ/mol
        -analytic 6.9026e+002 2.4309e-001 -1.6165e+004 -2.7989e+002 -2.7475e+002
#       -Range:  0-200

3.0000 H+ + 2.0000 HPO4--  =  H3P2O7- +1.0000 H2O
        -llnl_gamma           4.0   
        log_k           +14.4165
-delta_H 21.8112 kJ/mol # Calculated enthalpy of reaction H3P2O7-
# Enthalpy of formation: -544.1 kcal/mol
        -analytic 2.3157e+002 1.0161e-001 -4.3723e+003 -9.4050e+001 -6.8295e+001
#       -Range:  0-300

2.0000 HPO4-- + 1.0000 H+  =  HP2O7--- +1.0000 H2O
        -llnl_gamma           4.0   
        log_k           +5.4498
-delta_H 23.3326 kJ/mol # Calculated enthalpy of reaction HP2O7-3
# Enthalpy of formation: -2274.99 kJ/mol
        -analytic 3.9159e+002 1.5438e-001 -8.7071e+003 -1.6283e+002 -1.3598e+002
#       -Range:  0-300

 USER_GRAPH 1
 -headings "temperature"  "H3PO4" "H2PO4-" "H3P2O7-" "H4P2O7" "pH"
 -axis_titles "Temperature","Molality (m)" ,"pH"
 -chart_title  "Ions in 200degC H3PO4 solution"
      -axis_scale x_axis 0 auto 
      -axis_scale y_axis 1e-2 auto auto auto
      -axis_scale sy_axis -1 auto
 -start
 10 pH=-LA("H+")
 20 GRAPH_X TC
 30 GRAPH_Y MOL("H3PO4"), MOL("H2PO4-"), MOL("H3P2O7-") ,MOL("H4P2O7")
 40 GRAPH_SY pH
 -end
 -active                 true


[dlparkhurst]

First, you should not use LLNL_AQUEOUS_MODEL_PARAMETERS with pitzer.dat. I actually don't know what will be calculated using the Pitzer approach versus the llnl.dat approach in terms of the resulting activity coefficients. I thought there would be an error message, but I guess I need to add one.

Second, you did add Pitzer parameters for some species, but then you added additional species with no Pitzer parameters. I expect that you are doubly accounting for interactions between the Pitzer species with their parameters and the additional species that you added. If you are going to use the Pitzer approach, you should not add the additional species, which defeats your stated purpose. The purpose of the Pitzer approach is not to calculate speciation, but to calculate experimentally measurable quantities like mean activity coefficients and mineral solubility. Note also that the Pitzer parameters that you added have no temperature dependence, so the results at 200 C would be uncertain.

Finally, I am skeptical that the speciation calculated by any model at 200 C and high ionic strength will be meaningful. If you want to look at all multiple P species, I think you need to use an ion association model, perhaps llnl.dat. Note that the llnl.dat database contains these species H3PO4, H2PO4-, H4P2O7, H3P2O7-, H2P2O7-2, HPO4-2, HP2O7-3, P2O7-4, PO4-3.  Some have temperature dependence of equilibrium constants and some do not.       

[jungsdao]

Thank you for your reply.

I made mistake of mixing Pitzer parameter and LLNL parameters.
I'd like to simulate speciation between aqueous ions in 150 degC with phosphoric acid with high concentration.
Although concentration of phosphoric acid is highly, (around 85wt%), I presume ionic strength of phosphoric acid solution is not so high when compared to other strong acid solution. (H3PO4 is weak acid)
Phreeqc result also indicates ionic strength is lower than 1 mol/kgw throughout the temperature range.
With this low ionic strength, Do you think speciation calculated with LLNL model can be trustworthy?
I'd like to convinced that this calculated speciation is reliable around 150 degC @ high concentration as stated above.

Best regards,
Hyunwook

Code: [Select]

DATABASE c:\phreeqc\database\llnl.dat

USER_PUNCH
-head Spec_Cond
PUNCH SC
SELECTED_OUTPUT
-ionic_strength true
-pH
-water
-totals water
-charge_balance true
-saturation_indices
-equilibrium_phases

 # PITZER
 # # M. Cherif et al. / Fluid Phase Equilibria 175 (2000) 197–212
 # -B0
  # H2PO4- H+ 0.0227
 # -B1
  # H2PO4- H+ -3.4000
 # -LAMDA
  # H3PO4 H+  3.3600
 # -MU
  # H3PO4 H3PO4 H3PO4 0.0135
 
 
SOLUTION 1
pH 7 charge
temp 25
units mol/kgw
#density 1.1 calculate
water 0.015   # kg
P 57.826 # as P

SOLUTION_MASTER_SPECIES
#P(5) PO4-3 2.0 P 30.9738
P     HPO4--         2.0     P                  30.9738
#P(5) PO4-3 2.0 P 30.9738
#P        HPO4-2         2.0     P               30.9738
#P PO4-3 2.0 P 30.9738
#F        F-             0.0     F               18.9984
#O(0)     O2             0.0     O

REACTION_TEMPERATURE 1
25 200 in 50 steps

SOLUTION_SPECIES
HPO4-- =  HPO4--
-llnl_gamma 4.0000
log_k 0
-delta_H 0 kJ/mol # Calculated enthalpy of reaction HPO4-2
# Enthalpy of formation: -308.815 kcal/mol
2.0000 H+ + 1.0000 HPO4--  =  H3PO4
        -llnl_gamma           3.0   
        log_k           +9.3751
-delta_H 3.74468 kJ/mol # Calculated enthalpy of reaction H3PO4
# Enthalpy of formation: -307.92 kcal/mol
        -analytic 1.8380e+002 6.7320e-002 -3.7792e+003 -7.3463e+001 -5.9025e+001
#       -Range:  0-300

1.0000 HPO4-- + 1.0000 H+  =  H2PO4-
        -llnl_gamma           4.0   
        log_k           +7.2054
-delta_H -4.20492 kJ/mol # Calculated enthalpy of reaction H2PO4-
# Enthalpy of formation: -309.82 kcal/mol
        -analytic 8.2149e+001 3.4077e-002 -1.0431e+003 -3.2970e+001 -1.6301e+001
#       -Range:  0-300

1.0000 HPO4--  =  PO4--- +1.0000 H+
        -llnl_gamma           4.0   
        log_k           -12.3218
-delta_H 14.7068 kJ/mol # Calculated enthalpy of reaction PO4-3
# Enthalpy of formation: -305.3 kcal/mol
        -analytic -7.6170e+001 -3.3574e-002 1.3405e+002 2.9658e+001 2.1140e+000
#       -Range:  0-300

4.0000 H+ + 2.0000 HPO4--  =  H4P2O7 +1.0000 H2O
        -llnl_gamma           3.0   
        log_k           +15.9263
-delta_H 29.7226 kJ/mol # Calculated enthalpy of reaction H4P2O7
# Enthalpy of formation: -2268.6 kJ/mol
        -analytic 6.9026e+002 2.4309e-001 -1.6165e+004 -2.7989e+002 -2.7475e+002
#       -Range:  0-200

3.0000 H+ + 2.0000 HPO4--  =  H3P2O7- +1.0000 H2O
        -llnl_gamma           4.0   
        log_k           +14.4165
-delta_H 21.8112 kJ/mol # Calculated enthalpy of reaction H3P2O7-
# Enthalpy of formation: -544.1 kcal/mol
        -analytic 2.3157e+002 1.0161e-001 -4.3723e+003 -9.4050e+001 -6.8295e+001
#       -Range:  0-300

2.0000 HPO4-- + 1.0000 H+  =  HP2O7--- +1.0000 H2O
        -llnl_gamma           4.0   
        log_k           +5.4498
-delta_H 23.3326 kJ/mol # Calculated enthalpy of reaction HP2O7-3
# Enthalpy of formation: -2274.99 kJ/mol
        -analytic 3.9159e+002 1.5438e-001 -8.7071e+003 -1.6283e+002 -1.3598e+002
#       -Range:  0-300

 USER_GRAPH 1
 -headings "temperature"  "H3PO4" "H2PO4-" "H3P2O7-" "H4P2O7" "pH"
 -axis_titles "Temperature","Molality (m)" ,"pH"
 -chart_title  "Ions in 200degC H3PO4 solution"
      -axis_scale x_axis 0 auto 
      -axis_scale y_axis 1e-2 auto auto auto
      -axis_scale sy_axis -1 auto
 -start
 #10 pH=-LA("H+")
 10 GRAPH_X TC
 20 GRAPH_Y MOL("H3PO4"), MOL("H2PO4-"), MOL("H3P2O7-") ,MOL("H4P2O7")
 #40 GRAPH_SY pH
 -end
 -active                 true
 
 USER_GRAPH 2
-axis_titles "Temperature", "Ionic strength"
#-headings Addition "Ca(OH)2"
-chart_title  "Ionic strength"
-axis_scale x_axis  0 auto
-axis_scale y_axis  auto auto 1
-start
10 GRAPH_X TC
20 GRAPH_Y MU
-end
-active                 true


[dlparkhurst]

Sorry, I'm not convinced.

First note that you are using some species defined in llnl.dat. You defined 7 phosphate species in SOLUTION_SPECIES, but there are 9 in the output.

Code: [Select]

P(5)          4.462e+01
   H4P2O7          2.065e+01   2.065e+01     1.315     1.315     0.000     (0) 
   H3PO4           3.109e+00   3.109e+00     0.493     0.493     0.000     (0) 
   H3P2O7-         8.981e-02   5.942e-02    -1.047    -1.226    -0.179     (0) 
   H2PO4-          2.739e-02   1.812e-02    -1.562    -1.742    -0.179     (0) 
   H2P2O7-2        1.613e-03   2.971e-04    -2.792    -3.527    -0.735     (0) 
   HP2O7-3         2.300e-06   5.027e-08    -5.638    -7.299    -1.660     (0) 
   HPO4-2          1.457e-08   2.684e-09    -7.836    -8.571    -0.735     (0) 
   P2O7-4          1.160e-14   1.282e-17   -13.936   -16.892    -2.956     (0) 
   PO4-3           4.965e-19   1.085e-20   -18.304   -19.965    -1.660     (0) 

Even though the ionic strength is not that large, the concentrations are quite high, and I doubt the ion-association model accounts well for the interactions. Also, several of your reactions include H2O, so the activity of water is important. The ion-association model has a very simplistic formula for the activity of water.

If you are doing spectroscopic work and are actually measuring concentrations of species, I can understand you interest in speciation, but otherwise, I'm not sure  why you are interested.

[jungsdao]

Defined reaction in input stream is 7 but output has 9 phosphorus species because I just duplicated reactions in LLNL database.
Therefore, reaction definitions were actually redundant.
I'm interested in the concentration of ionic species because I suspect certain ionic species might more reactive toward surface dissolution.
At High temperature and high concentration, concentration of certain species dominates (according to "not so convinced" modeling data) over other species and might these species might play a key role when the H3PO4 solution is used for surface reaction.
Actually, relatively high concentration of certain species around 150degC and high conc. (around 85wt%) have been confirmed with several experimental literature(ex: raman & NMR): I just want to show corresponding result with Phreeqc modeling.
Hope this might have explained why I'm trying to calculate speciation of H3PO4 solution :)

[MichaelZ20]

Why don't you want to try using sit.dat database?

[jungsdao]

Thank you for your comment. 
As far as I know, SIT parameter for phosphoric acid is not well-parameterized...
I've look into the SIT database given by phreeqc, but I found insufficient parameter to model H3PO4 aqueous solution.
Any reason to suggest SIT model??

best regards,
Hyunwook

[MichaelZ20]

The only (except pitzer.dat) that works at high ionic strength.