H2S(g) solubility

[Helge]

Hello,
I am not a beginner of PHREEQC but meet unexptected difficulties with the calculation of the solubility of H2S(g). I already produced nice graphs matching experimental data for the solubility of H2S(g) as a function of p,T, m(NaCl). However, to use PHREEQC with an external fit-code I need PHREEQC to re-calculate experimental data at the point of interest. My strategy is to adjust total p, partial p(H2S), temperature and if necessary m(NaCl) to the published data and then calculate the molal solubility of H2S(g). I should also add that I managed this already for various p(H2S), p_tot<=1atm and several electrolytes. But with higher pressure, even slightly, this doesn't work, which tells me that there must be a profound misunderstanding on my side of how non-ideal gas behaviour has to be managed with PHREEQC.

These are some test data taken from
Suleimenov, O. M.; Krupp R. E. (1994): Solubility of hydrogen sulfide iп pure water and in NaCl solutions, from 20 to 320C and at saturation pressures. Geochim. Cosmochim. Acta (58,11), 2433-2444.
p(H2S) in this example was calculated as the difference between total pressure (published) and the p(H2O) at this temperature.

ID   p(H2S) / atm   +/-   T / K   +/-   m(H2S) / mol kg-1   +/-   p(tot) / atm
[SUL/KRU1994]_01   2.8832   0.0003   296.7   0.2   0.4179   0.0008   2.9114
[SUL/KRU1994]_02   4.7385   0.0005   333.4   0.2   0.4179   0.0008   4.9346
[SUL/KRU1994]_03   6.3263   0.0006   374.0   0.2   0.4179   0.0008   7.3427
[SUL/KRU1994]_06   3.8522   0.0004   295.6   0.2   0.5765   0.0012   3.8786
[SUL/KRU1994]_09   5.4306   0.0005   308.8   0.2   0.5975   0.0012   5.4873
[SUL/KRU1994]_15   4.1705   0.0004   294.0   0.2   0.7250   0.0014   4.1944
[SUL/KRU1994]_18   4.8689   0.0005   295.5   0.2   0.9246   0.0018   4.8951
[SUL/KRU1994]_25   7.9652   0.0008   359.0   0.2   0.7018   0.0014   8.5468
[SUL/KRU1994]_27   5.3122   0.0005   296.0   0.2   1.0100   0.0020   5.3393
[SUL/KRU1994]_30   4.2092   0.0004   294.5   0.2   0.7438   0.0015   4.2339
[SUL/KRU1994]_34   2.5320   0.0003   313.6   0.2   0.3517   0.0007   2.6055
[SUL/KRU1994]_36   4.3988   0.0004   313.2   0.2   0.5911   0.0012   4.4708
[SUL/KRU1994]_39   8.8012   0.0009   357.8   0.2   0.7274   0.0015   9.3560
[SUL/KRU1994]_41   2.5412   0.0003   294.6   0.2   0.4879   0.0010   2.5660
[SUL/KRU1994]_43   2.1615   0.0002   297.4   0.2   0.4002   0.0008   2.1910
[SUL/KRU1994]_46   7.7696   0.0008   296.0   0.2   1.5079   0.0030   7.7967

My (vain) attempts to recalculate these data are marked demo_01 and demo_02 in the script attached below. The calculation is done with pitzer.dat.

With demo_01 I observe that the target value for p(H2S) is well reproduced by not the total pressure, which seems to be systematically too high. The calculated m(H2S) is in bad agreement with the publshed data, and I am sure that pitzer.dat can perform better than that. With demo_02 no H2S at all is transferred to the system.

I am sure I am missing a very basic point here but I am stuck at the moment. Any hint would be appreciated.

Cheers,
Helge

### demo_01

SOLUTION 001; pH 3 charge; temp 23.5; SAVE SOLUTION 001; END;
SOLUTION 002; pH 3 charge; temp 60.2; SAVE SOLUTION 002; END;
SOLUTION 003; pH 3 charge; temp 100.8; SAVE SOLUTION 003; END;
SOLUTION 006; pH 3 charge; temp 22.4; SAVE SOLUTION 006; END;
SOLUTION 009; pH 3 charge; temp 35.6; SAVE SOLUTION 009; END;
SOLUTION 015; pH 3 charge; temp 20.8; SAVE SOLUTION 015; END;
SOLUTION 018; pH 3 charge; temp 22.3; SAVE SOLUTION 018; END;
SOLUTION 025; pH 3 charge; temp 85.8; SAVE SOLUTION 025; END;
SOLUTION 027; pH 3 charge; temp 22.8; SAVE SOLUTION 027; END;
SOLUTION 030; pH 3 charge; temp 21.3; SAVE SOLUTION 030; END;
SOLUTION 034; pH 3 charge; temp 40.4; SAVE SOLUTION 034; END;
SOLUTION 036; pH 3 charge; temp 40; SAVE SOLUTION 036; END;
SOLUTION 039; pH 3 charge; temp 84.6; SAVE SOLUTION 039; END;
SOLUTION 041; pH 3 charge; temp 21.4; SAVE SOLUTION 041; END;
SOLUTION 043; pH 3 charge; temp 24.2; SAVE SOLUTION 043; END;
SOLUTION 046; pH 3 charge; temp 22.8; SAVE SOLUTION 046; END;

