R/onelump_var.R

In mrke/NicheMapR: R implementation of Niche Mapper software for biophysical modelling

#' One-lump Transient Heat Budget (for use with deSolve package)
#'
#' Transient, 'one-lump', heat budget for computing rate of change of temperature
#' under environmental conditions that vary with time, using interpolation functions to
#' estimate environmental conditions at particular time intervals.
#' Michael Kearney, Raymond Huey and Warren Porter developed this R function and example in September 2017.
#' @param t time intervals (s) at which output is required
#' @param Tc_init = 5, initial temperature (°C) Organism shape, 0-5, Determines whether standard or custom shapes/surface area/volume relationships are used: 0=plate, 1=cyl, 2=ellips, 3=lizard (desert iguana), 4=frog (leopard frog), 5=custom (see details)
#' @param Ww_g = 500 weight (g)
#' @param rho_body = 1000, animal density (kg/m3)
#' @param q = 0, metabolic rate (W/m3)
#' @param c_body = 3073, specific heat of flesh (J/kg-C)
#' @param k_flesh = 0.5, conductivity of flesh (W/mK)
#' @param emis = 0.95, emissivity of skin (-)
#' @param alpha = 0.85, animal solar absorptivity (-)
#' @param geom = 2, organism shape, 0-5, Determines whether standard or custom shapes/surface area/volume relationships are used: 0=plate, 1=cyl, 2=ellips, 3=lizard (desert iguana), 4=frog (leopard frog), 5=custom (see parameter 'shape_coeffs')
#' @param shape_b = 1/5, proportionality factor (-) for going from volume to area, represents ratio of width:height for a plate, length:diameter for cylinder, b axis:a axis for ellipsoid
#' @param shape_c = 1/5, proportionality factor (-) for going from volume to area, represents ratio of length:height for a plate, c axis:a axis for ellipsoid
#' @param shape_coefs = c(10.4713,.688,0.425,0.85,3.798,.683,0.694,.743), Custom surface area coefficients. Operates if posture = 5, and consists of 4 pairs of values representing the parameters a and b of a relationship AREA=a*Ww_g^b, where AREA is in cm2 and Ww_g is in g. The first pair are a and b for total surface area, then a and b for ventral area, then for silhouette area normal to the sun, then silhouette area perpendicular to the sun
#' @param posture = 'n', pointing normal 'n' or parallel 'p' to the sun's rays, or average 'a'?
#' @param orient = 1, does the object orient toward the sun? (0,1)
#' @param fatosk = 0.4, solar configuration factor to sky (-)
#' @param fatosb = 0.4, solar configuration factor to substrate (-)
#' @param dyn_Q = 0, dynamic metabolic heat generation as a function of temperature, based on approxfun called qf (1) or constant based on parameter q (0)
#' @param alpha_sub = 0.2, substrate solar reflectivity, decimal percent
#' @param pdif = 0.1, proportion of solar energy that is diffuse (rather than direct beam)
#' @param fluid = 0, fluid type, air (0) or water (1)
#' @param Tairf air temperature function with time, generated by 'approxfun' (°C)
#' @param Tradf radiant temperature function with time, generated by 'approxfun'(°C), averaging ground and sky
#' @param velf wind speed function with time, generated by 'approxfun' (m/s)
#' @param Qsolf radiation function with time, generated by 'approxfun' (W/m2)
#' @param Zenf zenith angle of sun function with time, generated by 'approxfun' (90 is below horizon), degrees
#' @param press air pressure (Pa)
#' @return Tc Core temperature (°C)
#' @return Tcf Final (steady state) temperature (°C), if conditions remained constant indefinitely
#' @return tau Time constant (s)
#' @return dTc Rate of change of core temperature (°C/s)
#' @usage ode(y = Tb_init, t = times, func = onelump_var, parms = indata)
#' @examples
#' library(deSolve) # note due to some kind of bug in deSolve, it must be loaded before NicheMapR!
#' library(NicheMapR)
#'
#' # get microclimate data
#' loc <- c(133.8779, -23.6987) # Alice Springs, Australia
#' Usrhyt <- 0.05 # height of midpoint of animal, m
#' micro <- micro_global(loc = loc, Usrhyt = Usrhyt) # run the model with specified location and animal height
#' metout <- as.data.frame(micro$metout) # above ground microclimatic conditions, min shade
#' soil <- as.data.frame(micro$soil) # soil temperatures, minimum shade
#'
#' # append dummy dates
#' days <- rep(seq(1, 12), 24)
#' days <- days[order(days)]
#' dates <- days + metout$TIME / 60 / 24 - 1 # dates for hourly output
#' dates2 <- seq(1, 12, 1) # dates for daily output
#' metout <- cbind(dates, metout)
#' soil <- cbind(dates, soil)
#'
#' # combine relevant input fields
#' microclimate <- cbind(metout[, 1:5], metout[, 8], soil[, 4], metout[, 13:15], metout[, 6])
#' colnames(microclimate) <- c('dates', 'DOY', 'TIME', 'TALOC', 'TA1.2m', 'VLOC', 'TS', 'ZEN', 'SOLR', 'TSKYC', 'RHLOC')
#'
#' # define animal parameters - here simulating a 1000 g ellipsoid
#' c_body <- 3342 # specific heat of flesh, J/kg-C
#' rho_body <- 1000 # animal density, kg/m3
#' q <- 0 # metabolic rate, W/m3
#' k_flesh <- 0.5 # thermal conductivity of flesh, W/mK
#' geom <- 2 # shape, -
#' posture <- 'n' # pointing normal 'n' or parallel 'p' to the sun's rays, or average 'a'?
#' orient <- 1 # does the object orient toward the sun? (0,1)
#' shape_b <- 1/5 # shape coefficient a, -
#' shape_c <- 1/5 # shape coefficient b, -
#' shape_coefs <- c(10.4713, 0.688, 0.425, 0.85, 3.798, 0.683, 0.694, 0.743)
#' fatosk <- 0.4 # solar configuration factor to sky, -
#' fatosb <- 0.4 # solar configuration factor to substrate, -
#' alpha <- 0.9 # animal solar absorptivity, -
#' emis <- 0.95 # emissivity of skin, -
#' Ww_g <- 1000 # weight, g
#' alpha_sub <- 0.8 # substrate solar absorptivity, -
#' press <- 101325 # air pressure, Pa
#' pdif <- 0.1 # proportion of solar energy that is diffuse, -
#' dyn_q <- 0 # not dynamically varying q, -
#' fluid <- 0 # air, -
#'
#' # loop through middle day of each month
#' DOYs = c(15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349)
#' mons = c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December")
#'
#' for (i in 1:length(DOYs)) {
#'   simday = DOYs[i]
#'   microclim <- subset(microclimate, microclimate$DOY == simday)
#'
#'   # use approxfun to create interpolations for the required environmental variables
#'   time <- seq(0, 60 * 24, 60) #60 minute intervals from microclimate output
#'   time <- time * 60 # minutes to seconds
#'   hours <- time/3600 # seconds to hours
#'   Qsolf <- approxfun(time, c(microclim[, 9], (microclim[1, 9] + microclim[24, 9]) /
#'       2), rule = 2)
#'   # approximate radiant temperature as the average of sky and substrate temperature
#'   Tradf <- approxfun(time, rowMeans(cbind(c(microclim[, 7], (microclim[1, 7] + microclim[24, 7]) / 24), c(microclim[, 10], (microclim[1, 10] + microclim[24, 10]) / 24)), na.rm = TRUE), rule = 2)
#'   velf <- approxfun(time, c(microclim[, 6], (microclim[1, 6] + microclim[24, 6]) / 2), rule = 2)
#'   Tairf <- approxfun(time, c(microclim[, 4], (microclim[1, 4] + microclim[24, 4]) / 2), rule = 2)
#'   Zenf <- approxfun(time, c(microclim[, 8], 90), rule = 2)
#'
#'   t = seq(1, 3600 * 24, 60) # sequence of times for predictions (1 min intervals)
#'   indata<-list(alpha = alpha, emis = emis, alpha_sub = alpha_sub, press = press, Ww_g = Ww_g, c_body = c_body, rho_body = rho_body, q = q, k_flesh = k_flesh, geom = geom, posture = posture, shape_b = shape_b, shape_c = shape_c, shape_coefs = shape_coefs, pdif = pdif, fatosk = fatosk, fatosb = fatosb, fluid = fluid, dyn_q = dyn_q)
#'
#'   Tc_init<-Tairf(1) # set inital Tc as air temperature
#'
#'   Tbs_ode <- as.data.frame(ode(y = Tc_init, times = time, func = onelump_var, parms = indata))
#'   colnames(Tbs_ode) <- c('time', 'Tc', 'Tcf', 'tau', 'dTdt')
#'   Tbs_ode$time <- Tbs_ode$time / 3600 # convert to hours
#'
#'   with(Tbs_ode, plot(Tc ~ time, type = 'l', col = '1', ylim = c(-10, 80), xlim = c(0, 23), ylab='Temperature, °C',xlab = 'hour of day', main = paste0(round(loc[1], 1), round(loc[2], 1), ", ", mons[i], ", ", Ww_g," g")))
#'   with(Tbs_ode, points(Tcf ~ time, type = 'l', col = '2'))
#'   points(Tairf(time) ~ hours, type = 'l', col = 'blue', lty = 2)
#'   legend(0,70, c("Tc", "Tcf", "Tair"), lty = c(1, 1, 2), lwd = c(2.5, 2.5, 2.5), col = c("black", "red", "blue"))
#' }
#' @export
onelump_var <- function(t, y, indata) {
  with(as.list(c(indata, y)), {
    sigma <- 5.67e-8 # Stefan-Boltzman, W/(m.K)
    Tair <- Tairf(t) # interpolate air temp to current time step
    vel <- max(0,velf(t)) # interpolate wind speed to current time step
    Qsol <- max(0,Qsolf(t)) # interpolate solar to current time step
    Trad <- Tradf(t) # interpolate radiant temp to current time step
    Zen <- min(max(0,Zenf(t)),90) # interpolate zenith angle to current time step
    Zenith <- Zen * pi / 180 # zenith angle in radians
    Tc <- y
    Tskin <- Tc + 0.01
    vel[vel < 0.01] <- 0.01 # don't let wind speed go too low - always some free convection
    S2 <- 0.0001 # shape factor, arbitrary initialization, because not defined except for ellipsoid
    if(fluid == 0){
      DENSTY <- press / (287.04 * (Tair + 273.15)) # air density, kg/m3
      THCOND <- 0.02425 + (7.038 * 10 ^ -5 * Tair) # air thermal conductivity, W/(m.K)
      VISDYN <- (1.8325 * 10 ^ -5 * ((296.16 + 120) / ((Tair + 273.15) + 120))) * (((Tair + 273.15) / 296.16) ^ 1.5) # dynamic viscosity of air, kg/(m.s)
      DIFVPR <- 2.26e-05 * (((Tair + 273.15) / 273.15) ^ 1.81) * (1e+05 / press)
    }else{
      WATERPROP.out <- WATERPROP(Tair)
      DENSTY <- WATERPROP.out$DENSTY
      THCOND <- WATERPROP.out$THCOND
      VISDYN <- WATERPROP.out$VISDYN
      DIFVPR <- 0
    }
    VISKIN <- VISDYN / DENSTY
    # geometry section ############################################################
    m <- Ww_g / 1000 # convert weight to kg
    C <- m * c_body # thermal capacitance, J/K
    V <- m / rho_body # volume, m3
    if(dyn_q == 1){
      Q_gen <- qf(Tc) * V # total metabolic heat, J
    }else{
      Q_gen <- q * V
    }
    L <- V ^ (1 / 3) # characteristic dimension, m

    # FLAT PLATE geometry
    if (geom == 0) {
      AHEIT <- (V / (shape_b * shape_c)) ^ (1 / 3) # length, m
      AWIDTH <- AHEIT * shape_b # width, m
      ALENTH <- AHEIT * shape_c # height, m
      ATOT <- ALENTH * AWIDTH * 2 + ALENTH * AHEIT * 2 + AWIDTH * AHEIT * 2 # total area, m2
      ASILN <- ALENTH * AWIDTH # max silhouette area, m2
      ASILP <- AWIDTH * AHEIT # min silhouette area, m2
      L <- AHEIT # characteristic dimension, m
      if (AWIDTH <= ALENTH) {
        L <- AWIDTH
      } else{
        L <- ALENTH
      }
      R <- ALENTH / 2 # 'radius', m
    }

    # CYLINDER geometry
    if (geom == 1) {
      R1 <- (V / (pi * shape_b * 2)) ^ (1 / 3) # radius, m
      ALENTH <- 2 * R1 * shape_b # length, m
      ATOT <- 2 * pi * R1 ^ 2 + 2 * pi * R1 * ALENTH # total surface area, m2
      AWIDTH <- 2 * R1 # width, m
      ASILN <- AWIDTH * ALENTH # max silhouette area, m2
      ASILP <- pi * R1 ^ 2 # min silhouette area, m2
      L <- ALENTH # characteristic dimension, m
      R2 <- L / 2
      if (R1 > R2) {
        # choose shortest dimension as R
        R <- R2
      } else{
        R <- R1
      }
    }

    # Ellipsoid geometry
    if (geom == 2) {
      A1 <- ((3 / 4) * V / (pi * shape_b * shape_c)) ^ 0.333 # axis A, m
      B1 <- A1 * shape_b # axis B, m
      C1 <- A1 * shape_c # axis C, m
      P1 <- 1.6075 # a constant
      ATOT <- (4 * pi * (((A1 ^ P1 * B1 ^ P1 + A1 ^ P1 * C1 ^ P1 + B1 ^ P1 * C1 ^ P1)) / 3) ^ (1 / P1)) # total surface area, m2
      ASILN <- max(pi * A1 * C1, pi * B1 * C1) # max silhouette area, m2
      ASILP <- min(pi * A1 * C1, pi * B1 * C1) # min silhouette area, m2
      S2 <- (A1 ^ 2 * B1 ^ 2 * C1 ^ 2) / (A1 ^ 2 * B1 ^ 2 + A1 ^ 2 * C1 ^ 2 + B1 ^ 2 * C1 ^ 2) # fraction of semi-major and minor axes, see Porter and Kearney 2009 supp1
      #k_flesh <- 0.5 + 6.14 * B1 + 0.439 # thermal conductivity of flesh as a function of radius, see Porter and Kearney 2009
    }

    # Lizard geometry - DESERT IGUANA (PORTER ET AL. 1973 OECOLOGIA)
    if (geom == 3) {
      ATOT <- (10.4713 * Ww_g ^ .688) / 10000. # total surface area, m2
      AV <- (0.425 * Ww_g ^ .85) / 10000. # ventral surface area, m2
      # NORMAL AND POINTING @ SUN SILHOUETTE AREA: PORTER & TRACY 1984
      ASILN <- (3.798 * Ww_g ^ .683) / 10000. # Max. silhouette area (normal to the sun), m2
      ASILP <- (0.694 * Ww_g ^ .743) / 10000. # Min. silhouette area (pointing toward the sun), m2
      R <- L
    }

    # Frog geometry - LEOPARD FROG (C.R. TRACY 1976 ECOL. MONOG.)
    if (geom == 4) {
      ATOT <- (12.79 * Ww_g ^ 0.606) / 10000. # total surface area, m2
      AV <- (0.425 * Ww_g ^ 0.85) / 10000. # ventral surface area, m2
      # NORMAL AND POINTING @ SUN SILHOUETTE AREA: EQ'N 11 TRACY 1976
      ZEN <- 0
      PCTN <- 1.38171E-06 * ZEN ^ 4 - 1.93335E-04 * ZEN ^ 3 + 4.75761E-03 * ZEN ^ 2 - 0.167912 * ZEN + 45.8228
      ASILN <- PCTN * ATOT / 100 # Max. silhouette area (normal to the sun), m2
      ZEN <- 90
      PCTP <- 1.38171E-06 * ZEN ^ 4 - 1.93335E-04 * ZEN ^ 3 + 4.75761E-03 * ZEN ^ 2 - 0.167912 * ZEN + 45.8228
      ASILP <- PCTP * ATOT / 100 # Min. silhouette area (pointing toward the sun), m2
      R <- L
    }

    # user defined geometry
    if (geom == 5) {
      ATOT <- (shape_coefs[1] * Ww_g ^ shape_coefs[2]) / 10000. # total surface area, m2
      AV <- (shape_coefs[3] * Ww_g ^ shape_coefs[4]) / 10000 # ventral surface area, m2
      # NORMAL AND POINTING @ SUN SILHOUETTE AREA: PORTER & TRACY 1984
      # User must define Max. silhouette area (normal to the sun)
      ASILN <- (shape_coefs[5] * Ww_g ^ shape_coefs[6]) / 10000 # Max. silhouette area (normal to the sun), m2
      # User must define Min. silhouette area (pointing toward the sun)
      ASILP <- (shape_coefs[7] * Ww_g ^ shape_coefs[8]) / 10000 # Min. silhouette area (pointing toward the sun), m2
      R <- L
    }
    # end geometry section ############################################################

    if (max(Zen) >= 90) {
      Q_norm <- 0
    }else{
      if(orient == 1){
        Q_norm <- (Qsol / cos(Zenith))
      }else{
        Q_norm <- Qsol
      }
    }
    if (Q_norm > 1367) {
      Q_norm <- 1367 #making sure that low sun angles don't lead to solar values greater than the solar constant
    }
    if (posture == 'p') {
      Qabs <- (Q_norm * (1 - pdif) * ASILP + Qsol * pdif * fatosk * ATOT + Qsol * (1 - alpha_sub) * fatosb * ATOT) * alpha
    }
    if (posture == 'n') {
      Qabs <- (Q_norm * (1 - pdif) * ASILN + Qsol * pdif * fatosk * ATOT + Qsol * (1 - alpha_sub) * fatosb * ATOT) * alpha
    }
    if (posture == 'a') {
      Qabs <- (Q_norm * (1 - pdif) * (ASILN + ASILP) / 2 + Qsol * pdif * fatosk * ATOT + Qsol * (1 - alpha_sub) * fatosb * ATOT) * alpha
    }

    Re <- DENSTY * vel * L / VISDYN # Reynolds number
    PR <- 1005.7 * VISDYN / THCOND # Prandlt number
    if(fluid == 0){
      SC <- VISDYN / (DENSTY * DIFVPR) # Schmidt number
    }else{
      SC <- 1
    }

    if (geom == 0) {
      NUfor <- 0.102 * Re ^ 0.675 * PR ^ (1. / 3.)
    }
    if (geom == 3 | geom == 5) {
      NUfor <- 0.35 * Re ^ 0.6
    }
    if (geom == 1) {
      #       FORCED CONVECTION OF A CYLINDER
      #       ADJUSTING NU - RE CORRELATION FOR RE NUMBER (P. 260 MCADAMS,1954)
      if (Re < 4) {
        NUfor = 0.891 * Re ^ 0.33
      } else{
        if (Re < 40) {
          NUfor = 0.821 * Re ^ 0.385
        } else{
          if (Re < 4000) {
            NUfor = 0.615 * Re ^ 0.466
          } else{
            if (Re < 40000) {
              NUfor = 0.174 * Re ^ 0.618
            } else{
              if (Re < 400000) {
                NUfor = 0.0239 * Re ^ 0.805
              } else{
                NUfor = 0.0239 * Re ^ 0.805
              }
            }
          }
        }
      }
    }
    if (geom == 2 | geom == 4) {
      NUfor <- 0.35 * Re ^ (0.6) # Nusselt number, forced convection
    }
    h_conv_forced <- NUfor * THCOND / L # convection coefficent, forced

    GR <- abs(DENSTY ^ 2 * (1 / (Tair + 273.15)) * 9.80665 * L ^ 3 * (Tskin - Tair) / VISDYN ^ 2) # Grashof number
    Raylei <- GR * PR # Rayleigh number

    # get Nusselt for Free Convect
    if (geom == 0) {
      NUfre = 0.55 * Raylei ^ 0.25
    }
    if (geom == 1 | geom == 3 | geom == 5) {
      if (Raylei < 1.0e-05) {
        NUfre = 0.4
      } else{
        if (Raylei < 0.1) {
          NUfre = 0.976 * Raylei ^ 0.0784
        } else{
          if (Raylei < 100) {
            NUfre = 1.1173 * Raylei ^ 0.1344
          } else{
            if (Raylei < 10000.) {
              NUfre = 0.7455 * Raylei ^ 0.2167
            } else{
              if (Raylei < 1.0e+09) {
                NUfre = 0.5168 * Raylei ^ 0.2501
              } else{
                if (Raylei < 1.0e+12) {
                  NUfre = 0.5168 * Raylei ^ 0.2501
                } else{
                  NUfre = 0.5168 * Raylei ^ 0.2501
                }
              }
            }
          }
        }
      }
    }

    if (geom == 2 | geom == 4) {
      Raylei = (GR ^ 0.25) * (PR ^ 0.333)
      NUfre = 2 + 0.60 * Raylei
    }
    Nutotal <- (NUfre^3 + NUfor^3)^(1./3.)
    sh <- Nutotal * (SC / PR) ^ (1 / 3)
    h_conv <- Nutotal*(THCOND/L) # combined convection coefficient
    #R_conv <- 1 / (h_conv * ATOT) # convective resistance, eq. 3 of Kearney, Huey and Porter in prep. Appendix 1
    h_rad <- 4 * emis * sigma * ((Tc + Trad) / 2 + 273.15) ^ 3 # radiation heat transfer coefficient, eq. 46 of Transient Equations Derivation vignette

    if(geom == 2){ # ellipsoid
      j <- (Qabs + Q_gen + h_conv * ATOT * ((q * S2) / (2 * k_flesh) + Tair) + h_rad * ATOT * ((q * S2) / (2 * k_flesh) + Trad)) / C #based on eq. 52 of Kearney, Huey and Porter in prep. Appendix 1
    }else{ # assume cylinder
      j <- (Qabs + Q_gen + h_conv * ATOT * ((q * R ^ 2) / (4 * k_flesh) + Tair) + h_rad * ATOT * ((q * R ^ 2) / (2 * k_flesh) + Trad)) / C #based on eq. 52 of Kearney, Huey and Porter in prep. Appendix 1
    }
    theta_Tc <- ATOT * (Tc * h_conv + Tc * h_rad) / C #based on eq. 48 of Transient Equations Derivation vignette
    theta <- ATOT * (h_conv + h_rad) / C #based on eq. 48 of Transient Equations Derivation vignette
    Tcf <- j / theta # final Tc = j/theta, based on eq. 23 of Transient Equations Derivation vignette
    Tci <- Tc # initial temperature
    Tc <- (Tci - Tcf) * exp(-1 * theta * t) + Tcf # Tc at time t, based on eq. 33 of Transient Equations Derivation vignette
    tau <- 1 / theta # time constant
    dTc <- j - theta_Tc # rate of temperature change (°C/sec)
    list(y = dTc, x = Tcf, z = tau, zz = dTc)
  })
}