Problems calculating the saturation vapor pressure of aqueous solutions containing only water

[martinvoigt]

I ran into a problem where the solver is not able to find the correct or any solution for a simple case.
Essentially, it's about finding the water vapor pressure:

import reaktoro as rkt
db = rkt.PhreeqcDatabase('pitzer.dat')

solution = rkt.AqueousPhase(rkt.speciate("H O C"))
solution.setActivityModel(rkt.ActivityModelPitzerHMW())

gases = rkt.GaseousPhase("CO2(g) H2O(g)")
gases.setActivityModel(rkt.ActivityModelPengRobinson())

system = rkt.ChemicalSystem(db, solution, gases)

specs = rkt.EquilibriumSpecs(system)
specs.temperature()
specs.phaseAmount("GaseousPhase")

solver = rkt.EquilibriumSolver(specs)

state = rkt.ChemicalState(system)
state.set("H2O", 1.0, "kg")
#state.set("CO2", 0.05, "kg")

conditions = rkt.EquilibriumConditions(specs)
conditions.temperature(10, "celsius")
conditions.phaseAmount("GaseousPhase", 1, "umol")
conditions.setLowerBoundPressure(0.001, "bar")
conditions.setUpperBoundPressure(1000, "bar")

solver.solve(state, conditions)
print(state)

Without adding the CO2 (or adding very little), the pressure should be quite low at 10 °C, essentially the water vapor pressure, something like 0.01 bar. However, the solver doesn't converge and seems to go to the 100 bar upper bound.

I played around a bit with the bounds, as well as using other activity models, GasPhase amounts, etc, but didn't manage to get the right solution. It works fine with adding substantial amounts of CO2 (removing the # in the code above).

What's the problem here, and how do I solve this so that the code works robustly for both no/little amounts of CO2 as well as substantial amounts of CO2?

[allanleal]

Hi @martinvoigt , your code works after changing the upper bound pressure from 1000 bar to 2000 bar. It prints:

+-----------------+-------------+------+
| Property        |       Value | Unit |
+-----------------+-------------+------+
| Temperature     |    283.1500 |    K |
| Pressure        |   1302.1482 |  bar |
| Charge:         | -3.0000e-16 |  mol |
| Element Amount: |             |      |
| :: H            |  1.1101e+02 |  mol |
| :: C            |  4.0000e-16 |  mol |
| :: O            |  5.5506e+01 |  mol |
| Species Amount: |             |      |
| :: H+           |  8.0121e-08 |  mol |
| :: H2O          |  5.5506e+01 |  mol |
| :: CO3-2        |  1.0000e-16 |  mol |
| :: CO2          |  1.0000e-16 |  mol |
| :: HCO3-        |  1.0000e-16 |  mol |
| :: OH-          |  8.0121e-08 |  mol |
| :: CO2(g)       |  1.0000e-16 |  mol |
| :: H2O(g)       |  1.0000e-06 |  mol |
+-----------------+-------------+------+

However, if we compute the saturation pressure of water at 10 °C using the Wagner & Pruss (1995) EOS:

Psat = rkt.waterSaturationPressureWagnerPruss(10 + 273.15)
print(Psat)

we get 1228.11 bar. The difference is because Peng-Robinson is not as accurate for water vapor as Wagner-Pruss.

Something I tend to do when the calculation does not succeed is to print the log of the computation:

options = EquilibriumOptions()
options.optima.output.active = True
solver.setOptions(options)

You would most likely find out that the iterations got stuck on the upper bound 1000 bar, because pressure could not go beyond this value.

I see you are using ActivityModelPitzerHMW. Please switch to the new ActivityModelPitzer available in Reaktoro v2.3.0, which is a superior model and implementation. It is similar to the one in PHREEQC and it uses the same Pitzer interaction parameters available in the pitzer.dat file of PHREEQC v3.3.7.

[allanleal]

I just tried your code with the new Pitzer model. It works fine if the initial state is not pure water. Otherwise, a nan appears somewhere - most likely because of zero ionic strength in the denominator of some expression which I'll hunt down now! :)

[allanleal]

Thanks @allanleal for looking into this. I'm a bit confused - the water saturation/vapor pressure at 10 °C should be around 0.02 bar, not more than 1000 bar?!

Ops, 1228.11 is in Pa, not bar. I'm checking what is going on. I think it has to do with some recent performance improvements of the Wagner and Pruss EOS.

[allanleal]

If I set a small amount for CO2 as shown below:

state.set("CO2", 0.00005, "kg")

the calculation succeeds with pressure 0.033 bar. But decreasing it further, the calculated pressure jumps to ~1300 bar. I think I know the issue.

[martinvoigt]

