twolump: Two-lump Transient Heat Budget (for use with deSolve package)

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

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

Description

Transient, 'two-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 = c(Tc_init, Tsk_init, Ts_init), times = time, func = twolump, parms = indata)

Arguments

t	time intervals (s) at which output is required
Tc_init	= 5, initial core temperature (°C)
Tsk_init	= 5.1, initial shell (skin) temperature (°C)
Ts_init	= 5.2, initial surface temperature (°C)
Ww_g	= 500, animal weight (g)
rho_body	= 932, animal density (kg/m3)
x_shell	= 0.001, shell thickness, m
q	= 0, metabolic heat production rate W/m3
c_body_inner	= 3073, Specific heat of flesh J/(kg-K)
c_body_outer	= 3073, Specific heat of outer shell J/(kg-K)
k_inner	= 0.5, Thermal conductivity of inner shell (W/mK, range: 0.412-2.8)
k_outer	= 0.5, Thermal conductivity of outer shell (W/mK, range: 0.412-2.8)
emis	= 0.95, Emissivity of animal (0-1)
alpha	= 0.85, solar absorptivity, decimal percent
geom	= 2, Organism shape, 1-2, Determines surface area/volume relationships: 1=cyl, 2=ellips
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
posture	= 'n', pointing normal 'n', parallel 'p' to the sun's rays, or 'a' in between?
orient	= 1, does the object orient toward the sun? (0,1)
fatosk	= 0.4, Configuration factor to sky for infrared calculations (-)
fatosb	= 0.4, Configuration factor to substrate for infrared calculations (-)
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	= 101325, air pressure (Pa)

Value

Tc Core temperature (°C)

Tsk 'Skin' / shell temperature (°C)

Ts Outer surface temperature (°C)

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

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
c_body_inner <- c_body
c_body_outer <- c_body
rho_body <- 1000 # animal density, kg/m3
q <- 0 # metabolic rate, W/m3
k_inner <- 0.5 # thermal conductivity of core, W/mK
k_outer <- 0.5 # thermal conductivity of shell, W/mK
geom <- 2 # shape, -
x_shell <- 0.005 # thickness of outer shell (m)
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, -
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(Ww_g = Ww_g, x_shell = x_shell, geom = geom, k_inner = k_inner, k_outer = k_outer, q = q, c_body_inner = c_body_inner, c_body_outer = c_body_outer, emis = emis, rho_body = rho_body, alpha = alpha, shape_b = shape_b, shape_c = shape_c, posture = posture, orient = orient, fatosk = fatosk, fatosb = fatosb, alpha_sub = alpha_sub, pdif = pdif, press = press, fluid = fluid, dyn_q = dyn_q)

  Tc_init<-Tairf(1) # set inital Tc as air temperature
  Tsk_init <- Tc_init
  Ts_init <- Tc_init

  Tbs_ode <- as.data.frame(ode(y = c(Tc_init, Tsk_init, Ts_init), times = time, func = twolump, parms = indata))
  colnames(Tbs_ode) = c("time", "Tc", "Tsk", "Ts", "Tcf")
  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(loc, ", ", mons[i], ", ", Ww_g," g")))
  with(Tbs_ode, points(Ts ~ time, lty = 2, type = 'l', col = '1'))
  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", "Ts", "Tcf", "Tair"), lty = c(1, 2, 1, 2), lwd = c(2.5, 2.5, 2.5, 2.5), col = c("black", "black", "red", "blue"))
}

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