Endotherm Model Tutorial I: component functions

Michael Kearney

2019-10-08

Overview

This document provides an overview of an R implementation of the NicheMapR endotherm model, via the sub-functions of ‘endoR’. This implementation modularises the original version such that many of the internal subroutines and functions can be called separately. The idea is that the user has more flexibility in setting the model up for the behavioural and physiological responses of a particular animal. There is also a tutorial for using ‘endoR’ itself Endotherm Model Tutorial.

library(NicheMapR)
library(knitr)
library(plotrix)

Running the separate subroutines

IRPROP

This subroutine computes the feasible hair spacing and the three parameters needed for computing conduction and infrared radiation through the fur. It assumes a thermal conductivity of hair (keratin) of 0.209 W m{-1} K{-1} and that the thermal conductivity of air is 0.02425+(7.038e-5*TA).

# environmental input
TA <- 20  # air temperature, for calculation of conductivity of air (°C)

# shape input
SHAPE_B_MAX <- 2.7  # max possible ratio between long and short axis (-)
SHAPE_B_REF <- 2.7  # initial ratio between long and short axis (-)
SHAPE_B <- 2.7  # current ratio between long and short axis (-)
MAXPTVEN <- 0.3  # maxium fraction of surface area that is ventral (fractional, 0-1)

# fur properties
DHAIRD <- 3e-05  # hair diameter, dorsal (m)
DHAIRV <- 3e-05  # hair diameter, ventral (m)
LHAIRD <- 0.0239  # hair length, dorsal (m)  
LHAIRV <- 0.0239  # hair length, ventral (m)  
ZFURD <- 0.009  # fur depth, dorsal (m)
ZFURV <- 0.009  # fur depth, ventral (m)
ZFURCOMP <- ZFURV  # depth of compressed fur (for conduction) (m)
RHOD <- 39680000  # hair density, dorsal (1/m2) 
RHOV <- 27810000  # hair density, ventral (1/m2)
REFLD <- 0.301  # fur reflectivity dorsal (fractional, 0-1) 
REFLV <- 0.301  # fur reflectivity ventral (fractional, 0-1)

# call the subroutine
IRPROP.out <- IRPROP(TA, SHAPE_B_MAX, SHAPE_B_REF, SHAPE_B, DHAIRD, 
    DHAIRV, LHAIRD, LHAIRV, ZFURD, ZFURV, RHOD, RHOV, REFLD, 
    REFLV, MAXPTVEN, ZFURCOMP)

# output
E4B1 <- IRPROP.out[1]  # exponential integral, E4, of optical thickness, B1 (m), i.e. E4(B1)
KEFARA <- IRPROP.out[2:4]  # effective thermal conductivity of fur array, mean, dorsal, ventral (W/mK)
BETARA <- IRPROP.out[5:7]  # term involved in computing optical thickess
B1ARA <- IRPROP.out[8:10]  # optical thickness array, mean, dorsal, ventral (-)
DHAR <- IRPROP.out[11:13]  # fur diameter array, mean, dorsal, ventral (m)
LHAR <- IRPROP.out[14:16]  # fur length array, mean, dorsal, ventral (m)
RHOAR <- IRPROP.out[17:19]  # fur density array, mean, dorsal, ventral (fibers/m2)  
ZZFUR <- IRPROP.out[20:22]  # fur depth array, mean, dorsal, ventral (m)  
REFLFR <- IRPROP.out[23:25]  # fur reflectivity array, mean, dorsal, ventral (fractional, 0-1)
FURTST <- IRPROP.out[26]  # test of presence of fur (length x diamater x density x depth) (-)
KFURCMPRS <- IRPROP.out[27]  # effictive thermal conductivity of compressed ventral fur (W/mK)

IRPROP.lab <- c("E4B1", "KEFARA mean", "KEFARA dorsal", "KEFARA ventral", 
    "BETARA mean", "BETARA dorsal", "BETARA ventral", "B1ARA mean", 
    "B1ARA dorsal", "B1ARA ventral", "DHAR mean", "DHAR dorsal", 
    "DHAR ventral", "LHAR mean", "LHAR dorsal", "LHAR ventral", 
    "RHOAR mean", "RHOAR dorsal", "RHOAR ventral", "ZZFUR mean", 
    "ZZFUR dorsal", "ZZFUR ventral", "REFLFR mean", "REFLFR dorsal", 
    "REFLFR ventral", "FURTST", "KFURCMPRS")