Thanks - yeah that's similar to what I noticed - would be great to fix this :)

[allanleal]

I'm working on the fix for this.

[allanleal]

Hi Martin, I investigated this a bit more. You are using a PHREEQC database. Reaktoro will then resort to the same equation of state for water used in PHREEQC, which depends on first computing the density of water along the saturation pressure line and then interpolating to the desired pressure. I believe this is the reason why computing saturation pressure of pure water fails when PHREEQC databases are used. We would need to replace this simplified water EOS by the full Wagner & Pruss EOS.

The SUPCRT databases supported in Reaktoro will actually use Wagner & Pruss EOS. However, the calculation still fails. Now, for a different reason I believe. I recently improved the performance of this WP-EOS by creating a large table of densities for given T,P across a wide range. An approximate density can be quickly extracted from this table for a given (T, P) condition. This is then used as initial guess in the solution of the WP-EOS (i.e., the full EOS will refine the density and other properties). However, it's important to ensure that for the aqueous phase, initial guess for density reflects a liquid density. I think a gaseous density is being retrieved from this table and this, as an initial guess for the WP-EOS, ends up producing vapor-like properties instead of liquid-like ones.

I'll try to ensure that liquid densities are always produced for the aqueous solution (no matter the current value of P in the iterations) and then determine if this solves the issue (when using SUPCRT databases).

However, I think if water saturation properties are really needed, then it's best to just use the available methods that compute it rather quickly, e.g., waterSaturationPressureWagnerPruss.

In summary, I'd say that vapor pressure calculations for brines will work just fine with the setup above as long as the aqueous solution does not collapse into pure water.

[allanleal]

Hi Martin, I executed this PHREEQC script to compute the CO2 solubility in pure water at 50 °C and 50 atm:

SOLUTION
    units mol/kgw
    temp 50

GAS_PHASE
    -fixed_pressure
    -pressure 50

    CO2(g) 50
    H2O(g) 0

This produced 0.7636 moles of C(4) in the solution (which is the solubility of CO2 at 50 °C and 50 atm). This value is compatible with Reaktoro (using the same database) and experimental data.

Thus, we can say that 50 atm is the saturation pressure of an aqueous solution at 50 °C with 0.7636 moles of dissolved CO2. I then executed the following script to check if the calculated saturation pressure for this aqueous solution matches the expected 50 atm:

SOLUTION
    units mol/kgw
    temp 50

    C(4)    7.636e-01

GAS_PHASE
    -fixed_volume
    -volume 1e-8
    H2O(g) 0
    CO2(g) 0

I got a total pressure of 5.04 atm for the gas phase, which differs from the 50 atm.

If we use the following Reaktoro script to compute the saturation pressure for the solution:

import reaktoro as rkt

db = rkt.PhreeqcDatabase('phreeqc.dat')

solution = rkt.AqueousPhase(rkt.speciate("H O C"))
solution.setActivityModel(rkt.ActivityModelDavies())

gases = rkt.GaseousPhase("CO2(g) H2O(g)")
gases.setActivityModel(rkt.ActivityModelPengRobinson())

system = rkt.ChemicalSystem(db, solution, gases)

specs = rkt.EquilibriumSpecs(system)
specs.temperature()
specs.phaseAmount("GaseousPhase")

solver = rkt.EquilibriumSolver(specs)

state = rkt.ChemicalState(system)
state.set("H2O", 1.0, "kg")
state.set("CO2", 7.636e-01, "mol")

conditions = rkt.EquilibriumConditions(specs)
conditions.temperature(50, "celsius")
conditions.phaseAmount("GaseousPhase", 1, "umol")
conditions.setLowerBoundPressure(0.001, "bar")
conditions.setUpperBoundPressure(2000, "bar")

res = solver.solve(state, conditions)

assert res.succeeded()

print(state)

we get 49.7787 bar, which is close to the expected 50 atm.

I'm trying to understand if PHREEQC is actually calculating the saturation pressure for the aqueous solution or something else. I could be doing something completely wrong here, so please let me know if that's the case.

[martinvoigt]

PHREEQC can be a bit tricky in terms of how you have to express things.
For your second script, you would have to add the CO2 using REACTION. If you add it as C(4) in the SOLUTION, it add's only C / HCO3- or something similar, and things go wrong...

SOLUTION
    temp 50

REACTION
	CO2 7.636e-01

GAS_PHASE
    -fixed_volume
    -volume 1e-8
    H2O(g) 0
    CO2(g) 0

This gives something close to 50 bar.

[allanleal]

Indeed, your script produces 50 atm as expected!

