Computing Gibbs molar energy for water, hydrogen, and oxygen --- why is the result so strange?

[B-LIE]

I'm considering electrolysis of water into hydrogen and oxygen. As part of the study, I need Gibbs molar energy for each of these substances. I have implemented an ideal gas model assuming constant heat capacities, which is given as:

Here, decoration ~ means "molar". Superscript o in energies/entropy is "standard state", and subscript o in energies/entropy is "at reference temperature" $T^\circ$.

Standard state enthalpy at reference temperature = enthalpy of formation. I find enthalpy of formation + standard state/reference temperature entropy + heat capacities in Aylward and Findley (2002): SI Chemical Data, 5th edition.

Problem: My ideal gas expression gives wildly different results than Clapeyron (using PR). As an example, for Hydrogen:

julia> G_o(50+273.15,Ht_oo_H2,St_oo_H2,ctp_H2), gibbs_energy(mod_H2_pr,p_o,50+273.15,1)
(-42362.23012868961, -13568.639305566825)

I'm relatively sure that my ideal gas model is ok.

Question: Why this extreme discrepancy?

---

Julia code:

using Clapeyron
# 
# Data from Aylward and Findley, etc.
R = 8.314_462_618_153_24        # J/K mol
#
p_o = 1e5                       # Pa, Standard state
T_o = 25 + 273.15               # K, Reference temperature
#
Ht_oo_H2 = 0                    # J/mol, molar enthalpy at standard state/reference temperature = enthalpy of formation
Ht_oo_O2 = 0                    # J/mol
Ht_oo_H2O_g = -242e3            # J/mol
Ht_oo_H2O_L = -286e3            # J/mol
St_oo_H2 = 131                  # J/K mol, molar entropy at standard state/reference temperature
St_oo_O2 = 205                  # J/K mol
St_oo_H2O_g = 189               # J/K mol
St_oo_H2O_L = 70                # J/K mol
#
ctp_H2 = 29                     # J/K mol, molar heat capacity
ctp_O2 = 29                     # J/K mol
ctp_H2O_g = 34                  # J/K mol
ctp_H2O_L = 75                  # J/K mol
# 
# Standard state ideal gas molar Gibbs energy:
G_o(T,Ht_oo,St_oo,ctp) = Ht_oo - T*St_oo + ctp*(T-T_o - T*log(T/T_o))
# 
# Clapeyron models 
mod_H2O_pr = PR(["water"])
mod_H2_pr = PR(["hydrogen"])
mod_O2_pr = PR(["oxygen"])

[B-LIE]

Hm... If I add ideal gas:

mod_H2O_ig = BasicIdeal(["water"])
mod_H2_ig = BasicIdeal(["hydrogen"])
mod_O2_ig = BasicIdeal(["oxygen"])

the result is:

julia> G_o(50+273.15,Ht_oo_H2,St_oo_H2,ctp_H2), gibbs_energy(mod_H2_ig,p_o,50+273.15,1), gibbs_energy(mod_H2_pr,p_o,50+273.15,1)
(-42362.23012868961, -13569.40800266611, -13568.639305566825)

julia> G_o(50+273.15,Ht_oo_O2,St_oo_O2,ctp_O2), gibbs_energy(mod_O2_ig,p_o,50+273.15,1), gibbs_energy(mod_O2_pr,p_o,50+273.15,1)
(-66275.33012868962, -13569.40800266611, -13571.171789817974)

julia> G_o(50+273.15,Ht_oo_H2O_g,St_oo_H2O_g,ctp_H2O_g), gibbs_energy(mod_H2O_ig,p_o,50+273.15,1), gibbs_energy(mod_H2O_pr,p_o,50+273.15,1)
(-303110.0301508775, -13569.40800266611, -19505.96059134758)

NOTE:

• The ideal gas models give identical result for H2, O2, and H2O. This cannot possibly be right????
• PR gives only a tiny correction to the ideal gas model for H2 and O2, but a relatively substantial correction for H2O -- which makes sense.

Questions:

