From 396ba6c42252426a132be56b9f417786a1f276e3 Mon Sep 17 00:00:00 2001
From: Ranke Johannes <johannes.ranke@agroscope.admin.ch>
Date: Tue, 13 Aug 2024 10:58:16 +0200
Subject: Support units in drift, runoff und sediment PECs

---
 R/PEC_sw_drift.R          |  27 ++++--
 R/PEC_sw_exposit.R        | 210 ++++++++++++++++++++++++++++++++++++++++++++++
 R/PEC_sw_exposit_runoff.R | 188 -----------------------------------------
 R/PEC_sw_sed.R            |  22 +++--
 4 files changed, 243 insertions(+), 204 deletions(-)
 create mode 100644 R/PEC_sw_exposit.R
 delete mode 100644 R/PEC_sw_exposit_runoff.R

(limited to 'R')

diff --git a/R/PEC_sw_drift.R b/R/PEC_sw_drift.R
index 22aa012..f503c74 100644
--- a/R/PEC_sw_drift.R
+++ b/R/PEC_sw_drift.R
@@ -4,9 +4,16 @@
 #' concentration in surface water based on complete, instantaneous mixing
 #' with input via spray drift.
 #'
+#' It is recommened to specify the arguments `rate`, `water_depth` and
+#' `water_width` using [units::units] from the `units` package.
+#'
 #' @inheritParams drift_percentages_rautmann
+#' @importFrom units as_units set_units
 #' @seealso [drift_parameters_focus], [drift_percentages_rautmann]
-#' @param rate Application rate in units specified below
+#' @param rate Application rate in units specified below, or with units defined via the
+#' `units` package.
+#' @param rate_units Defaults to g/ha. For backwards compatibility, only used
+#' if the specified rate does not have [units::units]].
 #' @param drift_percentages Percentage drift values for which to calculate PECsw.
 #'   Overrides 'drift_data' and 'distances' if not NULL.
 #' @param drift_data Source of drift percentage data. If 'JKI', the [drift_data_JKI]
@@ -16,7 +23,6 @@
 #' @param crop_group_JKI When using the 'JKI' drift data, one of the German names
 #'   as used in [drift_parameters_focus].
 #' @param water_depth Depth of the water body in cm
-#' @param rate_units Defaults to g/ha
 #' @param PEC_units Requested units for the calculated PEC. Only µg/L currently supported
 #' @param water_width Width of the water body in cm
 #' @param side_angle The angle of the side of the water relative to the bottom which