OK, so clearly it should be possible to have Reaktoro producing the same saturation pressure when the aqueous solution collapses to pure water (while using PHREEQC databases). I'll need to check if I need to change a bit the Reaktoro implementation of the water EOS PHREEQC uses. There are some if-branches in the original implementation (and kept in Reaktoro) that affects automatic differentiation (and thus "wrong" Newton steps).

I think also that PHREEQC does not run into the same issue because it may iterate the solution in a different way. We'll see.

Thanks a lot for the interesting discussion.

[martinvoigt]

Thanks - keep me updated, it would be great to have this sort of capability when using PHREEQC databases, as for my applications these are sometimes the only suitable ones.

BTW Congrats on getting v2 out - that's a big milestone, looking forward to testing this more :)

[allanleal]

I'll keep you posted! Thanks once more.

[allanleal]

Just discovered that the issue lies on how Peng-Robinson EOS behaves. If you replace Peng-Robinson by Redlich-Kwong, you should get the desired saturation pressures when the solution is pure water:

gases.setActivityModel(rkt.ActivityModelRedlichKwong())

Executing the script below produces Psat = 0.0123 bar (WG-EOS produces 0.01228 bar) at 10 °C:

from reaktoro import *

db = PhreeqcDatabase('phreeqc.dat')

solution = AqueousPhase()
solution.setActivityModel(ActivityModelDavies())

gases = GaseousPhase("CO2(g) H2O(g)")
gases.set(ActivityModelRedlichKwong())

system = ChemicalSystem(db, solution, gases)

specs = EquilibriumSpecs(system)
specs.temperature()
specs.phaseAmount("GaseousPhase")

solver = EquilibriumSolver(specs)

state = ChemicalState(system)
state.set("H2O", 1.0, "kg")

conditions = EquilibriumConditions(specs)
conditions.temperature(10, "celsius")
conditions.phaseAmount("GaseousPhase", 1, "umol")
conditions.setLowerBoundPressure(0.001, "bar")
conditions.setUpperBoundPressure(2000, "bar")

res = solver.solve(state, conditions)

assert res.succeeded()

print(state)

At 50 °C, the script above produces 0.1215 bar (WG-EOS produces 0.12352 bar).

At least now I know that there is nothing fundamentally wrong.

What is happening is that Peng-Robinson is producing liquid-like properties for the gas phase because at these temperature and pressure conditions, the liquid state is more stable (perhaps by a tiny numerical margin under these conditions). So, I need to revisit this and ensure that thermodynamic properties are always gas-like.

[martinvoigt]

Aha, so the culprit is found :)

On a slightly related matter, does Reaktoro support the Akinfiev and Diamond (2003) thermodynamic model (e.g. for use of the aq17 database)?

[allanleal]

I'll ping @gdmiron here to comment on this! :)
If I'm not wrong, this is supported through the use of ThermoFun (developed by Dan Miron).

[gdmiron]

@martinvoigt if you use aq17 database, or your own version in thermofun format, when you define a reaktoro chemical system, the properties of some species are calculated using the Akinfiev-Diamond model and used in reaktoro to calculate the chemical equilibrium.

In thermofun format you can combine different methods for different species, for example hkf for some AD for others like in aq17.

You can of course just use aq17 or any other prepared database but if you are interested in making edits here are some examples on working with thermofun database, adding your own data, extending existing database, using reaction data, etc.

https://mybinder.org/v2/gh/thermohub/thermofun-jupyter/master?urlpath=lab/tree/gold2021workshop/

[martinvoigt]

This is great, many thanks @gdmiron !

[allanleal]

Hi @martinvoigt , just to let you know that in the next release (possibly tomorrow), we'll have a new ActivityModelPengRobinsonPhreeqc that is an equivalent implementation of PHREEQC's Peng-Robinson EOS (which contains binary interaction parameters for H2O-CO2, H2O-CH4, H2O-H2S, H2O-N2, H2O-C2H6, and H2O-C3H8).

ActivityModelPengRobinson is a plain Peng-Robinson model with zero binary interaction parameters. By replacing this with the new one, your H2O-CO2 gas/liquid/supercritical phase will be better handled thermodynamically. Stay tuned!

[allanleal]

This pull request here, resulting in Reaktoro v2.5.0, will give you new Peng-Robinson options:

• ActivityModelPengRobinsonPhreeqc
• ActivityModelPengRobinsonPhreeqcOriginal
• ActivityModelPengRobinsonSoreideWhitson

ActivityModelPengRobinsonPhreeqcOriginal is the exact implementation of PHREEQC.

By using ActivityModelPengRobinsonPhreeqc, I could reproduce the H2S solubility at high temperatures and pressures that PHREEQC produce:

which is comparable to this plot from https://www.hydrochemistry.eu/exmpls/other_gas.html: