diff --git a/DESCRIPTION b/DESCRIPTION index 2cb01fa..f39a277 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: LakeMetabolizer Title: Tools for the Analysis of Ecosystem Metabolism Maintainer: Luke Winslow -Version: 1.5.2 +Version: 1.5.3 Author: Luke Winslow, Jake Zwart, Ryan Batt, Jessica Corman, Hilary Dugan, Paul Hanson, Aline Jaimes, Jordan Read, Richard Woolway Description: A collection of tools for the calculation of freewater metabolism. diff --git a/NEWS.md b/NEWS.md index cc7c000..a4e64f5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ # Version 1.5.1 (2017-02-12) -* Fix `k.heiskanen` bug where it fails when trying to calculate lwnet \ No newline at end of file +* Fix `k.heiskanen` bug where it fails when trying to calculate lwnet + +# Version 1.5.3 (2020-04-23) + +* Fixes `k.read`, `k.read.soloviev`, and `k.macIntyre` bug where base functions fail when calculating kinematic viscosity diff --git a/R/k.macIntyre.R b/R/k.macIntyre.R index 7c125f1..99aabb7 100644 --- a/R/k.macIntyre.R +++ b/R/k.macIntyre.R @@ -1,79 +1,79 @@ -# ---Author: Hilary Dugan, 2013-10-20 --- +# ---Author: Hilary Dugan, 2013-10-20 --- # modified from K600 code by Richard Woolway and Jordan Read # INPUTS; # wnd.z: Height of anemometer # Kd: Diffuse attenuation coefficient -# atm.press: Atmospheric pressure in hPa or Mb. If unknown use '1013' for sea level +# atm.press: Atmospheric pressure in hPa or Mb. If unknown use '1013' for sea level # dateTime: date and time vector -# wtr: dataframe of water temperatures +# wtr: dataframe of water temperatures # depth: vector of water temperature depths # airT: vector of air temperatures in C # Uz: vector of wind speed in m/s # Rh: vector of relative humidity in % # sw: vector of shortwave radiation in W/m2 # lw: vector of longwave radiation. If missing, run calc.lw.net function first -# par: vector of par data in umol m-2 s-1 +# par: vector of par data in umol m-2 s-1 -# OUTPUT: returns the gas exchange velocity for O2 in units of m/(timeStep*min) (i.e. 30 minute sampling +# OUTPUT: returns the gas exchange velocity for O2 in units of m/(timeStep*min) (i.e. 30 minute sampling # interval will return kO2 in units of m/(1/48) - converts to fraction of day) #'@export k.macIntyre = function(ts.data, wnd.z, Kd, atm.press,params=c(1.2,0.4872,1.4784)){ - + S_B <- 5.67E-8 # Stefan-Boltzman constant (K is used) emiss <- 0.972 # emissivity; Kelvin = 273.15 #conversion from C to Kelvin - - # Get short wave radiation data - if(has.vars(ts.data, 'sw')){ + + # Get short wave radiation data + if(has.vars(ts.data, 'sw')){ sw <- get.vars(ts.data, 'sw') - + } else if (has.vars(ts.data, 'par')){ tmp.par = get.vars(ts.data, 'par') sw = par.to.sw(tmp.par) - } else { + } else { stop("Data must have PAR or SW column\n") } - + # Get water temperature data wtr <- get.vars(ts.data, 'wtr') Ts <- get.Ts(ts.data) - + airT <- get.vars(ts.data, 'airt') - + RH <- get.vars(ts.data, 'rh') - + # Get long wave radiation data - if(has.vars(ts.data, 'lwnet')){ + if(has.vars(ts.data, 'lwnet')){ lwnet <- get.vars(ts.data,'lwnet') - + } else if(has.vars(ts.data, 'lw')){ lw_in <- get.vars(ts.data, 'lw') # long wave in Tk <- Ts+Kelvin # water temperature in Kelvin LWo <- S_B*emiss*Tk^4 # long wave out lwnet <- lw_in[,2]-LWo - + lwnet = data.frame(datetime=lw_in$datetime, lwnet=lwnet) - - } else { + + } else { stop("no longwave radiation available") } - + wnd <- get.vars(ts.data, 'wnd') m.d = ts.meta.depths(wtr) - - k600 = k.macIntyre.base(wnd.z, Kd, atm.press, ts.data$datetime, Ts[,2], m.d$top, + + k600 = k.macIntyre.base(wnd.z, Kd, atm.press, ts.data$datetime, Ts[,2], m.d$top, airT[,2], wnd[,2], RH[,2], sw[,2], lwnet[,2],params) - + return(data.frame(datetime=ts.data$datetime, k600=k600)) - + } #'@export k.macIntyre.base <- function(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet, params=c(1.2,0.4872,1.4784)){ - + #Constants S_B <- 5.67E-8 # Stefan-Boltzman constant (K is used) emiss <- 0.972 # emissivity; @@ -95,14 +95,14 @@ k.macIntyre.base <- function(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wn C_D <- mm$C_D # drag coefficient for momentum E <- mm$alh # latent heat flux H <- mm$ash # sensible heat flux - + # calculate total heat flux dUdt <- sw*0.93 - E - H + lwnet Qo <- sw*(1-albedo_SW)*swRat - + # calculate water density rho_w <- water.density(Ts) - + # calculate u* if (wnd.z != 10) { e1 <- sqrt(C_D) @@ -110,23 +110,23 @@ k.macIntyre.base <- function(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wn }else{ u10 = wnd } - + rhoAir <- 1.2 # air density vonK <- 0.41 # von Karman constant tau <- C_D*u10^2*rhoAir uSt <- sqrt(tau/rho_w) - + # calculate the effective heat flux q1 <- 2-2*exp(z.aml*-Kd) q2 <- z.aml*Kd q3 <- exp(z.aml*-Kd) H_star <- dUdt-Qo*(q1/q2-q3) # Kim 1976 - - - # calculate the thermal expansion coefficient - thermalExpFromTemp <- function(Ts) { - V <- 1 - dT <- 0.001 + + + # calculate the thermal expansion coefficient + thermalExpFromTemp <- function(Ts) { + V <- 1 + dT <- 0.001 T1 <- water.density(Ts) T2 <- water.density(Ts+dT) V2 <- T1/T2 @@ -135,12 +135,12 @@ k.macIntyre.base <- function(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wn return (alpha) } tExp <- thermalExpFromTemp(Ts) - + B1 = H_star*tExp*g B2 = rho_w*C_w Bflx = B1/B2 - - + + # calculate kinematic viscosiy getKinematicVis <- function(Ts) { # from Mays 2005, Water Resources Engineering @@ -149,26 +149,26 @@ k.macIntyre.base <- function(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wn visTable <- c(1.792,1.519,1.308,1.141,1.007,0.897, 0.804,0.727,0.661,0.605,0.556,0.513,0.477,0.444, 0.415,0.39,0.367,0.347,0.328,0.311,0.296) - v <- data.frame(approx(tempTable,visTable,xout = Ts))[2] + v <- approx(tempTable,visTable,xout = Ts)$y v <- v*1e-6 return(v) } kinV <- getKinematicVis(Ts) KeNm = uSt^3 - + #SmE = 0.84*(-0.58*Bflx+1.76*KeNm/(vonK*z.aml)) SmE = params[1]*-Bflx+params[2]*KeNm/(vonK*z.aml) #change to two coefficients SmE[SmE<0] = 0 # set negative to 0 Sk = SmE*kinV Sk = Sk*100^4*3600^4 # Sally's K now in cm4/h4 - Sk600 = params[3]*600^(-0.5)*Sk^(1/4) # in cm/hr (Total) - + Sk600 = params[3]*600^(-0.5)*Sk^(1/4) # in cm/hr (Total) + k600 <- as.numeric(Sk600) # why is this not already numeric? k600 <- k600*24/100 #units in m d-1 return(k600) } -# -- References +# -- References #MACINTRYE, sALLY, ANDERS JONSSON, MATS JANSSON, JAN ABERG, DAMON E. TURNEY AND SCOTT D. MILLER. #2010. Buoyancy flux, turbulence, and the gas transfer coefficient in a stratified lake. -#Geophysical Research Letters. 37: L24604. doi:10.1029/2010GL044164 +#Geophysical Research Letters. 37: L24604. doi:10.1029/2010GL044164 diff --git a/R/k.read.soloviev.R b/R/k.read.soloviev.R index ac40956..453194d 100644 --- a/R/k.read.soloviev.R +++ b/R/k.read.soloviev.R @@ -1,63 +1,63 @@ #'@export k.read.soloviev = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ - + S_B <- 5.67E-8 # Stefan-Boltzman constant (K is used) emiss <- 0.972 # emissivity; Kelvin = 273.15 #conversion from C to Kelvin - + data = ts.data - # Get short wave radiation data - if(has.vars(data, 'sw')){ + # Get short wave radiation data + if(has.vars(data, 'sw')){ sw <- get.vars(data, 'sw') - + } else if (has.vars(data, 'par')){ tmp.par = get.vars(data, 'par') sw = par.to.sw(tmp.par) - } else { + } else { stop("Data must have PAR or SW column\n") } - + wtr <- get.vars(data,'wtr') Ts <- get.Ts(data) airT <- get.vars(data, 'airt') - + RH <- get.vars(data, 'rh') - + # Get long wave radiation data - if(has.vars(data, 'lwnet')){ + if(has.vars(data, 'lwnet')){ lwnet <- get.vars(data,'lwnet') - + } else if(has.vars(data, 'lw')){ lw_in <- get.vars(data, 'lw') # long wave in Tk <- Ts+Kelvin # water temperature in Kelvin LWo <- S_B*emiss*Tk^4 # long wave out lwnet <- lw_in[,2]-LWo - + lwnet = data.frame(datetime=lw_in$datetime, lwnet=lwnet) - - } else { + + } else { stop("no longwave radiation available") } - + wnd <- get.vars(data, 'wnd') - + m.d = ts.meta.depths(wtr) - - k600 = k.read.soloviev.base(wnd.z, Kd, lat, lake.area, atm.press, data$datetime, Ts[,2], m.d$top, + + k600 = k.read.soloviev.base(wnd.z, Kd, lat, lake.area, atm.press, data$datetime, Ts[,2], m.d$top, airT[,2], wnd[,2], RH[,2], sw[,2], lwnet[,2]) - + return(data.frame(datetime=data$datetime, k600=k600)) } #'@export -k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet){ - +k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet){ + # define constants used in function - Kelvin <- 273.15 # temp mod for deg K + Kelvin <- 273.15 # temp mod for deg K emiss <- 0.972 # emissivity; S_B <- 5.67E-8 # Stefan-Boltzman constant (?K is used) vonK <- 0.41 # von Karman constant @@ -70,30 +70,30 @@ k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, g <- 9.81 # gravity C_w <- 4186 # J kg-1 ?C-1 (Lenters et al. 2005) mnWnd <- 0.2 # minimum wind speed - + # impose limit on wind speed rpcI <- wnd < mnWnd wnd[rpcI] <- mnWnd - + # calculate sensible and latent heat fluxes mm <- calc.zeng(dateTime,Ts,airT,wnd,RH,atm.press,wnd.z) C_D <- mm$C_D # drag coefficient for momentum E <- mm$alh # latent heat flux H <- mm$ash # sensible heat flux - + # calculate total heat flux dUdt <- sw*0.93 - E - H + lwnet Qo <- sw*(1-albedo_SW)*swRat #PAR - + # calculate the effective heat flux q1 <- 2-2*exp(z.aml*-Kd) q2 <- z.aml*Kd q3 <- exp(z.aml*-Kd) H_star <- dUdt-Qo*(q1/q2-q3) #Effective surface heat flux Kim 1976 - + # calculate water density rho_w <- water.density(Ts) - + # calculate u* if (wnd.z != 10) { e1 <- sqrt(C_D) @@ -103,11 +103,11 @@ k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, tau <- C_D*wnd^2*rhoAir uSt <- sqrt(tau/rho_w) uSta <- sqrt(tau/rhoAir) #friction velocity in air - - # calculate the thermal expansion coefficient - thermalExpFromTemp <- function(Ts) { - V <- 1 - dT <- 0.001 + + # calculate the thermal expansion coefficient + thermalExpFromTemp <- function(Ts) { + V <- 1 + dT <- 0.001 T1 <- water.density(Ts) T2 <- water.density(Ts+dT) V2 <- T1/T2 @@ -116,14 +116,14 @@ k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, return (alpha) } tExp <- thermalExpFromTemp(Ts) - + # calculate buoyancy flux and w* - B1 <- H_star*tExp*g #Hstar * coefficient of thermal expansion * gravity + B1 <- H_star*tExp*g #Hstar * coefficient of thermal expansion * gravity B2 <- rho_w*C_w Bflx <- B1/B2 Bflx[Bflx>0] = 0 - wSt <- (-Bflx*z.aml)^1/3 + wSt <- (-Bflx*z.aml)^1/3 # calculate kinematic viscosiy getKinematicVis <- function(Ts) { @@ -133,13 +133,13 @@ k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, visTable <- c(1.792,1.519,1.308,1.141,1.007,0.897, 0.804,0.727,0.661,0.605,0.556,0.513,0.477,0.444, 0.415,0.39,0.367,0.347,0.328,0.311,0.296) - v <- data.frame(approx(tempTable,visTable,xout = Ts))[2] + v <- approx(tempTable,visTable,xout = Ts)$y v <- v*1e-6 return(v) } kinV <- getKinematicVis(Ts) kinVa <- getKinematicVis(airT) - + KeDe <- (kinV*g) KeNm <- uSt^3 Ke <- KeNm/KeDe @@ -147,28 +147,28 @@ k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, euPw <- (1+Ke/KeCrit) # tau_shear = tau/(1+Ke/Kecr) Ke is the Keulegan number # Could calculate KeCrit (critical Keulegan number) from wave age #KeCrit <- (kinVa/kinV)*((rhoAir/rho_w)^1.5)*(Rbcr/Aw) # Eq1.16-Soloviev et al(2007) - + tau_t <- tau/euPw # tau_t = tangential wind stress, tau = total wind stress - uTanS <- tau_t/rho_w + uTanS <- tau_t/rho_w uTanS <- uTanS^0.5 - + # calculate viscous sublayer Sv <- C1*kinV/uTanS # effective thickness of the aqueous viscous sublayer eu_N <- uTanS^3 # e_u(0) = (tau_t/rho)^1.5/(vonK*Sv) eu_D <- vonK*Sv # denominator eu_0 <- eu_N/eu_D # in m2/s3 ec_0 <- -1.0*Bflx # buoyancy flux, but only when outward - + #ewave_0 turbulence due to wave breaking lake.area <- lake.area/1e6 # convert surface area to km Fetch <- 2*sqrt(lake.area/pi) # fetch in km (assuming a conical lake) - Hs <- 0.0163*(Fetch^0.5)*wnd # significant wave height - Woolf (2005) + Hs <- 0.0163*(Fetch^0.5)*wnd # significant wave height - Woolf (2005) Aw <- (1/(2*pi))*(( (g*Hs*rhoAir)/(0.062*rho_w* uSt^2))^(2/3)) # wave age - eqn 1.11 Soloviev et al. (2007) W <- 3.8e-6*wnd^3.4 # simplified whitecap fraction (Fariall et al. 2000) - - + + Ap <- 2.45*W*((1/(W^0.25))-1) alphaW <- 100 # p. 185 - Soloviev et al. (2007) B <- 16.6 # p. 185 - Soloviev et al. (2007) @@ -178,23 +178,23 @@ k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, (((Ke/KeCrit)^1.5)/((1+Ke/KeCrit)^1.5)) * (uSt*g*kinV)/(0.062*vonK*cT*((2*pi*Aw)^1.5)) * (rhoAir/rho_w) - + #------------------------------------ e_0 <- ec_0+eu_0+ewave_0 # e(0) from Soloviev (w/o wave effects) Kc <- ec_0*kinV*100^4*3600^4 # convective component now in cm4/hr4 (Total) Ku <- eu_0*kinV*100^4*3600^4 # shear component now in cm4/hr4 (Total) Kwave <- ewave_0*kinV*100^4*3600^4 # wave component now in cm4/hr4 (Total) Kall <- e_0*kinV*100^4*3600^4 # turbulent kinetic energy now in cm4/hr4 (Total) - - #Schmidt number could be calculated as temperature dependent + + #Schmidt number could be calculated as temperature dependent #Sc <- 1568+(-86.04*Ts)+(2.142*Ts^2)+(-0.0216*Ts^3) k600org <- nu*600^(-0.5)*(Kc+Ku)^(1/4) # in cm/hr (Total) k600org <- k600org*24/100 #now in units in m d-1 - + k600 <- nu*600^(-0.5)*Kall^(1/4) # in cm/hr (Total) k600 <- k600*24/100 #now in units in m d-1 - - # ---Breaking Wave Component, Author: R I Woolway, 2014-11-13 --- + + # ---Breaking Wave Component, Author: R I Woolway, 2014-11-13 --- # bubble mediated component - Woolf 1997 kbi <- W*2450 beta_0 <- 2.71*1e-2 # Ostwald gas solubility (Emerson and Hedges, 2008) @@ -203,7 +203,7 @@ k.read.soloviev.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, kb <- kbi/((beta_0*kbiii)) kb <- kb*24/100 #units in m d-1 #---------------------------------------------------------------- - + k600b = k600+kb allks = data.frame(Ku,Kc,Kwave,kb,k600org,k600,k600b) colnames(allks) = c("shear","convective","wave","bubble",