Update e_s to a more accurate formulation and add e_i (based on Murphy and Koop, 2005)

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

docs/api/references.rst

@@ -141,6 +141,10 @@ References
            *Journal of Geodesy* **74**, 128-133, doi:`10.1007/s001900050278
            <https://doi.org/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 <https://doi.org/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
            <https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf>`_,

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
     )
 
 

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)