kable(cbind(IRPROP.lab, IRPROP.out[1:27]))
IRPROP.lab
E4B1 0.000802058267814305
KEFARA mean 0.029306880523346
KEFARA dorsal 0.0296911081471579
KEFARA ventral 0.0289267341979939
BETARA mean 613.67390356235
BETARA dorsal 674.176485875967
BETARA ventral 553.171321248733
B1ARA mean 5.52306513206115
B1ARA dorsal 6.0675883728837
B1ARA ventral 4.9785418912386
DHAR mean 3e-05
DHAR dorsal 3e-05
DHAR ventral 3e-05
LHAR mean 0.0239
LHAR dorsal 0.0239
LHAR ventral 0.0239
RHOAR mean 36119000
RHOAR dorsal 39680000
RHOAR ventral 32558000
ZZFUR mean 0.009
ZZFUR dorsal 0.009
ZZFUR ventral 0.009
REFLFR mean 0.301
REFLFR dorsal 0.301
REFLFR ventral 0.301
FURTST 0.25605504
KFURCMPRS 0.0289267341979939

Example analysis of fur diameter on effective thermal conductivity of fur.

Example analysis of fur density on effective thermal conductivity of fur.

Example analysis of fur depth on effective thermal conductivity of fur.

GEOM

This subroutine computes aspects of the animal’s geometry. Note that GEOM requires outputs from IRPROP on fur properties. Note also that this routine is called ‘allom’ in past versions of ‘Niche Mapper’.

# input
AMASS <- 10  # kg
ANDENS <- 1000  # kg/m3
SUBQFAT <- 1  # is subcutaneous fat present? (0 is no, 1 is yes)
FATPCT <- 20  # % body fat
SHAPE <- 4  # shape, 1 is cylinder, 2 is sphere, 3 is plate, 4 is ellipsoid
SHAPE_B_REF <- 2.7  # initial ratio between long and short axis (-)
SHAPE_B <- 2.7  # current ratio between long and short axis (-)
SHAPE_C <- SHAPE_B  # current ratio of length:height (plate)
DHARA <- DHAR[1]  # fur diameter, mean (m) (from IRPROP)
RHOARA <- RHOAR[1]  # hair density, mean (1/m2) (from IRPROP)
ZFUR <- ZZFUR[1]  # fur depth, mean (m) (from IRPROP)
PCOND <- 0.1  # fraction of body in contact with substrate (fractional, 0-1)
BIRD <- 0  # 
SAMODE <- 0  # if 1, uses bird skin surface area scaling from Walsberg and 1978. JEB Biology 76:185–189, if 2, uses mammal surface area scaling from Stahl (1967) J. of App. Physiology, 453–460.
ORIENT <- 0  # if 1, largest surface area normal to sun's ray's, if 2, largest surface parallel to sun's rays, if 0, average of normal/parallel posture, 

# call the subroutine
GEOM.out <- GEOM(AMASS, ANDENS, FATPCT, SHAPE, ZFUR, SUBQFAT, 
    SHAPE_B, SHAPE_B_REF, SHAPE_C, DHARA, RHOARA, PCOND, SAMODE, 
    ORIENT)

