onelump_var: One-lump Transient Heat Budget (for use with deSolve package)

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

One-lump Transient Heat Budget (for use with deSolve package)

Description

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.

Usage

ode(y = Tb_init, t = times, func = onelump_var, parms = indata)

Arguments

t	time intervals (s) at which output is required
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)
Ww_g	= 500 weight (g)
rho_body	= 1000, animal density (kg/m3)
q	= 0, metabolic rate (W/m3)
c_body	= 3073, specific heat of flesh (J/kg-C)
k_flesh	= 0.5, conductivity of flesh (W/mK)
emis	= 0.95, emissivity of skin (-)
alpha	= 0.85, animal solar absorptivity (-)
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')
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
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
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
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)
fatosk	= 0.4, solar configuration factor to sky (-)
fatosb	= 0.4, solar configuration factor to substrate (-)
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)
alpha_sub	= 0.2, substrate solar reflectivity, decimal percent
pdif	= 0.1, proportion of solar energy that is diffuse (rather than direct beam)
fluid	= 0, fluid type, air (0) or water (1)
Tairf	air temperature function with time, generated by 'approxfun' (°C)
Tradf	radiant temperature function with time, generated by 'approxfun'(°C), averaging ground and sky
velf	wind speed function with time, generated by 'approxfun' (m/s)
Qsolf	radiation function with time, generated by 'approxfun' (W/m2)
Zenf	zenith angle of sun function with time, generated by 'approxfun' (90 is below horizon), degrees
press	air pressure (Pa)

Value

Tc Core temperature (°C)

Tcf Final (steady state) temperature (°C), if conditions remained constant indefinitely

tau Time constant (s)

dTc Rate of change of core temperature (°C/s)

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"))
}

mrke/NicheMapR documentation built on April 3, 2024, 10:05 a.m.