EFSA PEC soil guidance from 2017

- Implement the new guidance as well as possible
- Maintenance work addressing CRAN checks

4 files changed, 192 insertions, 50 deletions

diff --git a/R/PEC_soil.R b/R/PEC_soil.R
index 15ccc90..1227f6a 100644
--- a/R/PEC_soil.R
+++ b/R/PEC_soil.R
@@ -17,23 +17,34 @@
 
 # Register global variables
 if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_type", "TP_identifier",
-                                                       "soil_scenario_data_EFSA_2015"))
+                                                       "soil_scenario_data_EFSA_2015",
+                                                       "soil_scenario_data_EFSA_2017", "bottom"))
 
 #' Calculate predicted environmental concentrations in soil
 #'
 #' This is a basic calculation of a contaminant concentration in bulk soil
-#' based on complete, instantaneous mixing. If an interval is given, an 
+#' based on complete, instantaneous mixing. If an interval is given, an
 #' attempt is made at calculating a long term maximum concentration using
-#' the concepts layed out for example in the PPR panel opinion (EFSA 2012).
-#' 
+#' the concepts layed out in the PPR panel opinion (EFSA PPR panel 2012
+#' and in the EFSA guidance on PEC soil calculations (EFSA, 2015, 2017).
+#'
 #' This assumes that the complete load to soil during the time specified by
 #' 'interval' (typically 365 days) is dosed at once. As in the PPR panel
-#' opinion cited below (PPR panel 2012), only temperature correction using the
-#' Arrhenius equation is performed.  
-#' 
+#' opinion cited below (EFSA PPR panel 2012), only temperature correction using the
+#' Arrhenius equation is performed.
+#'
 #' Total soil and porewater PEC values for the scenarios as defined in the EFSA
-#' guidance (2015, p. 13) can easily be calculated.
-#' 
+#' guidance (2017, p. 14/15) can easily be calculated.
+#' @note While time weighted average (TWA) concentrations given in the examples
+#' from the EFSA guidance from 2015 (p. 80) are be reproduced, this is not
+#' true for the TWA concentrations given for the same example in the EFSA guidance
+#' from 2017 (p. 92).
+#' @note According to the EFSA guidance (EFSA, 2017, p. 43), leaching should be
+#'   taken into account for the EFSA 2017 scenarios, using the evaluation depth
+#'   (here mixing depth) as the depth of the layer from which leaching takes
+#'   place.  However, as the amount leaching below the evaluation depth
+#'   (often 5 cm) will partly be mixed back during tillage, the default in this function 
+#'   is to use the tillage depth for the calculation of the leaching rate.
 #' @note If temperature information is available in the selected scenarios, as
 #'   e.g. in the EFSA scenarios, the DT50 for groundwater modelling
 #'   (destination 'PECgw') is taken from the chent object, otherwise the DT50
@@ -49,7 +60,17 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ
 #' @param PEC_units Requested units for the calculated PEC. Only mg/kg currently supported
 #' @param PEC_pw_units Only mg/L currently supported
 #' @param tillage_depth Periodic (see interval) deeper mixing in cm
-#' @param chent An optional chent object holding substance specific information. Can 
+#' @param leaching_depth EFSA (2017) uses the mixing depth (ecotoxicological
+#'   evaluation depth) to calculate leaching for annual crops where tillage
+#'   takes place. By default, losses from the layer down to the tillage 
+#'   depth are taken into account in this implementation.
+#' @param cultivation Does mechanical cultivation in the sense of EFSA (2017)
+#'   take place, i.e. twice a year to a depth of 5 cm? Ignored for scenarios
+#'   other than EFSA_2017
+#' @param crop Ignored for scenarios other than EFSA_2017. Only annual crops
+#' are supported when these scenarios are used. Only crops with a single cropping
+#' cycle per year are currently supported.
+#' @param chent An optional chent object holding substance specific information. Can
 #'   also be a name for the substance as a character string
 #' @param DT50 If specified, overrides soil DT50 endpoints from a chent object
 #'   If DT50 is not specified here and not available from the chent object, zero
@@ -65,6 +86,8 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ
 #'   'chents' is specified, the DegT50 with destination 'PECgw' will be used),
 #'   and corrected using an Arrhenius activation energy of 65.4 kJ/mol. Also
 #'   model and scenario adjustment factors from the EFSA guidance are used.
