Deviation of hydrogen solubilities from experimental values at higher pressures

[kennedyantwi1]

Code: [Select]

Dear Dr. Parkhurst,
I'm trying to compare hydrogen solubilities in pure water between numerical and experimental models. At higher pressures (eg., > 300 bars), the numerical values tend to deviate from the experimental. Please, is there a way i may fix these deviations up to about 600 bars? The attempted code is attached, thanks:



[SOLUTION 1
    temp  30.00
    pH    7.0 charge
    units  mol/kgw
    -water 1.0 #1.0 kg
    -redox pe
    pe 4.0

GAS_PHASE 1
    -fixed_pressure
    -pressure 592.154
    -volume 10.0
    -temperature 30.00
    Hdg(g) 1.0

REACTION_PRESSURE
0.1 592.154 in 10

USER_GRAPH 1
    -headings Pressure Mol(Hdg) Gas_molar_volume
    -chart_title "Solubility of Hydrogen in Water as a Function of Pressure"
    -axis_titles "Pressure, atm", "Hdg(aq), mol/kgw" "Gas molar volume, L/mol"
    -plot_concentration_vs x
    -start
10 GRAPH_X PRESSURE
20 GRAPH_Y MOL("Hdg")
30 GRAPH_SY GAS_VM
    -end
ENDcode]

[dlparkhurst]

For the databases pitzer.dat and phreeqc.dat, the parameters that affect the solubility of Hdg(g) are defined in PHASES and SOLUTION_SPECIES. By defining T_c, P_c, and Omega (critical T, critical P, and eccentricity), the Peng-Robinson equation of state is used. The values affect the relation between T, P, and V for the gas. The molar volume of Hdg(aq) (-Vm definition) affects the pressure dependence of the solubility. See the manual, and Appelo, Parkhurst and Post, 2014. Geochim. Cosmochim. Acta 125, 49?67 for more details. You would have to decide if it is reasonable to adjust any of these parameters.

Code: [Select]

PHASES
Hdg(g)
Hdg = Hdg
-analytic -9.3114e+0 4.6473e-3 -4.9335e+1 1.4341e+0 1.2815e+5
-T_c 33.2; -P_c 12.8; -Omega -0.225
SOLUTION_SPECIES
Hdg = Hdg # H2
-Vm 6.52 0.78 0.12 # supcrt
-dw 5.13e-9
END
Alternatively, you can use the ideal gas law instead of Peng-Robinson. The ideal gas law will be used if T_c, P_c, and Omega are not defined. Vm will still affect the pressure dependence, but could be eliminated. The Henry's Law constant will determine the solubility as follows:

Kh = act(Hdg(aq)) / P(Hdg(g))

[kennedyantwi1]

Thanks very much, Dr Parkhurst, it worked perfectly well. I'm closest to the Experimental values!
My second polite request: Can you guide me on writing the code (in User Graph as well) to determine the change in Ionic Strength as the Hydrogen Solubility changes.
Very much appreciated in advance.

Code: [Select]

Dear Dr. Parkhurst,
I'm trying to compare hydrogen solubilities in pure water between numerical and experimental models. At higher pressures (eg., > 300 bars), the numerical values tend to deviate from the experimental. Please, is there a way i may fix these deviations up to about 600 bars? The attempted code is attached, thanks:



[SOLUTION 1
    temp  30.00
    pH    7.0 charge
    units  mol/kgw
    -water 1.0 #1.0 kg
    -redox pe
    pe 4.0

GAS_PHASE 1
    -fixed_pressure
    -pressure 592.154
    -volume 10.0
    -temperature 30.00
    Hdg(g) 1.0

REACTION_PRESSURE
0.1 592.154 in 10

USER_GRAPH 1
    -headings Pressure Mol(Hdg) Gas_molar_volume
    -chart_title "Solubility of Hydrogen in Water as a Function of Pressure"
    -axis_titles "Pressure, atm", "Hdg(aq), mol/kgw" "Gas molar volume, L/mol"
    -plot_concentration_vs x
    -start
10 GRAPH_X PRESSURE
20 GRAPH_Y MOL("Hdg")
30 GRAPH_SY GAS_VM
    -end
ENDcode]

[dlparkhurst]

The ionic strength is derived only from H+ and OH-, and so it is quite small for all pressures--on the order of 1e-7.

You can plot the ionic strength with the Basic function MU. Look at the section The Basic Interpreter in the manual for version 3 for all Basic functions.

[kennedyantwi1]

Dear Mr. Parkhurst, thanks for your guidance.

[kennedyantwi1]

