Consistent and improved set of equations for saturation vapor pressur…

…e and derivative

- Previously, we had two functions `svapor` based on Hess 1959 that estimated saturation vapor pressure over water in [mmHg] and `slope_svp_to_t` that estimated the derivative based on Allen et al. 2005, but not of the implemented equation. The equation by Hess 1959 was increasingly (but slightly) underestimating saturation vapor pressure at temperatures > 20 C compared to Haynes 2017.

--> Replace functions with `svp` that calculates saturation vapor pressure over water and ice based on equations by Huang 2018 and estimates the gradient by taking the derivation of these equations.

- Updated unit tests

Haynes, W. M., editor. 2017. Vapor pressure and other saturation properties of water. Page 6.5-6.6 CRC Handbook of Chemistry and Physics. 97th Edition (Internet Version). CRC Press/Taylor & Francis, Boca Raton, FL.
Huang, J. 2018. A Simple Accurate Formula for Calculating Saturation Vapor Pressure of Water and Ice. Journal of Applied Meteorology and Climatology 57:1265–1272.

SW_Flow_lib_PET.c

@@ -899,24 +899,6 @@ double blackbody_radiation(double T) {
 
 
 
-/** @brief Slope of the saturation vapor pressure-temperature curve
-
-  Based on equation 13 (chapter 3) of Allen et al. (1998) @cite allen1998 and
-  equation 5 of Allen et al. (2005) @cite ASCE2005.
-
-  Prior to 2012-Oct-11, equation `vapor * 3010.21 / (kelvin * kelvin)`
-  was of unknown origin.
-  However, results are very similar to the current implementation.
-
-  @param es_at_tmean Saturation vapor pressure at average temperature [kPa]
-  @param tmean Average temperature for the day [°C]
-
-  @return slope of es:T at T=Ta [kPa / K]
-*/
-double slope_svp_to_t(double es_at_tmean, double tmean) {
-  return 4098. * es_at_tmean / squared(tmean + 237.3);
-}
-
 /** @brief Atmospheric pressure based on elevation
 
   Based on equation 7 (chapter 3) of Allen et al. (1998) @cite allen1998 and
@@ -947,26 +929,51 @@ double psychrometric_constant(double pressure) {
 }
 
 
-/** @brief Saturation vapor pressure of water.
 
-  Equations based on Hess 1959. @cite Hess1961
+/** @brief Saturation vapor pressure of water and ice and the
+           slope of the saturation vapor pressure-temperature curve
+
+  Empirical equation derived from integrating the Clausius–Clapeyron equation
+  by Huang 2018 @cite huang2018JAMC.
 
-  @author SLC
+  The slope of the svp-T curve is obtained by derivation for temperature.
 
-  @param  temp Average temperature for the day (°C).
+  @param[in] T Temperature [C]
+  @param[out] Slope of es:T [kPa / K]
 
-  @return svapor Saturation vapor pressure (mm of Hg/°F).
+  @return Saturation vapor pressure [kPa]
 */
-double svapor(double temp) {
-	double par1, par2;
+double svp(double T, double *slope_svp_to_t) {
+  double tmp, tmp0, tmp1, tmp2, tmp3, dp, svp;
+
+  if (T > 0) {
+    // Derivation for Huang 2018: eq. 17
+    tmp0 = T + 105.;
+    tmp1 = pow(tmp0, 1.57);
+    tmp = T + 237.1;
+    tmp2 = 4924.99 / squared(tmp);
+    tmp3 = exp(34.494 - 4924.99 / tmp);
+    dp = 1.57 * pow(tmp0, 0.57);
+
+  } else {
+    // Derivation for Huang 2018: eq. 18
+    tmp0 = T + 868.;
+    tmp1 = squared(tmp0);
+    tmp = T + 278.;
+    tmp2 = 6545.8 / squared(tmp);
+    tmp3 = exp(43.494 - 6545.8 / tmp);
+    dp = 2. * tmp0;
+  }
 
-	par1 = 1. / (temp + 273.);
-	par2 = log(6.11) + 5418.38 * (.00366 - par1); /*drs: par2 = ln(es [mbar]) = ln(es(at T = 273.15K) = 6.11 [mbar]) + (mean molecular mass of water vapor) * (latent heat of vaporization) / (specific gas constant) * (1/(273.15 [K]) - 1/(air temperature [K])) */
+  svp = tmp3 / tmp1;
+  *slope_svp_to_t = svp * tmp2 - tmp3 / squared(tmp1) * dp;
+  (*slope_svp_to_t) *= 1e-3;
 
-	return (exp(par2) * .75);
+  return 1e-3 * svp;
 }
 
 
+
 /**
 @brief Daily potential evapotranspiration
 

SW_Flow_lib_PET.h

@@ -50,8 +50,7 @@ double blackbody_radiation(double T);
 double petfunc(double H_g, double avgtemp, double elev,
   double reflec, double humid, double windsp, double cloudcov);
 
-double svapor(double temp);
-double slope_svp_to_t(double es_at_tmean, double tmean);
+double svp(double T, double *slope_svp_to_t);
 double atmospheric_pressure(double elev);
 double psychrometric_constant(double pressure);
 

test/test_SW_Flow_Lib_PET.cc

@@ -697,7 +697,7 @@ namespace
         {0.7, 0.7, 0.4, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.4},
       // Climate normals for Madison, WI
       // "WMO Climate Normals for MADISON/DANE CO REGIONAL ARPT, WI 1961–1990".
-      // National Oceanic and Atmospheric Administration. Retrieved March 10, 2014.
+      // National Oceanic and Atmospheric Administration. Retrieved Jul 3, 2020.
       // ftp://ftp.atdd.noaa.gov/pub/GCOS/WMO-Normals/TABLES/REG_IV/US/GROUP4/72641.TXT
       cloud_cover[12] =
         // Element 20:  Sky Cover (Cloud Cover)
@@ -745,24 +745,34 @@ namespace
 
 
 
-  // Test saturated vapor pressure function
-  TEST(SW2_PET_Test, svapor)
+  // Test saturation vapor pressure functions
+  TEST(SW2_PET_Test, svp)
   {
     int i;
     double
       // Temperature [C]
       temp_C[] = {-30, -20, -10, 0, 10, 20, 30, 40, 50, 60},
 
-      // Expected saturated vapor pressure
+      // Expected saturation vapor pressure [kPa]
+      check_svp,
       expected_svp[] = {
-        0.3889344, 0.9389376, 2.1197755,
-        4.5085235,
-        9.0911046, 17.4746454, 32.1712519, 56.9627354, 97.3531630, 161.1126950
+        0.0380009, 0.103226, 0.2598657, 0.6112912,
+        1.2281879, 2.3393207, 4.247004, 7.3849328, 12.3517837, 19.9461044
+      },
+
+      // Expected slope of svp - temperature curve [kPa / K]
+      check_svp_to_t,
+      expected_svp_to_T[] = {
+        0.0039537, 0.0099076, 0.0230775, 0.0503666,
+        0.0822986, 0.1449156, 0.2437929, 0.3937122, 0.6129093, 0.9231149
       };
 
     for (i = 0; i < 10; i++)
     {
-      EXPECT_NEAR(svapor(temp_C[i]), expected_svp[i], tol6);
+      check_svp = svp(temp_C[i], &check_svp_to_t);
+
+      EXPECT_NEAR(check_svp, expected_svp[i], tol6);
+      EXPECT_NEAR(check_svp_to_t, expected_svp_to_T[i], tol6);
     }
   }