# output
VOL <- GEOM.out[1]  # volume, m3
D <- GEOM.out[2]  # characteristic dimension for convection, m
MASFAT <- GEOM.out[3]  # mass body fat, kg
VOLFAT <- GEOM.out[4]  # volume body fat, m3
ALENTH <- GEOM.out[5]  # length, m
AWIDTH <- GEOM.out[6]  # width, m
AHEIT <- GEOM.out[7]  # height, m
ATOT <- GEOM.out[8]  # total area at fur/feathers-air interface, m2
ASIL <- GEOM.out[9]  # silhouette area to use in solar calcs, m2 may be normal, parallel or average set via ORIENT
ASILN <- GEOM.out[10]  # silhouette area normal to sun, m2
ASILP <- GEOM.out[11]  # silhouette area parallel to sun, m2
GMASS <- GEOM.out[12]  # mass, g
AREASKIN <- GEOM.out[13]  # area of skin, m2
FLSHVL <- GEOM.out[14]  # flesh volume, m3
FATTHK <- GEOM.out[15]  # fat layer thickness, m
ASEMAJ <- GEOM.out[16]  # semimajor axis length, m
BSEMIN <- GEOM.out[17]  # b semiminor axis length, m
CSEMIN <- GEOM.out[18]  # c semiminor axis length, m (currently only prolate spheroid)
CONVSK <- GEOM.out[19]  # area of skin for evaporation (total skin area - hair area), m2
CONVAR <- GEOM.out[20]  # area for convection (total area minus ventral area, as determined by PCOND), m2
R1 <- GEOM.out[21]  # shape-specific core-skin radius in shortest dimension, m
R2 <- GEOM.out[22]  # shape-specific core-fur radius in shortest dimension, m

# nest properties
NESTYP <- 3  # for nest calculations - 0 = none, 1 = flat, 2 = cup, 3 = cylinder, 4 = half cylinder, 5 = sphere, 6 = dome
D_NEST <- 0.35  # Outer diameter(m)
TKNEST <- 0.01  # Nest wall thickness (m) 
RONEST <- D_NEST/2
RINEST <- RONEST - TKNEST
if (RINEST < R2) {
    # NEST INNER DIAMETER GREATER THAN ANIMAL WITH OR WITHOUT FUR
    # OUTER DIAMETER. ENLARGE NEST
    RINEST <- R2
    RONEST <- RINEST + TKNEST
}
DENEST <- 1  # Density nest material (kg/m3) 
THCONW <- 0.071  # Nest wall (wood: 0.10-0.35;sheep wool:0.05) thermal conductivity (W/m-C) #!
ABSHEL <- 0.71  # Nest solar absorptivity (decimal: 1.0 = 100%)
EMISHEL <- 0.95  # NEest emissivity
SHELEN <- 0.3  # Length(m) #!

GEOM.lab <- c("VOL", "D", "MASFAT", "VOLFAT", "ALENTH", "AWIDTH", 
    "AHEIT", "ATOT", "ASIL", "ASILN", "ASILP", "GMASS", "AREASKIN", 
    "FLSHVL", "FATTHK", "ASEMAJ", "BSEMIN", "CSEMIN", "CONVSK", 
    "CONVAR", "R1", "R2")
kable(cbind(GEOM.lab, t(GEOM.out)))
GEOM.lab
VOL 0.01
D 0.215443459147026
MASFAT 2
VOLFAT 0.00221975582685905
ALENTH 0.529649771785497
AWIDTH 0.196166582142777
AHEIT 0.196166582142777
ATOT 0.308349671796026
ASIL 0.0615297646245973
ASILN 0.0884364061999146
ASILP 0.0346231230492801
GMASS 10000
AREASKIN 0.269773486865498
FLSHVL 0.00778024417314096
FATTHK 0.00980642236948426
ASEMAJ 0.264824885892749
BSEMIN 0.0980832910713883
CSEMIN 0.0980832910713883
CONVSK 0.262885895484012
CONVAR 0.277514704616423
R1 0.0980832910713883
R2 0.107083291071388

Plot of current animal configuration.

F_FACTOR

This subroutine computes animal configuration factors for the sky and overhead vegetation, given the shade level.

F_FACTOR.lab
FAVEG 0
FASKY 0.4
FAGRD 0.4
FANEST 0
C3 6.05620958918289e-09
C4 6.05620958918289e-09
C5 0
C6 0
C7 0

SOLAR

Compute the solar radiation absorbed by different parts of the organism and its environment.

# environmental variables
QSOLR <- QSOLR  # solar radiation, horizontal plane (W/m2) (defined above for F_FACTOR)
SHADE <- SHADE  # shade (fractional, 0-1) (defined above for F_FACTOR)
Z <- 20  # zenith angle of sun (degrees from overhead)
ABSSB <- 0.8  # solar absorptivity of substrate (fractional, 0-1)