GAS_PHASE 1
   -fixed_volume
   -volume 1
   H2Sg(g) 0
   H2O(g) 0

SELECTED_OUTPUT
   -file demo_01_pitzer-dat.calc

USER_PRINT
   -start
   10 print "SI(H2S) = ", SI("H2Sg(g)")
   20 print "Pressure / atm = ", PRESSURE
   -end

USE SOLUTION 001; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.4599; REACTION_PRESSURE 2.91142; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.41793)/0.41793; END; #[SUL/KRU1994]_01
USE SOLUTION 002; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.6756; REACTION_PRESSURE 4.93462; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.41793)/0.41793; END; #[SUL/KRU1994]_02
USE SOLUTION 003; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.8012; REACTION_PRESSURE 7.34271; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.41793)/0.41793; END; #[SUL/KRU1994]_03
USE SOLUTION 006; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.5857; REACTION_PRESSURE 3.87861; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.57651)/0.57651; END; #[SUL/KRU1994]_06
USE SOLUTION 009; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.7348; REACTION_PRESSURE 5.48729; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.59748)/0.59748; END; #[SUL/KRU1994]_09
USE SOLUTION 015; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.6202; REACTION_PRESSURE 4.19442; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.72499)/0.72499; END; #[SUL/KRU1994]_15
USE SOLUTION 018; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.6874; REACTION_PRESSURE 4.89514; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.92461)/0.92461; END; #[SUL/KRU1994]_18
USE SOLUTION 025; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.9012; REACTION_PRESSURE 8.54675; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.70175)/0.70175; END; #[SUL/KRU1994]_25
USE SOLUTION 027; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.7253; REACTION_PRESSURE 5.33925; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-1.01004)/1.01004; END; #[SUL/KRU1994]_27
USE SOLUTION 030; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.6242; REACTION_PRESSURE 4.2339; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.74384)/0.74384; END; #[SUL/KRU1994]_30
USE SOLUTION 034; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.4035; REACTION_PRESSURE 2.60548; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.35173)/0.35173; END; #[SUL/KRU1994]_34
USE SOLUTION 036; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.6433; REACTION_PRESSURE 4.47076; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.59115)/0.59115; END; #[SUL/KRU1994]_36
USE SOLUTION 039; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.9445; REACTION_PRESSURE 9.35603; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.72739)/0.72739; END; #[SUL/KRU1994]_39
USE SOLUTION 041; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.405; REACTION_PRESSURE 2.566; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.48793)/0.48793; END; #[SUL/KRU1994]_41
USE SOLUTION 043; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.3348; REACTION_PRESSURE 2.19097; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-0.40024)/0.40024; END; #[SUL/KRU1994]_43
USE SOLUTION 046; USE GAS_PHASE 1; EQUILIBRIUM_PHASES; H2Sg(g) 0.8904; REACTION_PRESSURE 7.79669; USER_PUNCH; -headings T_C p_tot_atm p_H2S_atm S_molal error_rel; 10 PUNCH TC; 15 PUNCH PRESSURE; 20 PUNCH pr_p("H2Sg(g)"); 30 PUNCH tot("Sg"); 40 PUNCH 100*(tot("Sg")-1.50786)/1.50786; END; #[SUL/KRU1994]_46

END

### demo_02

