Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Analytic Saturation Vapor Pressure #3726

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Analytic Saturation Vapor Pressure #3726

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9da7a65
Implement analytic satvap
dcamron Dec 25, 2024
dd12f6b
Update test values
dcamron Dec 25, 2024
6243bfe
Add satvap over ice
dcamron Jan 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 58 additions & 6 deletions src/metpy/calc/thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1545,7 +1545,7 @@ def saturation_vapor_pressure(temperature):
>>> from metpy.calc import saturation_vapor_pressure
>>> from metpy.units import units
>>> saturation_vapor_pressure(25 * units.degC).to('hPa')
<Quantity(31.6742944, 'hectopascal')>
<Quantity(31.623456, 'hectopascal')>

See Also
--------
Expand All @@ -1556,14 +1556,66 @@ def saturation_vapor_pressure(temperature):
Instead of temperature, dewpoint may be used in order to calculate
the actual (ambient) water vapor (partial) pressure.

The formula used is that from [Bolton1980]_ for T in degrees Celsius:
Implemented solution from [Ambaum2020]_, Eq. 13,
.. math:: e = e_{s0} \frac{T_0}{T}^{(c_{pl} - c_{pv}) / R_v} \exp{
\frac{L_0}{R_v T_0} - \frac{L}{R_v T}}

.. math:: 6.112 e^\frac{17.67T}{T + 243.5}
"""
latent_heat = water_latent_heat_vaporization._nounit(temperature)
heat_power = (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v) / mpconsts.nounit.Rv
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After reading through section 5, I don't think we should rely on our constants here, but use the tuned value for $C_{pl} - C_{pv}$ that's used in the paper, 2180 J / kg / K (eq. 14). Our values are essentially assuming an ideal gas for water vapor, and nominally for the triple point; the discussion in Ambaum 2020 explains how their value minimizes errors against empirical data from IAPWS95 in the 0-100 degC range.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm open to the idea that we allow users to override this value.

exp_term = ((mpconsts.nounit.Lv / mpconsts.nounit.T0 - latent_heat / temperature)
/ mpconsts.nounit.Rv)

return (
mpconsts.nounit.sat_pressure_0c
* (mpconsts.nounit.T0 / temperature) ** heat_power
* np.exp(exp_term)
)
Comment on lines +1564 to +1573
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like it matches the paper and I don't see any better ways to approach the array math, so I'd say we're good to go here.



@exporter.export
@preprocess_and_wrap(wrap_like='temperature')
@process_units({'temperature': '[temperature]'}, '[pressure]')
def saturation_vapor_pressure_ice(temperature):
r"""Calculate the saturation water vapor (partial) pressure over ice.

Parameters
----------
temperature : `pint.Quantity`
Air temperature

Returns
-------
`pint.Quantity`
Saturation water vapor (partial) pressure

Examples
--------
>>> from metpy.calc import saturation_vapor_pressure_ice
>>> from metpy.units import units
>>> saturation_vapor_pressure_ice(-25 * units.degC).to('hPa')
<Quantity(0.632434749, 'hectopascal')>

See Also
--------
saturation_vapor_pressure, vapor_pressure

Notes
-----
Implemented solution from [Ambaum2020]_, Eq. 17,
.. math:: e = e_{i0} \frac{T_0}{T}^{(c_{pi} - c_{pv}) / R_v} \exp{
\frac{L_{s0}}{R_v T_0} - \frac{L_s}{R_v T}}

"""
# Converted from original in terms of C to use kelvin.
return mpconsts.nounit.sat_pressure_0c * np.exp(
17.67 * (temperature - 273.15) / (temperature - 29.65)
latent_heat = water_latent_heat_sublimation._nounit(temperature)
heat_power = (mpconsts.nounit.Cp_i - mpconsts.nounit.Cp_v) / mpconsts.nounit.Rv
exp_term = ((mpconsts.nounit.Ls / mpconsts.nounit.T0 - latent_heat / temperature)
/ mpconsts.nounit.Rv)

return (
mpconsts.nounit.sat_pressure_0c
* (mpconsts.nounit.T0 / temperature) ** heat_power
* np.exp(exp_term)
)


Expand Down
10 changes: 5 additions & 5 deletions tests/calc/test_indices.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ def test_precipitable_water():
"""Test precipitable water with observed sounding."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
pw = precipitable_water(data['pressure'], data['dewpoint'], top=400 * units.hPa)
truth = 22.60430651 * units.millimeters
truth = 22.5880919 * units.millimeters
assert_array_almost_equal(pw, truth, 4)


Expand All @@ -32,7 +32,7 @@ def test_precipitable_water_no_bounds():
pressure = data['pressure']
inds = pressure >= 400 * units.hPa
pw = precipitable_water(pressure[inds], dewpoint[inds])
truth = 22.60430651 * units.millimeters
truth = 22.5880919 * units.millimeters
assert_array_almost_equal(pw, truth, 4)


Expand All @@ -43,7 +43,7 @@ def test_precipitable_water_bound_error():
dewpoint = np.array([25.5, 24.1, 23.1, 21.2, 21.1, 19.4, 19.2, 19.2, -87.1, -86.5, -86.5,
-86.5, -88.1]) * units.degC
pw = precipitable_water(pressure, dewpoint)
truth = 89.86846252697836 * units('millimeters')
truth = 89.7845131 * units('millimeters')
assert_almost_equal(pw, truth, 5)


Expand All @@ -59,7 +59,7 @@ def test_precipitable_water_nans():
-28.3, np.nan, -32.6, np.nan, -33.8, -35., -35.1, -38.1, -40.,
-43.3, -44.6, -46.4, -47., -49.2, -50.7]) * units.degC
pw = precipitable_water(pressure, dewpoint)
truth = 4.003660322395436 * units.mm
truth = 3.99738502 * units.mm
assert_almost_equal(pw, truth, 5)


Expand Down Expand Up @@ -161,7 +161,7 @@ def test_precipitable_water_xarray():
press = xr.DataArray(data['pressure'].m, attrs={'units': str(data['pressure'].units)})
dewp = xr.DataArray(data['dewpoint'], dims=('press',), coords=(press,))
pw = precipitable_water(press, dewp, top=400 * units.hPa)
truth = 22.60430651 * units.millimeters
truth = 22.5880919 * units.millimeters
assert_almost_equal(pw, truth)


Expand Down
Loading
Loading