+#' @param leaching Should leaching be taken into account? The default is FALSE,
+#'   except when the EFSA_2017 scenarios are used.
 #' @param porewater Should equilibrium porewater concentrations be estimated
 #'   based on Kom and the organic carbon fraction of the soil instead of total
 #'   soil concentrations?  Based on equation (7) given in the PPR panel opinion
@@ -76,23 +99,38 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ
 #'   selection and scenario parameterisation for predicting environmental
 #'   concentrations of plant protection products in soil. \emph{EFSA Journal}
 #'   \bold{10}(2) 2562, doi:10.2903/j.efsa.2012.2562
-#' 
+#'
+#'   EFSA (European Food Safety Authority) 2017) EFSA guidance document for
+#'   predicting environmental concentrations of active substances of plant
+#'   protection products and transformation products of these active substances
+#'   in soil. \emph{EFSA Journal} \bold{15}(10) 4982
+#'   doi:10.2903/j.efsa.2017.4982
+#'
 #'   EFSA (European Food Safety Authority) (2015) EFSA guidance document for
 #'   predicting environmental concentrations of active substances of plant
 #'   protection products and transformation products of these active substances
 #'   in soil. \emph{EFSA Journal} \bold{13}(4) 4093
 #'   doi:10.2903/j.efsa.2015.4093
+#'
 #' @author Johannes Ranke
 #' @export
 #' @examples
 #' PEC_soil(100, interception = 0.25)
 #' 
+#' # This is example 1 starting at p. 92 of the EFSA guidance (2017)
+#' # Note that TWA concentrations differ from the ones given in the guidance
+#' # for an unknown reason (the values from EFSA (2015) can be reproduced).
+#' PEC_soil(1000, interval = 365, DT50 = 250, t_avg = c(0, 21),
+#'                Kom = 1000, scenarios = "EFSA_2017")
+#' PEC_soil(1000, interval = 365, DT50 = 250, t_av = c(0, 21),
+#'                Kom = 1000, scenarios = "EFSA_2017", porewater = TRUE)
+#'
 #' # This is example 1 starting at p. 79 of the EFSA guidance (2015)
 #' PEC_soil(1000, interval = 365, DT50 = 250, t_avg = c(0, 21),
 #'                scenarios = "EFSA_2015")
 #' PEC_soil(1000, interval = 365, DT50 = 250, t_av = c(0, 21),
 #'                Kom = 1000, scenarios = "EFSA_2015", porewater = TRUE)
-#' 
+#'
 #' # The following is from example 4 starting at p. 85 of the EFSA guidance (2015)
 #' # Metabolite M2
 #' # Calculate total and porewater soil concentrations for tier 1 scenarios