1. Could it be that Clapeyron.jl -- since it probably is not designed for reaction equilibrium -- has not put a lot of emphasis into formation energies? For many cases without chemical reactions, formation energies are not important -- they get canceled out. (Not sure whether that is 100% true for all non-reacting cases, though). Anyways, lack of formation energies does not completely explain identical results for Gibbs energy for ideal gas; after all, the gases have different heat capacities.
2. Should I instead use my own function for ideal gas, and add residual energies (gibbs_energy_res)? I guess that should fix the problem?

[longemen3000]

Hello,

We are slowly getting there on terms of reactive equilibria. but first, let's talk about reference states, mainly to get this written somewhere.

Entropies, Enthalpies, Gibbs energies, Internal Energies and Helmholtz energies are defined in relation to an equation of state and a reference state, that is:

$$G = G_{EoS}(p,T) - G_0 (p,T)$$ $$H = H_{EoS}(p,T) - H_0 (p,T)$$ $$A = A_{EoS}(p,T) - A_0 (p,T)$$ $$S = S_{EoS}(p,T) - S_0 (p,T)$$ $$U = U_{EoS}(p,T) - U_0 (p,T)$$

If we work in a helmholtz or gibbs base, and after a little bit of math, we can express the functional form of $G_0$:

$$G = G_{EoS}(p,T)+ a_0 + a_1 T$$

On multicomponent equations of state:

$$G = G_{EoS}(p,T,n)+ \sum_{i}{n_i(a_{0,i} + a_{1,i}T)}$$

On a lot of thermodynamic properties, reference states are not used. that is because most algorithms are expressed as a difference between two states, and in equilibria, the temperature is the same, so the reference state terms cancel each other.

I'm intentionally omitting the fact that $a_0$ is an enthalpy and $a_1$ is an entropy, because in this formulation, we don't care about the value of $a_1$ and $a_2$, but we care about the result of evaluating the property at the reference state, that is:

$$S_0 = S(p0,T0) = S_{EoS}(p_0,T_0) + a_1$$ $$H_0 =H(p0,T0) = H_{EoS}(p_0,T_0) + a_0$$

If we have any equation of state and a pair of $S_0$ and $H_0$, we can find a set of constants $a_1$ and $a_2$ that satisfy the reference state conditions.

If we only care about normal equilibria, in pressure-temperature or pressure-volume specifications, then reference states are not used. However, we need
We added support for setting custom reference states in Clapeyron 0.5.11. We can access those in the following way:

julia> model = cPR(["isopentane","toluene"],idealmodel=ReidIdeal,reference_state = :stp)
PR{ReidIdeal, TwuAlpha, NoTranslation, vdW1fRule} with 2 components:
 "isopentane"
 "toluene"
Contains parameters: a, b, Tc, Pc, Mw
Reference state: stp

julia> reference_state(model)
ReferenceState(:stp) with 2 components:
 "isopentane" => 29230.6 - 136.942*T
 "toluene" => 42049.7 - 153.594*T

We created a model that follows the Standard Temperature and Pressure conditions (h = s = 0 at 1 bar, 0 °C fluid of the most stable phase):

julia> entropy(model,1e5,273.15,[1.0,0.0])
-0.0

julia> enthalpy(model,1e5,273.15,[1.0,0.0])
-4.980904577678296e-13

julia> enthalpy(model,1e5,273.15,[0.5,0.5]) #deviation from ideal behaviour
204.3629321343185

If you create an EoSModel without choosing a reference state, then the default is setting $a_0 = a_1 = 0$, that is, whatever reference state is hardcoded in the EoS is the state used.

On the technical implementation, the reference state struct is normally stored as a parameter of the ideal model:

julia> id_model = model.idealmodel
ReidIdeal with 2 components:
 "isopentane"
 "toluene"
Contains parameters: a, b, c, d, e, coeffs, reference_state, Mw
Reference state: stp

julia> id_model.params.reference_state
ReferenceState(:stp) with 2 components:
 "isopentane" => 29230.6 - 136.942*T
 "toluene" => 42049.7 - 153.594*T