SOLUTION 001; pressure 2.91142; temp 23.5; EQUILIBRIUM_PHASES H2Sg(g) 0.45988; END;
SOLUTION 002; pressure 4.93462; temp 60.2; EQUILIBRIUM_PHASES H2Sg(g) 0.67564; END;
SOLUTION 003; pressure 7.34271; temp 100.8; EQUILIBRIUM_PHASES H2Sg(g) 0.80115; END;
SOLUTION 006; pressure 3.87861; temp 22.4; EQUILIBRIUM_PHASES H2Sg(g) 0.58571; END;
SOLUTION 009; pressure 5.48729; temp 35.6; EQUILIBRIUM_PHASES H2Sg(g) 0.73485; END;
SOLUTION 015; pressure 4.19442; temp 20.8; EQUILIBRIUM_PHASES H2Sg(g) 0.62019; END;
SOLUTION 018; pressure 4.89514; temp 22.3; EQUILIBRIUM_PHASES H2Sg(g) 0.68743; END;
SOLUTION 025; pressure 8.54675; temp 85.8; EQUILIBRIUM_PHASES H2Sg(g) 0.9012; END;
SOLUTION 027; pressure 5.33925; temp 22.8; EQUILIBRIUM_PHASES H2Sg(g) 0.72527; END;
SOLUTION 030; pressure 4.2339; temp 21.3; EQUILIBRIUM_PHASES H2Sg(g) 0.6242; END;
SOLUTION 034; pressure 2.60548; temp 40.4; EQUILIBRIUM_PHASES H2Sg(g) 0.40346; END;
SOLUTION 036; pressure 4.47076; temp 40; EQUILIBRIUM_PHASES H2Sg(g) 0.64334; END;
SOLUTION 039; pressure 9.35603; temp 84.6; EQUILIBRIUM_PHASES H2Sg(g) 0.94454; END;
SOLUTION 041; pressure 2.566; temp 21.4; EQUILIBRIUM_PHASES H2Sg(g) 0.40503; END;
SOLUTION 043; pressure 2.19097; temp 24.2; EQUILIBRIUM_PHASES H2Sg(g) 0.33476; END;
SOLUTION 046; pressure 7.79669; temp 22.8; EQUILIBRIUM_PHASES H2Sg(g) 0.8904; END;

[dlparkhurst]

A few comments.

(1) Your initial solutions will be pH 7 after charge balance. Charge balance with an anion, like Cl, if you want pH 3 for the initial condition.

(2) In demo_01, you need a semicolon after REACTION_PRESSURE. The pressure needs to be listed on the next line.

(3) Similarly in demo_02, you need a semicolon after EQUILIBRIUM_PHASES.

I would probably do the calculation as follows:

Code: [Select]

SOLUTION 1
END
USE solution 1
GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    H2O(g)    0
    H2Sg(g)   0
REACTION_PRESSURE 1
    2.91142
REACTION_TEMPERATURE 1
    23.5
REACTION 1
    Sg         1
    10 moles
END

I will see if I can find some files that demonstrate the fit of H2S solubility with pitzer.dat, but it will have to be later.

[dlparkhurst]

Here are some files from Tony Appelo that demonstrate the fit to a couple of papers of experimental data. They are a little complicated in detail, but the comparisons to experiment are plotted by USER_GRAPH.

First, these two files for Jiang and others, 2020. A file H2S:

Code: [Select]

REACTION_TEMPERATURE 15
0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 70, 80, 90, 94
REACTION_TEMPERATURE 5
30, 60, 90, 120, 150
REACTION_TEMPERATURE 6
30, 60, 90, 120, 150, 180
REACTION_TEMPERATURE 7
30, 60, 90, 120, 150, 180, 210
REACTION_TEMPERATURE 8
30, 60, 90, 120, 150, 180, 210, 240
REACTION_TEMPERATURE 10
30, 60, 90, 120, 150, 180, 210, 240, 270, 300
REACTION_TEMPERATURE 101
30, 60, 90, 120, 150, 180, 210, 240, 270, 295

SOLUTION 2
pH 7 charge
USE solution none
END

SOLUTION 1; -water 1; pH 7 charge; Sg 9e3

SELECTED_OUTPUT 2; -reset false; -file H2S.pqi; -new_line false
USER_PUNCH 2
10    sc$ = chr$(59) : hs$ = ' ' + chr$(35) + ' '
20    data '  1    ', '  29.6 ', '  49.3 ', '  98.7 ', ' 148.0  ', ' 197.4 ', ' 246.7 ', ' 296.1 ', ' 345.4 ', ' 394.8 ', ' 444.1 ', ' 493.5 ', ' 542.8 ', ' 592.1 ', ' 641.5 ', ' 690.8 ', ' 740.2 ', ' 789.5 ', ' 888.2 ', ' 986.9 '
30    dim p$(20) : for i = 1 to 20 : read p$(i) : next i
40    data 15, 7, 8, 101, 10, 10, 10, 7, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5
50    dim t(20) : for i = 1 to 20 : read t(i) : next i
60    gp$ = 'GAS_PHASE 1' + sc$ + ' -fixed_pr' + sc$ + ' -vol 0.01' + sc$ + 'H2Sg(g) 0' + sc$ + 'H2O(g) 0' + sc$ + ' -pr '
62    gp2$ = 'GAS_PHASE 1' + sc$ + ' -fixed_pr' + sc$ + ' -vol 0.01' + sc$ + 'H2Sg(g) 0' + sc$ + 'H2O(g) 0' + sc$ + 'Ntg(g) 100' + sc$ + ' -pr '
70    for i = 1 to 20
80      punch eol$ + 'USE solution 1' + sc$ + ' USE reaction_temperature ' + str_f$(t(i), 1, 0)
90      if i < 5 then punch eol$ + gp$ + p$(i) else punch eol$ + gp$ + p$(i)
92      if i > 1 then punch eol$ + 'USER_GRAPH 1' + sc$ + ' -headings x' #+ p$(i)
100     punch eol$ + 'END'
110   next i
END