Dear Mr. Parkhurst,
Please refer to my enclosed input file.
The calculated H2 solubility value depicted on the Plot is 2.319e-05 mol/kgw.
This value however does not appear in the Output file.
Please, can you offer any explanation and your guidance on how to locate this value in the Output file.
Thanks.

Code: [Select]

SOLUTION_SPECIES
Hdg = Hdg # H2
#-Vm 6.52 0.78 0.12 # supcrt
      -Vm  3 0.78 0.12
-dw 5.13e-9
END

PHASES
Hdg(g)
Hdg = Hdg
-analytic -9.3114e+0 4.6473e-3 -4.9335e+1 1.4341e+0 1.2815e+5
-T_c 33.2; -P_c 12.8; -Omega -0.225
END

SOLUTION 1
    temp 53.35
    pH    7.0 charge
    units mol/kgw
    pe 4.0
    Na 5.0
    Cl 5.0

GAS_PHASE 1
    -fixed_pressure
    -volume 20.0
    Hdg(g) 1
REACTION_PRESSURE
5.595854922
REACTION_TEMPERATURE
53.35

USER_GRAPH 1
    -headings               Pressure 326.50K #53.35C
    -axis_titles            "Pressure, bar" "Hdg(aq), mol/kgw" ""
    -chart_title            "Solubility of Hydrogen in NaCl as a Function of Pressure"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
    -plot_tsv_file HdgWater.txt
  -start
10 GRAPH_X PRESSURE*1.01325 #convert pressures to Bars
20 GRAPH_Y MOL("Hdg") / (TOT("water")*1000/GFW("H2O"))
  -end
END


[dlparkhurst]

Let's start with the units. I don't know if I made the mistake or you did, but the function MOL("Hdg") returns the concentration of Hdg (which is equivalent to H2(aq)) in units of mol/kgw. So, line 20 should be simply 20 GRAPH_Y MOL("Hdg"), if the units to be plotted are mol/kgw.

If the units should be mass, like mg/kgw or mg/L, then a conversion is needed.

If you still have questions, please post the file with the experimental values in addition to the script.

[kennedyantwi1]

Dear Mr. Parkhurst,
Thanks for your response.
Please enclosed a snipnet of the calculated vrs. measured H2 solubilities.
The issue:
the calculated values match very well the experimental values; they do not however appear in the output file!!

Code: [Select]

SOLUTION_SPECIES
Hdg = Hdg # H2
#-Vm 6.52 0.78 0.12 # supcrt
      -Vm  3 0.78 0.12
-dw 5.13e-9
END

PHASES
Hdg(g)
Hdg = Hdg
-analytic -9.3114e+0 4.6473e-3 -4.9335e+1 1.4341e+0 1.2815e+5
-T_c 33.2; -P_c 12.8; -Omega -0.225
END

SOLUTION 1
    temp 74.75
    pH    7.0 charge
    units mol/kgw
    pe 4.0
    Na 1.0
    Cl 1.0

GAS_PHASE 1
    -fixed_pressure
    -volume 20.0
    Hdg(g) 1
REACTION_PRESSURE
42.3459166  79.6180607  123.862818  168.497409  213.377745
REACTION_TEMPERATURE
74.75

USER_GRAPH 1
    -headings               Pressure 347.9K #74.75C
    -axis_titles            "Pressure [atm]" "Hdg(aq) [mol/kgw]" ""
    -chart_title            "Solubility of H2 gas in 1.0M NaCl solution as a Function of Pressure"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
    -plot_tsv_file HdgWater.txt
  -start
10 GRAPH_X PRESSURE*1.01325 #convert pressures to Bars
20 GRAPH_Y MOL("Hdg") / (TOT("water")*1000/GFW("H2O"))
  -end
END


[dlparkhurst]

Please read my posts. Your units are wrong in the USER_GRAPH section.

As for not showing up in the output file, I don't know what you are talking about!! The results are definitely in the output file and are plotted in the graph. I have limited your output to one pressure and removed the experimental data from the plots. Here is the output that shows a concentration of  2.457e-02 in the Solution composition and the Distribution of species.

Code: [Select]

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Cl                1.000e+00   1.000e+00
Hdg               2.457e-02   2.457e-02
Na                1.000e+00   1.000e+00

----------------------------Description of solution----------------------------

                                       pH  =   6.297      Charge balance
                                       pe  =  -0.946      Adjusted to redox equilibrium
      Specific Conductance (uS/cm,  75oC)  = 167782
                          Density (g/cm3)  =   1.01389
                               Volume (L)  =   1.04400
                        Viscosity (mPa s)  =   0.42618
                        Activity of water  =   0.966
                 Ionic strength (mol/kgw)  =   1.000e+00
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =  -2.068e-16
                         Temperature (oC)  =  74.75
                           Pressure (atm)  = 42.35
                  Electrical balance (eq)  =   2.079e-13
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =   0.00
                               Iterations  =   8
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550622e+01