# traits
AREA <- ATOT  # surface area for solar exchange, m2 (from GEOM)
ABSAND <- 0.8  # solar absorptivity of dorsal fur (fractional, 0-1)
ABSANV <- 0.8  # solar absorptivity of ventral fur (fractional, 0-1)
ASIL <- ASIL  # silhouette area normal to sun, m2 (from GEOM)
PCTDIF <- 0.15  # proportion of solar radiation that is diffuse (fractional, 0-1)
FASKY <- FASKY  # configuration factor to sky (-) (from F_FACTOR)
FATOBJ <- FATOBJ  # configuration factor to object (-) (from F_FACTOR)
FAVEG <- FAVEG  # configuration factor to vegetation (-) (from F_FACTOR)

# solar radiation normal to sun's rays
ZEN <- pi/180 * Z  # convert degrees to radians
if (Z < 90) {
    # compute solar radiation on a surface normal to the direct
    # rays of the sun
    CZ = cos(ZEN)
    QNORM = QSOLR/CZ
} else {
    # diffuse skylight only
    QNORM = QSOLR
}

SOLAR.out <- SOLAR(AREA, ABSAND, ABSANV, ABSSB, ASIL, PCTDIF, 
    QNORM, SHADE, QSOLR, FASKY, FATOBJ, FAVEG)

QSOLAR <- SOLAR.out[1]  # total (global) solar radiation (W)
QSDIR <- SOLAR.out[2]  # direct solar radiaton (W)
QSOBJ <- SOLAR.out[3]  # lateral diffuse solar radiation (W)
QSSKY <- SOLAR.out[4]  # diffuse solar radiation from sky (W)
QSRSB <- SOLAR.out[5]  # diffuse solar radiation reflected from substrate (W)
QSDIFF <- SOLAR.out[6]  # total diffuse solar radiation (W)
QDORSL <- SOLAR.out[7]  # dorsal direct solar radiation (W)
QVENTR <- SOLAR.out[8]  # ventral diffuse solar radiaton (W)

SOLAR.lab <- c("QSOLAR", "QSDIR", "QSOBJ", "QSSKY", "QSRSB", 
    "QSDIFF", "QDORSL", "QVENTR")
kable(cbind(SOLAR.lab, t(SOLAR.out)))
SOLAR.lab
QSOLAR 88.9278060828641
QSDIR 44.5254533442364
QSOBJ 0
QSSKY 14.8007842462092
QSRSB 29.6015684924185
QSDIFF 44.4023527386277
QDORSL 59.3262375904456
QVENTR 29.6015684924185

CONV

Compute convective heat exchange.

CONV.lab
QCONV -38.3826769592207
HC 13.8308624086326
HCFREE 4.18906008479261
HCFOR 0
HD 0.0166500419458664
HDFREE 0.00389858046772071
HDFORC 0.0127514614781457
ANU 115.050184709447
RE 14271.7036555494
GR 14679651.4926261
PR 0.712515300549836
RA 10459476.2952353
SC 0.595551257906283
BP 101325

SEVAP

Compute evaporative heat exchange.

SEVAP.lab
QSEVAP 9.01497449278671
WEYES 4.23023290455338e-08
WCUTHF 0
WCUTF 3.63184136245377e-06
WCUT 3.63184136245377e-06
WTFUR 0
QFSEVAP 0

SIMULSOL

Setup and the SIMULSOL algorithm to SIMULtaneously SOLve for skin temperature \(T_skin\) and the temperature of the fur/feather interface \(T_fa\) given an environment, shape/posture and skin wetness. SIMULSOL is the core algorithm for computing the heat balance.

# environment
FLTYPE <- 0  # FLUID TYPE: 0 = AIR; 1 = FRESH WATER; 2 = SALT WATER - need's to be looked at - only invoked in main program when the dive table is set upTC <- 37 # core temperature (°C)
TAREF <- TA  # 1.2 m reference air temperature (°C)
TGRD <- TA  # ground temperature (°C)
TSKY <- TA  # sky temperature (°C)
TCONDSB <- TA  # surface temperature for conduction (°C)
TBUSH <- TA  # bush temperature (°C)
TVEG <- TAREF  # assume vegetation casting shade is at 1.2 m (reference) air temperature (°C)
SKYIR <- C3 * (TSKY + 273.15)^4  # sky infrared incoming (W)
VEGIR <- C6 * (TVEG + 273.15)^4  # vegetation infrared incomming (W)
SKYRAD <- SKYIR + VEGIR
SKYIN <- SKYRAD
GRDIN <- C4 * (TGRD + 273.15)^4  # note, MK put C4 here wherease before it was just SIG
TLOWER <- TGRD
RH <- 20  # relative humidity (%)
VEL <- 1  # wind speed (m s-1)
BP <- BP  # Pa, negative means altitude is used (from CONV)
ELEV <- 0  # m

