Data Analysis / Holiday Projects

Climate Change & Heat Stress Indicators: Exploring the relaxation of Tw

In this series of posts I am exploring the relaxation of the natural wet bulb temperature (Tw). This is the limiting step in the approximation of heat stress indicator wet bulb global temperature (WBGT), which is widely used to model the effects of temperature upon humans (Bernard and Pourmoghani (1999); Lemke and Kjellstron (2012)). In this blog post we explore the iterative model for the relaxation of Tw.

The iterative relaxation of the wet bulb temperature (Tw) is fundamental to the calculation of heat stress indicator WBGT. While Tw can be measured, the equipment is notoriously inaccurate and it is not commonly recorded. Hence, considerable effort has gone into building models and approximations. Currently, the most commonly accepted method is an iterative algorithm that takes two inputs: atmospheric temperature (Ta) and the dew point temperature (Td). Tw is initialised at the same value as the dew point and then iteratively refined. A very basic function for the relaxation of Tw is given below (written in R):

calcE <- function (t) 
{
    return (0.6106 * exp(17.27 * t / (237.3 + t)))
}
calcWBGT <- function (ta, td) 
{
    # tw counters used as stopping conditions
    tw.previous <- 10000
    tw.current <- 10000

    # initialise ed and tw (tw := td + 0.2)
    ed <- calcE(td)
    tw <- td + 0.2
    
    # Relax tw
    stop.point <- FALSE
    while (!stop.point) {    

        # stopping conditions    
        stop.point <- if (tw.previous > 0 & tw.current < 0) TRUE
                      else if (tw >= ta) TRUE
                      else FALSE
        
        ew <- calcE(tw)
        tw.previous <- tw.current
        tw.current <- (ed - ew) * (1556 - 1.484*tw) + 101 * (ta - tw)

        # step forward tw
        tw <- tw + 0.2
    }
    
    tw <- if (tw > (td + 0.3)) tw - 0.3
          else td
    
    wbgt <- 0.67 * tw + 0.33 * ta
    return (wbgt)
      
}

To explore the relaxation of Tw we made a slight adjustment to the function above and recorded the intermediate values of tw at each step. These values were kept in a vector called iter.tw and plotted as shown below:

Relaxation of Tw with inputs: Ta=33.3, Td=-26.1

Relaxation of Tw with inputs: Ta=33.3, Td=-26.1

We can see from the plot above that the relaxation of Tw is a smooth, continuous decay. We should be able to model this with a high level of certainty. However, through simulation we were able to show that the decay function is strongly dependent on the inputs as shown in the figure below:

Relaxation of Tw with simulated values for Ta and Td.

Relaxation of Tw with simulated values for Ta and Td.

On the surface, the relaxation function appears to be very similar for each set of inputs. But closer inspection will reveal that the axes are all very different scales. The differing relaxation profiles become clear when we fix the axis:

Simulated Analysis with fixed scales.

Simulated Analysis with fixed scales.

There is a fair amount of variation in the above relaxation profiles. We can clearly see that the y-intercept and rate of decay are strongly dependent upon the inputs Ta and Td. Less apparent in the above plots, it seems that the relative difference between Ta and Td has a strong effect on the rate of relaxation. We investigated the relationship between the inputs and the number of steps required to fully relax Tw:

# TMax - even steps
# Tdew - even steps with some randomness
sim.data <- data.frame(tmax = seq(-20, 50, length.out=1000),
                       tdew = seq(-50, 20, length.out=1000) + runif(1000, -40, 15),
                       rate_decay = rep(0, 1000))

for (i in 1:1000) {  
    iter.tw <- c()
    calcWBGT(sim.data[i, "tmax"], sim.data[i, "tdew"])
    sim.data$rate_decay[i] <- length(iter.tw)
}
sim.data$diff <- sim.data$tmax - sim.data$tdew
sim.data$ratio <- sim.data$tmax / sim.data$tdew

pairs(sim.data, upper.panel=panel.cor, diag.panel=panel.hist)
Investigating the relationship between inputs and number of steps to fully relax Tw.

Investigating the relationship between inputs and number of steps to fully relax Tw.

The simulation supports our theory, that the rate of decay is strongly dependant on the inputs. More specifically, it seems that the rate of decay is strongly correlated to the dew point temperature and the difference between Tmax and Tdew (Tmax – Tdew). Observationally, it seems that the rate of decay increases as Tdew decreases. And (Tmax – Tdew) is strongly correlated to Tdew which explains the correlation to the rate of decay.