SELECTED_OUTPUT 2; -active false
SELECTED_OUTPUT 1; -reset false; -file H2S.in.pqi; -new_line false
USER_PUNCH 1
-headings dif '' number.....tc.....P_calc..m_obs..m_calc
10  if get(1) = 0 then put(1, 1)

20  data 0.205, 0.174, 0.147, 0.127, 0.111, 0.097, 0.086, 0.076, 0.068, 0.054, 0.043, 0.031, 0.021, 0.011, 0.008,\
 1.932, 1.4483, 1.0793, 0.8461, 0.6825, 0.5223, 0.306,\
 1.9613, 2.1741, 1.842, 1.4887, 1.2722, 1.0691, 0.819, 0.4494,\
 1.9701, 2.2091, 2.5341, 2.753, 2.7876, 2.5869, 2.3883, 2.0201, 1.4239, 0.676,\
 2.0242, 2.2403, 2.6089, 3.1645, 3.6322, 3.8084, 3.8235, 3.6576, 3.2501, 2.5257,\
 2.0638, 2.3276, 2.7037, 3.5511, 3.5511, 3.5511, 5.3284, 5.5584, 5.4757, 4.874,\
 2.1094, 2.4236, 2.7872, 3.6007, 4.6227, 5.4793, 6.3812, 7.0389, 7.2365, 7.2987,\
 2.1546, 2.4581, 2.8515, 3.6984, 4.9065, 5.8992, 7.8033,\
  2.187, 2.4922, 2.9082, 3.7183, 5.2602, 6.5191,\
  2.2314, 2.5637, 2.986, 3.7918, 5.5007, 6.9683,\
 2.3262, 2.6256, 3.2248, 4.2306, 5.7343, 7.355,\
 2.3789, 2.6693, 3.2946, 4.3223, 5.9335, 7.705,\
 2.4056, 2.7116, 3.3599, 4.4433, 6.1232, 8.0869,\
 2.5025, 2.7524, 3.4157, 4.5237, 6.3934, 8.4285,\
 2.5239, 2.7848, 3.4541, 4.5971, 6.5064,\
 2.554, 2.8484, 3.5112, 4.7119, 6.6309,\
 2.6049, 2.8835, 3.5788, 4.7674, 6.8413,\
 2.6155, 2.931, 3.6326, 4.8544, 7.0222,\
 2.7046, 3.0362, 3.7248, 5.003, 7.4328,\
 2.7956, 3.08, 3.8067, 5.1491, 7.6829

30  no = get(1)
40  dim mo(150) : restore 20 : for i = 1 to no : read mo(i) : next i : i = no
50  pc = pressure : s = tot("Sg") : dif = (s / mo(i) - 1) * 1 #: if pressure < 190 then dif = 0
60  punch eol$ + str_e$(dif, 20, 13), ' ', str_f$(i, 1, 0) + str_f$(tc, 5, 0) + str_f$(pc, 8, 2) + str_f$(mo(i), 8, 3) + str_f$(s, 8, 3) : put(get(2) + dif^2, 2)
70  put(no + 1, 1)
160 if no = 143 then punch " phi =", get(2)
-end

#PRINT; -reset false; -echo_input false
USER_GRAPH 1 H2S solubility in water, data from Jiang et al., 2020
 -plot_csv_file H2S.dat
 -axis_titles "Temperature / Celsius" "H2S / (mol/kgw)"
 -axis_scale x_axis 0 310
 -axis_scale sy_axis 0 9
 # -initial_solutions true
 -connect_simulations false
 10 graph_x tc # tot("Na")
 20 y = tot("Sg") : if pressure < 2 then y = y * 10 : graph_sy y : end
 40 data 0, 0, 0, 4, 5, 6, 7, 8, 9, 0, 0, 0, 13, 0, 0, 0, 0, 18, 0, 0, 0, 22
 50 if sim_no > no then no = sim_no : dim p(30) : restore 40 : for i = 1 to no : read p(i) : next i
 60 if no = p(no) then graph_sy y
 -end
INCLUDE$ H2S.pqi
END