@@ -103,34 +141,97 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ
 #'                            Kom = 100, scenarios = "EFSA_2015", porewater = TRUE)
 
 PEC_soil <- function(rate, rate_units = "g/ha", interception = 0,
-                     mixing_depth = 5, 
+                     mixing_depth = 5,
                      PEC_units = "mg/kg", PEC_pw_units = "mg/L",
                      interval = NA, n_periods = Inf,
-                     tillage_depth = 20, 
-                     chent = NA, 
-                     DT50 = NA, 
+                     tillage_depth = 20,
+                     leaching_depth = tillage_depth,
+                     crop = "annual",
+                     cultivation = FALSE,
+                     chent = NA,
+                     DT50 = NA,
                      Koc = NA, Kom = Koc / 1.724,
                      t_avg = 0,
-                     scenarios = c("default", "EFSA_2015"),
+                     scenarios = c("default", "EFSA_2017", "EFSA_2015"),
+                     leaching = scenarios == "EFSA_2017",
                      porewater = FALSE)
 {
+  # Comments with equation numbers in parentheses refer to
+  # the numbering in the EFSA guidance from 2017, appendix A
   rate_to_soil = (1 - interception) * rate
   rate_units = match.arg(rate_units)
   PEC_units = match.arg(PEC_units)
   scenarios = match.arg(scenarios)
-  sce <- switch(scenarios, 
+  if (scenarios == "EFSA_2017") {
+    if (crop != "annual") stop("Only annual crops are currently supported")
+    if (cultivation) stop("Permanent crops with mechanical cultivation are currently not supported")
+  }
+  sce <- switch(scenarios,
     default = data.frame(rho = 1.5, T_arr = NA, theta_fc = 0.2, f_om = 1.724 * 0.02,
                          f_sce = 1, f_mod = 1, row.names = "default"),
-    EFSA_2015 = if (porewater) soil_scenario_data_EFSA_2015[4:6, ] 
-                else soil_scenario_data_EFSA_2015[1:3, ]
+    EFSA_2015 = if (porewater) soil_scenario_data_EFSA_2015[4:6, ]
+                else soil_scenario_data_EFSA_2015[1:3, ],
+    EFSA_2017 = if (porewater) soil_scenario_data_EFSA_2017[4:6, ]
+                else soil_scenario_data_EFSA_2017[1:3, ]
   )
   n_sce = nrow(sce)
 
   soil_volume = 100 * 100 * (mixing_depth/100)   # in m3
   soil_mass = soil_volume * sce$rho * 1000  # in kg
 
+  # In EFSA (2017), f_om is depth dependent for permanent crops
+  # For annual crops, the correction factor is 1 (uniform f_om is
+  # assumed)
+  mixing_depth_string <- paste(mixing_depth, "cm")
+  tillage_depth_string <- paste(tillage_depth, "cm")
+  if (scenarios == "EFSA_2017" & crop != "annual") {
+    # Correction factors f_f_om with depth according to EFSA 2017, p. 15
+    f_f_om_depth = data.frame(
+      depth = c("0-5", "5-10", "10-20", "20-30"),
+      bottom = c(5, 10, 20, 30),
+      thickness = c(5, 5, 10, 10),
+      f_f_om_no_cultivation = c(1.95, 1.30, 0.76, 0.62),
+      f_f_om_cultivation = c(1.50, 1.20, 0.90, 0.75))
+    # Averages for the 0-5 cm and 0-20 cm layers
+    f_f_om_layer = data.frame(
+      layer = c("0-5", "0-20"),
+      f_f_om_no_cultivation = c(1.95, (5 * 1.95 + 5 * 1.3 + 10 * 0.76)/20),
+      f_f_om_cultivation = c(1.50, (5 * 1.5 + 5 * 1.2 + 10 * 0.9)/20))
+    # The resulting mean value for 0-20 cm and no cultivation of 1.1925 is
+    # consistent with the value of 1.19 given in Table B.4 on p. 54 of the
+    # 2017 EFSA guidance
+
+    f_f_om_average <- function(depth, cultivation) {
+      rownames(f_f_om_layer) = paste(f_f_om_layer$layer, "cm")
+      if (depth %in% c(5, 20)) {
+        if (cultivation) {
+          return(f_f_om_layer[paste0("0-", depth, " cm"), "f_f_om_cultivation"])
+        } else {
+          return(f_f_om_layer[paste0("0-", depth, " cm"), "f_f_om_no_cultivation"])
+        }
+      } else {
+        stop("Depths other than 5 and 20 cm are not supported when using EFSA 2017 scenarios for permanent crops")
+      }
+    } 
+
+    # For the loss via leaching, the equilibrium and therefore the f_om at the
+    # bottom of the layer is probably most relevant. Unfortunately this is not
+    # clarified in the guidance.
+    f_f_om_bottom <- function(depth, cultivation) {
+      bottom_depth <- depth # rename to avoid confusion when subsetting
+      if (cultivation) {
+        f_f_om <- subset(f_f_om_depth, bottom == bottom_depth)$f_f_om_cultivation
+      } else {
+        f_f_om <- subset(f_f_om_depth, bottom == bottom_depth)$f_f_om_no_cultivation
+      }
+      return(f_f_om)
+    } 
+  } else {
+    f_f_om_average <- f_f_om_bottom <- function(depth, cultivation) 1
+  }
+
   # The following is C_T,ini from EFSA 2012, p. 22, but potentially with interception > 0
-  PEC_soil_ini = rate_to_soil * 1000 / soil_mass     # in mg/kg
+  PEC_soil_ini = rate_to_soil * 1000 / soil_mass     # in mg/kg (A1)
 
   # Decide which DT50 to take, or set degradation to zero if no DT50 available
   if (is.na(DT50) & is(chent, "chent")) {
@@ -142,26 +243,56 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0,
     if (length(DT50) > 1) stop("More than one PECsoil DT50 in chent object")
     if (length(DT50) == 0) DT50 <- Inf
   }
-  k = log(2)/DT50
+  k_ref = log(2)/DT50 # (A5)
 
   # Temperature correction of degradation (accumulation)
   if (all(is.na(sce$T_arr))) {   # No temperature correction
     f_T = 1
   } else {
     # Temperature correction as in EFSA 2012 p. 23
-    f_T = ifelse(sce$T_arr == 0, 
-                 0, 
-                 exp(- (65.4 / 0.008314) * (1/(sce$T_arr + 273.15) - 1/293.15)))
+    f_T = ifelse(sce$T_arr == 0,
+                 0, # (A4b)
+                 exp(- (65.4 / 0.008314) * (1/(sce$T_arr + 273.15) - 1/293.15))) # (A4a)
   }
 
-  # X is the fraction left after one period (EFSA guidance p. 23)
-  X = exp(- k * f_T * interval)
-  
+  # Define Kom if needed
+  if (leaching | porewater) {
+    # If Kom is not specified, try to get K(f)oc
+    if (is.na(Kom)) {
+      # If Koc not specified, try to get K(f)oc from chent
+      if (is.na(Koc) & is(chent, "chent")) {
+        Koc <- soil_Kfoc(chent)
+      } else {
+        stop("No Kom information specified")
+      }
+      Kom <- Koc / 1.724
+    }
+  }
+
+  if (leaching) {
+    leaching_depth_string <- paste(leaching_depth, "cm")
+    f_q <- c("1 cm" = 0.8, "2.5 cm" = 0.75, "5 cm" = 0.7, "20 cm" = 0.5) # EFSA 2017 p. 54
+    if (leaching_depth_string %in% names(f_q)) {
+      q_mm_year = f_q[leaching_depth_string] * sce$prec # Irrigation at tier 1? I have not found values for Tier 1
+      q_dm_day = q_mm_year / (100 * 365)
+      leaching_depth_dm <- leaching_depth / 10
+
+      k_leach = q_dm_day/(leaching_depth_dm * (sce$theta_fc + sce$rho * f_f_om_average(leaching_depth, cultivation) * sce$f_om * Kom))
+    } else {
+      stop("Leaching can not be calculated, because f_q for this leaching depth is undefined")
+    }
+  } else {
+    k_leach = 0
+  }
+
+  # X is the fraction left after one period (EFSA 2017 guidance p. 23)
+  X = exp(- (k_ref * f_T + k_leach) * interval) # (A3)
+
   # f_accu is the fraction left after n periods (X + X^2 + ...)
-  f_accu = 0 
+  f_accu = 0
   if (!is.na(interval)) {
     if (n_periods == Inf) {
-      f_accu = X/(1 - X)
+      f_accu = X/(1 - X) # part of (A2)
     } else {
       for (i in 1:n_periods) {
         f_accu = f_accu + X^i
@@ -171,25 +302,14 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0,
 
   f_tillage = mixing_depth / tillage_depth
 
-  PEC_background = f_accu * f_tillage * PEC_soil_ini
+  PEC_background = f_accu * f_tillage * PEC_soil_ini # (A2)
 
-  PEC_soil = (1 + f_accu * f_tillage) * PEC_soil_ini
+  PEC_soil = PEC_soil_ini + PEC_background # (A6)
 
   # Get porewater PEC if requested
   if (porewater) {
 
-    # If Kom is not specified, try to get K(f)oc
-    if (is.na(Kom)) {
-      # If Koc not specified, try to get K(f)oc from chent
-      if (is.na(Koc) & is(chent, "chent")) {
-        Koc <- soil_Kfoc(chent)
-      } 
-      Kom <- Koc / 1.724
-    }
-
-    if (is.na(Kom)) stop("No Kom information specified")
-
-    PEC_soil = PEC_soil/((sce$theta_fc/sce$rho) + sce$f_om * Kom)
+    PEC_soil = PEC_soil/((sce$theta_fc/sce$rho) + f_f_om_average(mixing_depth, cultivation) * sce$f_om * Kom) # (A7)
   }
 
   # Scenario adjustment factors
@@ -200,14 +320,15 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0,
 
   result <- matrix(NA, ncol = n_sce, nrow = length(t_avg),
                    dimnames = list(t_avg = t_avg, scenario = rownames(sce)))
-  
+
   result[1, ] <- PEC_soil_sce_mod
 
   for (i in seq_along(t_avg)) {
     t_av_i <- t_avg[i]
+    k_avg <- f_T * k_ref # Leaching not taken into account, EFSA 2017 p. 43
     if (t_av_i > 0) {
       # Equation 10 from p. 24 (EFSA 2015)
-      result[i, ] <- PEC_soil_sce_mod/(t_av_i * f_T * k) * (1 - exp(- f_T * k * t_av_i))
+      result[i, ] <- PEC_soil_sce_mod/(t_av_i * k_avg) * (1 - exp(- k_avg * t_av_i)) # (A8)
     }
   }
 
diff --git a/R/PEC_sw_exposit_runoff.R b/R/PEC_sw_exposit_runoff.R
index f56b8f8..733f621 100644
--- a/R/PEC_sw_exposit_runoff.R
+++ b/R/PEC_sw_exposit_runoff.R
@@ -6,6 +6,7 @@
 #' @format A data frame with percentage values for the dissolved fraction and the fraction
 #'   bound to eroding particles, with Koc classes used as row names
 #'   \describe{
+#'     \item{Koc_lower_bound}{The lower bound of the Koc class}
 #'     \item{dissolved}{The percentage of the applied substance transferred to an
 #'       adjacent water body in the dissolved phase}
 #'     \item{bound}{The percentage of the applied substance transferred to an
diff --git a/R/soil_scenario_data_EFSA_2015.R b/R/soil_scenario_data_EFSA_2015.R
index fb096bc..660cafe 100644
--- a/R/soil_scenario_data_EFSA_2015.R
+++ b/R/soil_scenario_data_EFSA_2015.R
@@ -10,10 +10,10 @@
 #'   e.g. CTN for the Northern scenario for the total concentration in soil. Columns are 
 #'   mostly self-explanatory. \code{rho} is the dry bulk density of the top soil.
 #' @source EFSA (European Food Safety Authority) (2015)
-#' EFSA guidance document for predicting environmental concentrations
-#' of active substances of plant protection products and transformation products of these
-#' active substances in soil. \emph{EFSA Journal} \bold{13}(4) 4093
-#' doi:10.2903/j.efsa.2015.4093
+#'   EFSA guidance document for predicting environmental concentrations
+#'   of active substances of plant protection products and transformation products of these
+#'   active substances in soil. \emph{EFSA Journal} \bold{13}(4) 4093
+#'   doi:10.2903/j.efsa.2015.4093
 #' @keywords datasets
 #' @examples
 #' \dontrun{
diff --git a/R/soil_scenario_data_EFSA_2017.R b/R/soil_scenario_data_EFSA_2017.R
new file mode 100644
index 0000000..79ee15f
--- /dev/null
+++ b/R/soil_scenario_data_EFSA_2017.R
@@ -0,0 +1,20 @@
+#' Properties of the predefined scenarios from the EFSA guidance from 2017
+#' 
+#' Properties of the predefined scenarios used at Tier 1, Tier 2A and Tier 3A for the 
+#' concentration in soil as given in the EFSA guidance (2017, p. 14/15). Also, the 
+#' scenario and model adjustment factors from p. 16 and p. 18 are included.
+#' 
+#' @name soil_scenario_data_EFSA_2017
+#' @docType data
+#' @format A data frame with one row for each scenario. Row names are the scenario codes, 
+#'   e.g. CTN for the Northern scenario for the total concentration in soil. Columns are 
+#'   mostly self-explanatory. \code{rho} is the dry bulk density of the top soil.
+#' @source EFSA (European Food Safety Authority) (2017)
+#'   EFSA guidance document for predicting environmental concentrations
+#'   of active substances of plant protection products and transformation products of these
+#'   active substances in soil. \emph{EFSA Journal} \bold{15}(10) 4982
+#'   doi:10.2903/j.efsa.2017.4982
+#' @keywords datasets
+#' @examples
+#' soil_scenario_data_EFSA_2017
+NULL