diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index f6eeeaf3724..99fd4c28eaa 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -35,6 +35,7 @@ Moist Thermodynamics dewpoint_from_relative_humidity dewpoint_from_specific_humidity equivalent_potential_temperature + ice_saturation_vapor_pressure mixing_ratio mixing_ratio_from_relative_humidity mixing_ratio_from_specific_humidity diff --git a/docs/api/references.rst b/docs/api/references.rst index 5c64c0603cb..9bfd417ade7 100644 --- a/docs/api/references.rst +++ b/docs/api/references.rst @@ -141,6 +141,10 @@ References *Journal of Geodesy* **74**, 128-133, doi:`10.1007/s001900050278 `_. +.. [MurphyKoop2005] Murphy, D.M. and T. Koop, 2005: Review of the vapour pressures of ice and + supercooled water for atmospheric applications. *Q.J.R. Meteorol. Soc.*, **131**, + 1539-1565, doi:`10.1256/qj.04.94 `_. + .. [NOAA1976] National Oceanic and Atmospheric Administration, National Aeronautics and Space Administration, and U. S. Air Force, 1976: `U. S. Standard Atmosphere 1976 `_, diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a85c4d6df4d..909cff99b8f 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1014,7 +1014,7 @@ def vapor_pressure(pressure, mixing_ratio): @preprocess_and_wrap(wrap_like='temperature') @process_units({'temperature': '[temperature]'}, '[pressure]') def saturation_vapor_pressure(temperature): - r"""Calculate the saturation water vapor (partial) pressure. + r"""Calculate the saturation (equilibrium) water vapor (partial) pressure. Parameters ---------- @@ -1035,14 +1035,51 @@ 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: + The formula used is from eq. 10 in [MurphyKoop2005]_ for T in degrees Kelvin: - .. math:: 6.112 e^\frac{17.67T}{T + 243.5} + .. math:: e^[54.842763 - \frac{6763.22}{T} - 4.210\:ln(T) + 0.000367T + + tanh(0.0415(T - 218.8))(53.878 - \frac{1331.22}{T} + - 9.44523\:ln(T) + 0.014025T)] """ - # 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) + # valid for 123 < T < 332 + return np.exp( + 54.842763 - 6763.22 / temperature - 4.210 * np.log(temperature) + + 0.000367 * temperature + np.tanh(0.0415 * (temperature - 218.8)) + * (53.878 - 1331.22 / temperature - 9.44523 * np.log(temperature) + + 0.014025 * temperature) + ) + + +@exporter.export +@preprocess_and_wrap(wrap_like='temperature') +@process_units({'temperature': '[temperature]'}, '[pressure]') +def ice_saturation_vapor_pressure(temperature): + r"""Calculate the ice saturation (equilibrium) water vapor (partial) pressure. + + Parameters + ---------- + temperature : `pint.Quantity` + Air temperature + + Returns + ------- + `pint.Quantity` + Ice saturation water vapor (partial) pressure + + See Also + -------- + vapor_pressure + + The formula used is from eq. 7 in [MurphyKoop2005]_ for T in degrees Kelvin: + + .. math:: e^[9.550426 - \frac{5723.265}{T} + 3.53068\:ln(T) - 0.00728332T] + + """ + # valid for T > 110 K + return np.exp( + 9.550426 - 5723.265 / temperature + 3.53068 * np.log(temperature) + - 0.00728332 * temperature ) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index cc8e5827bfc..c71f7df059d 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -25,7 +25,8 @@ relative_humidity_from_specific_humidity, relative_humidity_wet_psychrometric, saturation_equivalent_potential_temperature, saturation_mixing_ratio, - saturation_vapor_pressure, showalter_index, + saturation_vapor_pressure, ice_saturation_vapor_pressure, + showalter_index, specific_humidity_from_dewpoint, specific_humidity_from_mixing_ratio, static_stability, surface_based_cape_cin, temperature_from_potential_temperature, thickness_hydrostatic, @@ -295,6 +296,16 @@ def test_sat_vapor_pressure(): assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 2) +def test_ice_sat_vapor_pressure(): + """Test ice_saturation_vapor_pressure calculation.""" + temp = (np.arange(1, 38) * units.degC).to(units.degK) + assert np.all(saturation_vapor_pressure(temp) + < ice_saturation_vapor_pressure(temp)) + temp = (temp.to(units.degC) * (-1)).to(units.degK) + assert np.all(saturation_vapor_pressure(temp) + > ice_saturation_vapor_pressure(temp)) + + def test_sat_vapor_pressure_scalar(): """Test saturation_vapor_pressure handles scalar values.""" es = saturation_vapor_pressure(0 * units.degC)