aboutsummaryrefslogtreecommitdiff
path: root/R/PEC_sw_focus.R
blob: 4d5139e79ee19f27c0d70ed8fd258c088533da52 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
#' 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)
}

Contact - Imprint