Skip to contents

Sensitivity Analysis

Source: vignettes/articles/sensitivity.Rmd
sensitivity.Rmd

This report is to document the sensitivity of the TUV model, and the subsequent calculation of water quality guidelines for PAHs, to variation in certain TUV parameters.

Input parameters

The input parameters for the calculation of light attenuation through water for the TUV model are:

  • Vertical Attenuation Coefficient at wavelength wvl: \(kvat = k_de^{-Sk(ref\_wvl-wvl)}\)
    • Kd Light attenuation coefficient (set directly or calculated from DOC)
    • Sk (default 0.018)
    • ref_wvl Reference wavelength for Kd (default 305)
  • depth_m Water depth in m
  • lat Latitude in decimal degrees
  • lon Longitude in decimal degrees
  • elev_m Surface elevation, m above sea level
  • year Year
  • month Month
  • day Day
  • tzone Timezone Local Time - UTC (default 0; UTC)
  • tstart Start time, hours local time (default 0)
  • tstop Stop time, hours local time (default 23)
  • tsteps Number of time steps (default 24)
  • albedo Surface albedo (default 0.07)
  • o3_tc Ozone column, Dobson Units (DU) - Looked up in climatology (recommended default 300)
  • so2_tc SO2 column, DU (default 0)
  • no2_tc NO2 column, DU (default 0)
  • taucld Cloud optical depth (default 0)
  • zbase Cloud base, km (default 4)
  • ztop Cloud top, km (default 5)
  • tauaer Aerosol optical depth at 550 nm - Looked up in climatology (recommended default 0.235)
  • ssaaer Aerosol single scattering albedo (default 0.990)
  • alpha Aerosol Angstrom exponent(default 1.0)
  • wvl_start Starting wavelength, nm (default 279.5)
  • wvl_end End wavelength, nm (default 420.5)
  • wvl_steps Number of wavelength intervals (default 141)
  • nstr TUV run type; use -2 for fast, 4 for slightly more accurate (default -2)

Those currently under consideration for sensitivity analysis are:

  • Kd Light attenuation coefficient - includes use of DOC for calculating Kd
  • albedo Surface albedo (default 0.07)
  • o3_tc Ozone column, Dobson Units (DU) - Looked up in climatology (recommended default 300)
  • tauaer Aerosol optical depth at 550 nm - Looked up in climatology (recommended default 0.235)
  • ssaaer Aerosol single scattering albedo (default 0.990)
  • lat Latitude in decimal degrees
  • elev_m Surface elevation, m above sea level
  • depth_m Water depth in m
  • nstr TUV run type; use -2 for fast, 4 for slightly more accurate (default -2)

The analyses are conducted at three different sample lakes (Southern Interior, Vancouver Island, Northeast), using two PAHs (Anthracene and Benzo[a]pyrene).

For each analysis, we will test a range of reasonable values for the parameter of interest, and plot the PLC50 value calculated across that range. To look at the relative effect of that parameter on the photoxicity of the PAH, we also plot the ratio of PLC50:NLC50.

Initially, sensitivity analyses will be univariate (varying the input of interest while holding the others constant). If there are significant interactions expected between certain variables, these can be explored.

Basic usage of the pahwq package

To start, we will demonstrate a typical straightforward use case at a single site.

To calculate the acute phototoxic water quality guideline (PLC50) for Anthracene at 0.25 m depth in Okanagan Lake on June 21, 2023, with a measured DOC of 5 g/m^3, you would use the following code:

First, set up the options for the model run:

library(pahwq)

set_tuv_aq_params(
  depth_m = 0.25,
  lat = 49.601632,
  lon = -119.605862,
  elev_m = 342,
  DOC = 5,
  date = "2023-06-21",
  tzone = -8,
  albedo = 0.05
)

After setting them, you can view the options that will be used by TUV. Some are set via the function inputs, some are looked up (e.g., o3_tc, tauaer) or calculated (e.g., kd(305) is calculated from the input DOC).

view_tuv_aq_params()
#> a,b,c for: kvdom = a exp(-b(wvl-c)). a = kd(305), b = Sk, c = wavelength, wvl = 305: 20.11 0.018 305
#> ydepth, m: 0.25
#> lat, negative S of Equator: 49.601632
#> lon, negative W of Greenwich (zero) meridian: -119.605862
#> surface elevation, km above sea level: 0.342
#> timezone: Local Time - UTC: -8
#> iyear: 2023
#> imonth: 6
#> iday: 21
#> tstart, hours local time: 0
#> tstop, hours local time: 23
#> number of time steps: 24
#> surface albedo: 0.05
#> o3_tc  ozone column, Dobson Units (DU): 359.937
#> so2_tc SO2 column, DU: 0
#> no2_tc NO2 column, DU: 0
#> taucld - cloud optical depth: 0
#> zbase - cloud base, km: 4
#> ztop - cloud top, km: 5
#> tauaer - aerosol optical depth at 550 nm: 0.0641989811085006
#> ssaaer - aerosol single scattering albedo: 0.99
#> alpha - aerosol Angstrom exponent: 1
#> starting wavelength, nm: 279.5
#> end wavelength, nm: 700.5
#> number of wavelength intervals: 421
#> nstr, use -2 for fast, 4 for slightly more accurate: -2
#> out_irrad_y, T/F, planar spectral irradiance at ydepth: T
#> out_aflux_y, T/F, scalar spectral irradiance (actinic flux)  at depth: F
#> out_irrad_ave, T/F, planar irrad., averaged 0-ydepth: F
#> out_aflux_ave, T/F, scalar, ave 0-ydepth: F
#> out_irrad_atm, T/F, planar, in atmosphere: F
#> out_aflux_atm, T/F, scalar, in atmosphere: F

Run the TUV model and get the results. The results show the underwater irradiance at the specified depth, at each wavelength and time step.

run_tuv()

