Start of an Steps 1/2 calculator in R

6 files changed, 256 insertions, 0 deletions

diff --git a/R/FOCUS_Step_12_scenarios.R b/R/FOCUS_Step_12_scenarios.R
new file mode 100644
index 0000000..cb0654a
--- /dev/null
+++ b/R/FOCUS_Step_12_scenarios.R
@@ -0,0 +1,37 @@
+#' Step 1/2 scenario data as distributed with the FOCUS Step 1/2 calculator
+#'
+#' The data were extracted from the scenario.txt file using the R code shown below.
+#' The text file is not included in the package as its licence is not clear.
+#'
+#' @name FOCUS_Step_12_scenarios
+#' @docType data
+#' @format A list containing the scenario names in a character vector called 'names',
+#'   the drift percentiles in a matrix called 'drift', interception percentages in
+#'   a matrix called 'interception' and the runoff/drainage percentages for Step 2
+#'   calculations in a matrix called 'rd'.
+#' @keywords datasets
+#' @examples
+#'
+#' \dontrun{
+#'   # This is the code that was used to extract the data
+#'   scenario_path <- "inst/extdata/FOCUS_Step_12_scenarios.txt"
+#'   scenarios <- readLines(scenario_path)[9:38]
+#'   FOCUS_Step_12_scenarios <- list()
+#'   sce <- read.table(text = scenarios, sep = "\t", header = TRUE, check.names = FALSE,
+#'     stringsAsFactors = FALSE)
+#'   FOCUS_Step_12_scenarios$names = sce$Crop
+#'   rownames(sce) <- sce$Crop
+#'   FOCUS_Step_12_scenarios$drift = sce[, 3:11]
+#'   FOCUS_Step_12_scenarios$interception = sce[, 12:15]
+#'   sce_2 <- readLines(scenario_path)[41:46]
+#'   rd <- read.table(text = sce_2, sep = "\t")[1:2]
+#'   rd_mat <- matrix(rd$V2, nrow = 3, byrow = FALSE)
+#'   dimnames(rd_mat) = list(Time = c("Oct-Feb", "Mar-May", "Jun-Sep"),
+#'                           Region = c("North", "South"))
+#'   FOCUS_Step_12_scenarios$rd = rd_mat
+#'   save(FOCUS_Step_12_scenarios, file = "data/FOCUS_Step_12_scenarios.RData")
+#' }
+#'
+#' # And this is the resulting data
+#' FOCUS_Step_12_scenarios
+NULL
diff --git a/R/PEC_sw_focus.R b/R/PEC_sw_focus.R
new file mode 100644
index 0000000..5f3b3d2
--- /dev/null
+++ b/R/PEC_sw_focus.R
@@ -0,0 +1,109 @@
+#' Calculate FOCUS Step 1 PEC surface water
+#'
+#' This is an attempt to reimplement the FOCUS Step 1 and 2 calculator authored
+#' by Michael Klein. The Step 1 and 2 scenario assumes an area ratio of 10:1
+#' between field and waterbody, and a water depth of 30 cm.
+#' I did not (yet) implement the TWA formulas for times later than day 1, as I
+#' did not understand them right away.
+#' Also, Step 2 is not implemented (yet).
+#'
+#' @export
+#' @param parent A list containing substance specific parameters
+#' @param rate The application rate in g/ha. Overriden when
+#'   applications are given explicitly
+#' @param n The number of applications
+#' @param i The application interval
+#' @param applications A dataframe containing times and amounts of each application
+#' @param step At the moment, only Step 1 is implemented
+#' @examples
+#' dummy_1 <- chent_focus_sw(cwsat = 6000, DT50_ws = 6, Koc = 344.8)
+#' PEC_sw_focus(dummy_1, 3000, f_drift = 0)
+PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
+  applications = data.frame(time = seq(0, 0 + n * i, length.out = n),
+                            amount = rate),
+  met = NULL,
+  step = 1,
+  f_drift = 0.02759, f_rd = 0.1)
+{
+  if (is.null(met)) {
+    mw_ratio = 1
+    max_soil = 1
+    max_ws = 1
+    Koc = parent$Koc
+    DT50_ws = parent$DT50_ws
+  } else {
+    mw_ratio = met$mw / parent$mw
+    max_soil = met$max_soil
+    max_ws = met$max_ws
+    Koc = met$Koc
+    DT50_ws = met$DT50_ws
+  }
+
+  # Rates for a single application
+  eq_rate_drift_s = mw_ratio * max_ws * rate
+  eq_rate_rd_s = mw_ratio * max_soil * rate
+
+  # Drift input
+  input_drift_s = f_drift * eq_rate_drift_s / 10 # mg/m2
+  input_drift = n * input_drift_s
+
+  # Runoff/drainage input
+  ratio_field_wb = 10 # 10 m2 of field for each m2 of the waterbody
+  input_rd_s = f_rd * eq_rate_rd_s * ratio_field_wb / 10
+  input_rd = n * input_rd_s
+
+  # Fraction of compound entering the water phase via runoff/drainage
+  depth_sw  = 0.3  # m
+  depth_sed = 0.05 # m
+  depth_sed_eff = 0.01 # m, only relevant for sorption
+  rho_sed = 0.8 # kg/L
+  f_OC = 0.05 # 5% organic carbon in sediment
+  f_rd_sw = depth_sw / (depth_sw + (depth_sed_eff * rho_sed * f_OC * Koc))
+  f_rd_sed = 1 - f_rd_sw
+
+  # Initial PECs
+  PEC_sw_0 = (input_drift + input_rd * f_rd_sw) / depth_sw        # µg/L
+  PEC_sed_0 = (input_rd * f_rd_sed) / (depth_sed * rho_sed)       # µg/kg
+
+  # Initial PECs when assuming partitioning also of drift input
+  PEC_sw_0_part = (input_drift + input_rd) * f_rd_sw / depth_sw                  # µg/L
+  PEC_sed_0_part = (input_drift + input_rd) * f_rd_sed / (depth_sed * rho_sed)   # µg/kg
+
+  t_out = c(0, 1, 2, 4, 7, 14, 21, 28, 42, 50, 100)
+  PEC = matrix(NA, nrow = length(t_out), ncol = 4,
+               dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+
+  PEC["0", "PECsw"]  = PEC_sw_0
+  PEC["0", "PECsed"] = PEC_sed_0
+
+  # Degradation after partitioning according to Koc
+  k_ws = log(2)/DT50_ws
+  PEC[as.character(t_out[-1]), "PECsw"]  = PEC_sw_0_part  * exp( - k_ws * t_out[-1])
+  PEC[as.character(t_out[-1]), "PECsed"] = PEC_sed_0_part * exp( - k_ws * t_out[-1])
+
+  # TWA concentrations
+  PEC["1", "TWAECsw"] = (PEC_sw_0 + PEC["1", "PECsw"]) / 2
+  PEC["1", "TWAECsed"] = (PEC_sed_0 + PEC["1", "PECsed"]) / 2
+
+  list(eq_rate_drift_s = eq_rate_drift_s,
+       eq_rate_rd_s = eq_rate_rd_s,
+       input_drift_s = input_drift_s,
+       input_rd_s = input_rd_s,
+       f_rd_sw = f_rd_sw, f_rd_sed = f_rd_sed,
+       PEC = PEC)
+}
+
+#' Create an chemical compound object for FOCUS Step 1 and 2 calculations
+#'
+#' @export
+#' @param cwsat Water solubility in mg/L
+#' @param DT50_ws Half-life in water/sediment systems in days
+#' @param Koc Partition coefficient between organic carbon and water
+#'   in L/kg.
+#' @return A list with the substance specific properties
+chent_focus_sw <- function(Koc, DT50_ws, cwsat = 1000)
+{
+  list(Koc = Koc, DT50_ws = DT50_ws, cwsat = cwsat)
+}
+
+
diff --git a/data/FOCUS_Step_12_scenarios.RData b/data/FOCUS_Step_12_scenarios.RData
new file mode 100644
index 0000000..08b8f4d
--- /dev/null
+++ b/data/FOCUS_Step_12_scenarios.RData
diff --git a/man/PEC_sw_focus.Rd b/man/PEC_sw_focus.Rd
new file mode 100644
index 0000000..0f444ca
--- /dev/null
+++ b/man/PEC_sw_focus.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PEC_sw_focus.R
+\name{PEC_sw_focus}
+\alias{PEC_sw_focus}
+\title{Calculate FOCUS Step 1 PEC surface water}
+\usage{
+PEC_sw_focus(parent, rate, n = 1, i = NA, applications = data.frame(time =
+  seq(0, 0 + n * i, length.out = n), amount = rate), met = NULL, step = 1,
+  f_drift = 0.02759, f_rd = 0.1)
+}
+\arguments{
+\item{parent}{A list containing substance specific parameters}
+
+\item{rate}{The application rate in g/ha. Overriden when
+applications are given explicitly}
+
+\item{n}{The number of applications}
+
+\item{i}{The application interval}
+
+\item{applications}{A dataframe containing times and amounts of each application}
+
+\item{step}{At the moment, only Step 1 is implemented}
+}
+\description{
+This is an attempt to reimplement the FOCUS Step 1 and 2 calculator authored
+by Michael Klein. The Step 1 and 2 scenario assumes an area ratio of 10:1
+between field and waterbody, and a water depth of 30 cm.
+I did not (yet) implement the TWA formulas for times later than day 1, as I
+did not understand them right away.
+Also, Step 2 is not implemented (yet).
+}
+\examples{
+dummy_1 <- chent_focus_sw(cwsat = 6000, DT50_ws = 6, Koc = 344.8)
+PEC_sw_focus(dummy_1, 3000, f_drift = 0)
+}
diff --git a/man/chent_focus_sw.Rd b/man/chent_focus_sw.Rd
new file mode 100644
index 0000000..26742d0
--- /dev/null
+++ b/man/chent_focus_sw.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PEC_sw_focus.R
+\name{chent_focus_sw}
+\alias{chent_focus_sw}
+\title{Create an chemical compound object for FOCUS Step 1 and 2 calculations}
+\usage{
+chent_focus_sw(Koc, DT50_ws, cwsat = 1000)
+}
+\arguments{
+\item{Koc}{Partition coefficient between organic carbon and water
+in L/kg.}
+
+\item{DT50_ws}{Half-life in water/sediment systems in days}
+
+\item{cwsat}{Water solubility in mg/L}
+}
+\value{
+A list with the substance specific properties
+}
+\description{
+Create an chemical compound object for FOCUS Step 1 and 2 calculations
+}
diff --git a/tests/testthat/test_step_1.R b/tests/testthat/test_step_1.R
new file mode 100644
index 0000000..abef72c
--- /dev/null
+++ b/tests/testthat/test_step_1.R
@@ -0,0 +1,52 @@
+context("FOCUS Step 1 calculations")
+
+test_that("Results of Steps 1/2 calculator for Dummy 1 are reproduced", {
+  dummy_1 <- chent_focus_sw(cwsat = 6000, DT50_ws = 6, Koc = 344.8)
+  res_dummy_1 <- PEC_sw_focus(dummy_1, 3000, f_drift = 0)
+
+  t_out <- c(0, 1, 2, 4) # Checking the first four days should be sufficient for Step 1
+  PEC_orig = matrix(NA, nrow = length(t_out), ncol = 4,
+    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+
+  PEC_orig_1[, "PECsw"]  = c(685.06, 610.32, 543.73, 431.56)
+  PEC_orig_1["1", "TWAECsw"]  = 647.69 # Later TWAEC not implemented
+  PEC_orig_1[, "PECsed"]  = c(2.36, 2.1, 1.87, 1.49) * 1e3
+  PEC_orig_1["1", "TWAECsed"]  = 2.23e3 # Later TWAEC not implemented
+
+  expect_equal(res_dummy_1$PEC[1:4, c(1, 2)], PEC_orig_1[, c(1, 2)], tolerance = 0.01, scale = 1)
+  expect_equal(res_dummy_1$PEC[1:4, c(3, 4)], PEC_orig_1[, c(3, 4)], tolerance = 10, scale = 1)
+})
+
+test_that("Results of Steps 1/2 calculator for Dummy 2 are reproduced", {
+  dummy_2 <- chent_focus_sw(cwsat = 30, DT50_ws = 26, Koc = 110)
+  res_dummy_2 <- PEC_sw_focus(dummy_2, 1000)
+
+  t_out <- c(0, 1, 2, 4) # Checking the first four days should be sufficient for Step 1
+  PEC_orig_2 = matrix(NA, nrow = length(t_out), ncol = 4,
+    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+
+  PEC_orig_2[, "PECsw"]  = c(299.89, 290.86, 283.21, 268.50)
+  PEC_orig_2["1", "TWAECsw"]  = 295.38 # Later TWAEC not implemented
+  PEC_orig_2[, "PECsed"]  = c(319.77, 319.95, 311.53, 295.35)
+  PEC_orig_2["1", "TWAECsed"]  = 319.86 # Later TWAEC not implemented
+
+  expect_equal(res_dummy_2$PEC[1:4, c(1, 2)], PEC_orig_2[, c(1, 2)], tolerance = 0.01, scale = 1)
+  expect_equal(res_dummy_2$PEC[1:4, c(3, 4)], PEC_orig_2[, c(3, 4)], tolerance = 10, scale = 1)
+})
+
+test_that("Results of Steps 1/2 calculator for Dummy 2 are reproduced", {
+  dummy_2 <- chent_focus_sw(cwsat = 30, DT50_ws = 26, Koc = 110)
+  res_dummy_2 <- PEC_sw_focus(dummy_2, 1000)
+
+  t_out <- c(0, 1, 2, 4) # Checking the first four days should be sufficient for Step 1
+  PEC_orig_2 = matrix(NA, nrow = length(t_out), ncol = 4,
+    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+
+  PEC_orig_2[, "PECsw"]  = c(299.89, 290.86, 283.21, 268.50)
+  PEC_orig_2["1", "TWAECsw"]  = 295.38 # Later TWAEC not implemented
+  PEC_orig_2[, "PECsed"]  = c(319.77, 319.95, 311.53, 295.35)
+  PEC_orig_2["1", "TWAECsed"]  = 319.86 # Later TWAEC not implemented
+
+  expect_equal(res_dummy_2$PEC[1:4, c(1, 2)], PEC_orig_2[, c(1, 2)], tolerance = 0.01, scale = 1)
+  expect_equal(res_dummy_2$PEC[1:4, c(3, 4)], PEC_orig_2[, c(3, 4)], tolerance = 10, scale = 1)
+})