added CWD demo and extreme value distribution fitting

NAMESPACE

@@ -1,4 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(convert_et)
 export(cwd)
+export(simulate_snow)
 import(dplyr)

R/convert_et.R

@@ -0,0 +1,214 @@
+#' Water energy to mass conversion
+#'
+#' Convert the latent heat flux into evapotranspiration in (water) mass units.
+#' Based on the SPLASH model (Davis et al. 2017 http://www.geosci-model-dev.net/10/689/2017/)
+#'
+#' @param et_e Latent heat flux in energy units (W m-2; representative for the
+#' mean value over one time step)
+#' @param tc Air temperature in degrees Celsius
+#' @param patm Atmospheric pressure (Pa)
+#' @param return_df A logical specifying whether to return a data frame with a
+#' single column containing the water mass flux. Defaults to \code{FALSE}
+#'
+#' @details Returns the water mass flux in mm d-1 (equivalent to kg m-2)
+#'
+#' @export
+#'
+convert_et <- function(et_e, tc, patm, return_df = FALSE){
+
+  par_splash <- list(
+		kTkelvin = 273.15,  # freezing point in K (= 0 deg C)
+		kTo = 298.15,       # base temperature, K (from P-model)
+		kR  = 8.31446262,   # universal gas constant, J/mol/K (Allen, 1973)
+		kMv = 18.02,        # molecular weight of water vapor, g/mol (Tsilingiris, 2008)
+		kMa = 28.963,       # molecular weight of dry air, g/mol (Tsilingiris, 2008) XXX this was in SPLASH (WITH 1E-3 IN EQUATION) XXX
+		kfFEC = 2.04,       # from flux to energy conversion, umol/J (Meek et al., 1984)
+		kPo = 101325,       # standard atmosphere, Pa (Allen, 1973)
+		kL  = 0.0065,       # temperature lapse rate, K/m (Cavcar, 2000)
+		kG  = 9.80665,      # gravitational acceleration, m/s^2 (Allen, 1973)
+		k_karman = 0.41,    # Von Karman constant; from bigleaf R package
+		eps = 9.999e-6,     # numerical imprecision allowed in mass conservation tests
+		cp = 1004.834,      # specific heat of air for constant pressure (J K-1 kg-1); from bigleaf R package
+		Rd = 287.0586       # gas constant of dry air (J kg-1 K-1) (Foken 2008 p. 245; from bigleaf R package)
+		)
+
+  sat_slope <- calc_sat_slope(tc)
+  lv <- calc_enthalpy_vap(tc)
+  pw <- calc_density_h2o(tc, patm)
+  gamma <- calc_psychro(tc, patm, par_splash)
+  econ <- sat_slope / (lv * pw * (sat_slope + gamma))
+
+  ## J m-2 d-1 -> mm / d
+  et_mm <- et_e * econ * 60 * 60 * 24 * 1000
+
+  if (return_df){
+    return(tibble(et_mm = et_mm))
+  } else {
+    return(et_mm)
+  }
+
+}
+
+
+calc_patm <- function( elv, par ){
+  #----------------------------------------------------------------
+  # Calculates atmospheric pressure for a given elevation, assuming
+  # standard atmosphere at sea level (kPo)
+  # Ref:      Allen et al. (1998)
+  # This function is copied from SPLASH
+  #----------------------------------------------------------------
+  # use md_params_core, only: kPo, kL, kTo, kG, kMa, kR
+
+  # # arguments
+  # real, intent(in) :: elv # elevation above sea level, m
+
+  # # function return value
+  # real ::  patm ! atmospheric pressure (Pa)
+
+  patm <- par$kPo * (1.0 - par$kL * elv / par$kTo) ^ (par$kG * par$kMa * 1.e-3 / (par$kR * par$kL))
+
+  return(patm)
+
+}
+
+
+calc_sat_slope <- function( tc ){
+  #----------------------------------------------------------------
+  # Calculates the slope of the sat pressure temp curve, Pa/K
+  # Ref:      Eq. 13, Allen et al. (1998)
+  #----------------------------------------------------------------
+  # # arguments
+  # real, intent(in) :: tc # air temperature, degrees C
+
+  # # function return value
+  # real :: sat_slope  # slope of the sat pressure temp curve, Pa/K
+
+  sat_slope <- (17.269)*(237.3)*(610.78)*(exp(tc*17.269/(tc + 237.3))/((tc + 237.3)^2))
+
+	return( sat_slope )
+
+}
+
+
+calc_enthalpy_vap <- function( tc ){
+  #----------------------------------------------------------------
+  # Calculates the enthalpy of vaporization, J/kg
+  # Ref:      Eq. 8, Henderson-Sellers (1984)
+  #----------------------------------------------------------------
+  # # arguments
+  # real, intent(in) :: tc # air temperature, degrees C
+
+  # # function return value
+  # real ::  enthalpy_vap # enthalpy of vaporization, J/kg
+
+  enthalpy_vap <- 1.91846e6*((tc + 273.15)/(tc + 273.15 - 33.91))^2
+
+  return( enthalpy_vap )
+
+}
+
+
+calc_density_h2o <- function( tc, press ){
+  #----------------------------------------------------------------
+  # Calculates density of water at a given temperature and pressure
+  # Ref: Chen et al. (1977)
+  #----------------------------------------------------------------
+  # # arguments
+  # real, intent(in) :: tc     # air temperature (degrees C)
+  # real, intent(in) :: press  # atmospheric pressure (Pa)
+
+  # # local variables
+  # real :: po, ko, ca, cb
+  # real :: pbar               # atmospheric pressure (bar)
+
+  # # function return value
+  # real :: density_h2o  # density of water, kg/m^3
+
+  # Calculate density at 1 atm:
+  po <- 0.99983952
+	  + 6.788260e-5  *tc
+	  - 9.08659e-6   *tc*tc
+	  + 1.022130e-7  *tc*tc*tc
+	  - 1.35439e-9   *tc*tc*tc*tc
+	  + 1.471150e-11 *tc*tc*tc*tc*tc
+	  - 1.11663e-13  *tc*tc*tc*tc*tc*tc
+	  + 5.044070e-16 *tc*tc*tc*tc*tc*tc*tc
+	  - 1.00659e-18  *tc*tc*tc*tc*tc*tc*tc*tc
+
+  # Calculate bulk modulus at 1 atm:
+  ko <- 19652.17
+	  + 148.1830   *tc
+	  - 2.29995    *tc*tc
+	  + 0.01281    *tc*tc*tc
+	  - 4.91564e-5 *tc*tc*tc*tc
+	  + 1.035530e-7*tc*tc*tc*tc*tc
+
+  # Calculate temperature dependent coefficients:
+  ca <- 3.26138
+	  + 5.223e-4  *tc
+	  + 1.324e-4  *tc*tc
+	  - 7.655e-7  *tc*tc*tc
+	  + 8.584e-10 *tc*tc*tc*tc
+
+  cb <- 7.2061e-5
+	  - 5.8948e-6  *tc
+	  + 8.69900e-8 *tc*tc
+	  - 1.0100e-9  *tc*tc*tc
+	  + 4.3220e-12 *tc*tc*tc*tc
+
+  # Convert atmospheric pressure to bar (1 bar <- 100000 Pa)
+  pbar <- (1.0e-5)*press
+
+  density_h2o <- 1000.0*po*(ko + ca*pbar + cb*pbar^2.0)/(ko + ca*pbar + cb*pbar^2.0 - pbar)
+
+  return(density_h2o)
+
+}
+
+
+calc_psychro <- function( tc, press, par_splash ){
+  #----------------------------------------------------------------
+  # Calculates the psychrometric constant for a given temperature and pressure
+  # Ref: Allen et al. (1998); Tsilingiris (2008)
+  #----------------------------------------------------------------
+  # # arguments
+  # real, intent(in) :: tc     # air temperature, degrees C
+  # real, intent(in) :: press  # atmospheric pressure, Pa
+
+  # # local variables
+  # real :: lv  # latent heat of vaporization (J/kg)
+  # real :: cp
+
+  # # function return value
+  # real :: psychro  # psychrometric constant, Pa/K
+
+  # # local variables
+  # real :: my_tc    # adjusted temperature to avoid numerical blow-up
+
+  # Adopted temperature adjustment from SPLASH, Python version
+  my_tc <- tc
+
+	my_tc <- ifelse(tc > 100, 100, ifelse(tc < 0, 0, tc))
+
+  # Calculate the specific heat capacity of water, J/kg/K
+  # Eq. 47, Tsilingiris (2008)
+  cp <- 1.0e3*(1.0045714270
+			    + 2.050632750e-3  *my_tc
+				  - 1.631537093e-4  *my_tc*my_tc
+				  + 6.212300300e-6  *my_tc*my_tc*my_tc
+				  - 8.830478888e-8  *my_tc*my_tc*my_tc*my_tc
+				  + 5.071307038e-10 *my_tc*my_tc*my_tc*my_tc*my_tc
+          )
+
+  # Calculate latent heat of vaporization, J/kg
+  lv <- calc_enthalpy_vap(tc)
+
+  # Calculate psychrometric constant, Pa/K
+  # Eq. 8, Allen et al. (1998)
+  psychro <- cp * par_splash$kMa * press / (par_splash$kMv * lv)
+
+  return(psychro)
+
+}
+
+