# Get the results
tuv_res <- get_tuv_results(file = "out_irrad_y")
head(tuv_res)
#>    wl wavelength_start wavelength_end Kd_lambda t_00.00.00 t_01.00.00
#> 1 280            279.5          280.5      31.5          0          0
#> 2 281            280.5          281.5      31.0          0          0
#> 3 282            281.5          282.5      30.4          0          0
#> 4 283            282.5          283.5      29.9          0          0
#> 5 284            283.5          284.5      29.3          0          0
#> 6 285            284.5          285.5      28.8          0          0
#>   t_02.00.00 t_03.00.00 t_04.00.00 t_05.00.00 t_06.00.00 t_07.00.00 t_08.00.00
#> 1          0          0   3.16e-38   1.00e-37   1.90e-37   4.44e-37   7.97e-34
#> 2          0          0   3.21e-35   1.02e-34   1.94e-34   4.63e-34   6.16e-31
#> 3          0          0   3.74e-32   1.18e-31   2.27e-31   5.57e-31   5.39e-28
#> 4          0          0   1.84e-30   5.80e-30   1.12e-29   2.79e-29   2.25e-26
#> 5          0          0   2.19e-28   6.92e-28   1.35e-27   3.47e-27   2.16e-24
#> 6          0          0   1.75e-26   5.52e-26   1.09e-25   2.90e-25   1.36e-22
#>   t_09.00.00 t_10.00.00 t_11.00.00 t_12.00.00 t_13.00.00 t_14.00.00 t_15.00.00
#> 1   3.54e-29   2.17e-26   6.84e-25   2.05e-24   6.97e-25   2.26e-26   3.79e-29
#> 2   1.04e-26   3.65e-24   8.57e-23   2.34e-22   8.72e-23   3.79e-24   1.11e-26
#> 3   3.42e-24   6.74e-22   1.17e-20   2.90e-20   1.19e-20   6.97e-22   3.62e-24
#> 4   8.16e-23   1.16e-20   1.69e-19   3.97e-19   1.72e-19   1.20e-20   8.61e-23
#> 5   3.71e-21   3.41e-19   3.95e-18   8.62e-18   4.00e-18   3.51e-19   3.90e-21
#> 6   1.08e-19   6.33e-18   5.78e-17   1.17e-16   5.85e-17   6.49e-18   1.13e-19
#>   t_16.00.00 t_17.00.00 t_18.00.00 t_19.00.00 t_20.00.00 t_21.00.00 t_22.00.00
#> 1   8.93e-34   4.49e-37   1.91e-37   1.01e-37   3.22e-38          0          0
#> 2   6.83e-31   4.68e-34   1.95e-34   1.02e-34   3.27e-35          0          0
#> 3   5.92e-28   5.63e-31   2.28e-31   1.19e-31   3.81e-32          0          0
#> 4   2.46e-26   2.83e-29   1.13e-29   5.83e-30   1.87e-30          0          0
#> 5   2.34e-24   3.51e-27   1.36e-27   6.97e-28   2.23e-28          0          0
#> 6   1.46e-22   2.94e-25   1.09e-25   5.55e-26   1.78e-26          0          0
#>   t_23.00.00
#> 1          0
#> 2          0
#> 3          0
#> 4          0
#> 5          0
#> 6          0

We can inspect and verify the inputs that were used in the model run:

tuv_run_params(tuv_res)
#> a,b,c for: kvdom = a exp(-b(wvl-c)). a = kd(305), b = Sk, c = wavelength, wvl = 305 
#>                                                                   "20.11 0.018 305" 
#>                                                                           ydepth, m 
#>                                                                              "0.25" 
#>                                                          lat, negative S of Equator 
#>                                                                         "49.601632" 
#>                                        lon, negative W of Greenwich (zero) meridian 
#>                                                                       "-119.605862" 
#>                                               surface elevation, km above sea level 
#>                                                                             "0.342" 
#>                                                          timezone: Local Time - UTC 
#>                                                                                "-8" 
#>                                                                               iyear 
#>                                                                              "2023" 
#>                                                                              imonth 
#>                                                                                 "6" 
#>                                                                                iday 
#>                                                                                "21" 
#>                                                            tstart, hours local time 
#>                                                                                 "0" 
#>                                                             tstop, hours local time 
#>                                                                                "23" 
#>                                                                number of time steps 
#>                                                                                "24" 
#>                                                                      surface albedo 
#>                                                                              "0.05" 
#>                                              o3_tc  ozone column, Dobson Units (DU) 
#>                                                                           "359.937" 
#>                                                               so2_tc SO2 column, DU 
#>                                                                                 "0" 
#>                                                               no2_tc NO2 column, DU 
#>                                                                                 "0" 
#>                                                        taucld - cloud optical depth 
#>                                                                                 "0" 
#>                                                              zbase - cloud base, km 
#>                                                                                 "4" 
#>                                                                ztop - cloud top, km 
#>                                                                                 "5" 
#>                                            tauaer - aerosol optical depth at 550 nm 
#>                                                                "0.0641989811085006" 
#>                                           ssaaer - aerosol single scattering albedo 
#>                                                                              "0.99" 
#>                                                   alpha - aerosol Angstrom exponent 
#>                                                                                 "1" 
#>                                                             starting wavelength, nm 
#>                                                                             "279.5" 
#>                                                                  end wavelength, nm 
#>                                                                             "700.5" 
#>                                                      number of wavelength intervals 
#>                                                                               "421" 
#>                                 nstr, use -2 for fast, 4 for slightly more accurate 
#>                                                                                "-2" 
#>                              out_irrad_y, T/F, planar spectral irradiance at ydepth 
#>                                                                                 "T" 
#>               out_aflux_y, T/F, scalar spectral irradiance (actinic flux)  at depth 
#>                                                                                 "F" 
#>                                out_irrad_ave, T/F, planar irrad., averaged 0-ydepth 
#>                                                                                 "F" 
#>                                            out_aflux_ave, T/F, scalar, ave 0-ydepth 
#>                                                                                 "F" 
#>                                           out_irrad_atm, T/F, planar, in atmosphere 
#>                                                                                 "F" 
#>                                           out_aflux_atm, T/F, scalar, in atmosphere 
#>                                                                                 "F"

Next, calculate the site-specific light absorption (\(P_{abs}\)) for Anthracene from the TUV results. The p_abs() function uses a lookup table to get the molar absorption coefficient table for the specified PAH.

(Pabs <- p_abs(tuv_res, "Anthracene"))
#> [1] 448.6275

Finally, calculate PLC50 in 𝝁g/L, supplying the \(P_{abs}\) value.

plc50(Pabs, pah = "Anthracene")
#> [1] 2.133157

We can compare the PLC50 to the NLC50 to see the effect of the photoxicity of the PAH:

nlc50("Anthracene")
#> [1] 58.40685

A shortcut:

If you don’t need to inspect every step of the way, the above process can be completed in two function calls:

tuv_res <- tuv(
  depth_m = 0.25,
  lat = 49.601632,
  lon = -119.605862,
  elev_m = 342,
  DOC = 5,
  date = "2023-06-21",
  tzone = -8,
  albedo = 0.05
)

plc50(tuv_res, "Anthracene")
#> [1] 2.133157

Sensitivity Analysis: Setup

We’ll start by loading the packages we need and getting some data from B.C. EMS:

Show/Hide Code 
library(sf)
library(mapview)
library(rems)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(latex2exp)