So, we have almost all the pieces to add support for setting formation entropies and enthalpies as reference states. What is missing is:

• an API: how do we set the reference?
• a database: with S0,H0,T0,P0 and phase (some of those fields are redundant)

I haven't had the time to implement those last steps. any help in that regard is appreciated. a good api would probably be along the lines of:

model = cPR(["isopentane","toluene"],idealmodel=ReidIdeal,reference_state = :zero) #just to mark that we will initialize later
formation_reference_state!(model)

Now, as a design decision, if an ideal model is not selected, it defaults to BasicIdeal. That model corresponds to an ideal monoatomic gas (Cp = 5R/2), and it's fine for most calculations involving equilibria, excluding caloric and calculations where the reference state is relevant.

With no parameters, each instance of BasicIdeal is the same:

julia> BasicIdeal(["water"]) == BasicIdeal(["oxygen"]) == BasicIdeal()
true

That should explain the equality of gibbs free energies that you observed.

[B-LIE]

Thanks for detailed answer. Yes, I suspected you had used Cp/R = 5/2 or some similar standard value... Here is the molar Gibbs energies I got from BasicIdeal...

And here is what I got with my own implementation, which includes formation energies, etc.

I assume you are using the same expressions I found...

Entropy based on ideal gas, pure substance (tilde "~" on top means "molar"):

Entropy, mixture (the next-to-last term contains an entropy-of-mixing contribution):

Enthalpy, mixture:

I developed these expressions myself, but -- of course -- with a "cheat": I knew the result from one of the books of Prausnitz et al.

In the expressions, molar enthalpy at standard state and reference temperature equals molar enthalpy of formation, and is given in tables. Molar entropy at standard state and reference temperature is also given in many tables.

To extend this to internal energy and Helmholtz energy, molar volume at standard state and reference temperature is also needed...

• Would it be correct to assume ideal gas law for that, i.e,

? Or should a "real" value be used?

---

You seem to suggest constant heat capacity. Often, heat capacity of ideal gas is given as

Temperature dependent heat capacities are straightforward to integrate up in the expressions above.

[longemen3000]

for the gibbs energy, you can use the relation $G = H - TS$ to get the gibbs energy.

What you said about heat capacities is correct, for BasicIdeal we use a constant heat capacity, but other ideal models use different expressions for $C_p$. for example, The ReidIdeal model uses a 4th order polynomial, and AlyLeeIdeal model has terms of the form $A + B(CT⁻¹/sinh(CT⁻¹))² + D(ET⁻¹/cosh(ET⁻¹))²$. For all ideal models, we do the same integration procedure to calculate the ideal model contribution to the helmholtz energy (with one minor difference, we use volume instead of pressure). more details about the ideal helmholtz contribution exists in the CoolProp documentation.

On the reference state, the CRC handbook states:

The standard state of a pure gaseous substance is that of the substance as a (hypothetical) ideal gas at the standard state pressure (100 kPa)

So, $P_0 V_0 = R T_0$ seems valid to use.

[B-LIE]

I'm puzzled... I compared BasicIdeal to ReidIdeal. First, BasicIdeal:

# mod_H2O_pr = PR(["water"])
mod_H2O_ig = BasicIdeal(["water"])
mod_H2_ig = BasicIdeal(["hydrogen"])
mod_O2_ig = BasicIdeal(["oxygen"])
#
p_o = 1e5
G_ig(T;mod=mod_H2O_ig) = gibbs_energy(mod,p_o,T,1)

gives the following plot at standard state from 0 to 100 centigrades (repetition of plots above):

Next, ReidIdeal:

mod_H2O_igr = ReidIdeal(["water"])
mod_H2_igr = ReidIdeal(["hydrogen"])
mod_O2_igr = ReidIdeal(["oxygen"])

gives the following plot:

So -- almost mirrored slope. And "wrong" slope in ReidIdeal -- compared to my implementation of what I found in Prausnitz. Could this be due to the (default) choice of "reference state" :stp (standard pressure, I guess).

[B-LIE]

Hm. A little bit more of testing of ReidIdeal model with reference_state:

