aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/FOCUS_Step_12_scenarios.R37
-rw-r--r--R/PEC_sw_focus.R109
2 files changed, 146 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)
+}
+
+

Contact - Imprint