con <- connect_historic_db()
db <- attach_historic_data(con)
sites_sf <- db |> 
  filter(
    EMS_ID %in% c("0500236", "E207466", "0400390"),
    # grepl("(CHARLIE L)|(OKANAGAN)|(QUAMICHAN)", MONITORING_LOCATION), 
    COLLECTION_START > as.Date("2020-01-01"), PARAMETER_CODE == "1103"
  ) |> 
  collect() |> 
  st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = 4326, remove = FALSE) |> 
  filter(!is.na(RESULT)) |> 
  arrange(desc(COLLECTION_START)) |> 
  group_by(EMS_ID, MONITORING_LOCATION, LONGITUDE, LATITUDE) |> 
  mutate(doc_min = min(RESULT), doc_max = max(RESULT)) |> 
  slice(1) |> 
  select(emsid = EMS_ID, 
         name = MONITORING_LOCATION, 
         lon = LONGITUDE, 
         lat = LATITUDE, 
         date = COLLECTION_START, 
         DOC = RESULT,
         doc_min,
         doc_max) |> 
  left_join(
    tribble(
      ~ emsid, ~elev_m,
      "E207466", 25,
      "0400390", 693,
      "0500236", 342
    ), by = join_by("emsid")
  ) |> 
  ungroup() |> 
  relocate(elev_m, .after = lat)

disconnect_historic_db(con)

Locations

Show/Hide Code 
mapview(sites_sf, zcol = "name", legend = FALSE)

Okanagan Lake, Okanagan

  • EMS ID: 0500236
  • Region: Okanagan
  • Latitude: 49.8614 N
  • Longitude: 119.5134 W
  • Site Depth: 70.1 m
  • Maximum Lake Depth: 232 m
  • Lake Elevation: 342 m
  • Lake Surface Area: 350.08 sq km

Charlie Lake, Peace

  • EMS ID: 0400390
  • Region: Peace
  • Latitude: 56.3125 N
  • Longitude: 120.9642 W
  • Site Depth: 13 m
  • Maximum Lake Depth: 13 m
  • Lake Elevation: 693 m
  • Lake Surface Area: 17.56 sq km

Quamichan Lake, Vancouver Island

  • EMS ID: E207466
  • Region: Vancouver Island
  • Latitude: 48.8003 N
  • Longitude: 123.6625 W
  • Site Depth: 7 m
  • Maximum Lake Depth: 8 m
  • Lake Elevation: 25 m
  • Lake Surface Area: 2.88 sq km

Prepare the data for comparisons, including setting the date to be the same

for all sites:

Show/Hide Code 
sites <- sites_sf |> 
  st_drop_geometry() |> 
  mutate(date = "2023-08-01")
gt(sites)
emsid name lon lat elev_m date DOC doc_min doc_max
0400390 CHARLIE L DEEP STATION 1.2 KM EAST OF PARK -120.9642 56.3125 693 2023-08-01 14.00 0.96 15.40
0500236 OKANAGAN L D/S KELOWNA STP (DEEP) -119.5134 49.8614 342 2023-08-01 4.26 4.06 5.17
E207466 QUAMICHAN LAKE; CENTRE -123.6625 48.8003 25 2023-08-01 6.61 6.37 11.80

Basic analysis using defaults

To perform the sensitivity analysis, we need to define a function, multi_plc50() that allows us to do repeated runs of the TUV model and PLC50 calculation using a range of inputs at multiple sites:

Show/Hide Code 
multi_plc50 <- function(df, site = "name", pah, varying, vals = NULL, ...) {
  if (!varying %in% union(names(tuv_aq_defaults()), names(formals(set_tuv_aq_params)))) {
    stop(varying, " is not a valid argument for `set_tuv_aq_params()`")
  }
  
  if (!is.null(vals)) {
    var_df <- data.frame(vals)
    names(var_df) <- varying
    
    df <- df |> 
      select(!any_of(varying)) |> 
      dplyr::cross_join(var_df)
  }
  
  df |> 
    rowwise() |> 
    mutate(
      tuv_res = list(calc_tuv(
          date = date, 
          lat = lat, 
          lon = lon,
          elev_m = elev_m,
          varying = .data[[varying]],
          vary_var = varying,
          ...
        ))
    ) |> 
    cross_join(data.frame(PAH = pah)) |> 
    mutate(
      timing = attr(tuv_res, "timing"),
      Pabs = p_abs(tuv_res, PAH),
      plc50 = plc50(Pabs, pah = PAH),
      nlc50 = nlc50(PAH),
      plc_nlc_ratio = plc50 / nlc50
    ) |> 
    ungroup()
}

calc_tuv <- function(date, lat, lon, elev_m, varying, vary_var, ...) {
  args <- c(
    varying, 
    list(
      depth_m = 0.25,
      date = as.Date(date),
      lat = lat,
      lon = lon,
      elev_m = elev_m
    ),
    ...
  )
  names(args)[1] <- vary_var
  
  # allow overriding of one of the core args by one supplied in 'varying', 
  # this will keep the first of duplicated argument names, which will
  # be the one in 'varying'
  args <- args[unique(names(args))]
  
  do.call("set_tuv_aq_params", args)
  
  t <- system.time(run_tuv(quiet = TRUE))
  res <- get_tuv_results(file = "out_irrad_y")
  attr(res, "timing") <- unname(t["elapsed"])
  res
}

We then use the multi_plc50() function to calculate \(P_{abs}\), PLC50, and ratio of PLC50:NLC50 for Anthracene using recorded DOC values at the three sites. This also uses the utility of pahwq to look up ozone column and aerosol optical depth from climatologies based on latitude, longitude, and month.

Show/Hide Code 
diff_DOC <- multi_plc50(sites, pah = c("Anthracene", "Naphthalene"), varying = "DOC")
emsid name lon lat elev_m date DOC PAH timing Pabs plc50 nlc50 plc_nlc_ratio
0400390 CHARLIE L DEEP STATION 1.2 KM EAST OF PARK -120.9642 56.3125 693 2023-08-01 14.00 Anthracene 0.299 4.678927e+00 12.228584 58.40685 0.20936901
0400390 CHARLIE L DEEP STATION 1.2 KM EAST OF PARK -120.9642 56.3125 693 2023-08-01 14.00 Naphthalene 0.299 3.197590e-07 538.305842 540.10746 0.99666433
0500236 OKANAGAN L D/S KELOWNA STP (DEEP) -119.5134 49.8614 342 2023-08-01 4.26 Anthracene 0.303 5.303841e+02 1.991333 58.40685 0.03409417
0500236 OKANAGAN L D/S KELOWNA STP (DEEP) -119.5134 49.8614 342 2023-08-01 4.26 Naphthalene 0.303 7.389978e-02 328.314614 540.10746 0.60786906
E207466 QUAMICHAN LAKE; CENTRE -123.6625 48.8003 25 2023-08-01 6.61 Anthracene 0.297 1.766459e+02 3.117431 58.40685 0.05337440
E207466 QUAMICHAN LAKE; CENTRE -123.6625 48.8003 25 2023-08-01 6.61 Naphthalene 0.297 4.721068e-03 450.143041 540.10746 0.83343237