USE solution 2
REACTION_TEMPERATURE 1; 25      25      50      100     150     200.2   250     300     25      50      100     150     200.2   250     300#    25      50      100     150     200.2   250     300
REACTION_PRESSURE 1;    9.86    197.24  197.34  197.24  197.24  197.14  197.14  197.24  345.17  345.17  345.17  345.17  345.07  345.17  345.07# 1e3     1e3     1e3     1e3     1e3     1e3     1e3
REACTION 1; H2Sg 1;     0.232   0.247   0.249   0.251   0.252   0.254   0.266   0.26    0.246   0.249   0.25    0.252   0.253   0.267   0.259#  0.246   0.249   0.25    0.252   0.253   0.267   0.259

# SELECTED_OUTPUT 1; -reset false; -file in.prn; -new_line false
USER_PUNCH 1
# -headings dif '' number.....tc.....P_calc..m_obs..m_calc
10  if get(1) = 0 then put(1, 1)
20  data  34.8, 34.8, 36, 39, 43, 49.2, 58.3, 76.6, 35, 35.8, 38.4, 42.5, 47.7, 55.2, 67.9
30  no = step_no
40  dim mo(150) : restore 20 : for i = 1 to no : read mo(i) : next i : i = i - 1
50  pc = pressure : v = (vm("H2Sg") * mol("H2Sg") + vm("(H2Sg)2") * mol("(H2Sg)2") + vm("HSg-") * mol("HSg-")) / tot("Sg") : dif = (v / mo(i) - 1) * 1
52  dif = dif * 2.5
60  punch eol$ + str_e$(dif, 20, 13), ' ', str_f$(i, 1, 0) + str_f$(tc, 5, 0) + str_f$(pc, 8, 2) + str_f$(mo(i), 8, 3) + str_f$(v, 8, 3) : put(get(2) + dif^2, 2)
70  put(no + 1, 1)

160 if no = 15 then punch " phi =", get(2)

USER_GRAPH 1; -active false
USER_GRAPH 2
-axis_titles "temperature / K" "Overall Vm H2S species / (cm3/mol)"
-connect_simulations false
10 v = (vm("H2Sg") * mol("H2Sg") + vm("(H2Sg)2") * mol("(H2Sg)2") + vm("HSg-") * mol("HSg-")) / tot("Sg")
20 plot_xy tk, v, line_width = 0
# 20 plot_xy tk, vm("CO2"), line_width = 0
30 data 34.8, 34.8, 36, 39, 43, 49.2, 58.3, 76.6, 35, 35.8, 38.4, 42.5, 47.7, 55.2, 67.9
40 dim v(30) : restore 30 : for i = 1 to step_no : read v(i) : next i
50 plot_xy tk, v(i - 1), color = Black, line_width = 0, symbol_size = 12
# 60 data 8.77, 4.31, -17.1,  0,  -0.51, -1e-2
# 60 data 6.4, 15, -27, 0, -0.5, 0, 0, 0
# 60 data 6.4, 2.2, -27, 0, -0.5, 0, 0, 0
# 60 data 1.30620,  -4.59330, -30.00025,   0.00000,  -0.50000, 0
# 60 data -0.92986 ,  3.37217,  53.99450,   0.00000,  -0.50000, 0
# 70 dim v1(30) : restore 60 : for i = 1 to 6 : read v1(i) : next i
# 72 pb_s = 2600. + pressure * 1.01325 : TK_s = tc + 45.15
# 80 v2 = 41.84 * (v1(1) * 0.1 + v1(2) * 100 / (pb_s)  + v1(3) / (Tk_s) + v1(4) * 1e4 / (pb_s) / (Tk_s)) - v1(5) * 1e5 * QBrn + v1(6) * (pressure - 1)
# 90 plot_xy tk+step_no, v2, color = Magenta, symbol_size = 10, line_width = 0
END
and the data file H2S.dat:

Code: [Select]