using Clapeyron

mod_H2O_igref = ReidIdeal(["water"], reference_state=:stp)
mod_H2_igref = ReidIdeal(["hydrogen"], reference_state=:stp)
mod_O2_igref = ReidIdeal(["oxygen"], reference_state=:stp)

When I do:

mod_H2O_igref.params.reference_state

etc., I get:

# "water" => 831.289 - 41.9727*T
# "hydrogen" => 713.122 - 41.5586*T
# "oxygen" => 728.643 - 41.613*T

Question:

1. I understand that the constant term is enthalpy of formation (Ht_oo below), while the temperature dependent term is entropy at the reference state (St_oo below), right?
2. What is the reference state? I assume :stp means standard pressure, which normally is 1 bar (or sometimes: 1 atm). But what is the reference temperature? 25 C?

The reference values extracted above are very different from what I find in Aylward and Findlay:

Ht_oo_H2 = 0                    # J/mol
Ht_oo_O2 = 0                    # J/mol
Ht_oo_H2O_g = -242e3            # J/mol
Ht_oo_H2O_L = -286e3            # J/mol
St_oo_H2 = 131                  # J/K mol
St_oo_O2 = 205                  # J/K mol
St_oo_H2O_g = 189               # J/K mol

Question:
3. What database does Clapeyron.jl extract "Ht_oo" and "St_oo" from?
4. How could I change these values? (I assume that is what was meant by how an API should look like...)
5. The a,b,c,d,e parameters are parameters for the cp model in Reid, Prausnitz, and Poling, right? The one with cp/R = a + bt + ct^2 + dt^3 + et^(-2), and t = T/1000???
6. Perhaps it would be useful for an API to change a,b,c,d,e, too? [I haven't checked the numbers.]

[longemen3000]

The default values depend on the implementation, but we use the CoolProp reference states (plus some others). Here is the list from the ReferenceState docs:

• :ashrae: h = s = 0 at -40 °C saturated liquid
• :iir: h = 200.0 kJ·kg⁻¹, s=1.0 kJ·kg⁻¹·K⁻¹ at 0 °C saturated liquid
• :nbp: h = s = 0 at 1 atm saturated liquid
• :stp: h = s = 0 at 1 bar, 0 °C fluid of the most stable phase
• :stp_old: h = s = 0 at 1 atm, 0 °C fluid of the most stable phase
• :ntp: h = s = 0 at 1 atm, 20 °C fluid of the most stable phase

We took those reference states from what CoolProp uses: https://coolprop.org/coolprop/HighLevelAPI.html#reference-states .

In addition to those predefined states, we support passing directly the values of p0,T0,H0,S0 to the reference state constructor. For example reference_state = :stp is equal to reference_state = ReferenceState(:volume,H0 = [0.0],S0 = [0.0],T0 = 273.15, P0 = 1e5).

On the values from ReidIdeal, by default, we use the coefficients from Poling, Bruce E. The Properties of Gases and Liquids, 5th edition, if i remember correctly. The coefficients in ReidIdeal use the form $C_p(T) = a + b T + c T^2 + d T^3 + eT^4$. We have a model that uses the form $C_p(T) = a + b T + c T^2 + d T^3 + \frac{e}{T^2}$ under the name ShomateIdeal, We are missing a big enough database for that model.

You can provide your own values of a,b,c,d,d,e in the following way (as shown in the ReidIdeal documentation):

  idealmodel = ReidIdeal(["water","butane"];
              userlocations = (a = [32.24, 9.487],
                          b = [0.00192, 0.3313],
                          c = [1.06e-5, -0.0001108],
                          d = [-3.6e-9, -2.822e-9],
                          e = [0.0, 0.0])
                          )

On the gibbs energy of the reid ideal model, It seems like the reference state is the problem, (you can check the slopes of ReidIdeal("water", reference_state = :stp) for comparison)

[B-LIE]

Thanks again for very informative response. On the :stp reference state...

What I have read is that a normal reference state (as in Aylward and Findlay) is 1 bar and 25 C. Also, in this literature, it is stated that enthalpy and Gibbs energy of formation for stable atoms/molecules from the periodic table is 0, while entropy for these are not zero. (Probably, entropy is zero at a lower temperature??).

This is at odds with what ReidIdeal gives with reference_state=:std:

# "water" => 831.289 - 41.9727*T
# "hydrogen" => 713.122 - 41.5586*T
# "oxygen" => 728.643 - 41.613*T

Here, according to Aylward and Findlay, the constant term for hydrogen and oxygen (but not for water!) should be zero.

Another problem: In the above numbers, the enthalpy term is in the order of 20 times larger than the entropy term. In Aylward and Findlay, the enthalpy term is in the order of 1000 times larger than the entropy term.

---

OK -- I can accept that entropy and enthalpy are relative values, and that there is freedom in choosing reference values.
But the reaction properties (enthalpy, entropy, Gibbs energy) must be the same whichever reference state is used -- if not, something is wrong, isn't it?

For the H2O <-> H2 + 0.5 O2 reaction, Gibbs energy of reaction is:

DrG = -G_H2O + G_H2 + 0.5G_O2

So I must be able to get the same DrG no matter what reference state is used. OK -- here is a comparison:

In the figure I show Gibbs molar energy of reaction for H2O <-> H2 + 0.5*O2

• "ig" is my own ideal gas implementation, using data from Aylward and Findlay
• "pr res, Reid, :stp" is adding residual Gibbs energies from Clapeyron, PR with ReidIdeal and :stp, added to my ideal gas model
• "ig, Reid, :stp": using Clapeyron with ReidIdeal and reference_state = :stp.
• "pr, Reid, :stp": using Clapeyron with idealmodel=ReidIdeal" and reference_state = :stp`.

In conclusion:

• The problem appears to be the ideal gas model in Clapeyron. In particular, I strongly suspect that the reference entropy and enthalpies are wrong. The residue appears to be ok.

I compared the parameters in ReidIdeal with those in Poling et al. (see below). Here is what I find for water:

julia> mod_H2O_igr.params.coeffs
SingleParam{NTuple{5, Float64}}("Reid Coefficients") with 1 component:
 "water" => (36.54206320678348, -0.03480434051958945, 0.000116818199785053, -1.3003819534791665e-7, 5.2547403746728466e-11)

These numbers are identical to what I get by taking a0-a4 in Poling et al., and multiply with the gas constant. So: Poling gives parameters for cp/R, while ReidIdeal gives parameters for cp.

When it comes to the formation energies:

julia> mod_H2O_igr.params.reference_state
ReferenceState(:stp) with 1 component:
 "water" => 831.289 - 41.9727*T

From Poling et al.:

T_o = 25 + 273.15
julia> Ht_oo_H2O_pol = -241.81e3
julia> Gt_oo_H2O_pol = -228.42e3
julia> St_oo_H2O_pol = (Ht_oo_H2O_pol - Gt_oo_H2O_pol)/T_o
-44.9102800603723

Again: the reference state values used in ReidIdeal are very different from what Poling specifies. OK -- it is not 100% clear to me what "water" => 831.289 - 41.9727*T means. I am assuming that it means that the reference state molar enthalpy is 831.289 J/mol while the reference state molar entropy is 41.9727 J/mol.

So -- I realize that ReidIdeal seems to use a different standard state than Poling et al. uses. BUT the molar Gibbs energy of reaction must be the same in both cases. And I get a wrong Gibbs energy of reaction if I use ReidIdeal and the :stp.

How do I know that ReidIdeal gives the wrong result, and my ideal gas model/number is correct? Because my values are compatible with the published and experimental corresponding Nernst voltage at the reference state.

---

OK -- I tried your trick wrt. changing the reference state:

NOTE: I got an error message when I did:

reference_state = ReferenceState(:volume,H0 = [Ht_oo_H2O_g],S0 = [St_oo_H2O_g],T0 = T_o, P0 = p_o)

The error message disappeared when I instead used scalars for H0 and S0:

reference_state = ReferenceState(:volume,H0 = Ht_oo_H2O_g,S0 = St_oo_H2O_g,T0 = T_o, P0 = p_o)

---

One more double-peculiarity that I notice...

julia> my_ref_state = ReferenceState(:volume;H0=242e3,S0=189,T0=25+273.15,P0=1e5)
julia> mod_H2O_igr_r = ReidIdeal(["water"]; reference_state=my_ref_state)

julia> mod_H2O_igr_r.params.reference_state
ReferenceState(:volume) with 1 component:
 "water" => 241995.0 - 228.043*T

• Why does the values for H0 and S0 change from the values I specified???

And...

julia> mod_H2O_pr = PR(["water"]; idealmodel=ReidIdeal, reference_state=my_ref_state)

julia> propertynames(mod_H2O_pr_r)
(:components, :alpha, :mixing, :translation, :params, :idealmodel, :references)

julia> propertynames.idealmodel.params.reference_state
ReferenceState(:volume) with 1 component:
 "water" => 287655.0 - 351.118*T

Why is the change in H0 and S0 (which is really only used in the ideal gas model) changed even differently when I embed the ReidIdeal within PR?

---

Poling, Prausnitz, O'Connell (2001): The Properties of Gases and Liquids, 5th ed.
-- Provides enthalpy and Gibbs energy of formation, but not entropy, at reference state 25 C and 1 atm (not 1 bar!). They provide parameters for the Cp model of the form in ReidIdeal. (It is not clear to me whether the parameters a-e in ReidIdeal are for cp/R, or for cp.)

Elliott, Diky, Knotts IV, and Wilding (2023). The Properties of Gases and Liquids, 6th ed.
-- Provides enthalpy and Gibbs energy of formation + entropy at reference pressure. The reference state is still 25 C and 1 bar according to their info in Appendix A (but 25 C and 1 atm according to their Chapter 4... perhaps not completely updated text).

When it comes to cp, they provide a mixture of 3 models: the classical with polynomial in T up to fourth order (parameters a-e), one with sinh and cosh terms (parameters a-e), and one with exp terms (parameters a-g).

[longemen3000]

On the ReidIdeal model. our current implementation does have some constant terms (there is an integration from a temperature T0 = 298 to T). The constants printed by the ReferenceState are calculated in a way that allows the equation of state object to match the enthalpy and entropy specifications at the reference state conditions. They do not have a physical meaning, in the sense that they don't represent a meaninful value of entropy (or enthalpy). After all, our approach is to consider inside the evaluation of the a_ideal function (reduced ideal helmholtz energy) as a black box.

Take for example the IAPWS95 EoS. their reference state is the triple point, but we can "translate" the reference state to s.t.p. conditions if we add a pair of constants a0 and a1. That is the role of those values, they work as a translation from the reference state stored in the EoS and any reference state we want.

The error provided vector is genuine, but i believe it is fixed in the last release. I also added what could be seen as a start of the support for ideal gas formation reference states:
If you check the master branch:

julia> ref = Clapeyron.IGFormReferenceState(["water"])
ReferenceState(:ideal_gas) (not set)

julia> model = ReidIdeal(["water"],reference_state = ref)
ReidIdeal with 1 component:
 "water"
Contains parameters: a, b, c, d, e, coeffs, reference_state, Mw
Reference state: ideal_gas

julia> enthalpy(model,1e5,298.15)
-241800.0

julia> entropy(model,1e5,298.15)
188.8

julia> model.params.reference_state.H0
1-element Vector{Float64}:
 -241800.0

julia> model.params.reference_state.S0
1-element Vector{Float64}:
 188.8

julia> model.params.reference_state
ReferenceState(:ideal_gas) with 1 component:
 "water" => -241805.0 - 227.843*T

note that, while, at 1 bar, 298.15 K, the values of entropy and enthalpy evaluated by the EoS correspond to the formation values, the constants a0 and a1 do not match.

I'm starting to think that displaying a0 and a1 wasn't the best decision (the original motivation was to see if the reference state calculation was done correctly). The display method for ReferenceState needs to be clearer in that regard.