To compare the sites using the same input parameters other than location (lat, lon, elevation), we modify the data to set a constant [DOC] (DOC = 5), and calculate \(P_{abs}\), PLC50, and ratio of PLC50:NLC50 for Anthracene using constant DOC = 5:

Show/Hide Code 
same_DOC <- sites |>
  mutate(DOC = 5) |>
  multi_plc50(pah = "Anthracene", varying = "DOC")
emsid name lon lat elev_m date DOC PAH timing Pabs plc50 nlc50 plc_nlc_ratio
0400390 CHARLIE L DEEP STATION 1.2 KM EAST OF PARK -120.9642 56.3125 693 2023-08-01 5 Anthracene 0.297 349.4598 2.362989 58.40685 0.04045740
0500236 OKANAGAN L D/S KELOWNA STP (DEEP) -119.5134 49.8614 342 2023-08-01 5 Anthracene 0.296 375.3538 2.294903 58.40685 0.03929168
E207466 QUAMICHAN LAKE; CENTRE -123.6625 48.8003 25 2023-08-01 5 Anthracene 0.296 376.1315 2.292960 58.40685 0.03925842

Light Attenuation Coefficicent (kd)

The TUV model calculates the light attenuation coefficient \(k_d(λ)\) for each wavelength, based on a reference \(k_d(λ_{ref})\), where \(\lambda_{ref}\) is the reference wavelength.

\(k_d(λ_{ref})\) at 305nm (\(k_{d,305}\)) can be estimated from Dissolved Organic Carbon concentration [DOC] using the equation:

\(k_{d,305} = a_{305}[DOC]^{b_{305}} + 0.13\); \(a_{305}\) = 2.76 and \(b_{305}\) = 1.2 from Morris et al. (1995); Equation 4-1a in Jourabchi (2023)

For the sensitivity analysis for \(k_d\), it is more interpretable to supply a range of \(k_d\) values to the TUV model directly, even though in practice it is more likely to use [DOC] to estimate \(k_d\).

We will test the sensitivity to \(k_d\) using sample values from Table 3 (Jourabchi (2023)), from 0.08 to 275:

Show/Hide Code 
kd_vals <- seq(0.08, 275, length.out = 20)

Effect of varying kd on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
kd_test_a <- multi_plc50(sites, 
                       pah = "Anthracene", 
                       varying = "Kd_ref", 
                       vals = kd_vals, 
                       Kd_wvl = 305,
                       o3_tc = 300,
                       tauaer = 0.235)