tc 0.1.MPa.(m*10) 3.MPa 5.MPa 10.MPa 15.MPa 20.MPa 25.MPa 45.MPa 70.MPa 100.MPa
0 2.05
5 1.74
10 1.47
15 1.27
20 1.11
25 0.97
30 0.86
35 0.76
40 0.68
50 0.54
60 0.43
70 0.31
80 0.21
90 0.11
94 0.08
30 1.932
60 1.4483
90 1.0793
120 0.8461
150 0.6825
180 0.5223
210 0.306
30 1.9613
60 2.1741
90 1.842
120 1.4887
150 1.2722
180 1.0691
210 0.819
240 0.4494
30 1.9701
60 2.2091
90 2.5341
120 2.753
150 2.7876
180 2.5869
210 2.3883
240 2.0201
270 1.4239
300 0.5554
30 2.0242
60 2.2403
90 2.6089
120 3.1645
150 3.6322
180 3.8084
210 3.8235
240 3.6576
270 3.2501
300 2.5257
30 2.0638
60 2.3276
90 2.7037
120 3.5511
150 4.2857
180 4.8725
210 5.3284
240 5.5584
270 5.4757
300 4.874
30 2.1094
60 2.4236
90 2.7872
120 3.6007
150 4.6227
180 5.4793
210 6.3812
240 7.0389
270 7.2365
300 7.2987
30 2.3262
60 2.6256
90 3.2248
120 4.2306
150 5.7343
180 7.355
30 2.554
60 2.8484
90 3.5112
120 4.7119
150 6.6309
30 2.7956
60 3.08
90 3.8067
120 5.1491
150 7.6829

And here is a comparison to data from Xia and others, 2000. There is the PHREEQC script, H2S_NaCl_Na2SO4, followed by data in H2S_NaCl.dat and H2S_Na2SO4.dat.

Code: [Select]

SOLUTION 1; pH 7 charge; Sg 1.7e3

SELECTED_OUTPUT 1; -reset false; -file H2S_NaCl.pqi; -new_line false
USER_PUNCH 1
10 sc$ = chr$(59) : hs$ = ' ' + chr$(35) + ' '

20 data 3.38, 3.59, 8.02, 10.31, 12.45, 17.34, 20.84, 24.34, 4.24, 9.62, 13.02, 25.15, 25.99, 40.07, 5.76, 5.77, 16.15, 30.00, 44.66, 59.46, 6.07, 19.30, 19.38, 40.46, 65.71, 82.19, 95.73, 2.95, 4.75, 8.16, 9.27, 16.01, 17.39, 17.51, 23.29, 26.56, 4.96, 7.42, 15.39, 17.99, 29.47, 38.02, 3.29, 10.69, 17.86, 17.95, 28.86, 38.08, 58.55, 2.45, 10.79, 22.39, 22.41, 34.60, 44.69, 73.01, 94.85,            3.27, 10.00, 16.31, 16.97, 17.33, 19.54, 27.85, 4.52, 11.70, 12.41, 23.34, 30.44, 31.14, 35.82, 38.11, 8.42, 13.87, 32.61, 39.20, 52.99, 9.26, 17.83, 35.09, 52.76, 83.54, 5.97, 6.35, 6.20, 8.14, 10.74, 11.85, 12.38, 13.79, 17.42, 20.89, 24.93, 7.55, 9.77, 10.99, 14.45, 22.31, 25.30, 26.21, 36.74, 40.47, 5.62, 18.58, 30.56, 36.83, 41.03, 45.39, 55.43, 18.32, 21.68, 22.46, 45.02, 64.29, 65.86, 90.68, 96.56
 
30 no = 57 + 60 : dim p(150) : restore 20 : for i = 1 to no : read p(i) : next i

40  sol$ = eol$ + 'USE solution 1' + sc$
50  gp$ = ' GAS_PHASE 1' + sc$ + ' -fixed_pr' + sc$ + ' -vol 0.01' + sc$ + 'H2Sg(g) 0' + sc$ + 'H2O(g) 0' + sc$ + ' -pr '
60  gp0$ = ' GAS_PHASE 1' + sc$ + ' -fixed_vol' + sc$ + ' -vol 0.01' + sc$ + 'H2Sg(g) 0' + sc$ + 'H2O(g) 0'
70  r0$ = ' REACTION 1' + sc$ + ' H2Sg 1' + sc$ + ' -1.69 in 10' + sc$
80  r4$ = ' REACTION 1' + sc$ + ' NaCl 4.007' + sc$
90  r6$ = ' REACTION 1' + sc$ + ' NaCl 5.953' + sc$
100 t40$ = ' REACTION_TEMPERATURE 1' + sc$ + ' 40'   + sc$
110 t120$ = ' REACTION_TEMPERATURE 1' + sc$ + ' 120' + sc$
120 punch hs$ + ' 0 M NaCl...' + eol$ + sol$ + t40$ + r0$ + gp0$ + eol$ + 'END' + eol$
130 for i = 1 to 32
140   punch eol$ + hs$ + ' i = ' + str_f$(i, 1, 0)
150   if i > 17 then 200 REM 120.C # ==>
160     if i < 9 then punch sol$ + t40$ + r4$ + gp$ + str_f$(p(i), 1, 2) else punch sol$ + t40$ + r6$ + gp$ + str_f$(p(i + 19), 1, 2) # 8 and 9 samples
170     if i = 1 then punch eol$ + 'USER_GRAPH 1' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Green, symbol= None' + eol$
180     if i = 9 then punch eol$ + 'USER_GRAPH 1' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Blue, symbol= None' + eol$
190     GOTO 260 punch eol$ + ' END'  # ==>
200   REM 120.C # <==
210     if i = 18 then punch hs$ + ' 0 M NaCl...' + eol$ + sol$ + t120$ + r0$ + gp0$ + eol$
220     if i = 18 then punch eol$ + 'USER_GRAPH 1' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Orange, symbol= None' + eol$ + 'END' + eol$
230     if i < 25 then punch sol$ + t120$ + r4$ + gp$ + str_f$(p(i + 3), 1, 2) else punch sol$ + t120$ + r6$ + gp$ + str_f$(p(i + 25), 1, 2) # 7 and 8 samples
240     if i = 18 then punch eol$ + 'USER_GRAPH 1' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Magenta, symbol= None'    + eol$
250     if i = 25 then punch eol$ + 'USER_GRAPH 1' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Black , symbol= None' + eol$
260   punch eol$ + ' END' # <==
270 next i
END