----------------------------Distribution of species----------------------------

                                               Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma   cm3/mol

   OH-             7.244e-07   3.998e-07    -6.140    -6.398    -0.258     -2.02
   H+              6.958e-07   5.047e-07    -6.157    -6.297    -0.139      0.00
   H2O             5.551e+01   9.656e-01     1.744    -0.015     0.000     18.44
Cl            1.000e+00
   Cl-             1.000e+00   5.807e-01    -0.000    -0.236    -0.236     18.90
   HCl             2.860e-08   7.621e-08    -7.544    -7.118     0.426     (0) 
H(0)          1.400e-14
   H2              7.002e-15   8.815e-15   -14.155   -14.055     0.100     28.56
Hdg           2.457e-02
   Hdg             2.457e-02   3.094e-02    -1.610    -1.510     0.100     13.83
Na            1.000e+00
   Na+             1.000e+00   7.000e-01     0.000    -0.155    -0.155      1.05
O(0)          0.000e+00
   O2              0.000e+00   0.000e+00   -50.675   -50.575     0.100     32.86

When using the corrected USER_GRAPH, the graph plots the same concentration.

Code: [Select]

SOLUTION_SPECIES
Hdg = Hdg # H2
#-Vm 6.52 0.78 0.12 # supcrt
      -Vm  3 0.78 0.12
-dw 5.13e-9
END

PHASES
Hdg(g)
Hdg = Hdg
-analytic -9.3114e+0 4.6473e-3 -4.9335e+1 1.4341e+0 1.2815e+5
-T_c 33.2; -P_c 12.8; -Omega -0.225
END

SOLUTION 1
    temp 74.75
    pH    7.0 charge
    units mol/kgw
    pe 4.0
    Na 1.0
    Cl 1.0

GAS_PHASE 1
    -fixed_pressure
    -volume 20.0
    Hdg(g) 1
REACTION_PRESSURE
42.3459166  #79.6180607  123.862818  168.497409  213.377745
REACTION_TEMPERATURE
74.75

USER_GRAPH 1
    -headings               Pressure 347.9K #74.75C
    -axis_titles            "Pressure [atm]" "Hdg(aq) [mol/kgw]" ""
    -chart_title            "Solubility of H2 gas in 1.0M NaCl solution as a Function of Pressure"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
    #-plot_tsv_file HdgWater.txt
  -start
10 GRAPH_X PRESSURE*1.01325 #convert pressures to Bars
20 GRAPH_Y MOL("Hdg")
  -end
END

If you want to add more data to the output file, then you can use USER_PRINT as shown.

Code: [Select]

----------------------------------User print-----------------------------------

 Temperature    Pressure  H2(aq),mol/kgw
  7.4750e+01   4.2346e+01   2.4574e-02

If you are talking about putting data in the selected output file, then you need to use SELECTED_OUTPUT and possibly USER_PUNCH.

[kennedyantwi1]

Thanks, Mr. Parkhurst for your response.
Please below are the calculated H2 solubilities after running your enclosed/corrected input file.
My question is:
Are these calculated solubilities supposed to appear in the Output file?
Because, I'm NOT seeing them anywhere in the Output file.
Thanks.

         x                347.9K   
  4.2907e+01     4.4272e-04   
  8.0673e+01     8.2702e-04   
  1.2550e+02     1.2779e-03   
  1.7073e+02     1.7283e-03   
  2.1621e+02     2.1775e-03   

[dlparkhurst]

The output file contains the molalities of total Hdg and the aqueous species Hdg, which is H2(aq). The concentrations are moles per kilogram water. I gave excerpts from the output file for the first pressure; you should be getting the same output. Here is line 172 of the output file, is your output the same?

Code: [Select]

Hdg               2.457e-02   2.457e-02

What are the units of the two columns that you list? They don't seem to be mol/kgw, since I calculate 2.457e-02 mol/kgw for the first value from my output file (above) and the pressure(first column?) is not quite bars or atmospheres.

Code: [Select]

         x                347.9K   
  4.2907e+01     4.4272e-04   
  8.0673e+01     8.2702e-04   
  1.2550e+02     1.2779e-03   
  1.7073e+02     1.7283e-03   
  2.1621e+02     2.1775e-03 


[kennedyantwi1]

Good morning Dr. Parkhurst,
The x-values are Pressures in bars.
Please, if you may want to see the generatded graph after running the input file.
I suspect those points on the graph are the calculated solubilities in mol/kgw.
I compared them with the experimental solubilities from the literature and they matched.
Thanks.