R/cwd.R

@@ -35,7 +35,7 @@
 #'
 cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_drop = 0.9){
 
-  ## corresponds to mct2.R
+  # corresponds to mct2.R
 
   if (thresh_terminate > thresh_drop){
     rlang::warn("Aborting. thresh_terminate must be smaller or equal thresh_drop. Setting it equal.")
@@ -51,12 +51,12 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d
     # dplyr::select(date, !!varname_wbal) %>%
     mutate(iinst = NA, dday = NA, deficit = 0)
 
-  ## search all dates
+  # search all dates
   while (idx <= (nrow(df)-1)){
     idx <- idx + 1
 
-    ## if the water balance (deficit = prec - et) is negative, start accumulating deficit
-    ## cumulate negative water balances (deficits)
+    # if the water balance (deficit = prec - et) is negative, start accumulating deficit
+    # cumulative negative water balances (deficits)
     if (df[[ varname_wbal ]][idx]<0){
 
       dday <- 0
@@ -65,33 +65,24 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d
       iidx <- idx
       done_finding_dropday <- FALSE
 
-      # ## xxx debug
-      # newlargeinst <- TRUE
-
-      ## continue accumulating deficit as long as the deficit has not fallen below (thresh_terminate) times the maximum deficit attained in this event
+      # continue accumulating deficit as long as the deficit has not fallen below (thresh_terminate) times the maximum deficit attained in this event
       while (iidx <= (nrow(df)-1) && (deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit)){
         dday <- dday + 1
         deficit <- deficit - df[[ varname_wbal ]][iidx]
 
-        # ## xxx debug
-        # if (newlargeinst && deficit>50){
-        #   print("entering large deficit...")
-        #   newlargeinst <- FALSE
-        # }
-
-        ## record the maximum deficit attained in this event
+        # record the maximum deficit attained in this event
         if (deficit > max_deficit){
           max_deficit <- deficit
           done_finding_dropday <- FALSE
         }
 
-        ## record the day when deficit falls below (thresh_drop) times the current maximum deficit
+        # record the day when deficit falls below (thresh_drop) times the current maximum deficit
         if (deficit < (max_deficit * thresh_drop) && !done_finding_dropday){
           iidx_drop <- iidx
           done_finding_dropday <- TRUE
         }
 
-        ## once, deficit has fallen below threshold, all subsequent dates are dropped (dday set to NA)
+        # once, deficit has fallen below threshold, all subsequent dates are dropped (dday set to NA)
         if (done_finding_dropday){
           df$iinst[iidx] <- NA
           df$dday[iidx]  <- NA
@@ -101,7 +92,7 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d
           iidx_drop <- iidx
         }
 
-        # ## (even before "drop-day" is found), drop data of days after rain, i.e., where current CWD is below the maximum (previously) attained in the same event
+        # # (even before "drop-day" is found), drop data of days after rain, i.e., where current CWD is below the maximum (previously) attained in the same event
         # if (deficit < max_deficit){
         #   df$iinst[iidx] <- NA
         #   df$dday[iidx]  <- NA
@@ -118,16 +109,11 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d
 
       }
 
-      # ## xxx debug
-      # if (!newlargeinst){
-      #   print("crucial moment")
-      # }
-
-      ## record instance
+      # record instance
       this_inst <- tibble( idx_start = idx, len = iidx_drop-idx, iinst = iinst, date_start=df[[varname_date]][idx], date_end = df[[varname_date]][iidx_drop-1], deficit = max_deficit )
       inst <- inst %>% bind_rows(this_inst)
 
-      ## update
+      # update
       iinst <- iinst + 1
       dday <- 0
       idx <- iidx

R/simulate_snow.R

@@ -0,0 +1,63 @@
+#' Simulate snow mass
+#'
+#' Simulate snow mass accumulation and melt based on Orth et al. (2013)
+#' https://www.jstor.org/stable/24914341
+#'
+#' @param df A data frame containing columns for air temperature (deg. C),
+#' precipitation in liquid form (rain, mm d-1), and precipitation in solid form
+#' (snow water equivalents, mm d-1). The column names of the respective
+#' variables are provided by the other arguments.
+#' @param varnam_temp A character string specifying the variable name for air
+#' temperature.
+#' @param varnam_prec A character string specifying the variable name for rain.
+#' @param varnam_snow A character string specifying the variable name for snow.
+#'
+#' @details Returns a data frame with two added columns: (1) \code{liquid_to_soil}
+#' is the rain plus snow melt in mm d-1; (2) \code{snow_pool} is the snow mass
+#' in water equivalents (mm) for each day.
+#'
+#' @export
+#'
+simulate_snow <- function(df, varnam_temp, varnam_prec, varnam_snow){
+
+  temp <- df |> dplyr::pull(!!varnam_temp)
+  prec <- df |> dplyr::pull(!!varnam_prec)
+  snow <- df |> dplyr::pull(!!varnam_snow)
+
+  ## fixed parameters
+  temp_threshold <- 1.0
+  maxmeltrate <- 1.0
+
+  snow_pool <- 0
+  liquid_to_soil <- rep(NA, length(prec))
+  snow_pool_out <- rep(NA, length(prec))
+
+  ## spinup 1 year
+  for (doy in 1:365){
+    if ( snow_pool > 0.0 && temp[doy] > temp_threshold ){
+      melt  <- min( snow_pool, maxmeltrate * ( temp[doy] - temp_threshold ) )
+    } else {
+      melt <- 0.0
+    }
+    snow_pool <- snow_pool + snow[doy] - melt
+  }
+
+  ## transient forward
+  for (doy in 1:length(prec)){
+    if ( snow_pool > 0.0 && temp[doy] > temp_threshold ){
+      melt  <- min( snow_pool, maxmeltrate * ( temp[doy] - temp_threshold ) )
+    } else {
+      melt <- 0.0
+    }
+    snow_pool <- snow_pool + snow[doy] - melt
+    liquid_to_soil[doy] <- prec[doy] + melt
+    snow_pool_out[doy] <- snow_pool
+  }
+
+  ## complement
+  df$liquid_to_soil <- liquid_to_soil
+  df$snow_pool <- snow_pool_out
+
+  return(df)
+
+}

cummulative_water_deficit.Rproj

@@ -18,3 +18,4 @@ StripTrailingWhitespace: Yes
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace

data-raw/get_data.R

@@ -0,0 +1,11 @@
+library(readr)
+library(dplyr)
+library(here)
+library(lubridate)
+
+# get fluxnet data for FR-Pue site
+df <- read_csv(paste0(here(), "/inst/extdata/FLX_CH-Lae_FLUXDATAKIT_FULLSET_DD_2004_2014_2-3.csv")) |>
+  select(TIMESTAMP, P_F, TA_F_MDS, PA_F, LE_F_MDS) |>
+  filter(year(TIMESTAMP) > 2004)
+
+saveRDS(df, file = paste0(here(), "/data/df_ch-lae.rds"))