SELECTED_OUTPUT 1; -active false
#PRINT; -reset false; -echo_input false
USER_GRAPH 1 H2S solubility in NaCl, data from Xia et al., 2000
 -plot_csv_file H2S_NaCl.dat
 -axis_titles "H2S / (mol / kgw)" "Pressure / atm"
 -axis_scale x_axis 0 1.4
 -axis_scale y_axis 0 100
 10 plot_xy tot("Sg"), pressure, symbol = None
INCLUDE$ H2S_NaCl.pqi
END

SOLUTION 1; pH 7 charge; Sg 2.2e3

USER_GRAPH 1; -active false

SELECTED_OUTPUT 2; -reset false; -file H2S_Na2SO4.pqi; -new_line false
USER_PUNCH 2
10 sc$ = chr$(59) : hs$ = ' ' + chr$(35) + ' '

20 data 3.27, 10.00, 16.31, 16.97, 17.33, 19.54, 27.85, 4.52, 11.70, 12.41, 23.34, 30.44, 31.14, 35.82, 38.11, 8.42, 13.87, 32.61, 39.20, 52.99, 9.26, 17.83, 35.09, 52.76, 83.54,        5.97, 6.35, 6.20, 8.14, 10.74, 11.85, 12.38, 13.79, 17.42, 20.89, 24.93, 7.55, 9.77, 10.99, 14.45, 22.31, 25.30, 26.21, 36.74, 40.47, 5.62, 18.58, 30.56, 36.83, 41.03, 45.39, 55.43, 8.32, 21.68, 22.46, 45.02, 64.29, 65.86, 90.68, 96.56
 
30 no = 25 + 35 : dim p(150) : restore 20 : for i = 1 to no : read p(i) : next i

40  sol$ = eol$ + 'USE solution 1' + sc$
50  gp$ = ' GAS_PHASE 1' + sc$ + ' -fixed_pr' + sc$ + ' -vol 0.01' + sc$ + 'H2Sg(g) 0' + sc$ + 'H2O(g) 0' + sc$ + ' -pr '
60  gp0$ = ' GAS_PHASE 1' + sc$ + ' -fixed_vol' + sc$ + ' -vol 0.01' + sc$ + 'H2Sg(g) 0' + sc$ + 'H2O(g) 0'
70  r0$ = ' REACTION 1' + sc$ + ' H2Sg 1' + sc$ + ' -2.19 in 10' + sc$
80  r4$ = ' REACTION 1' + sc$ + ' Na2SO4 0.492' + sc$
90  r6$ = ' REACTION 1' + sc$ + ' Na2SO4 0.95' + sc$
92  r62$ = ' REACTION 1' + sc$ + ' Na2SO4 0.953' + sc$
100 t40$ = ' REACTION_TEMPERATURE 1' + sc$ + ' 40'   + sc$
110 t120$ = ' REACTION_TEMPERATURE 1' + sc$ + ' 120' + sc$
120 punch hs$ + ' 0 M Na2SO4...' + eol$ + sol$ + t40$ + r0$ + gp0$ + eol$ + 'END' + eol$
130 for i = 1 to 31
140   punch eol$ + hs$ + ' i = ' + str_f$(i, 1, 0)
150   if i > 18 then 200 REM 120.C # ==>
160     if i < 8 then punch sol$ + t40$ + r4$ + gp$ + str_f$(p(i), 1, 2) else punch sol$ + t40$ + r6$ + gp$ + str_f$(p(i + 18), 1, 2) # 7 and 11 samples
170     if i = 1 then punch eol$ + 'USER_GRAPH 2' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Green, symbol= None' + eol$
180     if i = 8 then punch eol$ + 'USER_GRAPH 2' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Blue, symbol= None' + eol$
190     GOTO 260 punch eol$ + ' END'  # ==>
200   REM 120.C # <==
210     if i = 19 then punch hs$ + ' 0 M Na2SO4...' + eol$ + sol$ + t120$ + r0$ + gp0$ + eol$
220     if i = 19 then punch eol$ + 'USER_GRAPH 2' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Orange, symbol= None' + eol$ + 'END' + eol$
230     if i < 24 then punch sol$ + t120$ + r4$ + gp$ + str_f$(p(i + 2), 1, 2) else punch sol$ + t120$ + r6$ + gp$ + str_f$(p(i + 29), 1, 2) # 5 and 8 samples
240     if i = 19 then punch eol$ + 'USER_GRAPH 2' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Magenta, symbol= None'    + eol$
250     if i = 24 then punch eol$ + 'USER_GRAPH 2' + sc$ + ' 10 plot_xy tot("Sg"), pressure, color = Black , symbol= None' + eol$
260   punch eol$ + ' END' # <==
270 next i
END

