Endotherm Model Tutorial II: example simulations

Michael Kearney

2019-10-08

Overview

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.

Single calculation example

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)

treg
	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 (-)

morph
	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 (-)

enbal
	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)

masbal
	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)

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

# 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

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

# 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))

Endotherm Model Tutorial II: example simulations

Michael Kearney

2019-10-08

Overview

Single calculation example

Metabolic chamber example

Microclimate model integration example