# physiology and morphology
SKINW <- 0  # part of the skin surface that is wet (%)
AK1 <- 0.9  # initial thermal conductivity of flesh (0.412 - 2.8 W/mK)
AK2 <- 0.23  # conductivity of fat (W/mK)
BAREVAP <- 0  # is evaporation partly from bare skin? (0 = no, 1 = yes, % defined with PCTSKINEVAP)
PCTBAREVAP <- 0  # surface area for evaporation that is skin, e.g. licking paws (%)
PCTEYES <- 0.03  # surface area made up by the eye (%) - make zero if sleeping
# behaviour
FLYHR <- 0  # is flight occuring this hour? (imposes forced evaporative loss)

# configuration factors
FATOBJ <- 0  # configuration factor to nearby object
FAGRD <- 0.5  # configuration factor to ground
FASKY <- 0.5  # configuration factor to sky
FAVEG <- 0  # this is for overhead veg (at TAREF)
FABUSH <- 0  # this is for veg below/around animal (at TALOC)

# reference configuration factors
FABUSHREF <- FABUSH  # nearby bush
FATOBJREF <- FATOBJ  # nearby object
FASKYREF <- FASKY  # sky
FAGRDREF <- FAGRD  # ground
FAVEGREF <- FAVEG  # vegetation

FURTHRMK <- 0  # user-specified fur thermal conductivity (W/mK), not used if 0
SHAPE <- 4  # shape, 1 is cylinder, 2 is sphere, 3 is plate, 4 is ellipsoid

# Initial values
TS <- TC  # CURRENT GUESS OF OBJECT SURFACE TEMPERATURE
TFA <- TA  # current guess of fur/air interface temperature

DIFTOL <- 0.001  # tolerance for SIMULSOL
SIMULSOL.out <- matrix(data = 0, nrow = 2, ncol = 15)  # vector to hold the SIMULSOL results for dorsal and ventral side

