R/egg_water.R

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

#' Tracy et al. 1978 egg water budget model (requires the deSolve package)
#'
#' Function for computing the exchange of gaseous and liquid water between an egg and soil
#' under environmental conditions that vary with time, using interpolation functions to
#' estimate environmental conditions at particular time intervals.
#' Michael Kearney developed this R function based on Tracy, C. R., Packard, G. C.,
#' and Packard, M. J. (1978). Water relations of chelonian eggs. Physiological
#' Zoology 51, 378–387. Parameter examples come from this paper.
#' @param t time intervals (s) at which output is required
#' @param m_init = 5.8 / 1000, mass of freshly laid egg (kg)
#' @param psi_e_init = -707, water potential of freshly laid egg (J/kg)
#' @param shape = 0.6, ratio of minor axes major axis of ellipsoid (-)
#' @param f_air = 0.5, fraction of egg surface exposed to air
#' @param K_e = 1.12 * 60 * 24 / 1e6 / (3600 * 24 * 10), hydraulic conductance of egg (kg m-2 s-1 (J/kg)-1), converted from original μg cm-2 min-1 bar-1
#' @param spec_hyd = 0.5, water potential-specific hydration (m3 m-3 (J/kg)-1)
#' @param pct_wet = 0.24, percent of surface area acting as a free-water surface for evaporation (\%)
#' @param P_e = -0.5052209, soil air entry potential (J/kg), derived from soil texture data via the function 'pedotransfer' and used to get soil hydraulic conductivity
#' @param b = 1.41005, Campbell's b parameter (-), derived from soil texture data via the function 'pedotransfer' and used to get soil hydraulic conductivity
#' @param K_sat = 0.003733307, soil saturated hydraulic conductivity (kg m-1 s-1 (J/kg)-1) derived from soil texture data via the function 'pedotransfer' and used to get soil hydraulic conductivity
#' @param elev = 0, elevation (m), needed for properties of air
#' @param vel = 0.001, wind speed (m/s)
#' @param Tsoilf soil temperature function with time, generated by 'approxfun' (°C)
#' @param Tbf body temperature function with time, generated by 'approxfun'(°C)
#' @param RHsoilf relative humidity function with time, generated by 'approxfun' (fractional)
#' @param PSIsoilf soil water potential with time, generated by 'approxfun' (J/kg)
#' @return m Egg mass (kg)
#' @return psi_e Egg water potential (J/kg)
#' @usage ode(y = c(m_init, psi_e_init), times = t, func = egg_water, parms = indata)
#' @examples
#' library(NicheMapR)
#'
#' # egg functional traits
#' m_init <- 5.8 / 1000 # kg, mass of freshly laid egg
#' psi_e_init <- -707 # J/kg, water potential of freshly laid egg
#' shape <- 0.6 # -, ratio of minor axes major axis of ellipsoid (-)
#' f_air <- 0.5 # -, fraction of egg surface exposed to air
#' K_e <- 1.12 * 60 * 24 / 1e6 / (3600 * 24 * 10) # kg m-2 s-1 (J/kg)-1, hydraulic conductance of egg
#' spec_hyd <- 0.0304 / 100 # m3 m-3 (J/kg)-1, water potential-specific hydration, empirically determined from linear model
#' pct_wet <- 0.25 # %, percent of surface area acting as a free-water surface for evaporation
#'
#' # get soil properties based on texture
#' pct_sand <- 99 # %
#' pct_silt <- 0 # %
#' pct_clay <- 100 - pct_sand # %
#' bulkdens <- 1.3 # Mg/m3
#' soilpro <- matrix(data = c(0, bulkdens, pct_clay, pct_silt, pct_sand), nrow = 1, ncol = 5)
#' colnames(soilpro) <- c('depth', 'blkdens', 'clay', 'silt', 'sand')
#'
#' #Now get hydraulic properties for this soil using Cosby et al. 1984 pedotransfer functions.
#' soil.hydro <- pedotransfer(soilpro = as.data.frame(soilpro), DEP = 0, model = 2)
#' P_e <- soil.hydro$PE # J/kg
#' b <- soil.hydro$BB # -
#' K_sat <- soil.hydro$KS # kg s/m3
#'
#' t <- seq(0, 24 * 3600 * 150, 24 * 3600) # seconds
#'
#' elev <- 0 # m, need elevation for air properties
#' vel <- 0.001 # m/s, wind speed
#' T_soil <- 29 # deg C, soil temperature
#' T_egg <- T_soil + 0.15 # deg C, egg temperature
#' RH_soil <- 0.99 # fractional, relative humidity
#' PSI_soil <- -100 # J/kg, soil water potential
#'
#' # use approxfun to create interpolations for the required environmental variables
#' Tsoilf <- approxfun(t, rep(T_soil, length(t)), rule = 2) # deg C, soil temperature
#' Tbf <- approxfun(t, rep(T_egg, length(t)), rule = 2) # deg C, body temperature in soil
#' RHsoilf <- approxfun(t, rep(RH_soil, length(t)), rule = 2) # fractional, relative humidity
#' PSIsoilf <- approxfun(t, rep(PSI_soil, length(t)), rule = 2) # J/kg, soil water potential
#'
#' # input parameters
#' indata <- list(shape = shape, K_e = K_e, spec_hyd = spec_hyd, f_air = f_air, pct_wet = pct_wet, elev = elev, vel = vel, K_sat = K_sat, P_e = P_e, b = b)
#'
#' # solve
#' egg_ode <- as.data.frame(deSolve::ode(y = c(m_init, psi_e_init), times = t, func = egg_water, parms = indata))
#' colnames(egg_ode) = c("time", "m", "psi_e")
#'
#' # plot results
#' par(mfrow = c(2, 1))
#' with(egg_ode, plot(t / (24 * 3600), m * 1000, type = 'l', ylab = 'egg mass, g', xlab = 'time, days', col = 'black'))
#' with(egg_ode, plot(t / (24 * 3600), psi_e, type = 'l', ylab = 'egg water potential, J / kg', xlab = 'time, days', col = 'black', ylim = c(psi_e_init, PSI_soil)))
#' abline(h = PSI_soil, lty = 2, col = 'blue')
#' @export
egg_water <- function(t, y, indata) {
  with(as.list(c(indata, y)), {

    m <- as.numeric(y[1]) # kg, mass
    psi_e <- as.numeric(y[2]) # J/kg, water potential of egg


    # area calculation for ellipsoid
    get_area <- function(a = 1, b = b, mass = mass){ # from GEOM.f in ectotherm model
      ratio <- a / b
      volume <- mass / 1e6 # cm3 to m3
      b <- ((3 * volume) / (4 * 3.14159 * ratio)) ^ 0.333 # semi-minor axis 1
      c <- b # semi-minor axis 2
      a <- b * ratio # semi-major axis
      # formula for surface area of an ellipsoid
      Ecc <- sqrt(a ^ 2 - c ^ 2) / a # eccentricity
      A_tot <- (2 * pi * b ^ 2 + 2 * pi * ((a * b) / Ecc) * asin(Ecc)) # m2
    }

    # calculate areas
    A_tot <- get_area(b = shape, mass = m * 1000)
    A_e <- A_tot * f_air # m2, surface area exposed to air
    A_s <- A_tot * (1 - f_air) # m2, surface area in contact with soil

    # evaporative exchange sub function
    egg_evap <- function(T_air, T_egg, RH, vel, elev, m, A_e, pct_wet){
      # constants
      G <- 9.80665	# acceleration from gravity, m/s
      C_P <- 1.01E+03	# specific heat dry air, J/(kg K)
      BETA <-	1 / (T_air + 273.15) # coeff thermal expansion at constant density, 1/K
      rho_flesh <- 1000 # density of flesh (kg / m3)

      # air properties
      dryair <- DRYAIR(db = T_air, alt = elev)
      bp <- dryair$patmos # barometric pressure, Pa
      rho_dry <- dryair$densty # density, kg/m3
      mu <- dryair$visdyn # dynamic viscosity, kg/(m s)
      k_air <- dryair$thcond # thermal conductivity, W/(m K)
      D_H2O <- dryair$difvpr # diffusivity of water vapour in air, m2/s

      # free convection
      VOL <- m / rho_flesh # animal volume
      D <- VOL ^ (1 / 3) # m, characteristic dimension for convection
      PR_free <- C_P * mu / k_air # Prandlt number for free convection, -
      SC_free <- mu / (rho_dry * D_H2O) # Schmidt number for free convection, -
      delta_T <- T_egg - T_air # temperature difference between skin and air, °C
      delta_T[abs(delta_T) < 0.01] <- 0.01 * sign(delta_T[abs(delta_T) < 0.01]) # impose lower limit of 0.01 °C
      GR <- abs(((rho_dry ^ 2) * BETA * G * (D ^ 3) * delta_T) / (mu ^ 2)) # Grashof number, -
      RE <- rho_dry * vel * D / mu # Reynold's number, -
      RA <- (GR ^ (1/4)) * (PR_free ^ (1/3)) # Rayleigh number, - (formula for sphere or lizard)
      NU_free <- 2 + 0.60 * RA # Nusslet number for free convection, -
      HC_free <- (NU_free * k_air) / D	# Heat transfer coeff for free convection, W/(m2 K)
      HC_free[HC_free < 5] <- 5 # don't let it go below 5, as this is unrealistic
      NU_free <- HC_free * D / k_air # recalcualting NU_free in case HC_free had gone below lower limit of 5
      SH_free <- NU_free * (SC_free / PR_free) ^ (1/3) # Sherwood number, -
      HD_free <- SH_free * D_H2O / D # mass transfer coeff for free convection, m/s
      HD <- HD_free # Final mass transfer coefficient, m/s

      # vapour densities
      rho_H2O_air <- WETAIR(db = T_air, rh = RH)$vd # vapour density of air, kg/m3
      MW <- 0.018 # molar mass of water, kg/mol
      R <- 8.314 # gas constant, J/mol/K
      RH_egg <- exp(psi_e / (R / MW * (T_egg + 273.15))) * 100 # %, relative humidity
      rho_H2O_skin <- WETAIR(db = T_egg, rh = RH_egg)$vd # vapour density of egg, kg/m3

      # area wet
      A_eff <- A_e * pct_wet / 100 # get effective wet area

      # 'cutaneous' water loss rate
      E_cut <- A_eff * HD * (rho_H2O_skin - rho_H2O_air) # water loss rate, kg/s
    }

    # evaporative water exchange calculations
    T_air <- Tsoilf(t) # deg C, interpolate soil temperature to current time step
    T_egg <- Tbf(t) # deg C, Tb prediction including metabolic heat
    RH <- RHsoilf(t) * 100 # %, relative humidity
    if(RH > 100){RH <- 100}
    if(RH < 0){RH <- 0}
    m_a <- egg_evap(T_air, T_egg, RH, vel, elev, m, A_e, pct_wet) # kg/s, evaporative water exchange

    # liquid water exchange calculations
    psi_s <- min(0, PSIsoilf(t)) # interpolate water potential to current time, J/kg
    k_s <- K_sat * (P_e / psi_s) ^ (2 + 3 / b) # equation 9.2 from Campbell and Norman 2001
    m_s <- (psi_s - psi_e) / (1 / (A_s * K_e) + 1 / ((2 * pi * A_s) ^ (1 / 2) * k_s)) # kg/s, liquid water exchange (Darcy's law)

    # combined exchanges and egg water potential
    d_m <- m_s - m_a # kg s-1, change in egg mass
    d_psi_e <- d_m / (spec_hyd * m) # J/kg s-1, change in egg water potential

    list(y = c(d_m[1], d_psi_e[1]))
  })
}