Water fugacity

Flow laws of hydrous rocks and the solubility of water in nominally anhydrous minerals depend on the fugacity of water, $f_{{\mathrm{H}}_2\mathrm{O}}$. GeoParams provides the following routines; all take pressure in Pa and temperature in K and return the fugacity in Pa, or accept Unitful quantities.

Pitzer & Sterner (1994)

compute_water_fugacity solves the equation of state of Pitzer & Sterner (1994) for the molar density $\rho$ of H2O at given $P$ and $T$,

$$\frac{P}{RT}=\rho +c_1{\rho }^2-{\rho }^2\frac{c_3+2c_4\rho +3c_5{\rho }^2+4c_6{\rho }^3}{(c_2+c_3\rho +c_4{\rho }^2+c_5{\rho }^3+c_6{\rho }^4{)}^2}+c_7{\rho }^2{e}^{-c_8\rho }+c_9{\rho }^2{e}^{-c_{10}\rho },$$

where each coefficient is a function of temperature, $c_i=c_{i1}{T}^{-4}+c_{i2}{T}^{-2}+c_{i3}{T}^{-1}+c_{i4}+c_{i5}T+c_{i6}{T}^2$. The fugacity follows from the residual Helmholtz energy $A_{\mathrm{res}}$:

$$\ln f_{{\mathrm{H}}_2\mathrm{O}}=\ln (\rho RT)+\frac{A_{\mathrm{res}}}{RT}+\frac{P}{\rho RT}-1,$$$$\frac{A_{\mathrm{res}}}{RT}=c_1\rho +\frac{1}{c_2+c_3\rho +c_4{\rho }^2+c_5{\rho }^3+c_6{\rho }^4}-\frac{1}{c_2}-\frac{c_7}{c_8}(e^{-c_8\rho }-1)-\frac{c_9}{c_{10}}(e^{-c_{10}\rho }-1).$$

The implementation reproduces the online calculator of A. Withers to better than 0.4% for 0.1-2.8 GPa and 100-1500 °C.

Crustal approximation

compute_water_fugacity_Shinevar uses the closed-form approximation of Shinevar et al. (2015),

$$f_{{\mathrm{H}}_2\mathrm{O}}\approx A_1\exp (-\frac{A_2+P{A}_3}{RT}),$$

which is within a factor 1.5 of the Pitzer & Sterner value for 0.2-1 GPa and 600-1000 K.

Water content of olivine

compute_water_content_olivine gives the hydrogen solubility in olivine of Zhao et al. (2004),

$$C_{\mathrm{OH}}=A_w\exp (-\frac{E_w+P{V}_w}{RT})f_{{\mathrm{H}}_2\mathrm{O}},$$

in H/10⁶ Si (with $f_{{\mathrm{H}}_2\mathrm{O}}$ in MPa in this expression; the routine takes it in Pa).

Example

julia> using GeoParams

julia> fH2O = compute_water_fugacity(1GPa, 800C);    # Pitzer & Sterner (1994), in Pa

julia> round(MPa, fH2O, digits = 1)
1642.0 MPa

julia> round(MPa, compute_water_fugacity_Shinevar(1GPa, 800C), digits = 1)   # crustal approximation
1575.3 MPa

julia> round(compute_water_content_olivine(fH2O, 1GPa, 800C))   # H/10⁶ Si, water-saturated olivine
157.0

Routines

GeoParams.WaterFugacity.compute_water_fugacity Function

fH2O = compute_water_fugacity(P, T)

Fugacity of pure H2O [Pa] at pressure P [Pa] and temperature T [K], computed with the equation of state of Pitzer & Sterner (1994). The density of H2O is found by solving the equation of state with a bracketed root finder; the fugacity then follows from the residual Helmholtz energy $A_{\mathrm{res}}$:

$$\ln f_{{\mathrm{H}}_2\mathrm{O}}=\ln (\rho RT)+\frac{A_{\mathrm{res}}}{RT}+\frac{P}{\rho RT}-1$$

The equation of state is calibrated for 0-10 GPa and 0-2000 °C. P and T can also be given as Unitful quantities, in which case the fugacity is returned in Pa. For non-positive pressure (e.g. near the free surface of a numerical model) the fugacity is zero, its physical limit; a wet creep law then gives zero strain rate.

Reference: Pitzer, K.S. & Sterner, S.M. (1994), Equations of state valid continuously from zero to extreme pressures for H2O and CO2, J. Chem. Phys. 101, 3111-3116.

Example

julia> round(compute_water_fugacity(1.0e9, 1073.15) / 1.0e6, digits = 1)   # MPa
1642.0

julia> round(MPa, compute_water_fugacity(1GPa, 800C), digits = 1)
1642.0 MPa

source

GeoParams.WaterFugacity.compute_water_fugacity_Shinevar Function

fH2O = compute_water_fugacity_Shinevar(P, T)

Approximate fugacity of H2O [Pa] at pressure P [Pa] and temperature T [K] for crustal conditions, following Shinevar, Behn & Hirth (2015):

$$f_{{\mathrm{H}}_2\mathrm{O}}=A_1\exp (-\frac{A_2+P{A}_3}{RT})$$

with $A_1=5.521$ GPa, $A_2=31.28$ kJ/mol and $A_3=-2.009\times {10}^{-5}$ m³/mol (parameters as quoted in eq. 5.7 of Becker & Faccenna's geodynamics textbook). For 0.2-1 GPa and 600-1000 K this is within a factor 1.5 of compute_water_fugacity, the deviation grows for colder conditions at high pressure (factor ~2.3 at 1 GPa and 500 K). P and T can also be given as Unitful quantities.

Reference: Shinevar, W.J., Behn, M.D. & Hirth, G. (2015), Compositional dependence of lower crustal viscosity, Geophys. Res. Lett. 42, 8333-8340.

Example

julia> round(compute_water_fugacity_Shinevar(1.0e9, 800.0) / 1.0e6, digits = 1)   # MPa
1026.6

source

GeoParams.WaterFugacity.compute_water_content_olivine Function

C_OH = compute_water_content_olivine(fH2O, P, T)

Solubility of water (hydrogen) in olivine, $C_{\mathrm{OH}}$ in H/10⁶ Si, as function of water fugacity fH2O [Pa], pressure P [Pa] and temperature T [K] after Zhao, Ginsberg & Kohlstedt (2004):

$$C_{\mathrm{OH}}=A_w\exp (-\frac{E_w+P{V}_w}{RT})f_{{\mathrm{H}}_2\mathrm{O}}$$

with $A_w=26$ H/10⁶ Si/MPa, $E_w=40$ kJ/mol and $V_w=10\times {10}^{-6}$ m³/mol (as quoted in eq. 5.9 of Becker & Faccenna's geodynamics textbook); in this expression $f_{{\mathrm{H}}_2\mathrm{O}}$ is in MPa, which is taken care of internally. The arguments can also be given as Unitful quantities; the result is always a plain number in H/10⁶ Si.

For water-saturated conditions, combine it with compute_water_fugacity.

Reference: Zhao, Y.-H., Ginsberg, S.B. & Kohlstedt, D.L. (2004), Solubility of hydrogen in olivine: dependence on temperature and iron content, Contrib. Mineral. Petrol. 147, 155-161.

Example

julia> P, T = 2.0e9, 1573.15;

julia> round(compute_water_content_olivine(compute_water_fugacity(P, T), P, T))   # H/10⁶ Si
3018.0

source