It may be worth adding a wet-bulb temperature to the available variables in standard netCDF files produced by aircraft projects. Reasons include these:
Below cloud, the wet-bulb temperature is a useful indicator of the potential for rain from the cloud to produce downdrafts by evaporative cooling.
The wet-bulb temperature is an important indicator of human comfort and the ability of the human body to control its temperature in hot conditions. Wet-bulb temperatures in excess of \(90^\circ\)F are significant health hazards, and there will be an extreme hazard if the wet-bulb temperature exceeds normal body temperature because perspiration cannot control the body temperature. If global warming increases the likelihood of these extreme temperatures that can lead to significant new hazards. The close association between wet-bulb temperature and human comfort and safety makes it is surprising that wet-bulb temperature is almost never mentioned in connection with weather events. Instead, a substitute and unquantitative comment like “high temperatures are less comfortable when the humidity is high” is often heard without specific reference to the relevant quantitative measure of the hazard.
The calculation of the wet-bulb temperature is best done iteratively, so users can’t just apply an analytical formula to the measured temperature, pressure, and humidity. Therefore it might be useful for users to have this calculation already done.
Because it can also be useful to have the calculation available, a new function “wetbulbT()” has been added to the Ranadu functions. The code is shown below.
Consider a mixture of air (\(p\), \(T\), \(e\)) and liquid drops in an amount not sufficient to affect the specific heat of the mixture. If the vapor pressure is below the equilibrium value for the temperature \(T\), the water will evaporate, requiring balance between the latent heat consumed (\(L_vdr\)) by the change in mixing ratio \(r\) and the change in energy of the air caused by its change in temperature \(dT\), \((c_{pd}+rc_{pv})dT)\)). Both are expressed in change per unit mass of dry air, leading to \[L_v(T)dr=-c_p(r)dT\] Because \(L_v(T)\) and \(c_p(r)\) are variable depending on \(T\) and \(r\), appropriate values of each should be used. The liquid surface will be at the wet-bulb temperature \(T_{wb}\), so \(L_v^*=L_v(T_{wb})\) is the appropriate value to use for the latent heat of vaporization. Similarly, the air in contact with the liquid surface will have vapor pressure in equilibrium with that surface and so will have specific heat appropriate for that mixture: \(c_p^*=c_{pd}+r_s(T_{wb})c_{pv}\) where \(r_s(T)\approx 0.622e_s(T)/(p-e_s(T))\). The variation of \(c_p\) with \(r\) will affect the temperature profile near the water surface, but the value \(c_p^*\) will be the appropriate one to use to describe the balance at the surface of the water. These functions are known: \(L_v^*=L_0+L_1T_{wb}\) and \(c_p^*=c_{pd}+r_s(T_{wb})c_{pv}\) with \(L_0=2.501\times 10^6\) and \(L_1=2370\) if \(T\) is expressed in degrees Celsius.
This leads to the fundamental equation determining the wet-bulb temperature \(T_{wb}\): \[T_{wb}=T+\frac{L_v^*}{c_p^*}(r-r_s(T_{wb}))\] Because \(L_v^*\) and \(c_p^*\) depend on the wet-bulb temperature, an analytical solution is not available and the wet-bulb temperature must be found iteratively as in the following code chunk. Here the R routine “nleqslv()” is used to find the zero value of the twb() function:
wetbulbT <- function(P, AT, DPT) { # pressure [hPa], temperature and dewpoint [deg C]
twb <- function(tw) { # This function depends on the additional
# variables P, AT, and DPT being in the
# calling environment. The supplied argument tw is deg.C.
e <- MurphyKoop(DPT) # omitting dependence on p, as is conventional
e2 <- MurphyKoop(tw)
r <- MixingRatio(e / P)
rs <- MixingRatio(e2 / P)
Lv <- 2.501e6 - 2370 * tw
# Lv <- 2.501e6
cp <-SpecificHeats(e2 / P)[1]
# cp <- SpecificHeats()[1]
return(tw - AT - Lv * (r - rs) / cp) # for adjusting to zero return
}
## Avoid errors for missing values: This may be useful for processing vectors
# P <- zoo::na.approx (as.vector(P), na.rm=FALSE, rule = 2)
# AT <- zoo::na.approx (as.vector(AT), na.rm=FALSE, rule = 2)
# DPT <- zoo::na.approx (as.vector(DPT), na.rm=FALSE, rule = 2)
T1 <- 0.5 * (AT + DPT) # The first guess, input to nleqslv
X <- nleqslv(T1, twb) # X$x is the solution for tw
return(X$x)
}
wetbulbT(800, 30, 20) ## an example## [1] 22.45461The task is to find the value of tw for which twb(tw, …) is zero. Newton’s method is an iterative approach to the solution, where each successive approximation is given by \[x_{i+1}=x_i-\frac{f(x_i)}{f^\prime(x_i)}\] Because the function is not a simple function of the arguments, the derivative needed for Newton’s method is best determined numerically: e.g., for some small value of \(\delta\), \(twb^\prime(tw)=(twb(tw+\delta)-twb(tw-\delta))/(2\delta)\). Then \[tw_{i+1}=tw_i-twb(tw_i)\frac{2\delta}{twb(tw+\delta)-twb(tw-\delta)}\] If the preceding example is used along with the same first guess, successive iterations give this sequence of estimates:
twb <- function(tw) { # This function depends on calling-environment Tdry, p, and dp
# (i.e., temperature [degC], pressure [hPa], dewpoint [degC]
e <- MurphyKoop(DPT) # omitting dependence on p, as is conventional
e2 <- MurphyKoop(tw)
r <- MixingRatio(e / P)
rs <- MixingRatio(e2 / P)
Lv <- 2.501e6 - 2370 * tw
# Lv <- 2.501e6
cp <-SpecificHeats(e2 / P)[1]
# cp <- SpecificHeats()[1]
return(tw - AT - Lv * (r - rs) / cp) # for adjusting to zero return
}
AT <- 30
P <- 800
DPT <- 20
t1 <- 0.5 * (AT + DPT) ## first guess
d <- 0.0001
t <- rep(t1,6)
for (i in 1:5) {
t[i+1] <- t[i] - twb(t[i]) * 2 * d / (twb(t[i]+d)-twb(t[i]-d))
}
options(digits=9)
print(t)## [1] 25.0000000 22.5841506 22.4549537 22.4546108 22.4546108 22.4546108After only three iterations, the solution converges acceptably. This is typical. Therefore the following function is suggested for routine processing (after conversion to equivalent C++ code):
wetb <- function(P, AT, DPT) {
t2 <- 0.5 * (AT + DPT) # The first guess
d <- 0.0001
twb <- function(tw) { # This function depends on calling-environment AT, p, and DPT
# (i.e., temperature [degC], pressure [hPa], dewpoint [degC]
e <- MurphyKoop(DPT) # omitting dependence on p, as is conventional
e2 <- MurphyKoop(tw)
r <- MixingRatio(e / P)
rs <- MixingRatio(e2 / P)
Lv <- 2.501e6 - 2370 * tw
cp <-SpecificHeats(e2 / P)[1]
return(tw - AT - Lv * (r - rs) / cp) # for adjusting to zero return
}
knt <- 0
repeat {
t1 <- t2
t2 <- t1 - twb(t1) * 2*d / (twb(t1+d)-twb(t1-d))
# t2 <- 0.5 * (t1 + t2) ## might try something like this to improve stability
# print (sprintf('t1, t2 = %.7f %.7f', t1, t2)) ## comment this for normal usage
knt <- knt + 1 ## knt is a safeguard against runaway interation -- probably not needed.
if (knt > 10) {
print ('too many iterations in wetb; step test not met')
break
}
if (abs(t2-t1) < 0.0001) {break}
}
return(t2)
}
## Examples:
options(digits=5)
wetb(800,30,20)## [1] 22.455wetb(300, -40, -50)## [1] -40.613wetb(500,0,-10)## [1] -4.6282wetb(700,10,9)## [1] 9.3621wetb(1000,35,31)## [1] 31.769