@@ -49,7 +55,7 @@
 #' PEC_sw_drift(100, drift_data = "RF", formula = "FOCUS", side_angle = 45, water_width = 200)
 PEC_sw_drift <- function(rate,
   applications = 1,
-  water_depth = 30,
+  water_depth = as_units("30 cm"),
   drift_percentages = NULL,
   drift_data = c("JKI", "RF"),
   crop_group_JKI = c("Ackerbau", 
@@ -59,22 +65,27 @@ PEC_sw_drift <- function(rate,
     "fruit, late", "fruit, early", "aerial"),
   distances = c(1, 5, 10, 20),
   formula = c("Rautmann", "FOCUS"),
-  water_width = 100,
+  water_width = as_units("100 cm"),
   side_angle = 90,
   rate_units = "g/ha",
   PEC_units = "\u00B5g/L")
 {
   rate_units <- match.arg(rate_units)
   PEC_units <- match.arg(PEC_units)
+  # Set default units if not specified
+  if (!inherits(rate, "units")) rate <- set_units(rate, rate_units, mode = "symbolic")
+  if (!inherits(water_width, "units")) water_width <- set_units(water_width, "cm")
+  if (!inherits(water_depth, "units")) water_depth <- set_units(water_depth, "cm")
   drift_data <- match.arg(drift_data)
   crop_group_JKI <- match.arg(crop_group_JKI)
   crop_group_focus <- match.arg(crop_group_focus)
   formula <- match.arg(formula)
   if (side_angle < 0 | side_angle > 90) stop("The side anglemust be between 0 and 90 degrees")
-  mean_water_width <- if (side_angle == 90) water_width
+  mean_water_width <- if (side_angle == 90) water_width # Mean water width over waterbody depth
   else water_width - (water_depth / tanpi(side_angle/180))
-  water_volume <- 100 * mean_water_width * (water_depth/100) * 1000   # in L (for 1 ha)
-  PEC_sw_overspray <- rate * 1e6 / water_volume          # in µg/L
+  if (as.numeric(mean_water_width) < 0) stop("Undefined geometry")
+  relative_mean_water_width <- mean_water_width / water_width # Always <= 1
+  PEC_sw_overspray <- set_units(rate / (relative_mean_water_width * water_depth), PEC_units, mode = "symbolic")
   dist_index <- as.character(distances)
 
   if (is.null(drift_percentages)) {
@@ -82,7 +93,7 @@ PEC_sw_drift <- function(rate,
       JKI = pfm::drift_data_JKI[[applications]][dist_index, crop_group_JKI],
       RF =  drift_percentages_rautmann(distances, applications, 
         formula = formula,
-        crop_group_focus, widths = water_width/100)
+        crop_group_focus, widths = as.numeric(set_units(water_width, "m")))
     )
     names(drift_percentages) <- paste(dist_index, "m")
   } else {
diff --git a/R/PEC_sw_exposit.R b/R/PEC_sw_exposit.R
new file mode 100644
index 0000000..a47ca18
--- /dev/null
+++ b/R/PEC_sw_exposit.R
@@ -0,0 +1,210 @@
+#' Runoff loss percentages as used in Exposit 3
+#'
+#' A table of the loss percentages used in Exposit 3 for the twelve different Koc classes
+#'
+#' @name perc_runoff_exposit
+#' @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
+#'       adjacent water body bound to eroding particles}
+#'   }
+#' @source Excel 3.02 spreadsheet available from
+#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
+#' @docType data
+#' @examples
+#' print(perc_runoff_exposit)
+"perc_runoff_exposit"
+
+#' Runoff reduction percentages as used in Exposit
+#'
+#' A table of the runoff reduction percentages used in Exposit 3 for different vegetated buffer widths
+#'
+#' @name perc_runoff_reduction_exposit
+#' @format A named list of data frames with reduction percentage values for the
+#' dissolved fraction and the fraction bound to eroding particles, with
+#' vegetated buffer widths as row names. The names of the list items are the Exposit versions
+#' from which the values were taken.
+#'   \describe{
+#'     \item{dissolved}{The reduction percentage for the dissolved phase}
+#'     \item{bound}{The reduction percentage for the particulate phase}
+#'   }
+#' @source Excel 3.02 spreadsheet available from
+#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
+#'
+#'   Agroscope version 3.01a with additional runoff factors for 3 m and 6 m buffer zones received from Muris Korkaric (not published).
+#'   The variant 3.01a2 was introduced for consistency with previous calculations performed by Agroscope for a 3 m buffer zone.
+#' @docType data
+#' @examples
+#' print(perc_runoff_reduction_exposit)
+"perc_runoff_reduction_exposit"
+
+#' Calculate PEC surface water due to runoff and erosion as in Exposit 3
+#'
+#' This is a reimplementation of the calculation described in the Exposit 3.02 spreadsheet file,
+#' in the worksheet "Konzept Runoff".
+#'
+#' It is recommened to specify the arguments `rate`, `Koc`, `DT50`, `t_runoff`, `V_ditch` and `V_event`
+#' using [units::units] from the `units` package.
+#'
+#' @importFrom units as_units set_units drop_units
+#' @importFrom dplyr across mutate
+#' @param rate The application rate in g/ha
+#' @param interception The fraction intercepted by the crop
+#' @param Koc The sorption coefficient to soil organic carbon
+#' @param DT50 The soil half-life in days
+#' @param t_runoff The time between application and the runoff event, where degradation occurs, in days
+#' @param exposit_reduction_version The version of the reduction factors to be used. "3.02" is the current
+#'   version used in Germany, "3.01a" is the version with additional percentages for 3 m and 6 m buffer
+#'   zones used in Switzerland. "3.01a2" is a version introduced for consistency with previous calculations
+#'   performed for a 3 m buffer zone in Switzerland, with the same reduction being applied to the dissolved
+#'   and the bound fraction.
+#' @param V_ditch The volume of the ditch is assumed to be 1 m * 100 m * 30 cm = 30 m3
+#' @param V_event The unreduced runoff volume, equivalent to 10 mm precipitation on 1 ha
+#' @param dilution The dilution factor
+#' @return A list containing the following components
+#'   \describe{
+#'     \item{perc_runoff}{The runoff percentages for dissolved and bound substance}
+#'     \item{runoff}{A matrix containing dissolved and bound input for the different distances}
+#'     \item{PEC_sw_runoff}{A matrix containing PEC values for dissolved and bound substance
+#'       for the different distances. If the rate was given in g/ha, the PECsw are in microg/L.}
+#'   }
+#' @export
+#' @source Excel 3.02 spreadsheet available from
+#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
+#' @seealso \code{\link{perc_runoff_exposit}} for runoff loss percentages and \code{\link{perc_runoff_reduction_exposit}} for runoff reduction percentages used
+#' @examples
+#'   PEC_sw_exposit_runoff(500, Koc = 150)
+#'   PEC_sw_exposit_runoff(600, Koc = 10000, DT50 = 195, exposit = "3.01a")
+PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, DT50 = set_units(Inf, "d"), 
+  t_runoff = set_units(3, "days"),
+  exposit_reduction_version = c("3.02", "3.01a", "3.01a2", "2.0"),
+  V_ditch = set_units(30, "m3"), V_event = set_units(100, "m3"), dilution = 2)
+{
+  # Set default units if not specified
+  if (!inherits(rate, "units")) rate <- set_units(rate, "g/ha")
+  if (!inherits(Koc, "units")) Koc <- set_units(Koc, "L/kg")
+  if (!inherits(DT50, "units")) DT50 <- set_units(DT50, "d")
+  if (!inherits(t_runoff, "units")) t_runoff <- set_units(t_runoff, "d")
+  if (!inherits(V_ditch, "units")) V_ditch <- set_units(V_ditch, "m3")
+  if (!inherits(V_event, "units")) V_event <- set_units(V_event, "m3")
+
+  k_deg <- log(2)/DT50
+  
+  # The input is calculated for an area of 1 ha
+  input <- rate * as_units(1, "ha") * (1 - interception) * exp(as.numeric(-k_deg * t_runoff))
+  input_units <- units(input)
+  input_numeric <- drop_units(input)
+
+  if (length(Koc) > 1) stop("Only one compound at a time supported")
+
+  exposit_reduction_version <- match.arg(exposit_reduction_version)
+  reduction_runoff <- pfm::perc_runoff_reduction_exposit[[exposit_reduction_version]] / 100
+        
+  transfer_runoff <- 1 - reduction_runoff
+
+  V_runoff <- V_event * (1 - reduction_runoff[["dissolved"]])
+  V_flowing_ditch_runoff <- dilution * (V_ditch + V_runoff)
+  
+  f_runoff_exposit <- function(Koc) {
+    Koc_breaks <- c(pfm::perc_runoff_exposit$Koc_lower_bound, set_units(Inf, "L/kg"))
+    Koc_classes <- as.character(cut(Koc, Koc_breaks, labels = rownames(pfm::perc_runoff_exposit)))
+    perc_runoff <- pfm::perc_runoff_exposit[Koc_classes, c("dissolved", "bound")]
+    if (identical(Koc, 0)) perc_runoff <- c(dissolved = 0, bound = 0)
+    return(unlist(perc_runoff) / 100)
+  }
+  f_runoff <- f_runoff_exposit(Koc)
+
+  runoff_dissolved <- input_numeric * f_runoff["dissolved"] * transfer_runoff[, "dissolved"]
+  runoff_bound <- input_numeric * f_runoff["bound"] * transfer_runoff[, "bound"]
+  runoff_input <- cbind(dissolved = runoff_dissolved, bound = runoff_bound, 
+    total = runoff_dissolved + runoff_bound)
+  rownames(runoff_input) <- rownames(reduction_runoff)
+  units(runoff_input) <- input_units
+
+  dn <- dimnames(runoff_input)
+  PEC_sw_runoff <- set_units(runoff_input / V_flowing_ditch_runoff, "\u00B5g/L")
+  dimnames(PEC_sw_runoff) <- dn
+
+  result <- list(
+    perc_runoff = 100 * f_runoff,
+    runoff = as.data.frame(runoff_input),
+    PEC_sw_runoff = as.data.frame(PEC_sw_runoff))
+  return(result)
+}
+
+#' Calculate PEC surface water due to drainage as in Exposit 3
+#'
+#' This is a reimplementation of the calculation described in the Exposit 3.02 spreadsheet file,
+#' in the worksheet "Konzept Drainage". Although there are four groups of
+#' compounds ("Gefährdungsgruppen"), only one distinction is made in the
+#' calculations, between compounds with low mobility (group 1) and compounds
+#' with modest to high mobility (groups 2, 3 and 4). In this implementation,
+#' the group is derived only from the Koc, if not given explicitly. For
+#' details, see the discussion of the function arguments below.
+#'
+#' @param rate The application rate in g/ha
+#' @param interception The fraction intercepted by the crop
+#' @param Koc The sorption coefficient to soil organic carbon used to determine the mobility. A trigger
+#'   value of 550 L/kg is used in order to decide if Koc >> 500.
+#' @param mobility Overrides what is determined from the Koc.
+#' @param DT50 The soil half-life in days
+#' @param t_drainage The time between application and the drainage event, where degradation occurs, in days
+#' @param V_ditch The volume of the ditch is assumed to be 1 m * 100 m * 30 cm = 30 m3
+#' @param V_drainage The drainage volume, equivalent to 1 mm precipitation on 1 ha for spring/summer or 10 mm for
+#'   autumn/winter/early spring.
+#' @param dilution The dilution factor
+#' @return A list containing the following components
+#'   \describe{
+#'     \item{perc_runoff}{The runoff percentages for dissolved and bound substance}
+#'     \item{runoff}{A matrix containing dissolved and bound input for the different distances}
+#'     \item{PEC_sw_runoff}{A matrix containing PEC values for dissolved and bound substance
+#'       for the different distances. If the rate was given in g/ha, the PECsw are in microg/L.}
+#'   }
+#' @export
+#' @source Excel 3.02 spreadsheet available from
+#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
+#' @seealso \code{\link{perc_runoff_exposit}} for runoff loss percentages and \code{\link{perc_runoff_reduction_exposit}} for runoff reduction percentages used
+#' @examples
+#'   PEC_sw_exposit_drainage(500, Koc = 150)
+PEC_sw_exposit_drainage <- function(rate, interception = 0, Koc = NA, mobility = c(NA, "low", "high"), DT50 = Inf, t_drainage = 3,
+  V_ditch = 30, V_drainage = c(spring = 10, autumn = 100), dilution = 2)
+{
+  # Rückstand zum Zeitpunkt des Niederschlagsereignisses (residue at the time of the drainage event)
+  k_deg <- log(2)/DT50
+  residue <- rate * (1 - interception) * 1 * exp(-k_deg * t_drainage) # assumes 1 ha treated area
+
+  mobility <- match.arg(mobility)
+  if (is.na(mobility)) {
+    if (is.na(Koc)) stop("Koc is needed if the mobility is not specified")
+    else {
+      if (Koc > 550) mobility = "low"
+      else mobility = "high"
+    }
+  }
+
+  V_ditch_drainage <- V_ditch + V_drainage
+  V_flowing_ditch_drainage <- dilution * V_ditch_drainage
+
+  # Gesamtaustrag (total fraction of the residue drained)
+  if (mobility == "low") {
+    f_drainage_total <- c(spring = 0.01 * 1e-2,
+                          autumn = 0.05 * 1e-2)
+  } else {
+    f_drainage_total <- c(spring = 0.2 * 1e-2,
+                          autumn = 1.0 * 1e-2)
+  }
+
+  f_peak = c(spring = 0.125, autumn = 0.25) # Stoßbelastung (fraction drained at event)
+
+  PEC_sw_drainage <- 1000 * residue * f_drainage_total * f_peak / V_flowing_ditch_drainage
+
+  result <- list(
+    perc_drainage_total = 100 * f_drainage_total,
+    perc_peak = 100 * f_peak,
+    PEC_sw_drainage = PEC_sw_drainage)
+  return(result)
+}
diff --git a/R/PEC_sw_exposit_runoff.R b/R/PEC_sw_exposit_runoff.R
deleted file mode 100644
index 8b89cd9..0000000
--- a/R/PEC_sw_exposit_runoff.R
+++ /dev/null
@@ -1,188 +0,0 @@
-#' Runoff loss percentages as used in Exposit 3
-#'
-#' A table of the loss percentages used in Exposit 3 for the twelve different Koc classes
-#'
-#' @name perc_runoff_exposit
-#' @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
-#'       adjacent water body bound to eroding particles}
-#'   }
-#' @source Excel 3.02 spreadsheet available from
-#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
-#' @docType data
-#' @examples
-#' print(perc_runoff_exposit)
-"perc_runoff_exposit"
-
-#' Runoff reduction percentages as used in Exposit
-#'
-#' A table of the runoff reduction percentages used in Exposit 3 for different vegetated buffer widths
-#'
-#' @name perc_runoff_reduction_exposit
-#' @format A named list of data frames with reduction percentage values for the
-#' dissolved fraction and the fraction bound to eroding particles, with
-#' vegetated buffer widths as row names. The names of the list items are the Exposit versions
-#' from which the values were taken.
-#'   \describe{
-#'     \item{dissolved}{The reduction percentage for the dissolved phase}
-#'     \item{bound}{The reduction percentage for the particulate phase}
-#'   }
-#' @source Excel 3.02 spreadsheet available from
-#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
-#'
-#'   Agroscope version 3.01a with additional runoff factors for 3 m and 6 m buffer zones received from Muris Korkaric (not published).
-#'   The variant 3.01a2 was introduced for consistency with previous calculations performed by Agroscope for a 3 m buffer zone.
-#' @docType data
-#' @examples
-#' print(perc_runoff_reduction_exposit)
-"perc_runoff_reduction_exposit"
-
-#' Calculate PEC surface water due to runoff and erosion as in Exposit 3
-#'
-#' This is a reimplementation of the calculation described in the Exposit 3.02 spreadsheet file,
-#' in the worksheet "Konzept Runoff".
-#'
-#' @param rate The application rate in g/ha
-#' @param interception The fraction intercepted by the crop
-#' @param Koc The sorption coefficient to soil organic carbon
-#' @param DT50 The soil half-life in days
-#' @param t_runoff The time between application and the runoff event, where degradation occurs, in days
-#' @param exposit_reduction_version The version of the reduction factors to be used. "3.02" is the current
-#'   version used in Germany, "3.01a" is the version with additional percentages for 3 m and 6 m buffer
-#'   zones used in Switzerland. "3.01a2" is a version introduced for consistency with previous calculations
-#'   performed for a 3 m buffer zone in Switzerland, with the same reduction being applied to the dissolved
-#'   and the bound fraction.
-#' @param V_ditch The volume of the ditch is assumed to be 1 m * 100 m * 30 cm = 30 m3
-#' @param V_event The unreduced runoff volume, equivalent to 10 mm precipitation on 1 ha
-#' @param dilution The dilution factor
-#' @return A list containing the following components
-#'   \describe{
-#'     \item{perc_runoff}{The runoff percentages for dissolved and bound substance}
-#'     \item{runoff}{A matrix containing dissolved and bound input for the different distances}
-#'     \item{PEC_sw_runoff}{A matrix containing PEC values for dissolved and bound substance
-#'       for the different distances. If the rate was given in g/ha, the PECsw are in microg/L.}
-#'   }
-#' @export
-#' @source Excel 3.02 spreadsheet available from
-#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
-#' @seealso \code{\link{perc_runoff_exposit}} for runoff loss percentages and \code{\link{perc_runoff_reduction_exposit}} for runoff reduction percentages used
-#' @examples
-#'   PEC_sw_exposit_runoff(500, Koc = 150)
-#'   PEC_sw_exposit_runoff(600, Koc = 10000, DT50 = 195, exposit = "3.01a")
-PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, DT50 = Inf, t_runoff = 3,
-  exposit_reduction_version = c("3.02", "3.01a", "3.01a2", "2.0"),
-  V_ditch = 30, V_event = 100, dilution = 2)
-{
-  k_deg <- log(2)/DT50
-  input <- rate * (1 - interception) * 1 * exp(-k_deg * t_runoff) # assumes 1 ha treated area
-
-  if (length(Koc) > 1) stop("Only one compound at a time supported")
-
-  exposit_reduction_version <- match.arg(exposit_reduction_version)
-  red_water <- pfm::perc_runoff_reduction_exposit[[exposit_reduction_version]]["dissolved"] / 100
-  red_bound <- pfm::perc_runoff_reduction_exposit[[exposit_reduction_version]]["bound"] / 100
-  reduction_runoff <- pfm::perc_runoff_reduction_exposit[[exposit_reduction_version]] / 100
-  transfer_runoff <- 1 - reduction_runoff
-
-  V_runoff <- V_event * (1 - reduction_runoff[["dissolved"]]) # m3
-  V_flowing_ditch_runoff <- dilution * (V_ditch + V_runoff)
-
-  f_runoff_exposit <- function(Koc) {
-    Koc_breaks <- c(pfm::perc_runoff_exposit$Koc_lower_bound, Inf)
-    Koc_classes <- as.character(cut(Koc, Koc_breaks, labels = rownames(pfm::perc_runoff_exposit)))
-    perc_runoff <- pfm::perc_runoff_exposit[Koc_classes, c("dissolved", "bound")]
-    if (identical(Koc, 0)) perc_runoff <- c(dissolved = 0, bound = 0)
-    return(unlist(perc_runoff) / 100)
-  }
-  f_runoff <- f_runoff_exposit(Koc)
-  runoff_dissolved <- input * f_runoff["dissolved"] * transfer_runoff["dissolved"]
-  runoff_bound <- input * f_runoff["bound"] * transfer_runoff["bound"]
-  runoff_input <- cbind(runoff_dissolved, runoff_bound)
-  runoff_input$total <- runoff_input$dissolved + runoff_input$bound
-
-  PEC_sw_runoff <- 1000 * runoff_input / V_flowing_ditch_runoff
-
-  result <- list(
-    perc_runoff = 100 * f_runoff,
-    runoff = runoff_input,
-    PEC_sw_runoff = PEC_sw_runoff)
-  return(result)
-}
-
-#' Calculate PEC surface water due to drainage as in Exposit 3
-#'
-#' This is a reimplementation of the calculation described in the Exposit 3.02 spreadsheet file,
-#' in the worksheet "Konzept Drainage". Although there are four groups of
-#' compounds ("Gefährdungsgruppen"), only one distinction is made in the
-#' calculations, between compounds with low mobility (group 1) and compounds
-#' with modest to high mobility (groups 2, 3 and 4). In this implementation,
-#' the group is derived only from the Koc, if not given explicitly. For
-#' details, see the discussion of the function arguments below.
-#'
-#' @param rate The application rate in g/ha
-#' @param interception The fraction intercepted by the crop
-#' @param Koc The sorption coefficient to soil organic carbon used to determine the mobility. A trigger
-#'   value of 550 L/kg is used in order to decide if Koc >> 500.
-#' @param mobility Overrides what is determined from the Koc.
-#' @param DT50 The soil half-life in days
-#' @param t_drainage The time between application and the drainage event, where degradation occurs, in days
-#' @param V_ditch The volume of the ditch is assumed to be 1 m * 100 m * 30 cm = 30 m3
-#' @param V_drainage The drainage volume, equivalent to 1 mm precipitation on 1 ha for spring/summer or 10 mm for
-#'   autumn/winter/early spring.
-#' @param dilution The dilution factor
-#' @return A list containing the following components
-#'   \describe{
-#'     \item{perc_runoff}{The runoff percentages for dissolved and bound substance}
-#'     \item{runoff}{A matrix containing dissolved and bound input for the different distances}
-#'     \item{PEC_sw_runoff}{A matrix containing PEC values for dissolved and bound substance
-#'       for the different distances. If the rate was given in g/ha, the PECsw are in microg/L.}
-#'   }
-#' @export
-#' @source Excel 3.02 spreadsheet available from
-#'   \url{https://www.bvl.bund.de/SharedDocs/Downloads/04_Pflanzenschutzmittel/zul_umwelt_exposit.html}
-#' @seealso \code{\link{perc_runoff_exposit}} for runoff loss percentages and \code{\link{perc_runoff_reduction_exposit}} for runoff reduction percentages used
-#' @examples
-#'   PEC_sw_exposit_drainage(500, Koc = 150)
-PEC_sw_exposit_drainage <- function(rate, interception = 0, Koc = NA, mobility = c(NA, "low", "high"), DT50 = Inf, t_drainage = 3,
-  V_ditch = 30, V_drainage = c(spring = 10, autumn = 100), dilution = 2)
-{
-  # Rückstand zum Zeitpunkt des Niederschlagsereignisses (residue at the time of the drainage event)
-  k_deg <- log(2)/DT50
-  residue <- rate * (1 - interception) * 1 * exp(-k_deg * t_drainage) # assumes 1 ha treated area
-
-  mobility <- match.arg(mobility)
-  if (is.na(mobility)) {
-    if (is.na(Koc)) stop("Koc is needed if the mobility is not specified")
-    else {
-      if (Koc > 550) mobility = "low"
-      else mobility = "high"
-    }
-  }
-
-  V_ditch_drainage <- V_ditch + V_drainage
-  V_flowing_ditch_drainage <- dilution * V_ditch_drainage
-
-  # Gesamtaustrag (total fraction of the residue drained)
-  if (mobility == "low") {
-    f_drainage_total <- c(spring = 0.01 * 1e-2,
-                          autumn = 0.05 * 1e-2)
-  } else {
-    f_drainage_total <- c(spring = 0.2 * 1e-2,
-                          autumn = 1.0 * 1e-2)
-  }
-
-  f_peak = c(spring = 0.125, autumn = 0.25) # Stoßbelastung (fraction drained at event)
-
-  PEC_sw_drainage <- 1000 * residue * f_drainage_total * f_peak / V_flowing_ditch_drainage
-
-  result <- list(
-    perc_drainage_total = 100 * f_drainage_total,
-    perc_peak = 100 * f_peak,
-    PEC_sw_drainage = PEC_sw_drainage)
-  return(result)
-}
diff --git a/R/PEC_sw_sed.R b/R/PEC_sw_sed.R
index 2865ab7..28e12d8 100644
--- a/R/PEC_sw_sed.R
+++ b/R/PEC_sw_sed.R
@@ -10,24 +10,30 @@
 #' @param method The method used for the calculation
 #' @param sediment_depth Depth of the sediment layer
 #' @param water_depth Depth of the water body in cm