ggplot(kd_test_a, aes(x = Kd_ref, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = TeX(r"(PLC50 for Anthracene at 3 sites in B.C. for a range of $k_d(305)$)"), 
    x = TeX(r"($k_d(305)$)"),
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(kd_test_a, aes(x = Kd_ref, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = TeX(r"(Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of $k_d(305)$)"), 
    x = TeX(r"($k_d(305)$)"),
    y = "Ratio of PLC50:NLC50"
  )

For Anthracene, the ratio of PLC50:NLC50 approaches 1 when \(k_d(305)\) is ~ 250 — in other words at that level of \(k_d\), PLC50 is nearly the same as NLC50, indicating that the light is attenuated to such an extent that the phototoxicity of Anthracene is not activated at the light levels when \(k_d(305)\) >= 250.

According to the above equation for estimating \(k_d\) from [DOC]:

\(250 = 2.76[DOC]^{1.23} + 0.13\)

[DOC] ~= 40

Therefore \(k_d(305)\) of ~250 is representative of light attenuation when [DOC] is ~40. This is a significantly higher DOC concentration than observed in our three sample lakes, and is outside the range for which the above equation is recommended (it is recommended for [DOC] between 0.2 and 23 mg/L).

Effect of varying kd on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

We will use the same input values as we did for Anthracene:

Show/Hide Code 
kd_test_b <- multi_plc50(sites, 
                       pah = "Benzo[a]pyrene", 
                       varying = "Kd_ref", 
                       vals = kd_vals,
                       Kd_wvl = 305,
                       o3_tc = 300,
                       tauaer = 0.235)

ggplot(kd_test_b, aes(x = Kd_ref, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = TeX(r"(PLC50 for Benzo\[a\]pyrene at 3 sites in B.C. for a range of $k_d(305)$)"), 
    x = TeX(r"($k_d(305)$)"),
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(kd_test_b, aes(x = Kd_ref, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = TeX(r"(Ratio of PLC50:NLC50 for Benzo\[a\]pyrene at 3 sites in B.C. for a range of $k_d(305)$)"), 
    x = TeX(r"($k_d(305)$)"),
    y = "Ratio of PLC50:NLC50"
  )

Surface Albedo (albedo)

We test surface albedo values from 0.05 to 0.1 (typical for water), and include the suggested default of 0.7:

albedo_vals <- c(seq(0.05, 0.1, length.out = 10), 0.07)

For this and the remaining analyses, we will set DOC = 5 for all sites to ensure a constant \(k_d\) value. We can then test the range of surface albedo values at all three sites for Anthracene and Benzo[a]pyrene.

Effect of varying Surface Albedo on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
albedo_test_a <- multi_plc50(sites, 
                           pah = "Anthracene", 
                           varying = "albedo", 
                           vals = albedo_vals,
                           DOC = 5,
                           o3_tc = 300,
                           tauaer = 0.235)

ggplot(albedo_test_a, aes(x = albedo, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(albedo_test_a, albedo == 0.07),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(albedo_test_a, albedo == 0.07, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.07",
    nudge_x = 0.001,
    nudge_y = 0.05
  ) +
  labs(
    title = "PLC50 for Anthracene at 3 sites in B.C. for a range of Surface Albedo", 
    x = "Surface Albedo",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(albedo_test_a, aes(x = albedo, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(albedo_test_a, albedo == 0.07),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(albedo_test_a, albedo == 0.07, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.07",
    nudge_x = 0.001,
    nudge_y = 0.0005
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of\nSurface Albedo", 
    x = "Surface Albedo",
    y = "Ratio of PLC50:NLC50"
  )

Effect of varying Surface Albedo on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

Show/Hide Code 
albedo_test_b <- multi_plc50(sites, 
                           pah = "Benzo[a]pyrene", 
                           varying = "albedo", 
                           vals = albedo_vals,
                           DOC = 5,
                           o3_tc = 300,
                           tauaer = 0.235)

ggplot(albedo_test_b, aes(x = albedo, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(albedo_test_b, albedo == 0.07),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(albedo_test_b, albedo == 0.07, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.07",
    nudge_x = 0.001,
    nudge_y = 0.0005
  ) +
  labs(
    title = "PLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of Surface Albedo", 
    x = "Surface Albedo",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(albedo_test_b, aes(x = albedo, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(albedo_test_b, albedo == 0.07),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(albedo_test_b, albedo == 0.07, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.07",
    nudge_x = 0.001,
    nudge_y = 0.0003
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nSurface Albedo", 
    x = "Surface Albedo",
    y = "Ratio of PLC50:NLC50"
  )

Ozone Column (o3_tc)

The pahwq package currently uses average monthly ozone column data from 1980-1991 from Fortuin and Kelder (1998), which is bundled with the TUV model. Based on latitude and longitude and month, the ozone column value is looked up and supplied to the TUV model. This default behaviour can be overridden by supplying a value to the o3_tc argument in the set_tuv_aq_params function. In the absence of climatological data, the recommended default value is 300 DU (Jourabchi (2023)).

To test the sensitivity of PLC50 to ozone column, a range of values is tested from 280-420 DU, which are typical values in middle to Northern latitudes. See NASA’s ‘Ozone Watch’ page for more information on atmospheric ozone.

Show/Hide Code 
ozone_vals <- seq(280, 420, by = 10)

Effect of varying Ozone Total Column on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
ozone_test_a <- multi_plc50(sites, 
                          pah = "Anthracene", 
                          varying = "o3_tc", 
                          vals = ozone_vals,
                          DOC = 5,
                          tauaer = 0.235)

ggplot(ozone_test_a, aes(x = o3_tc, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ozone_test_a, o3_tc == 300),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ozone_test_a, o3_tc == 300, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 300",
    nudge_x = 5, 
    nudge_y = 0.03
  ) +
  labs(
    title = "PLC50 for Anthracene at 3 sites in B.C. for a range of Ozone Total Column", 
    x = "Ozone Total Column",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(ozone_test_a, aes(x = o3_tc, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ozone_test_a, o3_tc == 300),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ozone_test_a, o3_tc == 300, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 300",
    nudge_x = 5, 
    nudge_y = 0.0005
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of\nOzone Total Column", 
    x = "Ozone Total Column",
    y = "Ratio of PLC50:NLC50"
  )

Effect of varying Ozone Total Column on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

Show/Hide Code 
ozone_test_b <- multi_plc50(sites, 
                          pah = "Benzo[a]pyrene", 
                          varying = "o3_tc", 
                          vals = ozone_vals,
                          DOC = 5,
                          tauaer = 0.235)

ggplot(ozone_test_b, aes(x = o3_tc, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ozone_test_b, o3_tc == 300),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ozone_test_b, o3_tc == 300, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 300",
    nudge_x = 5, 
    nudge_y = 0.0005
  ) +
  labs(
    title = "PLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of Ozone Total Column", 
    x = "Ozone Total Column",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(ozone_test_b, aes(x = o3_tc, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ozone_test_b, o3_tc == 300),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ozone_test_b, o3_tc == 300, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 300",
    nudge_x = 5, 
    nudge_y = 0.0003
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nOzone Total Column", 
    x = "Ozone Total Column",
    y = "Ratio of PLC50:NLC50"
  )

Aerosol Optical Depth (tauaer)

The pahwq package currently uses average monthly aerosol optical depth data from 2002 to 2023, obtained from MODIS/Aqua satellite data. Based on latitude and longitude and month, the AOD value is looked up and supplied to the TUV model. This default behaviour can be overridden by supplying a value to the tauaer argument in the set_tuv_aq_params function. In the absence of climatological data, the recommended default value is 0.235 (Jourabchi (2023)).

To test the sensitivity of PLC50 to AOD, a range of values is tested from 0.1 (clear skies with maximum visibility) to 1.0 (very hazy skies). See the NASA Earth Observatory page on aerosols.

Show/Hide Code 
aod_vals <- c(seq(0.1, 1.0, length.out = 10), 0.235)

Effect of Varying Aerosol Optical Depth (AOD) on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
aod_test_a <- multi_plc50(sites, 
                        pah = "Anthracene", 
                        varying = "tauaer", 
                        vals = aod_vals,
                        DOC = 5,
                        o3_tc = 300)

ggplot(aod_test_a, aes(x = tauaer, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(aod_test_a, tauaer == 0.235),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(aod_test_a, tauaer == 0.235, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.235",
    nudge_y = 0.1
  ) +
  labs(
    title = "PLC50 for Anthracene at 3 sites in B.C. for a range of Aeorosol Optical Depth", 
    x = "Aeorosol Optical Depth",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(aod_test_a, aes(x = tauaer, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(aod_test_a, tauaer == 0.235),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(aod_test_a, tauaer == 0.235, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.235",
    nudge_y = 0.002
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of\nAeorosol Optical Depth", 
    x = "Aeorosol Optical Depth",
    y = "Ratio of PLC50:NLC50"
  )

Effect of varying Aerosol Optical Depth (AOD) on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

Show/Hide Code 
aod_test_b <- multi_plc50(sites, 
                        pah = "Benzo[a]pyrene", 
                        varying = "tauaer", 
                        vals = aod_vals,
                        DOC = 5,
                        o3_tc = 300)

ggplot(aod_test_b, aes(x = tauaer, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(aod_test_b, tauaer == 0.235),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(aod_test_b, tauaer == 0.235, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.235",
    nudge_y = 0.002
  ) +
  labs(
    title = "PLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nAeorosol Optical Depth", 
    x = "Aeorosol Optical Depth",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(aod_test_b, aes(x = tauaer, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(aod_test_b, tauaer == 0.235),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(aod_test_b, tauaer == 0.235, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.235",
    nudge_y = 0.001
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nAeorosol Optical Depth", 
    x = "Aeorosol Optical Depth",
    y = "Ratio of PLC50:NLC50"
  )

Aerosol Single Scattering Albedo (ssaaer)

Aerosol Single Scattering Albedo is currently set as a default constant value of 0.99. This can be overridden by setting the ssaaer parameter of the set_tuv_aq_params() function to a different value.

Here we test a range for values from 0.80 to 0.99, including the recommended default of 0.99 (Jourabchi (2023)).

Show/Hide Code 
ssa_vals <- seq(0.8, 0.99, length.out = 10)

Effect of varying Aerosol Single Scattering Albedo on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
ssa_test_a <- multi_plc50(sites, 
                        pah = "Anthracene", 
                        varying = "ssaaer", 
                        vals = ssa_vals,
                        DOC = 5,
                        o3_tc = 300,
                        tauaer = 0.235)

ggplot(ssa_test_a, aes(x = ssaaer, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ssa_test_a, ssaaer == 0.99),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ssa_test_a, ssaaer == 0.99, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.99",
    nudge_y = .06
  ) +
  labs(
    title = "PLC50 for Anthracene at 3 sites in B.C. for a range of\nAeorosol Single Scattering Albedo", 
    x = "Aeorosol Single Scattering Albedo",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(ssa_test_a, aes(x = ssaaer, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ssa_test_a, ssaaer == 0.99),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ssa_test_a, ssaaer == 0.99, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.99",
    nudge_y = .001
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of\nAeorosol Single Scattering Albedo", 
    x = "Aeorosol Single Scattering Albedo",
    y = "Ratio of PLC50:NLC50"
  )

Effect of varying Aerosol Single Scattering Albedo on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

Show/Hide Code 
ssa_test_b <- multi_plc50(sites, 
                        pah = "Benzo[a]pyrene", 
                        varying = "ssaaer", 
                        vals = ssa_vals,
                        DOC = 5,
                        o3_tc = 300,
                        tauaer = 0.235)

ggplot(ssa_test_b, aes(x = ssaaer, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ssa_test_b, ssaaer == 0.99),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ssa_test_b, ssaaer == 0.99, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.99",
    nudge_y = .001
  ) +
  labs(
    title = "PLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nAeorosol Single Scattering Albedo", 
    x = "Aeorosol Single Scattering Albedo",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(ssa_test_b, aes(x = ssaaer, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  geom_point(
    data = filter(ssa_test_b, ssaaer == 0.99),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(ssa_test_b, ssaaer == 0.99, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.99",
    nudge_y = .0005
  ) +
  labs(
    title = "Ratio of PLC50:NLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nAeorosol Single Scattering Albedo", 
    x = "Aeorosol Single Scattering Albedo",
    y = "Ratio of PLC50:NLC50"
  )

Latitude (lat)

Latitude has a strong effect on the angle of the sun, which will in turn affect the light penetration through water, thus it is important to investigate the sensitivity of the calculation of PLC50 to variation in latitude.

For this analysis, we will maintain the longitudes associated with each site, and make “mock” sites by creating a range of latitudes from 48 to 72 degrees N and pair those with the longitudes of the original sites.

lat_vals <- seq(48, 72, by = 2)

Effect of varying Latitude on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
lat_test_a <- multi_plc50(sites, 
                        pah = "Anthracene", 
                        varying = "lat", 
                        vals = lat_vals,
                        DOC = 5,
                        o3_tc = 300,
                        tauaer = 0.235)

ggplot(lat_test_a, aes(x = lat, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(lon), ncol = 1, labeller = as_labeller(\(x) paste("Longitude: ", x))) +
  labs(
    title = "PLC50 for Anthracene at 3 sites in B.C. for a range of latitudes", 
    x = "Latitude",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(lat_test_a, aes(x = lat, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(lon), ncol = 1, labeller = as_labeller(\(x) paste("Longitude: ", x))) +
  labs(
    title = "Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of latitudes",
    x = "Latitude",
    y = "Ratio of PLC50:NLC50"
  )

Effect of varying Latitude on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

Show/Hide Code 
lat_test_b <- multi_plc50(sites, 
                        pah = "Benzo[a]pyrene", 
                        varying = "lat", 
                        vals = lat_vals,
                        DOC = 5,
                        o3_tc = 300,
                        tauaer = 0.235)

ggplot(lat_test_b, aes(x = lat, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(lon), ncol = 1, labeller = as_labeller(\(x) paste("Longitude: ", x))) +
  labs(
    title = "PLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of latitudes", 
    x = "Latitude",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(lat_test_b, aes(x = lat, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(lon), ncol = 1, labeller = as_labeller(\(x) paste("Longitude: ", x))) +
  labs(
    title = "Ratio of PLC50:NLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nlatitudes",
    x = "Latitude",
    y = "Ratio of PLC50:NLC50"
  )

Elevation (elev_m)

We can test the effect of elevation on PLC50 by varying elevation from sea level to 1000m above sea level. We will keep the latitude and longitude associated with each site.

elevation_vals <- seq(0, 1000, by = 100)

Effect of varying Elevation on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
elev_test_a <- multi_plc50(sites, 
                         pah = "Anthracene", 
                         varying = "elev_m", 
                         vals = elevation_vals,
                         DOC = 5,
                         o3_tc = 300,
                         tauaer = 0.235)

ggplot(elev_test_a, aes(x = elev_m, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "PLC50 for Anthracene at 3 sites in B.C. for a range of elevations", 
    x = "Elevation (m)",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(elev_test_a, aes(x = elev_m, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of\nelevations",
    x = "Elevation (m)",
    y = "Ratio of PLC50:NLC50"
  )

Effect of varying Elevation on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

Show/Hide Code 
elev_test_b <- multi_plc50(sites, 
                         pah = "Benzo[a]pyrene", 
                         varying = "elev_m", 
                         vals = elevation_vals,
                         DOC = 5,
                         o3_tc = 300,
                         tauaer = 0.235)

ggplot(elev_test_b, aes(x = elev_m, y = plc50)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "PLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of elevations", 
    x = "Elevation (m)",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(elev_test_b, aes(x = elev_m, y = plc_nlc_ratio)) +
  geom_point() + 
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "Ratio of PLC50:NLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nelevations",
    x = "Elevation (m)",
    y = "Ratio of PLC50:NLC50"
  )

Water depth (depth_m)

The depth at which the photoxicity of a PAH is determined should be a conservative one - i.e., a depth at which sensitive species/life stages exist and thus would be exposed. The default is set at 0.25m.

We can test the effect of depth on PLC50 by varying elevation from water level to 1m below the surface:

depth_vals <- c(seq(0, 1, by = 0.1), 0.25)

Effect of varying Water Depth on PLC50 and ratio of PLC50:NLC50 for Anthracene

Show/Hide Code 
depth_test_a <- multi_plc50(sites, 
                          pah = "Anthracene", 
                          varying = "depth_m", 
                          vals = depth_vals,
                          DOC = 5,
                          o3_tc = 300,
                          tauaer = 0.235)

ggplot(depth_test_a, aes(x = depth_m, y = plc50)) +
  geom_point() + 
  geom_point(
    data = filter(depth_test_a, depth_m == 0.25),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(depth_test_a, depth_m == 0.25, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.25m",
    nudge_x = -0.05,
    nudge_y = 5
  ) +
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "PLC50 for Anthracene at 3 sites in B.C. for a range of water depths", 
    x = "Water Depth (m)",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(depth_test_a, aes(x = depth_m, y = plc_nlc_ratio)) +
  geom_point() + 
  geom_point(
    data = filter(depth_test_a, depth_m == 0.25),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(depth_test_a, depth_m == 0.25, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.25m",
    nudge_x = -0.05,
    nudge_y = 0.1
  ) +
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "Ratio of PLC50:NLC50 for Anthracene at 3 sites in B.C. for a range of\nwater depths",
    x = "Water Depth (m)",
    y = "Ratio of PLC50:NLC50"
  )

Effect of varying Water Depth on PLC50 and ratio of PLC50:NLC50 for Benzo[a]pyrene

Show/Hide Code 
depth_test_b <- multi_plc50(sites, 
                          pah = "Benzo[a]pyrene", 
                          varying = "depth_m", 
                          vals = depth_vals,
                          DOC = 5,
                          o3_tc = 300,
                          tauaer = 0.235)

ggplot(depth_test_b, aes(x = depth_m, y = plc50)) +
  geom_point() + 
  geom_point(
    data = filter(depth_test_b, depth_m == 0.25),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(depth_test_b, depth_m == 0.25, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.25m",
    nudge_x = -0.05,
    nudge_y = 0.1
  ) +
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "PLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of water depths", 
    x = "Water Depth (m)",
    y = "PLC50 (μg/L)"
  )
Show/Hide Code 

ggplot(depth_test_b, aes(x = depth_m, y = plc_nlc_ratio)) +
  geom_point() + 
  geom_point(
    data = filter(depth_test_b, depth_m == 0.25),
    colour = "red"
  ) + 
  geom_text_repel(
    data = filter(depth_test_b, depth_m == 0.25, name == "QUAMICHAN LAKE; CENTRE"),
    colour = "red", 
    label = "Recommended default = 0.25m",
    nudge_x = -0.05,
    nudge_y = 0.05
  ) +
  facet_wrap(vars(name), ncol = 1) +
  labs(
    title = "Ratio of PLC50:NLC50 for Benzo[a]pyrene at 3 sites in B.C. for a range of\nwater depths",
    x = "Water Depth (m)",
    y = "Ratio of PLC50:NLC50"
  )

Radiative transfer scheme (nstr)

Number of streams for radiative transfer calculations:

  • If nstr < 2, uses 2-stream delta-Eddington (faster)
  • if nstr >= 2, uses n-stream discrete ordinates (more accurate: must be even number, maximum = 32)
Show/Hide Code 
out <- multi_plc50(sites, 
          pah = "Anthracene", 
          varying = "nstr", 
          vals = c(-2, 4), 
          DOC = 5,
          o3_tc = 300,
          tauaer = 0.235) |> 
  mutate(nstr = factor(nstr))

ggplot(out, aes(x = name, y = plc50, colour = nstr)) + 
  geom_point() + 
  labs(
    title = "Calculated PLC50 using two different TUV methods for radiative transfer", 
    x = "Site"
  ) + 
  theme(axis.text.x = element_text(angle = 355))
Show/Hide Code 

ggplot(out, aes(x = name, y = timing, colour = nstr)) + 
  geom_point() + 
  labs(
    title = "Execution time of TUV model runs with different methods for radiative transfer", 
    x = "Site",
    y = "Execution Time (s)"
  ) + 
  theme(axis.text.x = element_text(angle = 355))

We can see that there are small differences in the calculated PLC50 when using different values of nstr, with the values when calculated using nstr = -2 being on average 0.61% lower than when using nstr = 4. Using nstr = 4 however, is about 8 times slower than using nstr = -2. In most cases, unless high precision is required, it is likely that using the much faster nstr = -2 should be the default.

Methylated Polycyclic Aromatic Hydrocarbons

Absorption spectra for methylated PAHs are not all available from published sources, so it would be beneficial to be able to approximate them for calculating water quality guidelines. In this section we will look at the feasibilty of using the absorbance spectra from a methylated PAH’s parent compound as a proxy for the specific spectra. We will do this by calculating the PLC50 for various methylated PAHs for which we have absorbance spectra, and comparing these values to those calculated using the absorbance spectra of the parent PAH.

For these comparisons we will limit it to just one site in the Okanagan.

okanagan <- filter(sites, emsid == "0500236")
Show/Hide Code 
# Define some functions:
ma_plot <- function(chemicals) {
  pahwq:::molar_absorption |> 
    filter(chemical %in% tolower(chemicals)) |> 
    mutate(
      chemical = factor(chemical, levels = tolower(chemicals))
    ) |> 
    ggplot(aes(x = wavelength, y = molar_absorption, colour = chemical)) + 
    geom_point() + 
    labs(
      title = paste0("Molar absorption of ", chemicals[1], " and methylated derivatives"), 
      y = "Molar absorption (L/mol/cm)"
    )
}

nlc50_plot <- function(chemicals) {
  pahwq:::nlc50_lookup |> 
    filter(chemical %in% tolower(chemicals)) |> 
    mutate(
      nlc50 = vapply(chemical, nlc50, FUN.VALUE = 1), 
      chemical = factor(chemical, levels = tolower(chemicals))
    ) |> 
    ggplot(aes(x = log_kow, y = nlc50, colour = chemical)) + 
    geom_point() + 
    labs(
      title = paste0("NLC50 of ", chemicals[1], " and its methylated derivatives"),
      y = "NLC50 (ug/L)"
    )
}

plc50_surrogates <- function(chemicals) {
  multi_plc50(
  okanagan,
  site = "name",
  pah = chemicals,
  varying = "Kd_ref",
  vals = c(1, 150),
  depth_m = 0.25
) |> 
  group_by(Kd_ref) |> 
  rowwise() |> 
  mutate(
    Pabs_parent = p_abs(tuv_res, tolower(chemicals[1])),
    plc50_parent = plc50(Pabs_parent, PAH)
  ) |> 
  ungroup() |> 
  mutate(
    PAH = factor(PAH, levels = chemicals)
  ) |> 
  pivot_longer(cols = c(plc50, plc50_parent), 
               names_to = "abs_spectra", 
               values_to = "plc50") |> 
  mutate(
    plc_nlc_ratio = plc50 / nlc50,
    abs_spectra = case_when(
      abs_spectra == "plc50" ~ "specific",
      abs_spectra == "plc50_parent" ~ "parent"
    )
  )
}

plc50_plot <- function(df) {

ggplot(df, aes(x = PAH, y = plc50, colour = abs_spectra)) + 
  geom_point() + 
  facet_wrap(vars(Kd_ref), ncol = 1, scales = "free_y", 
             labeller = as_labeller(function(x) paste0("Kd(305) = ", x))) + 
  geom_text_repel(
    data = filter(df, PAH != chemicals[1]),
    mapping = aes(label = round(plc_nlc_ratio, 3)),  
    show.legend = FALSE
  ) +
  labs(
    title = paste0("Comparison of PLC50 of methylated ", chemicals[1], ", using\nspecific absorption spectra vs absorption spectra of parent compound"),
    y = "PLC50 (ug/L)",
    colour = "Absorption spectra used",
    caption = "*Text labels indicate the ratio of PLC50:NLC50"
  )
}

Napthalene

To start, we can examine the molar absorption specra of three methylated Naphthalene derivatives, as compared to unmethylated Naphthalene:

Show/Hide Code 
chemicals <- c("Naphthalene", "1-Methylnaphthalene", "2,6-Dimethylnaphthalene", "1,6,7-Trimethylnaphthalene")

ma_plot(chemicals)

We can see in general that more highly methylated compounds have greater absorption, though it is not consistent across the spectrum.

We can also examine the properties and narcotic toxicity (NLC50) of Naphthalene and its methylated variants:

Show/Hide Code 
nlc50_plot(chemicals)

We can see here that 1,6,7-Trimethylnaphthalene has the highest log(KoW) and lowest NLC50.

Finally, we can compare the phototoxicity (PLC50) of Naphthalene and its methylated derivatives at the Okanagan site, at the two values of \(k_{d,305}\) (1 and 150).

Show/Hide Code 
df <- plc50_surrogates(chemicals)
plc50_plot(df) + 
  theme(axis.text.x = element_text(angle = 355))

Here we can see that a lower PLC50 (higher phototoxicity) is calculated for all methylated naphthalenes when the absorption spectra of the methylated compound is used, as compared to when that of the parent compound (Naphthalene) is used. This difference is much greater when there is higher light penetration (lower Kd) - which is to be expected, as when there is less light penetration the PLC50 approaches the NLC50 so the absorption spectra of the chemical becomes irrelevant.

Phenanthrene

As above, we can visualize the absorption spectra of Phenanthrene and 9-Methylphenanthrene:

Show/Hide Code 
chemicals <- c("Phenanthrene", "9-Methylphenanthrene")

ma_plot(chemicals)

Here we see that the parent compound Phenanthrene has much higher absorption than 9-Methylphenanthrene at the low end of the light spectrum.

Show/Hide Code 
nlc50_plot(chemicals)

9-Methylphenanthrene has higher narcotic toxicity (lower NLC50) than its parent Phenanthrene.

Show/Hide Code 
plc50_df_phenanthrene <- plc50_surrogates(chemicals)
plc50_plot(plc50_df_phenanthrene)

Using the parent molar absorption spectra for 9-Methylphenanthrene gives a more conservative PLC50 than using the specific spectra for the methylated compound. This is different than the case for Naphthalene.

Similar to Naphthalene, this only appears to be substantially different at high light penetration (low Kd(305)).

Anthracene

As above, we can visualize the absorption spectra of Anthracene and 9-Methylanthracene:

Show/Hide Code 
chemicals <- c("Anthracene", "9-Methylanthracene")
ma_plot(chemicals)

The absorption spectra for both Anthracene and 9-Methylanthracene is very broad compared to that of the other PAH compounds.

Show/Hide Code 
nlc50_plot(chemicals)

9-Methylanthracene has higher narcotic toxicity (lower NLC50) than its parent Anthracene.

Show/Hide Code 
plc50_df_anthracene <- plc50_surrogates(chemicals)
plc50_plot(plc50_df_anthracene) + 
  theme(axis.text.x = element_text(angle = 355))

Using the parent absorption spectra for 9-Methylanthracene results in a less conservative PLC50 than if the specific absorption spectra is used. Contrary to that seen in the other PAHs above, this effect is still quite pronounced at high light attenuation (Kd(305) = 150), presumably because of the broad absorption spectra exhibited by both Anthracene and 9-Methylanthracene.

Pyrene

Show/Hide Code 
chemicals <- c("Pyrene", "1-Methylpyrene")

ma_plot(chemicals)

Show/Hide Code 
nlc50_plot(chemicals)

9-Methylpyrene has higher narcotic toxicity (lower NLC50) than its parent Pyrene.

Show/Hide Code 
plc50_df_pyrene <- plc50_surrogates(chemicals)
plc50_plot(plc50_df_pyrene)

Using the parent absorption spectra for 9-Methylpyrene results in a less conservative PLC50 than if the specific absorption spectra is used.

Chrysene

Show/Hide Code 
chemicals <- c("Chrysene", "1-Methylchrysene")
ma_plot(chemicals)

Show/Hide Code 
nlc50_plot(chemicals)

Similar to the other methylated PAHs, 1-Methylchrysene has a lower NLC50 (higher narcotic toxicity) than its parent Chrysene.

Show/Hide Code 
plc50_df_chrysene <- plc50_surrogates(chemicals)
plc50_plot(plc50_df_chrysene)

The results for Chrysene and 1-Methylchrysene show that the effect of using the absorbance spectra of a parent compound can be different depending on the light intensity/attenuation. At low light attenuation (Kd(305) = 1), the specific spectra gives a more conservative PLC50 for 1-Methylchrysene than that obtained using the spectra from Chrysene, while the opposite is true at high light attenuation (Kd(305) = 100).

These results show that substituting the absorption spectra of a parent compound when that of a specific methylated PAH is unknown is possible, but also somewhat complicated. Differences in absorption spectra between methylated PAHs and that of their parent was quite variable across compounds. The difference in the resulting calculated PLC50 is generally quite small and thus it is likely reasonable to use the absorption spectra of a parent when it is not available for a methylated PAH. However, the magnitude and even the sign of the difference is highly variable. In all cases examined, the NLC50 of the methylated compounds was lower than that of the parent compounds (i.e., they are more toxic), however the final phototoxicity (PLC50) depended on the differences in the absorption spectra and, in some cases, the level of light attenuation.

References

Fortuin, J. Paul F., and Hennie Kelder. 1998. “An Ozone Climatology Based on Ozonesonde and Satellite Measurements.” Journal of Geophysical Research: Atmospheres 103 (D24): 31709–34. https://doi.org/https://doi.org/10.1029/1998JD200008.
Jourabchi, P. 2023. “Report on Estimating Light Attenuation in Support of the b.c. Water Quality Guideline Development for Phototoxic PAHs.” Vancouver, B.C.: B.C. Ministry of Water, Land; Resource Stewardship, prepared by ARIS Environmental Ltd.
Morris, Donald P., Horatio Zagarese, Craig E. Williamson, Esteban G. Balseiro, Bruce R. Hargreaves, Beatriz Modenutti, Robert Moeller, and Claudia Queimalinos. 1995. “The Attenuation of Solar UV Radiation in Lakes and the Role of Dissolved Organic Carbon.” Limnology and Oceanography 40 (8): 1381–91. https://doi.org/https://doi.org/10.4319/lo.1995.40.8.1381.