This document shows how to use the endoR.R
function of NicheMapR to compute endotherm heat and mass balances. This function collates all the different steps in Endotherm Components Tutorial into a single function, ‘endoR’, that can be applied to vectors of data. The latter calls a FORTRAN routine call ‘SOLVENDO’ (high speed calculation). An R equivalent of this, ‘endoR_devel’ (MUCH slower), can be used for the development of algorithms for organism-specific sequences of thermoregulatory responses, and then a version of SOLVENDO can be modified accordingly.
Here is how you run the calculation with all default values, specifying a local air temperature of 0 deg C.
library(NicheMapR)
library(knitr)
endo.out <- endoR(TA = 0)
treg <- as.data.frame(endo.out$treg)
morph <- as.data.frame(endo.out$morph)
enbal <- as.data.frame(endo.out$enbal)
masbal <- as.data.frame(endo.out$masbal)
value | description and units | |
---|---|---|
TC | 37 | core temperature (°C) |
TLUNG | 34.0428 | lung temperature (°C) |
TSKIN_D | 31.0855 | dorsal skin temperature (°C) |
TSKIN_V | 31.0855 | ventral skin temperature (°C) |
TFA_D | 15.7187 | dorsal fur/feather-air interface temperature (°C) |
TFA_V | 15.7187 | ventral fur-air interface temperature (°C) |
SHAPE_B | 3 | current ratio between long and short axis due to postural change (-) |
PANT | 1 | breathing rate multiplier (-) |
SKINWET | 0.5 | part of the skin surface that is wet (%) |
K_FLESH | 1.4 | thermal conductivity of flesh (W/mC) |
K_FUR | 0.0332 | thermal conductivity of fur/feathers (W/mC) |
K_FUR_D | 0.0332 | thermal conductivity of dorsal fur (W/mC) |
K_FUR_V | 0.0332 | thermal conductivity of ventral fur/feathers (W/mC) |
K_COMPFUR | 0.0332 | thermal conductivity of compressed fur/feathers (W/mC) |
Q10 | 1 | Q10 multiplier on metabolic rate (-) |
value | description and units | |
---|---|---|
AREA | 0.0611 | total outer surface area (m2) |
VOLUME | 0.001 | total volume (m3) |
CHAR_DIM | 0.09 | characteristic dimension for convection (m) |
MASS_FAT | 0.2 | fat mass (kg) |
FAT_THICK | 0 | thickness of fat layer (m) |
FLESH_VOL | 0.001 | flesh volume (m3) |
LENGTH | 0.2581 | length (m) |
WIDTH | 0.086 | width (m) |
HEIGHT | 0.086 | height (m) |
DIAM_FLESH | 0.043 | diameter, core to skin (m) |
DIAM_FUR | 0.045 | diameter, core to fur/feathers (m) |
AREA_SIL | 0.0124 | silhouette area (m2) |
AREA_SILN | 0.0185 | silhouette area normal to sun’s rays (m2) |
AREA_ASILP | 0.0064 | silhouette area parallel to sun’s rays (m2) |
AREA_SKIN | 0.0572 | total skin area (m2) |
AREA_SKIN_EVAP | 0.0559 | skin area available for evaporation (m2) |
AREA_CONV | 0.0611 | area for convection (m2) |
AREA_COND | 0 | area for conduction (m2) |
F_SKY | 0.5 | configuration factor to sky (-) |
F_GROUND | 0.5 | configuration factor to ground (-) |
value | description and units | |
---|---|---|
QSOL | 0 | solar radiation absorbed (W) |
QIRIN | 22.286 | longwave (infra-red) radiation absorbed (W) |
QMET | 12.3092 | characteristic dimension for convection (W) |
QEVAP | 0.2872 | evaporation (W) |
QIROUT | 27.4852 | longwave (infra-red) radiation lost (W) |
QCONV | 6.829 | convection (W) |
QCOND | 0 | conduction (W) |
ENB | -0.0062 | energy balance (W) |
NTRY | 1 | iterations required for a solution (-) |
SUCCESS | 1 | was a solution found (0=no, 1=yes) |
value | description and units | |
---|---|---|
AIR_L | 56.158 | breating rate (L/h) |
O2_L | 2.353 | oxgyen consumption rate (L/h) |
H2OResp_g | 0.249 | respiratory water loss (g/h) |
H2OCut_g | 0.1724 | cutaneous water loss (g/h) |
O2_mol_in | 0.5252 | oxygen inhaled (mol/h) |
O2_mol_out | 0.4202 | oxygen expelled (mol/h) |
N2_mol_in | 1.9811 | nitrogen inhaled (mol/h) |
N2_mol_out | 1.9811 | nitrogen expelled (mol/h) |
AIR_mol_in | 2.5071 | air inhaled (mol/h) |
AIR_mol_out | 2.406 | air expelled (mol/h) |
Now try a sequence of air temperatures for a specified core temperature, mass, fur/feather depth and base skin wetness, with the option to sweat at 0.25 % increments to cool down if needed. Note that the default settings make the environment a blackbody one (all surface and air temperatures the same) with low (5%) relative humidity and low wind speed (0.1 m/s), as is often the case in metabolic chamber experiments.
library(NicheMapR)
# environment
TAs <- seq(0, 50, 1) # air temperature (deg C)
VEL <- 0.002 # wind speed (m/s)
vd <- WETAIR(rh = 30, db = 40)$vd # Weather and Schoenbaechler had 16.7 mm Hg above 40 deg C = 30% RH at 40 deg C
vd_sat <- WETAIR(rh = 100, db = TAs)$vd # Weather and Schoenbaechler had 16.7 mm Hg above 40 deg C = 30% RH at 40 deg C
exp_rh <- vd/vd_sat * 100
exp_rh[exp_rh > 100] <- 100
exp_rh[TAs < 30] <- 15
hum <- exp_rh #rep(humidity,96)
# core temperature
TC <- 38 # core temperature (deg C)
TCMAX <- 43 # maximum core temperature (deg C)
RAISETC <- 0.25 # increment by which TC is elevated (deg C)
# size and shape
AMASS <- 0.0337 # mass (kg)
SHAPE_B_REF <- 1.1 # start off near to a sphere (-)
SHAPE_B_MAX <- 5 # maximum ratio of length to width/depth
UNCURL <- 0.1 # allows the animal to uncurl to SHAPE_B_MAX, the value being the increment SHAPE_B is increased per iteration
SHAPE <- 4 # use ellipsoid geometry
SAMODE <- 1 # use bird surface area relations (2 is mammal, 0 is based on shape specified in GEOM)
# feather properties
DHAIRD <- 3e-05 # hair diameter, dorsal (m)
DHAIRV <- 3e-05 # hair diameter, ventral (m)
LHAIRD <- 0.0231 # hair length, dorsal (m)
LHAIRV <- 0.0227 # hair length, ventral (m)
ZFURD <- 0.0058 # fur depth, dorsal (m)
ZFURV <- 0.0056 # fur depth, ventral (m)
RHOD <- 8e+07 # hair density, dorsal (1/m2)
RHOV <- 8e+07 # hair density, ventral (1/m2)
REFLD <- 0.248 # fur reflectivity dorsal (fractional, 0-1)
REFLV <- 0.351 # fur reflectivity ventral (fractional, 0-1)
# physiological responses
SKINW <- 0.5 # base skin wetness (%)
MXWET <- 5 # maximum skin wetness (%)
SWEAT <- 0.25 # intervals by which skin wetness is increased (%)
Q10s <- rep(1, length(TAs))
Q10s[TAs >= TCMAX] <- 2 # assuming Q10 effect kicks in only after air temp rises above TCMAX
QBASAL <- 10^(-1.461 + 0.669 * log10(AMASS * 1000)) # basal heat generation (W)
DELTAR <- 5 # offset between air temeprature and breath (°C)
EXTREF <- 25 # O2 extraction efficiency (%)
PANTING <- 0.1 # turns on panting, the value being the increment by which the panting multiplier is increased up to the maximum value, PANTMAX
PANTMAX <- 15 # maximum panting rate - multiplier on air flow through the lungs above that determined by metabolic rate
ptm <- proc.time() # start timing
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = VEL, TC = TC, TCMAX = TCMAX, RH = hum[x],
AMASS = AMASS, SHAPE = SHAPE, SHAPE_B_REF = SHAPE_B_REF,
SHAPE_B_MAX = SHAPE_B_MAX, SKINW = SKINW, SWEAT = SWEAT,
MXWET = MXWET, Q10 = Q10s[x], QBASAL = QBASAL, DELTAR = DELTAR,
DHAIRD = DHAIRD, DHAIRV = DHAIRV, LHAIRD = LHAIRD, LHAIRV = LHAIRV,
ZFURD = ZFURD, ZFURV = ZFURV, RHOD = RHOD, RHOV = RHOV,
REFLD = REFLD, RAISETC = RAISETC, PANTING = PANTING,
PANTMAX = PANTMAX, EXTREF = EXTREF, UNCURL = UNCURL,
SAMODE = SAMODE)
}) # run endoR across environments
proc.time() - ptm # stop timing
## user system elapsed
## 0.02 0.02 0.37
# extract the output
endo.out1 <- do.call("rbind", lapply(endo.out, data.frame))
# thermoregulation output
treg <- endo.out1[, grep(pattern = "treg", colnames(endo.out1))]
colnames(treg) <- gsub(colnames(treg), pattern = "treg.", replacement = "")
# morphometric output
morph <- endo.out1[, grep(pattern = "morph", colnames(endo.out1))]
colnames(morph) <- gsub(colnames(morph), pattern = "morph.",
replacement = "")
# heat balance
enbal <- endo.out1[, grep(pattern = "enbal", colnames(endo.out1))]
colnames(enbal) <- gsub(colnames(enbal), pattern = "enbal.",
replacement = "")
# mass aspects
masbal <- endo.out1[, grep(pattern = "masbal", colnames(endo.out1))]
colnames(masbal) <- gsub(colnames(masbal), pattern = "masbal.",
replacement = "")
QGEN <- enbal$QMET # metabolic rate (W)
H2O <- masbal$H2OResp_g + masbal$H2OCut_g # g/h water evaporated
TFA_D <- treg$TFA_D # dorsal fur surface temperature
TFA_V <- treg$TFA_V # ventral fur surface temperature
TskinD <- treg$TSKIN_D # dorsal skin temperature
TskinV <- treg$TSKIN_V # ventral skin temperature
TCs <- treg$TC # core temperature
par(mfrow = c(2, 2))
par(oma = c(2, 1, 2, 2) + 0.1)
par(mar = c(3, 3, 1.5, 1) + 0.1)
par(mgp = c(2, 1, 0))
plot(QGEN ~ TAs, type = "l", ylab = "metabolic rate, W", xlab = "air temperature, deg C",
ylim = c(0.2, 1.2))
points(Weathers1976Fig1$Tair, Weathers1976Fig1$mlO2gh * 20.1/3600 *
(AMASS * 1000), pch = 16, col = 2)
legend(x = 30, y = 1.2, legend = c("observed", "predicted"),
col = c("red", "black"), lty = c(NA, 1), bty = "n", pch = c(16,
NA))
plot(H2O ~ TAs, type = "l", ylab = "water loss, g/h", xlab = "air temperature, deg C",
ylim = c(0, 1.5))
points(masbal$H2OResp_g ~ TAs, type = "l", lty = 2)
points(masbal$H2OCut_g ~ TAs, type = "l", lty = 2, col = "blue")
legend(x = 3, y = 1.5, legend = c("total (obs)", "total (pred)",
"respiratory", "cutaneous"), col = c("red", "black", "black",
"blue"), lty = c(NA, 1, 2, 2), bty = "n", pch = c(16, NA,
NA, NA))
points(Weathers1976Fig3$Tair, Weathers1976Fig3$mgH2Ogh * AMASS,
pch = 16, col = 2)
plot(TFA_D ~ TAs, type = "l", col = "grey", ylab = "feather, skin and core temperature, deg C",
xlab = "air temperature, deg C", ylim = c(10, 50))
points(TFA_V ~ TAs, type = "l", col = "grey", lty = 2)
points(TskinD ~ TAs, type = "l", col = "orange")
points(TskinV ~ TAs, type = "l", col = "orange", lty = 2)
points(TCs ~ TAs, type = "l", col = "red")
legend(x = 30, y = 33, legend = c("core (obs)", "core (pred)",
"skin dorsal", "skin ventral", "feathers dorsal", "feathers ventral"),
col = c("red", "red", "orange", "orange", "grey", "grey"),
lty = c(NA, 1, 1, 2, 1, 2), bty = "n", pch = c(16, NA, NA,
NA, NA, NA))
points(Weathers1976Fig2$Tair, Weathers1976Fig2$Tb, pch = 16,
col = 2)
plot(masbal$AIR_L * 1000/60 ~ TAs, ylim = c(0, 250), lty = 1,
xlim = c(-5, 50), ylab = "ml / min", xlab = paste("air temperature (deg C)"),
type = "l")
legend(x = 0, y = 250, legend = c("observed", "predicted"), col = c("red",
"black"), lty = c(NA, 1), bty = "n", pch = c(16, NA))
points(Weathers1976Fig5$breaths_min * (13.2 * AMASS^1.08) * ((Weathers1976Fig5$Tair +
273.15)/273.15) ~ Weathers1976Fig5$Tair, col = "red", pch = 16) # tidal volume allometry from Lasiewski, R. C., and W. A. Calder. 1971, correcting volume according to PV = nRT equation, where V_2 = T_2 * V_1 / T_2, and T_1 is at STP, so 0 deg C
In this final example, the microclimate model is run for a location in central Australia, under unshaded conditions, and the then the relevant outputs are sent into the endotherm model along with the same animal parameters as above.
library(NicheMapR)
loc <- c(131.05, -22.75)
Usrhyt <- 0.07
maxshade <- 100
micro <- micro_global(loc = loc, Usrhyt = Usrhyt, maxshade = maxshade)
metout <- as.data.frame(micro$metout) # unshaded above-ground conditions
soil <- as.data.frame(micro$soil) # unshaded below-ground soil temperatures
shadmet <- as.data.frame(micro$shadmet) # shaded above-ground conditions
shadsoil <- as.data.frame(micro$shadsoil) # shaded below-ground soil temperatures
dates <- micro$dates
TAs <- metout$TALOC # air temperatures at height of animal (deg C)
TAREFs <- metout$TAREF # air temperatures at reference height (deg C)
TSKYs <- metout$TSKYC # sky temperatures (deg C)
TGRDs <- soil$D0cm # surface temperatures (deg C)
VELs <- metout$VLOC # wind speeds at animal height (m/s)
RHs <- metout$RHLOC # relative humidity at animal height (%)
QSOLRs <- metout$SOLR # solar radiation (W/m2)
Zs <- metout$ZEN # zenith angle of the sun (degrees)
ELEV <- micro$elev # elevation (m)
ABSSB <- 1 - micro$REFL # substrate solar absorptivity (%)
# core temperature
TC <- 38 # core temperature (deg C)
TCMAX <- 43 # maximum core temperature (deg C)
RAISETC <- 0.25 # increment by which TC is elevated (deg C)
# size and shape
AMASS <- 0.0337 # mass (kg)
SHAPE <- 4 # use ellipsoid geometry
SHAPE_B_REF <- 1.1 # start off near to a sphere (-)
SHAPE_B_MAX <- 5 # maximum ratio of length to width/depth
UNCURL <- 0.1 # allows the animal to uncurl to SHAPE_B_MAX, the value being the increment SHAPE_B is increased per iteration
# feather properties
DHAIRD = 3e-05 # hair diameter, dorsal (m)
DHAIRV = 3e-05 # hair diameter, ventral (m)
LHAIRD = 0.0231 # hair length, dorsal (m)
LHAIRV = 0.0227 # hair length, ventral (m)
ZFURD = 0.0058 # fur depth, dorsal (m)
ZFURV = 0.0056 # fur depth, ventral (m)
RHOD = 8e+07 # hair density, dorsal (1/m2)
RHOV = 8e+07 # hair density, ventral (1/m2)
REFLD = 0.248 # fur reflectivity dorsal (fractional, 0-1)
REFLV = 0.351 # fur reflectivity ventral (fractional, 0-1)
# physiological responses
SKINW <- 0.1 # base skin wetness (%)
MXWET <- 0.1 # maximum skin wetness (%)
SWEAT <- 0.25 # intervals by which skin wetness is increased (%)
Q10 <- 1 # Q10 effect of body temperature on metabolic rate (-)
QBASAL <- 10^(-1.461 + 0.669 * log10(AMASS * 1000)) # basal heat generation (W)
DELTAR <- 5 # offset between air temeprature and breath (°C)
EXTREF <- 25 # O2 extraction efficiency (%)
PANTING <- 0.1 # turns on panting, the value being the increment by which the panting multiplier is increased up to the maximum value, PANTMAX
PANTMAX <- 15 # maximum panting rate - multiplier on air flow through the lungs above that determined by metabolic rate
# run the model
ptm <- proc.time() # start timing
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], TAREF = TAREFs[x], TSKY = TSKYs[x], TGRD = TGRDs[x],
VEL = VELs[x], RH = RHs[x], QSOLR = QSOLRs[x], Z = Zs[x],
ELEV = ELEV, ABSSB = ABSSB, TC = TC, TCMAX = TCMAX, AMASS = AMASS,
SHAPE = SHAPE, SHAPE_B_REF = SHAPE_B_REF, SHAPE_B_MAX = SHAPE_B_MAX,
SKINW = SKINW, SWEAT = SWEAT, Q10 = Q10, QBASAL = QBASAL,
DELTAR = DELTAR, DHAIRD = DHAIRD, DHAIRV = DHAIRV, LHAIRD = LHAIRD,
LHAIRV = LHAIRV, ZFURD = ZFURD, ZFURV = ZFURV, RHOD = RHOD,
RHOV = RHOV, REFLD = REFLD, RAISETC = RAISETC, PANTING = PANTING,
PANTMAX = PANTMAX, EXTREF = EXTREF, UNCURL = UNCURL,
SAMODE = SAMODE, SHADE = 0)
})
proc.time() - ptm # end timing
## user system elapsed
## 0.38 0.25 2.08
# extract the output
endo.out1 <- do.call("rbind", lapply(endo.out, data.frame))
# thermoregulation output
treg <- endo.out1[, grep(pattern = "treg", colnames(endo.out1))]
colnames(treg) <- gsub(colnames(treg), pattern = "treg.", replacement = "")
# morphometric output
morph <- endo.out1[, grep(pattern = "morph", colnames(endo.out1))]
colnames(morph) <- gsub(colnames(morph), pattern = "morph.",
replacement = "")
# heat balance
enbal <- endo.out1[, grep(pattern = "enbal", colnames(endo.out1))]
colnames(enbal) <- gsub(colnames(enbal), pattern = "enbal.",
replacement = "")
# mass aspects
masbal <- endo.out1[, grep(pattern = "masbal", colnames(endo.out1))]
colnames(masbal) <- gsub(colnames(masbal), pattern = "masbal.",
replacement = "")
# extract variables for plotting
QGEN <- enbal$QMET # metabolic rate (W)
H2O <- masbal$H2OResp_g + masbal$H2OCut_g # g/h water evaporated
TFA_D <- treg$TFA_D # dorsal fur surface temperature
TFA_V <- treg$TFA_V # ventral fur surface temperature
TskinD <- treg$TSKIN_D # dorsal skin temperature
TskinV <- treg$TSKIN_V # ventral skin temperature
TCs <- treg$TC # core temperature
SkinW <- treg$SKINWET # skin wetness (%)
Pant <- treg$PANT # panting multiplier (-)
# plot the results
par(mfrow = c(2, 2))
par(oma = c(2, 1, 2, 2) + 0.1)
par(mar = c(3, 3, 1.5, 1) + 0.1)
par(mgp = c(2, 1, 0))
plot(QGEN ~ dates, type = "l", ylab = "metabolic rate, W", xlab = "time",
ylim = c(0, QBASAL * 5))
plot(H2O ~ dates, type = "l", ylab = "water loss, g / h", xlab = "time",
ylim = c(0, 2))
points(masbal$H2OResp_g ~ dates, type = "l", lty = 2)
points(masbal$H2OCut_g ~ dates, type = "l", lty = 2, col = "blue")
legend(x = 4, y = 2, legend = c("total", "respiratory", "cutaneous"),
col = c("black", "black", "blue"), lty = c(1, 2, 2), bty = "n")
plot(TFA_D ~ dates, type = "l", col = "grey", ylab = "temperature, deg C",
xlab = "time", ylim = c(0, 60))
points(TFA_V ~ dates, type = "l", col = "grey", lty = 2)
points(TskinD ~ dates, type = "l", col = "orange")
points(TskinV ~ dates, type = "l", col = "orange", lty = 2)
points(TCs ~ dates, type = "l", col = "red")
legend(x = 1, y = 65, legend = c("core", "skin dorsal", "skin ventral",
"feathers dorsal", "feathers ventral"), col = c("red", "orange",
"orange", "grey", "grey"), lty = c(1, 1, 2, 1, 2), bty = "n",
ncol = 2)
plot((H2O/(AMASS * 1000)) * 1000 ~ dates, type = "l", ylab = "H2O loss, % body mass / h",
xlab = "time", ylim = c(0, 50))