-#' @param sediment_density The density of the sediment in L/kg (equivalent to
+#' @param sediment_density The density of the sediment in kg/L (equivalent to
 #'   g/cm3)
 #' @param PEC_sed_units The units of the estimated sediment PEC value
 #' @return The predicted concentration in sediment
 #' @export
 #' @author Johannes Ranke
 #' @examples
+#' library(pfm)
+#' library(units)
 #' PEC_sw_sed(PEC_sw_drift(100, distances = 1), percentage = 50)
 PEC_sw_sed <- function(PEC_sw, percentage = 100, method = "percentage", 
-                       sediment_depth = 5, water_depth = 30,
-                       sediment_density = 1.3,
+                       sediment_depth = set_units(5, "cm"),
+                       water_depth = set_units(30, "cm"),
+                       sediment_density = set_units(1.3, "kg/L"),
                        PEC_sed_units = c("\u00B5g/kg", "mg/kg"))
 {
-  method = match.arg(method)
-  PEC_sed_units = match.arg(PEC_sed_units)
+  if (!inherits(PEC_sw, "units")) PEC_sw <- set_units(PEC_sw, "\u00B5g/L")
+  if (!inherits(sediment_depth, "units")) PEC_sw <- set_units(sediment_depth, "cm")
+  if (!inherits(water_depth, "units")) PEC_sw <- set_units(water_depth, "cm")
+  if (!inherits(sediment_density, "units")) PEC_sw <- set_units(sediment_density, "cm")
+  method <- match.arg(method)
+  PEC_sed_units <- match.arg(PEC_sed_units)
   if (method == "percentage") {
-    PEC_sed = PEC_sw * (percentage/100) * (water_depth / sediment_depth) * (1 / sediment_density)
-    if (PEC_sed_units == "mg/kg") PEC_sed <- PEC_sed / 1000
+    PEC_sed <- PEC_sw * (percentage/100) * as.numeric((water_depth / sediment_depth)) * (1 / sediment_density)
   }
-  return(PEC_sed)
+  return(set_units(PEC_sed, PEC_sed_units, mode = "symbolic"))
 }
-- 
cgit v1.2.1