SELECTED_OUTPUT 2; -active false
#PRINT; -reset false; -echo_input false
USER_GRAPH 2 H2S solubility in Na2SO4, data from Xia et al., 2000
 -plot_csv_file H2S_Na2SO4.dat
 -axis_titles "H2S / (mol / kgw)" "Pressure / atm"
 -axis_scale x_axis 0 2
 -axis_scale y_axis 0 100
 10 plot_xy tot("Sg"), pressure, symbol = None
INCLUDE$ H2S_Na2SO4.pqi
END


Code: [Select]

mH2S 40.C.==>0.0M ...==>4M ...==>6M 120.C.==>0.0M ...==>4M ...==>6M # M_NaCl
-1 -1
#
-1 -1
#
0.129 3.38
0.137 3.59
0.329 8.02
0.419 10.31
0.509 12.45
0.712 17.34
0.857 20.84
0.982 24.34
#
0.101 2.95
0.163 4.75
0.296 8.16
0.329 9.27
0.552 16.01
0.584 17.39
0.602 17.51
0.783 23.29
0.872 26.56
#
0.076 6.07
0.306 19.30
0.31 19.38
0.667 40.46
1.056 65.71
1.26 82.19
1.391 95.73
#
0.016 2.45
0.141 10.79
0.312 22.39
0.312 22.41
0.48 34.60
0.616 44.69
0.938 73.01
1.129 94.85


Code: [Select]

mH2S 40.C.==>0.0M ...==>0.49M ...==>0.95M 120.C.==>0.0M ...==>0.49M ...==>0.95M # M_Na2SO4
-1 -1
#
-1 -1
#
0.172 3.27
0.529 10.00
0.87 16.31
0.903 16.97
0.916 17.33
1.042 19.54
1.446 27.85
#
0.219 5.97
0.228 6.35
0.232 6.20
0.309 8.14
0.423 10.74
0.475 11.85
0.493 12.38
0.55 13.79
0.688 17.42
0.823 20.89
0.989 24.93
#
0.174 9.26
0.38 17.83
0.809 35.09
1.203 52.76
1.822 83.54
#
0.131 8.32
0.389 21.68
0.406 22.46
0.832 45.02
1.161 64.29
1.189 65.86
1.525 90.68
1.574 96.56


[Helge]

Thank you David. I am still busy swallowing the scripts.
If I understand correctly, the boundary conditions fix the total pressure and add an arbitrary mass of H2S (or Sg). The best agreement between set and calculated p(tot) and p(H2S) is with demo_02, which reads for the first line
SOLUTION 001; pressure 2.91142; temp 23.5; EQUILIBRIUM_PHASES; H2Sg(g) 0.45988; END;

I am not sure what the difference is between setting pressure in the solution block or setting reaction_pressure. I am also not sure if I have to add H2O(g) explicitly in the equlibrium_phases block or not. I'll be back next week.

[dlparkhurst]

Yes, these two calculations give slightly different results, but I think they should be close. I think calculations with the GAS_PHASE keyword are slightly better defined, and there are some assumptions about gasses for initial solution calculations that are more approximate.

Code: [Select]


SOLUTION 1
pH 3 charge;
pressure 2.91142
temp 23.5
Sg 1  H2Sg(g) 0.45988
END

SOLUTION 1
END
USE solution 1
GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    H2O(g)    0
    H2Sg(g)   0
REACTION_PRESSURE 1
    2.91142
REACTION_TEMPERATURE 1
    23.5
REACTION 1
    Sg         1
    10 moles
END