diff options
| -rw-r--r-- | R/FOCUS_Step_12_scenarios.R | 37 | ||||
| -rw-r--r-- | R/PEC_sw_focus.R | 109 | ||||
| -rw-r--r-- | data/FOCUS_Step_12_scenarios.RData | bin | 0 -> 1319 bytes | |||
| -rw-r--r-- | man/PEC_sw_focus.Rd | 36 | ||||
| -rw-r--r-- | man/chent_focus_sw.Rd | 22 | ||||
| -rw-r--r-- | tests/testthat/test_step_1.R | 52 |
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 Binary files differnew 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) +}) |