#' Calculate PEC surface water at FOCUS Step 1
#'
#' This is a reimplementation of the FOCUS Step 1 and 2 calculator version 3.2,
#' authored by Michael Klein, in R. Note that results for multiple
#' applications should be compared to the corresponding results for a
#' single application. At current, this is not done automatically in
#' this implementation. Only Step 1 PECs are calculated. However,
#' input files are generated that are suitable as input also for Step 2
#' to be used with the FOCUS calculator.
#'
#' @importFrom utils read.table
#' @references FOCUS (2014) Generic guidance for Surface Water Scenarios (version 1.4).
#' FOrum for the Co-ordination of pesticde fate models and their USe.
#' http://esdac.jrc.ec.europa.eu/public_path/projects_data/focus/sw/docs/Generic%20FOCUS_SWS_vc1.4.pdf
#' @references Website of the Steps 1 and 2 calculator at the Joint Research
#' Center of the European Union:
#' http://esdac.jrc.ec.europa.eu/projects/stepsonetwo
#' @note The formulas for input to the waterbody via runoff/drainage of the
#' parent and subsequent formation of the metabolite in water is not
#' documented in the model description coming with the calculator. As one would
#' expect, this appears to be (as we get the same results) calculated by
#' multiplying the application rate with the molar weight
#' correction and the formation fraction in water/sediment systems.
#' @note Step 2 is not implemented.
#' @export
#' @param parent A list containing substance specific parameters, e.g.
#' conveniently generated by [chent_focus_sw].
#' @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 comment A comment for the input file
#' @param met A list containing metabolite specific parameters. e.g.
#' conveniently generated by [chent_focus_sw]. If not NULL,
#' the PEC is calculated for this compound, not the parent.
#' @param f_drift The fraction of the application rate reaching the waterbody
#' via drift. If NA, this is derived from the scenario name and the number
#' of applications via the drift data defined by the
#' [FOCUS_Step_12_scenarios]
#' @param f_rd The fraction of the amount applied reaching the waterbody via
#' runoff/drainage. At Step 1, it is assumed to be 10%, be it the
#' parent or a metabolite
#' @param scenario The name of the scenario. Must be one of the scenario
#' names given in [FOCUS_Step_12_scenarios]
#' @param region 'n' for Northern Europe or 's' for Southern Europe. If NA, only
#' Step 1 PECsw are calculated
#' @param season 'of' for October to February, 'mm' for March to May, and 'js'
#' for June to September. If NA, only step 1 PECsw are calculated
#' @param interception One of 'no interception' (default), 'minimal crop cover',
#' 'average crop cover' or 'full canopy'
#' @param met_form_water Should the metabolite formation in water be taken into
#' account? This can be switched off to check the influence and to compare
#' with previous versions of the Steps 12 calculator
#' @param txt_file the name, and potentially the full path to the
#' Steps.12 input text file to which the specification of the run(s)
#' should be written
#' @param overwrite Should an existing file a the location specified in
#' \code{txt_file} be overwritten? Only takes effect if append is FALSE.
#' @param append Should the input text file be appended?
#' @examples
#' # Parent only
#' dummy_1 <- chent_focus_sw("Dummy 1", cwsat = 6000, DT50_ws = 6, Koc = 344.8)
#' PEC_sw_focus(dummy_1, 3000, f_drift = 0, overwrite = TRUE, append = FALSE)
#'
#' # Metabolite
#' new_dummy <- chent_focus_sw("New Dummy", mw = 250, Koc = 100)
#' M1 <- chent_focus_sw("M1", mw = 100, cwsat = 100, DT50_ws = 100, Koc = 50,
#' max_ws = 0, max_soil = 0.5)
#' PEC_sw_focus(new_dummy, 1000, scenario = "cereals, winter", met = M1)
PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
comment = "",
met = NULL,
f_drift = NA, f_rd = 0.1,
scenario = FOCUS_Step_12_scenarios$names,
region = c("n", "s"),
season = c(NA, "of", "mm", "js"),
interception = c("no interception", "minimal crop cover",
"average crop cover", "full canopy"),
met_form_water = TRUE,
txt_file = "pesticide.txt", overwrite = FALSE, append = TRUE)
{
if (n > 1 & is.na(i)) stop("Please specify the interval i if n > 1")
scenario = match.arg(scenario)
if (is.na(f_drift)) {
f_drift = FOCUS_Step_12_scenarios$drift[scenario, "1"] / 100
# For Step 2 we would select the reduced percentiles for multiple apps:
if (n <= 8) {
f_drift_2 = FOCUS_Step_12_scenarios$drift[scenario, as.character(n)] / 100
} else {
f_drift_2 = FOCUS_Step_12_scenarios$drift[scenario, ">8"] / 100
}
}
interception = match.arg(interception)
# Write to txt file if requested
header <- c("Active Substance", "Compound", "Comment",
"Mol mass a.i.", "Mol mass met.",
"Water solubility",
"KOC assessed compound", "KOC parent compound",
"DT50",
"Max. in Water",
"Max. in Soil asessed compound", # we reproduce the typo...
"App. Rate", "Number of App.", "Time between app.", "App. Type",
"DT50 soil parent compound", "DT50 soil",
"DT50 water", "DT50 sediment",
"Region / Season",
"Interception class")
add_line <- function(x) {
cat(paste0(x, "\r\n"), file = txt, append = TRUE)
}
if (file.exists(txt_file)) {
if (append) {
txt <- file(txt_file, "a")
} else {
if (overwrite) {
txt <- file(txt_file, "w")
add_line(paste(header, collapse = "\t"))
} else {
stop("The file", txt_file, "already exists, and you did not request",
"appending or overwriting it")
}
}
} else {
txt <- file(txt_file, "w")
add_line(paste(header, collapse = "\t"))
}
on.exit(close(txt))
region = match.arg(region)
season = match.arg(season)
interception = match.arg(interception)
if (is.null(met)) {
compound = parent$name
cwsat = parent$cwsat
mw_ratio = 1
max_soil = 1
max_ws = 1
Koc_assessed = parent$Koc
DT50_ws = parent$DT50_ws
DT50_soil = parent$DT50_soil
} else {
compound = met$name
cwsat = met$cwsat
mw_ratio = met$mw / parent$mw
max_soil = met$max_soil
max_ws = met$max_ws
Koc_assessed = met$Koc
DT50_ws = met$DT50_ws
DT50_soil = met$DT50_soil
}
# In the Step 12 input file, the application type is the index of the drift
# scenarios, with numbering starting at 0
app_type = which(scenario == rownames(FOCUS_Step_12_scenarios$drift)) - 1
# The interception class is a number from 1.00 (no interception) to 4.00 (full canopy)
int_class = which(interception == colnames(FOCUS_Step_12_scenarios$interception))
# We calculate Step 2 if a season is specified
# Interception and region/season code
if (is.na(season[1])) {
f_int = 0
reg_sea = 0
} else {
f_int = FOCUS_Step_12_scenarios$interception[scenario, interception]
# the code for region / season is the sum of the region code
reg_code = switch(region,
"n" = 0,
"s" = 3)
# and the season code
sea_code = switch(season,
"of" = 0,
"mm" = 1,
"js" = 2)
reg_sea = reg_code + sea_code
}
scenario_safe <- sub(" /", " or", scenario)
if (is.null(met)) {
name_input <- paste(parent$name, scenario_safe, region, season)
} else {
name_input <- paste(met$name, scenario_safe, region, season)
}
if (comment != "") name_input = paste(name_input, comment)
run_txt <- c(
ai = name_input, compound = name_input,
comment = comment)
run_numeric <- c(
mw_ai = NA, mw_met = NA,
cwsat = cwsat, Koc_assessed = Koc_assessed, Koc_parent = 0,
DT50_ws = DT50_ws,
max_ws = 0, max_soil = 0,
rate = rate, n = n, i = if(is.na(i)) 0 else i,
app_type = app_type,
DT50_soil_parent = 0, DT50_soil = 0, DT50_water = 0, DT50_sediment = 0,
reg_sea = reg_sea, int_class = int_class)
if (is.null(met)) {
run_numeric["DT50_soil"] = parent$DT50_soil
run_numeric["DT50_water"] = parent$DT50_water
run_numeric["DT50_sediment"] = parent$DT50_sediment
} else {
run_numeric["mw_ai"] = parent$mw
run_numeric["mw_met"] = met$mw
run_numeric["max_soil"] = 100 * met$max_soil
run_numeric["max_ws"] = 100 * met$max_ws
run_numeric["Koc_parent"] = parent$Koc
run_numeric["DT50_soil"] = met$DT50_soil
run_numeric["DT50_soil_parent"] = parent$DT50_soil
run_numeric["DT50_water"] = met$DT50_water
run_numeric["DT50_sediment"] = met$DT50_sediment
}
print_numeric <- function(x) {
ifelse(is.na(x), " -99.00",
ifelse(x < 0.01,
sprintf("%8.2E", x),
sprintf("%8.2f", x)))
}
run_line <- paste(c(run_txt, print_numeric(run_numeric)), collapse = "\t")
add_line(run_line)
# Rates for a single application (suffix _s)
eq_rate_drift_s = mw_ratio * max_ws * rate
# Parent only, or metabolite formed in soil:
eq_rate_rd_s = mw_ratio * max_soil * rate
# Metabolite formed in water (this part is not documented in the Help files
# of the Steps 1/2 calculator):
if (!is.null(met)) {
if (met_form_water) {
eq_rate_rd_parent_s = mw_ratio * max_ws * rate
eq_rate_rd_s_tot = eq_rate_rd_s + eq_rate_rd_parent_s
} else {
eq_rate_rd_parent_s = NA
eq_rate_rd_s_tot = eq_rate_rd_s
}
} else {
eq_rate_rd_parent_s = NA
eq_rate_rd_s_tot = eq_rate_rd_s
}
# Drift input
input_drift_s = f_drift * eq_rate_drift_s / 10 # mg/m2
# 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_tot * ratio_field_wb / 10
input_rd = n * input_rd_s
# Step 1 accumulation in the aquatic system
# No accumulation between multiple applications if 3 * DT50 < i
if (n > 1 && (3 * DT50_ws < i)) {
input_drift = input_drift_s
input_rd = input_rd_s
} else {
input_drift = n * input_drift_s
input_rd = n * input_rd_s
}
# Step 2 accumulation in soil
# Use interception and degradation in soil for Step 2
k_soil = log(2)/DT50_soil
f_deg_soil = exp(- 4 * k_soil) * # Four days of degradation
(1 - exp(-n * i * k_soil)) / (1 - exp(-i * k_soil)) # multiple applications
# 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_assessed))
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_sw_1 = PEC["1", "PECsw"]
PEC_sed_1 = PEC["1", "PECsed"]
TWAEC_sw_1 = (PEC_sw_0 + PEC_sw_1) / 2
TWAEC_sed_1 = (PEC_sed_0 + PEC_sed_1) / 2
TWAEC_after = function(TWAEC_1, PEC_1, t) {
(TWAEC_1 + PEC_1 * (1 - exp(- k_ws * (t - 1))) / k_ws) / t
}
PEC["1", "TWAECsw"] = TWAEC_sw_1
PEC["1", "TWAECsed"] = TWAEC_sed_1
PEC[as.character(t_out[-1]), "TWAECsw"] = TWAEC_after(TWAEC_sw_1, PEC_sw_1, t_out[-1])
PEC[as.character(t_out[-1]), "TWAECsed"] = TWAEC_after(TWAEC_sed_1, PEC_sed_1, t_out[-1])
# Check if PEC_sw_max is above water solubility
PEC_sw_max = max(PEC[, "PECsw"])
if (!is.na(PEC_sw_max) & PEC_sw_max > 1000 * cwsat) {
warning("The maximum PEC surface water exceeds the water solubility")
}
PEC_sed_max = max(PEC[, "PECsed"])
list(f_drift = f_drift,
eq_rate_drift_s = eq_rate_drift_s,
eq_rate_rd_s = eq_rate_rd_s,
eq_rate_rd_parent_s = eq_rate_rd_parent_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,
PEC_sw_max = PEC_sw_max,
PEC_sed_max = PEC_sed_max
)
}
#' Create a chemical compound object for FOCUS Step 1 calculations
#'
#' @export
#' @param name Length one character vector containing the name
#' @param cwsat Water solubility in mg/L
#' @param DT50_ws Half-life in water/sediment systems in days
#' @param DT50_soil Half-life in soil in days
#' @param DT50_water Half-life in water in days (Step 2)
#' @param DT50_sediment Half-life in sediment in days (Step 2)
#' @param Koc Partition coefficient between organic carbon and water
#' in L/kg.
#' @param mw Molar weight in g/mol.
#' @param max_soil Maximum observed fraction (dimensionless) in soil
#' @param max_ws Maximum observed fraction (dimensionless) in water/sediment
#' systems
#' @return A list with the substance specific properties
chent_focus_sw <- function(name, Koc, DT50_ws = NA, DT50_soil = NA,
DT50_water = NA, DT50_sediment = NA,
cwsat = 1000, mw = NA, max_soil = 1, max_ws = 1)
{
list(name = name, Koc = Koc,
DT50_ws = DT50_ws, DT50_soil = DT50_soil, DT50_water = DT50_water,
DT50_sediment = DT50_sediment,
cwsat = cwsat, mw = mw, max_soil = max_soil, max_ws = max_ws)
}