# repeat for each side, dorsal and ventral, of the animal
for (S in 1:2) {
    
    # Calculating solar intensity entering fur. This will depend
    # on whether we are calculating the fur temperature for the
    # dorsal side or the ventral side. The dorsal side will have
    # solar inputs from the direct beam hitting the silhouette
    # area as well as diffuse solar scattered from the sky and
    # objects. The ventral side will have diffuse solar scattered
    # off the substrate.
    
    # Setting config factors and solar depending on whether the
    # dorsal side (S=1) or ventral side (S=2) is being estimated.
    if (QSOLAR > 0) {
        if (S == 1) {
            FASKY <- FASKYREF/(FASKYREF + FATOBJREF + FAVEGREF)
            FATOBJ <- FATOBJREF/(FASKYREF + FATOBJREF + FAVEGREF)
            FAVEG <- FAVEGREF/(FASKYREF + FATOBJREF + FAVEGREF)
            FAGRD <- 0
            FABUSH <- 0
            if (FATOBJ == 0) {
                QSLR <- 2 * QSDIR + ((QSSKY/FASKYREF) * FASKY)
            } else {
                QSLR <- 2 * QSDIR + ((QSSKY/FASKYREF) * FASKY) + 
                  ((QSOBJ/FATOBJREF) * FATOBJ)
            }
        } else {
            FASKY <- 0
            FATOBJ <- 0
            FAVEG <- 0
            FAGRD <- FAGRDREF/(1 - FAGRDREF - FATOBJREF - FAVEGREF)
            FABUSH <- FABUSHREF/(1 - FAGRDREF - FATOBJREF - FAVEGREF)
            QSLR <- QVENTR/(1 - FASKYREF - FATOBJREF - FAVEGREF)
        }
    } else {
        QSLR <- 0
        if (S == 1) {
            FASKY <- FASKYREF/(FASKYREF + FATOBJREF + FAVEGREF)
            FATOBJ <- FATOBJREF/(FASKYREF + FATOBJREF + FAVEGREF)
            FAVEG <- FAVEGREF/(FASKYREF + FATOBJREF + FAVEGREF)
            FAGRD <- 0
            FABUSH <- 0
        } else {
            FASKY <- 0
            FATOBJ <- 0
            FAVEG <- 0
            FAGRD <- FAGRDREF/(1 - FAGRDREF - FATOBJREF - FAVEGREF)
            FABUSH <- FABUSHREF/(1 - FAGRDREF - FATOBJREF - FAVEGREF)
        }
    }
    # set fur depth and conductivity index for KEFARA, the
    # conductivity, is the average (1), front/dorsal (2),
    # back/ventral(3) of the body part
    if (QSOLR > 0 | ZFURD != ZFURV) {
        if (S == 1) {
            ZL <- ZFURD
            KEFF <- KEFARA[2]
        } else {
            ZL <- ZFURV
            KEFF <- KEFARA[3]
        }
    } else {
        ZL <- ZFUR
        KEFF <- KEFARA[1]
    }
    
    RDXDEP <- 1  # not used yet - relates to radiation through fur
    XR <- RDXDEP  # not used yet - relates to radiation through fur
    X <- RDXDEP  # not used yet - relates to radiation through fur
    RSKIN <- R1  # body radius (including fat), m
    RFLESH <- R1 - FATTHK  # body radius flesh only (no fat), m
    RFUR <- R1 + ZL  # body radius including fur, m
    D <- 2 * RFUR  # diameter, m
    RRAD <- RSKIN + (XR * ZL)  # effective radiation radius, m
    LEN <- ALENTH  # length, m
    
    # Correcting volume to account for subcutaneous fat
    if (SUBQFAT == 1 & FATTHK > 0) {
        VOL <- FLSHVL
    }
    
    # Getting compressed fur thermal conductivity
    AREACND <- ATOT * PCOND
    CD <- AREACND * ((KFURCMPRS/ZFURCOMP))
    
    # package up inputs
    FURVARS <- c(LEN, ZFUR, FURTHRMK, KEFF, BETARA, FURTST, ZL)
    GEOMVARS <- c(SHAPE, SUBQFAT, CONVAR, VOL, D, CONVAR, CONVSK, 
        RFUR, RFLESH, RSKIN, XR, RRAD, ASEMAJ, BSEMIN, CSEMIN, 
        CD)
    ENVVARS <- c(FLTYPE, TA, TS, TBUSH, TVEG, TLOWER, TSKY, TCONDSB, 
        RH, VEL, BP, ELEV, FASKY, FABUSH, FAVEG, FAGRD, QSLR)
    TRAITS <- c(TC, AK1, AK2, EMISAN, FATTHK, FLYHR, FURWET, 
        PCTBAREVAP, PCTEYES)
    
    # set IPT, the geometry assumed in SIMULSOL: 1 = cylinder, 2
    # = sphere, 3 = ellipsoid
    if (SHAPE %in% c(1, 3, 5)) {
        IPT <- 1
    }
    if (SHAPE == 2) {
        IPT <- 2
    }
    if (SHAPE == 4) {
        IPT <- 3
    }
    
    # call SIMULSOL
    SIMULSOL.out[S, ] <- SIMULSOL(DIFTOL, IPT, FURVARS, GEOMVARS, 
        ENVVARS, TRAITS, TFA, SKINW, TS)
}

SIMULSOL.out <- cbind(c(1, 2), SIMULSOL.out)
colnames(SIMULSOL.out) <- c("SIDE", "TFA", "TSKIN", "QCONV", 
    "QCOND", "QGENNET", "QSEVAP", "QRAD", "QSLR", "QRSKY", "QRBSH", 
    "QRVEG", "QRGRD", "QFSEVAP", "NTRY", "SUCCESS")
