Consolidate methods and standardize units for calculating saturated vapor pressure deficit

CHANGELOG.md

@@ -11,7 +11,8 @@ section for the next release.
 
 * Add function `clip_and_save_raster_file()` for subsetting rasters to match a polygon of interest (#3537).
 * Add CH4 and N2O to standard_vars in PEcAn.utils
-
+* New function `sat_vapor_pressure()` added for computing saturation vapor pressure from temperature using various methods.
+* 
 ## [1.9.0] - 2025-05-25
 
 ### Added

modules/data.atmosphere/NAMESPACE

@@ -94,6 +94,7 @@ export(process_gridded_noaa_download)
 export(qair2rh)
 export(read.register)
 export(rh2qair)
+export(sat_vapor_pressure)
 export(save.betas)
 export(save.model)
 export(site.lst)

modules/data.atmosphere/NEWS.md

@@ -2,7 +2,7 @@
 
 * `ERA5_met_process()` can now process ensemble data efficiently in parallel using new option `n_cores`
 * Dependency `ggplot2` is now suggested rather than required. It is used in two vignettes and for optional diagnostic plots from `debias_met_regression`.
-
+* New function `sat_vapor_pressure()` added for computing saturation vapor pressure from temperature using various methods.
 
 # PEcAn.data.atmosphere 1.9.0
 

modules/data.atmosphere/R/metutils.R

@@ -69,32 +69,94 @@ get.vpd <- function(rh, temp) {
   ## calculate saturation vapor pressure
   es <- get.es(temp)
   ## calculate vapor pressure deficit
-  return(((100 - rh)/100) * es)
+  return(((100 - rh) / 100) * es)
 } # get.vpd
 
-##' Calculate saturation vapor pressure
-##'
-##' @title get es
-##' @param temp temperature in degrees C
-##' @return saturation vapor pressure in mb
-##' @export
-##' @author David LeBauer
-##' @examples
-##' temp <- -30:30
-##' plot(temp, get.es(temp))
-get.es <- function(temp) {
-  return(6.11 * exp((2500000/461) * (1/273 - 1/(273 + temp))))
-} # get.es
+#' Saturation vapor pressure (t2es)
+#'
+#' Compute saturation vapor pressure from temperature using one of the
+#' following methods:
+#' - (Default) Clausius–Clapeyron (FAO-56 style) — Recommended for most applications.
+#'   Commonly used approximation for terrestrial ecosystem models, consistent with Penman-Monteith
+#'   and FAO-56 (Allen et al, 1998).
+#' - Magnus — More accurate in the range −40 to +50 C. Coefficients as in Alduchov & Eskridge (1996).
+#' - Goff–Gratch - Highest accuracy; use when following WMO-style recommendations. Goff–Gratch 1946; WMO, 2014.
+#'
+#' Each method uses different units internally, users can specify units
+#' for both inputs and outputs, with defaults "degC" and "kPa", respectively.
+#' 
+#' @param temp numeric vector of temperatures
+#' @param method one of "Magnus","ClausiusClapeyron" (default), or "GoffGratch".
+#' See details for references.
+#' @param temp_units input temperature units ("degC","K","degF"), default "degC"
+#' @param out_units output pressure units ("kPa","hPa","Pa","mb"), default "kPa"
+#' @return numeric vector in `out_units`
+#' @aliases t2es
+#'
+#' @references
+#' Alduchov, O. A., & Eskridge, R. E. (1996). Improved Magnus Form Approximation of Saturation Vapor Pressure. J. Appl. Meteor.*, 35(4), 601–609. <doi:10.1175/1520-0450(1996)035<0601:IMFAOS>2.0.CO;2>
+#'
+#' Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998). **Crop evapotranspiration – Guidelines for computing crop water requirements.** FAO Irrigation and Drainage Paper 56.
+#'
+#' Goff, J. A., & Gratch, S. (1946). Low-pressure properties of water from −160 to 212F. Trans. ASHVE, 52, 95–122.
+#'
+#' WMO (2014) Guide to Instruments and Methods of Observation (WMO-No. 8), ch. 4.
+#' @md
+#' @author David LeBauer
+#' @export
+sat_vapor_pressure <- function(
+    temp,
+    temp_units = "degC",
+    out_units = "kPa",
+    method = c("ClausiusClapeyron", "Magnus", "GoffGratch")) {
+  method <- match.arg(method)
+  # normalize common alias
+  if (tolower(out_units) == "mb") out_units <- "hPa"
+
+  if (method == "Magnus") {
+    # canonical temp: degC; canonical pressure: kPa
+    Tc <- units::ud_convert(temp, temp_units, "degC")
+    es_kPa <- 0.61078 * exp((17.27 * Tc) / (Tc + 237.3))
+    return(units::ud_convert(es_kPa, "kPa", out_units))
+  }
+
+  if (method == "ClausiusClapeyron") {
+    # canonical temp: degC; canonical pressure: hPa
+    Tc <- units::ud_convert(temp, temp_units, "degC")
+    L <- 2.5e6 # J kg^-1
+    Rv <- 461 # J kg^-1 K^-1
+    es_hPa <- 6.11 * exp((L / Rv) * (1 / 273 - 1 / (273 + Tc)))
+    return(units::ud_convert(es_hPa, "hPa", out_units))
+  }
 
-## TODO: merge SatVapPress with get.es; add option to choose method
-SatVapPres <- function(T) {
-  # /estimates saturation vapor pressure (kPa) Goff-Gratch 1946 /input: T = absolute temperature
-  T_st <- 373.15  ##steam temperature (K)
-  e_st <- 1013.25  ##/saturation vapor pressure at steam temp (hPa)
-  return(0.1 * exp(-7.90298 * (T_st/T - 1) + 5.02808 * log(T_st/T) - 1.3816e-07 * (10^(11.344 * (1 - T/T_st)) -
-    1) + 0.0081328 * (10^(-3.49149 * (T_st/T - 1)) - 1) + log(e_st)))
-} # SatVapPres
+  if (method == "GoffGratch") {
+    # canonical temp: K; canonical pressure: hPa
+    Tk <- units::ud_convert(temp, temp_units, "K")
+    Tst <- 373.15 # K
+    est <- 1013.246 # hPa at steam point
+    lg10 <- function(z) log10(z)
+    log10_es <- -7.90298 * (Tst / Tk - 1) +
+      5.02808 * lg10(Tst / Tk) -
+      1.3816e-7 * (10^(11.344 * (1 - Tk / Tst)) - 1) +
+      8.1328e-3 * (10^(-3.49149 * (Tst / Tk - 1)) - 1) +
+      lg10(est)
+    es_hPa <- 10^log10_es
+    return(units::ud_convert(es_hPa, "hPa", out_units))
+  }
+}
 
+# ---- Aliases for backward-compatibility ----
+
+#' @rdname sat_vapor_pressure
+#' @export
+get.es <- function(temp) {
+  sat_vapor_pressure(
+    temp = temp,
+    method = "ClausiusClapeyron",
+    temp_units = "degC",
+    out_units = "hPa"
+  )
+}
 
 ##' Calculate RH from temperature and dewpoint
 ##'
@@ -201,7 +263,7 @@ sw2ppfd <- function(sw) {
 ##' Campbell and Norman (1998). Introduction to Environmental Biophysics. pg 151 'the energy content of solar radiation in the PAR waveband is 2.35 x 10^5 J/mol'
 ##' See also the chapter radiation basics (10)
 ##' Here the input is the total solar radiation so to obtain in the PAR spectrum need to multiply by 0.486
-##' This last value 0.486 is based on the approximation that PAR is 0.45-0.50 of the total radiation
+##' This is based on the approximation that PAR is 0.45-0.50 of the total radiation
 ##' This means that 1e6 / (2.35e6) * 0.486 = 2.07
 ##' 1e6 converts from mol to mu mol
 ##' 1/3600 divides the values in hours to seconds

modules/data.atmosphere/tests/testthat/test.metutils.R

@@ -31,4 +31,52 @@ test_that("get.rh RH from dewpoint",{
   expect_equal(getrhtest(25, 10), 38.82, tolerance = 0.2)
   expect_equal(getrhtest(0, -5), 69, tolerance = 0.2)
 })
+
+test_that("different methods of sat_vapor_pressure work correctly", {
+  expect_equal(
+    sat_vapor_pressure(c(-10, 10), method = "Magnus"),
+    c(0.286, 1.228),
+    tolerance = 0.001
+  )
+  expect_equal(
+    sat_vapor_pressure(c(-10, 10), method = "ClausiusClapeyron"),
+    c(0.287, 1.233), 
+    tolerance = 0.001
+  )
+  expect_equal(
+    sat_vapor_pressure(c(-10, 10), method = "GoffGratch"),
+    c(0.286, 1.227),
+    tolerance = 0.001
+  )
+})
+
+test_that("sat_vapor_pressure works with different units", {
+  expect_equal(
+    sat_vapor_pressure(283.15,
+      method = "GoffGratch",
+      temp_units = "K",
+      out_units = "mb"
+    ),
+    12.27,
+    tolerance = 0.01
+  )
 
+    expect_equal(
+      sat_vapor_pressure(283.15,
+        method = "ClausiusClapeyron",
+        temp_units = "K",
+        out_units = "kPa"
+      ),
+      1.227,
+      tolerance = 0.01
+    )
+    expect_equal(
+      sat_vapor_pressure(283.15,
+        method = "Magnus",
+        temp_units = "K",
+        out_units = "Pa"
+      ),
+      1227,
+      tolerance = 1
+    )
+})