Interestingly the simulated data suggests that the rate of decay and (Tmax – Tdew) are not strongly correlated to Tmax. Suggesting that Tdew has a greater influence.

The ratio of Tmax to Tdew does not seems to be an important factor here.

Generalised Linear Regression

We used generalised linear regression to investigate the observations above. The rate of decay (i.e. the number of steps required to fully relax Tw) was modeled as count data assuming the following:

  • the rate of decay follows a poisson distribution
  • the observations are independent
  • the residuals are approximately normally distributed with constant variance
  • that the sample size is sufficient

Using stepwise regression an initial model (model.step) was obtained:

decay \; = \; tmax \; + \; tdew \; + \; ratio

However analysis of deviance showed that the ratio term was insignificant (p-value = 0.1358). This term was removed and the model recreated with only Tmax and Tdew (model.inputs). There was no significant difference between the two AIC values which were determined to be 9119.105 and 9119.33 for model.step and model.inputs respectively.

An interaction model (decay ~ tmax*tdew) (interaction.one) was built and compared to the additive model. Analysis of deviance confirmed that these two models were significantly different and that the interaction term was significant (p-value << 0.05). The resulting AIC of the interaction model was 7962.14.

A full additive mode (decay ~ tmax + tdew + diff + ratio) was investigated, and the diff term was found to be insignificant due to strong correlation to tdew and tmax. However, we are primarily interested in the rate of decay which the pairwise plot shown previously suggested was strongly correlated to the difference between tdew and tmax. Subsequently we built an interaction model with the diff variable in place of tmax: decay ~ tdew*diff (interaction.two). Analysis of deviance confirmed that the interaction term was significant (p-value << 0.05). We compared the AIC of the two interaction models and determined them to be 7962.14 and 7360.523 for interaction.one and interaction.two respectively. Thus, we conclude that interaction.two is more suitable in this case.

Finally, we considered the model accuracy of both interaction.one (residual deviance of 1248.9 on 996 residual degrees of freedom) and interaction.two (residual deviance of 647.26 on 996 residual degrees of freedom) and found that interaction.two provided an adequate fit to the data.

Standard residual analysis plots are shown below:

Residual vs. Fitted Values
Normal Q-Q

Normal Q-Q

The plot of Residuals vs. Fitted values shows strong signs of overdispersion and the Normal Q-Q plot indicates that the data deviates slightly from the expected normal distribution at the top end. These would give us concern for any predictions made from this model, and suggest that the confidence intervals may be wide for some sets of inputs. However, as we are not using this model for prediction then I am comfortable with making inferences as to the observed effects.

Conclusions

The iterative relaxation of Tw is a significant performance problem when applied to large data sets. However, the relaxation of Tw appears to be a smooth, continuous exponential function which we might be able to describe accurately.

Initial simulations suggest that the rate of decay (and therefore, overall performance) is srongly dependent on Tdew which is used as the initial estimate for Tw. Additionally, the difference between Tmax and Tdew is strongly correlated to the performance of the iterative algorithm. We might be able to exploit this to optimise the choice of initial value for Tw. However, this would still leave us with an iterative approach and potential problems at scale.

In future posts we will explore the relaxation more closely. Specifically, we will re-examine the empirial formula published by Stull (2011) and investigate nonlinear regression as a potential solution for modeling the decay.

References

Bernard, T.E., Pourmoghani, M. (1999). Prediction of workplace wet bulb global temperature. Applied occupational and environmental hygiene, 14(2), 126-134.

Lemke, B., & Kjellstrom, T. (2012). Calculating workplace WBGT from meteorological data: a tool for climate change assessment. Industrial health, 50(4), 267-278.

Stull, R. (2011). Wet-bulb temperature from relative humidity and air temperature. Journal of Applied Meteorology and climatology, 50(11), 2267-2269.

Advertisements

2 thoughts on “Climate Change & Heat Stress Indicators: Exploring the relaxation of Tw

  1. Pingback: Heat Stress Indicators & Climate Change: Exploring the relaxation of Wet Bulb Temperature | adatastory

  2. Pingback: WBGT & the relaxation of Tw | :: NickBurns

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s