Average light environment negative in certain case

[itowers1]

Hi team!

Finding certain cases where average light environment is less than 0, likely due to the interpolator, but unclear to me why this would be.

Checking out the simple_water_use_model branch, and running the following code, reproduces the error. Alternatively, repex_canopy_issue.R in example_vignettes can also be used (has the below code).

devtools::load_all(".")

library(tidyverse)

p0 <- scm_base_parameters("FF16w")

p_50_2_K_s <- function(p50){
  1.638999*(p50/2)^(-1.38)
}

calc_vul_b <- function(p_50, c){
  num <- p_50
  den <- (-log(1-50/100))^(1/c)
  num/den
}

calc_psi_crit <- function(b,c) {
  b*(log(1/0.05))^(1/c)
}

# set default values 
run_FF16w_model <- function(inputs){

  ctrl = scm_base_control()

  for(v in names(ctrl)){
    if(length(inputs[[v]]) > 0){
      if(!is.na(inputs[[v]])){
        ctrl[[v]] <- inputs[[v]]
      }
    }
  }

  p0 <- scm_base_parameters("FF16w")
  p0$max_patch_lifetime <- inputs$max_patch_lifetime

  p50_1 <- inputs$p50_1
  p50_2 <- inputs$p50_2

  if(is.na(inputs$epsilon_leaf)){
    epsilon_leaf = 0.00001
  } else{
    epsilon_leaf = inputs$epsilon_leaf
  }

  ## increase birth rates to 10 or 50

  p1 <- expand_parameters(trait_matrix(c(p50_1, p50_2), "p_50"), p0, mutant=FALSE, birth_rate_list = c(1,1))

  p1$strategies[[1]]$K_s <- p_50_2_K_s(p50_1)
  p1$strategies[[1]]$b <- calc_vul_b(p50_1, p1$strategies[[1]]$c)
  p1$strategies[[1]]$psi_crit <- calc_psi_crit(p1$strategies[[1]]$b, p1$strategies[[1]]$c)
  p1$strategies[[1]]$epsilon_leaf <- epsilon_leaf
  p1$strategies[[1]]$recruitment_decay <- inputs$recruitment_decay

  p1$strategies[[2]]$K_s <- p_50_2_K_s(p50_2)
  p1$strategies[[2]]$b <- calc_vul_b(p50_2, p1$strategies[[2]]$c)
  p1$strategies[[2]]$psi_crit <- calc_psi_crit(p1$strategies[[2]]$b, p1$strategies[[2]]$c)
  p1$strategies[[2]]$epsilon_leaf <- epsilon_leaf
  p1$strategies[[2]]$epsilon_leaf <- inputs$recruitment_decay

  env <- make_environment("FF16w", soil_initial_state = rep(inputs$initial_theta, 1), rainfall = inputs$rainfall)
  time1 <- system.time(result <- build_schedule(p1, env, ctrl))
  time2 <- system.time(result <- run_scm_collect(result, env, ctrl, collect_auxiliary_variables = TRUE))
  return(list(result, time1, time2))
}

poss_run_plant = safely(.f = run_FF16w_model, otherwise = "Error")

inputs <- expand_grid(rainfall = 2.93, p50_1 = 5.99, p50_2 = 8.46, initial_theta = c(0.1), max_patch_lifetime = 11, schedule_eps = NA, ode_tol_abs = NA, ode_tol_rel = NA, epsilon_leaf = NA, recruitment_decay = c(1)) %>%
  slice(sample(1:n()))

run_and_save_plant_water_model <- function(..., save = FALSE, dir = NULL){
  input_parameters <- tibble(...)

  if(save == TRUE){
  if(!file.exists(paste("examples_vignettes/outputs/", dir, sep = ""))){
  browser()

    dir.create(paste("examples_vignettes/outputs/", dir, sep = ""))

  }

  dir <- paste("examples_vignettes/outputs/", dir, "/", sep = "")

  filename <-paste("rainfall_", round(input_parameters$rainfall,2), 
                   "p501_", round(input_parameters$p50_1, 2), 
                   "p502_", round(input_parameters$p50_2, 2), 
                   "lifetime_", input_parameters$max_patch_lifetime,
                   "schedule_eps", round(input_parameters$schedule_eps, 6),
                   "ode_tol_abs", round(input_parameters$ode_tol_abs, 6),
                   "ode_tol_rel", round(input_parameters$ode_tol_rel, 6),
                   "epsilon_leaf", round(input_parameters$epsilon_leaf, 6),
                   "theta", input_parameters$initial_theta,
                   "recruitment_decay", input_parameters$recruitment_decay,
                   ".RDS", sep="")

  path <- paste(dir, filename, sep="")

  if(!file.exists(path)){
  result <- poss_run_plant(input_parameters)
  result<-list(result, input_parameters)
  saveRDS(result, path)
  }
  else {
    print(paste("skipping", path, "already done", sep = " "))
  }
  } else {

  result <- poss_run_plant(input_parameters)
  result<-list(result, input_parameters)

  return(result)
  }

  return
}

result <- pmap(inputs, run_and_save_plant_water_model, save = FALSE, dir = "recruitment_dynamics")

[itowers1]

Note that I've only recently come across this issue after setting recruitment_decay to 1, an equivalent example with recruitment_decay as 0 runs successfully.

[dfalster]

Confirming that I can reproduce the error with the code above.

I've looked more into the cause, the negative light numbers are coming about because there is a high density of leaf area just above this plant (i.e. just above 4m), causing canopy openness to drop from around 0.3 to near zero in a short height interval. Recall light environment is modelled with a spline. This spline is struggling with the transition to a value near zero. It wiggles around ever so much.

We've worried about this previously, but never so much do change anything.

Several possible avenues forward here

• When using the light environment, always insist it's > 0
• Ensure photosynthesis algorithm is robust against PPFD values close to zero, even if slightly negative
• A more conservative approach would be to setup the spline interpolator based on total leaf area above a height rather than canopy openness. (see Canopy.h