tSIMULSOL.out <- t(SIMULSOL.out)
colnames(tSIMULSOL.out) <- c("DORSAL", "VENTRAL")
kable(tSIMULSOL.out)
DORSAL VENTRAL
SIDE 1.0000000 2.0000000
TFA 39.8685652 31.0081465
TSKIN 37.9519603 34.9371225
QCONV 76.8841311 42.3778612
QCOND 1.7816783 1.4824179
QGENNET -1.8782209 4.0700645
QSEVAP 0.1434575 0.1169764
QRAD 37.9687017 19.2990675
QSLR 118.6524752 59.2031370
QRSKY 37.9687017 0.0000000
QRBSH 0.0000000 0.0000000
QRVEG 0.0000000 0.0000000
QRGRD 0.0000000 19.2990675
QFSEVAP 0.0000000 0.0000000
NTRY 1.0000000 1.0000000
SUCCESS 1.0000000 1.0000000

ZBRENT and RESPFUN

Finally, find a value of QRESP (heat lost by respiration) to balance the heat budget by solving RESPFUN with ZBRENT.

# define basal metabolic rate
QBASAL <- (70 * AMASS^0.75) * (4.185/(24 * 3.6))  # heat generation (W)
DELTAR <- 0  # offset between air temeprature and breath (°C)
O2GAS <- 20.95  # oxygen concentration of air (%)
N2GAS <- 79.02  # nitrogen concetration of air (%)
CO2GAS <- 0.03  # carbon dioxide concentration of air (%)
RQ <- 0.8  # respiratory quotient (fractional, typically 0.7 (fats) to 1 (carbs), with 0.8 typical for protein)
EXTREF <- 20  # O2 extraction efficiency (%)
RELXIT <- 100  # relative humidity of exhaled air (%)
TIMACT <- 1  # multiplier on metabolic rate for activity costs
PANT <- 1  # multiplier on breathing rate (-)

# Now compute a weighted mean heat generation for all the
# parts/components = (dorsal value
# *(FASKY+FAVEG+FATOBJ))+(ventral value*FAGRD)
GEND <- SIMULSOL.out[1, 6]
GENV <- SIMULSOL.out[2, 6]
DMULT <- FASKYREF + FAVEGREF + FATOBJ
VMULT <- 1 - DMULT  # Assume that reflectivity of veg below = ref of soil so VMULT left as 1 - DMULT
X <- GEND * DMULT + GENV * VMULT  # weighted estimate of metabolic heat generation

# lung temperature and temperature of exhaled air
TLUNG <- (TC + (SIMULSOL.out[1, 3] + SIMULSOL.out[1, 3]) * 0.5) * 
    0.5  # average of skin and core
TAEXIT <- min(TA + DELTAR, TLUNG)  # temperature of exhaled air, °C

QMIN <- QBASAL
QM1 <- X - (5 * QMIN)
QM2 <- X + (10 * QMIN)
QSUM <- X
TOL <- AMASS * 0.01

ZBRENT.in <- c(TA, O2GAS, N2GAS, CO2GAS, BP, QMIN, RQ, TLUNG, 
    GMASS, EXTREF, RH, RELXIT, TIMACT, TAEXIT, QSUM, PANT)

# call ZBRENT subroutine which calls RESPFUN
ZBRENT.out <- ZBRENT(QM1, QM2, TOL, ZBRENT.in)

colnames(ZBRENT.out) <- c("RESPFN", "QRESP", "GEVAP", "PCTO2", 
    "PCTN2", "PCTCO2", "RESPGEN", "O2STP", "O2MOL1", "N2MOL1", 
    "AIRML1", "O2MOL2", "N2MOL2", "AIRML2", "AIRVOL")
tZBRENT.out <- t(ZBRENT.out)
colnames(tZBRENT.out) <- c("OUTPUT")
kable(tZBRENT.out)
OUTPUT
RESPFN -0.0994386
QRESP 0.8442388
GEVAP 0.0003500
PCTO2 0.2095000
PCTN2 0.7902000
PCTCO2 0.0003000
RESPGEN 1.8407220
O2STP 0.0010124
O2MOL1 0.0002260
N2MOL1 0.0008524
AIRML1 0.0010787
O2MOL2 0.0001808
N2MOL2 0.0008524
AIRML2 0.0010353
AIRVOL 0.0241634