From ff7e67a4d3415419dd3f712ef1af7467ebf65508 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 21 Sep 2018 19:39:44 +0200 Subject: Support FOMC in PEC_soil --- ChangeLog | 6 +++++ DESCRIPTION | 5 ++-- R/PEC_soil.R | 74 +++++++++++++++++++++++++++++++++++++++++++++++++-------- _pkgdown.yml | 1 + build.log | 2 +- man/PEC_soil.Rd | 16 +++++++++---- 6 files changed, 85 insertions(+), 19 deletions(-) diff --git a/ChangeLog b/ChangeLog index 5e265c0..69d5e6b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,9 @@ +commit 03c3035ca01c66b6a1352f7e509753fe2d057af2 +Author: Johannes Ranke +Date: 2018-07-11 03:41:43 +0200 + + Improve PELMO tests + commit 22b36c824fe5e1561868a649216fe079c6fbfb85 Author: Johannes Ranke Date: 2018-07-10 18:06:29 +0200 diff --git a/DESCRIPTION b/DESCRIPTION index e795eb1..932a4c8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: pfm Type: Package Title: Utilities for Pesticide Fate Modelling -Version: 0.5.1 -Date: 2018-07-11 +Version: 0.5.2 +Date: 2018-09-20 Authors@R: person("Johannes Ranke", email = "jranke@uni-bremen.de", role = c("aut", "cre", "cph"), comment = c(ORCID = "0000-0003-4371-6538")) @@ -20,7 +20,6 @@ Suggests: testthat, chents, magrittr, - R.utils, PELMO.installeR License: GPL LazyLoad: yes diff --git a/R/PEC_soil.R b/R/PEC_soil.R index 1227f6a..513136e 100644 --- a/R/PEC_soil.R +++ b/R/PEC_soil.R @@ -43,7 +43,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ #' taken into account for the EFSA 2017 scenarios, using the evaluation depth #' (here mixing depth) as the depth of the layer from which leaching takes #' place. However, as the amount leaching below the evaluation depth -#' (often 5 cm) will partly be mixed back during tillage, the default in this function +#' (often 5 cm) will partly be mixed back during tillage, the default in this function #' is to use the tillage depth for the calculation of the leaching rate. #' @note If temperature information is available in the selected scenarios, as #' e.g. in the EFSA scenarios, the DT50 for groundwater modelling @@ -54,15 +54,16 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ #' @param rate_units Defaults to g/ha #' @param interception The fraction of the application rate that does not reach the soil #' @param mixing_depth Mixing depth in cm -#' @param interval Period of the deeper mixing, defaults to 365, which is a year if -#' degradation rate units are in days +#' @param interval Period of the deeper mixing. The default is NA, i.e. no +#' deeper mixing. For annual deeper mixing, set this to 365 when degradation +#' units are in days #' @param n_periods Number of periods to be considered for long term PEC calculations #' @param PEC_units Requested units for the calculated PEC. Only mg/kg currently supported #' @param PEC_pw_units Only mg/L currently supported #' @param tillage_depth Periodic (see interval) deeper mixing in cm #' @param leaching_depth EFSA (2017) uses the mixing depth (ecotoxicological #' evaluation depth) to calculate leaching for annual crops where tillage -#' takes place. By default, losses from the layer down to the tillage +#' takes place. By default, losses from the layer down to the tillage #' depth are taken into account in this implementation. #' @param cultivation Does mechanical cultivation in the sense of EFSA (2017) #' take place, i.e. twice a year to a depth of 5 cm? Ignored for scenarios @@ -75,6 +76,10 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ #' @param DT50 If specified, overrides soil DT50 endpoints from a chent object #' If DT50 is not specified here and not available from the chent object, zero #' degradation is assumed +#' @param FOMC If specified, it should be a named numeric vector containing +#' the FOMC parameters alpha and beta. This overrides any other degradation +#' endpoints, and the degradation during the interval and after the maximum PEC +#' is calculated using these parameters without temperature correction #' @param Koc If specified, overrides Koc endpoints from a chent object #' @param Kom Calculated from Koc by default, but can explicitly be specified #' as Kom here @@ -116,7 +121,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ #' @export #' @examples #' PEC_soil(100, interception = 0.25) -#' +#' #' # This is example 1 starting at p. 92 of the EFSA guidance (2017) #' # Note that TWA concentrations differ from the ones given in the guidance #' # for an unknown reason (the values from EFSA (2015) can be reproduced). @@ -150,6 +155,7 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0, cultivation = FALSE, chent = NA, DT50 = NA, + FOMC = NA, Koc = NA, Kom = Koc / 1.724, t_avg = 0, scenarios = c("default", "EFSA_2017", "EFSA_2015"), @@ -212,7 +218,7 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0, } else { stop("Depths other than 5 and 20 cm are not supported when using EFSA 2017 scenarios for permanent crops") } - } + } # For the loss via leaching, the equilibrium and therefore the f_om at the # bottom of the layer is probably most relevant. Unfortunately this is not @@ -225,7 +231,7 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0, f_f_om <- subset(f_f_om_depth, bottom == bottom_depth)$f_f_om_no_cultivation } return(f_f_om) - } + } } else { f_f_om_average <- f_f_om_bottom <- function(depth, cultivation) 1 } @@ -292,10 +298,22 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0, f_accu = 0 if (!is.na(interval)) { if (n_periods == Inf) { - f_accu = X/(1 - X) # part of (A2) + if (!is.na(FOMC[1])) { + f_accu = get_vertex(x = 8:10, y = PEC_FOMC_accu_rel(n = 10, interval, FOMC)[8:10])$yv + warning("The long term calculation is based on a pseudo-plateau constructed\n ", + "by fitting a parabola through the residues after 8, 9 and 10 years.\n ", + "This is the method used by the ESCAPE software tool.\n ", + "Please check the validity by specifying e.g. n_periods = 20 or 50.") + } else { + f_accu = X/(1 - X) # part of (A2) + } } else { for (i in 1:n_periods) { - f_accu = f_accu + X^i + if (!is.na(FOMC[1])) { + f_accu = f_accu + 1 / (((i * interval)/FOMC[["beta"]]) + 1)^FOMC[["alpha"]] + } else { + f_accu = f_accu + X^i + } } } } @@ -328,9 +346,45 @@ PEC_soil <- function(rate, rate_units = "g/ha", interception = 0, k_avg <- f_T * k_ref # Leaching not taken into account, EFSA 2017 p. 43 if (t_av_i > 0) { # Equation 10 from p. 24 (EFSA 2015) - result[i, ] <- PEC_soil_sce_mod/(t_av_i * k_avg) * (1 - exp(- k_avg * t_av_i)) # (A8) + if (!is.na(FOMC[1])) { + result[i, ] <- PEC_soil_sce_mod * (FOMC[["beta"]])/(t_av_i * (1 - FOMC[["alpha"]])) * + ((t_av_i/FOMC[["beta"]] + 1)^(1 - FOMC[["alpha"]]) - 1) + } else { + result[i, ] <- PEC_soil_sce_mod/(t_av_i * k_avg) * (1 - exp(- k_avg * t_av_i)) # (A8) + } } } return(result) } + +#' Get the relative accumulation of an FOMC model over multiples of an interval +#' @param n number of applications +#' @param interval Time between applications +#' @param FOMC Named numeric vector containing the FOMC parameters alpha and beta +#' @return A numeric vector containing all n accumulation factors for the n applications +#' @export +PEC_FOMC_accu_rel <- function(n, interval, FOMC) { + PEC_accu_rel <- 0 + for (i in 1:(n - 1)) { + PEC_accu_rel[i + 1] <- PEC_accu_rel[i] + mkin::FOMC.solution(i * interval, 1, + FOMC[["alpha"]], FOMC[["beta"]]) + } + return(PEC_accu_rel) +} + +#' Fit a parabola through three points +#' +#' This was inspired by an answer on stackoverflow +#' https://stackoverflow.com/a/717791 +get_vertex <- function(x, y) { + m <- cbind(x^2, x, 1) + m_inv <- solve(m) + res <- m_inv %*% y + A <- res[1] + B <- res[2] + C <- res[3] + xv = -B / (2*A) + yv = C - B**2 / (4*A) + return(list(xv = xv, yv = yv, A = A, B = B, C = C)) +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 7252624..e5b5705 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -21,6 +21,7 @@ reference: - PEC_soil - soil_scenario_data_EFSA_2015 - soil_scenario_data_EFSA_2017 + - PEC_FOMC_accu - title: Predicted environmental concentrations in groundwater contents: - FOCUS_GW_scenarios_2012 diff --git a/build.log b/build.log index c13de11..eb169b1 100644 --- a/build.log +++ b/build.log @@ -7,5 +7,5 @@ Removed empty directory ‘pfm/inst/testdata/SwashProjects/Project_1/MACRO’ Removed empty directory ‘pfm/inst/testdata/SwashProjects/Project_1’ Removed empty directory ‘pfm/inst/testdata/SwashProjects’ * looking to see if a ‘data/datalist’ file should be added -* building ‘pfm_0.5.1.tar.gz’ +* building ‘pfm_0.5.2.tar.gz’ diff --git a/man/PEC_soil.Rd b/man/PEC_soil.Rd index f9e82e8..c8c93ac 100644 --- a/man/PEC_soil.Rd +++ b/man/PEC_soil.Rd @@ -8,9 +8,9 @@ PEC_soil(rate, rate_units = "g/ha", interception = 0, mixing_depth = 5, PEC_units = "mg/kg", PEC_pw_units = "mg/L", interval = NA, n_periods = Inf, tillage_depth = 20, leaching_depth = tillage_depth, crop = "annual", cultivation = FALSE, chent = NA, DT50 = NA, - Koc = NA, Kom = Koc/1.724, t_avg = 0, scenarios = c("default", - "EFSA_2017", "EFSA_2015"), leaching = scenarios == "EFSA_2017", - porewater = FALSE) + FOMC = NA, Koc = NA, Kom = Koc/1.724, t_avg = 0, + scenarios = c("default", "EFSA_2017", "EFSA_2015"), leaching = scenarios + == "EFSA_2017", porewater = FALSE) } \arguments{ \item{rate}{Application rate in units specified below} @@ -25,8 +25,9 @@ PEC_soil(rate, rate_units = "g/ha", interception = 0, mixing_depth = 5, \item{PEC_pw_units}{Only mg/L currently supported} -\item{interval}{Period of the deeper mixing, defaults to 365, which is a year if -degradation rate units are in days} +\item{interval}{Period of the deeper mixing. The default is NA, i.e. no +deeper mixing. For annual deeper mixing, set this to 365 when degradation +units are in days} \item{n_periods}{Number of periods to be considered for long term PEC calculations} @@ -52,6 +53,11 @@ also be a name for the substance as a character string} If DT50 is not specified here and not available from the chent object, zero degradation is assumed} +\item{FOMC}{If specified, it should be a named numeric vector containing +the FOMC parameters alpha and beta. This overrides any other degradation +endpoints, and the degradation during the interval and after the maximum PEC +is calculated using these parameters without temperature correction} + \item{Koc}{If specified, overrides Koc endpoints from a chent object} \item{Kom}{Calculated from Koc by default, but can explicitly be specified -- cgit v1.2.1