Humid Heat Metrics

Humid Heat Metrics#

Overview#

Humid heat metrics index temperature and humidity, which is often more useful than comparing either variable alone when considering perceived temperature, as they together affect the human body’s ability to cool itself down. Humid heat metrics are especially important for the safety of outside workers, the elderly, or otherwise high-risk, high-exposure individuals[1].

Wet Bulb Globe Temperature (\(WBGT\)) is a measure of heat stress. The equations for outdoor (\(WBGT_{od}\)) and indoor/shaded (\(WBGT_{id}\)) WBGTs are[1]:

\(WBGT_{od} = 0.7*T_{nwb} + 0.2*T_g + 0.1*T_a\)

\(WBGT_{id} = 0.7*T_{nwb} + 0.3*T_g\)

where \(T_a\) refers to Dry Bulb Ambient Temperatue, \(T_{nwb}\) is the Natural Wet Bulb Temperature with exposure to wind and sun, and \(T_g\) is the Globe Temperature taken from inside a copper globe painted black and exposed to the sun[1].

However, this formula is complicated by the reality that Natural Wet Bulb Temperature and Globe Temperature are not always readily available variables from weather stations or atmospheric models.

In this notebook, we will demonstrate the Australian Bureau of Meteorology (ABM) and Bernard methods of predicting wet bulb global temperature with a focus on the July 1995 Chicago heatwave.

For our analysis, we have ERA5 reanalysis data for the lower contiguous United States (50°N, 24°S, -66°E, -125°W) from July 1995 with the variables: 2-meter temperature, 2-meter dew point temperature, surface pressure, and u/v wind components.

ERA5 is a reanalysis, a global weather/climate dataset that combines model output and observations with physical understanding to create a spatially and temporally consistent historic dataset, spanning from 1940 to today (updated every 5 days)[2].

import geocat.datafiles as gdf
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from scipy import optimize

Downloading file 'registry.txt' from 'https://github.com/NCAR/GeoCAT-datafiles/raw/main/registry.txt' to '/home/runner/.cache/geocat'.

era5 = xr.open_dataset(gdf.get("netcdf_files/era5_1995-07-14T12.nc"))
era5.u10

Downloading file 'netcdf_files/era5_1995-07-14T12.nc' from 'https://github.com/NCAR/GeoCAT-datafiles/raw/main/netcdf_files/era5_1995-07-14T12.nc' to '/home/runner/.cache/geocat'.

<xarray.DataArray 'u10' (latitude: 105, longitude: 237)> Size: 199kB
[24885 values with dtype=float64]
Coordinates:
  * longitude  (longitude) float32 948B -125.0 -124.8 -124.5 ... -66.25 -66.0
  * latitude   (latitude) float32 420B 50.0 49.75 49.5 49.25 ... 24.5 24.25 24.0
    time       datetime64[ns] 8B ...
Attributes:
    units:      m s**-1
    long_name:  10 metre U wind component
# Convert from Kelvin to Celsius

era5['t2m_C'] = era5['t2m'] - 273.15
era5['d2m_C'] = era5['d2m'] - 273.15

# Chicago coordinates
lat_chicago = 41.8781
lon_chicago = -87.6298

Australian Bureau of Meteorology (ABM)#

The Australian Bureau of Meteorology’s method of estimating Wet Bulb Global Temperature is attractive due to its simplicity. It only requires temperature and relative humidity[3].

Below is a chart of WBGT from relative humidity and temperature[3]:

This method tends to overpredict WBGT compared to other models and assumes full sunlight and light breeze.

def calc_abm_wbgt(t_a, rh):
    p = (
        (rh / 100) * 6.105 * np.exp(17.27 * t_a / (237.7 + t_a))
    )  # water vapor pressure [hPa]
    wbgt = (0.567 * t_a) + (0.393 * p) + 3.94
    return wbgt

To use our ERA5 data in this equation, we need to first calculate relative humidity (the ratio of vapor pressure to saturation pressure) from temperature and dewpoint. To do this we use the Magnus-Tetens Approximation for vapor pressure[4]:

\(e = 6.11 \exp {\left( \frac{17.625 \times t}{t + 243.04} \right)}\)

where \(e\) is vapor pressure and \(t\) is temperature in Kelvin.

def _calc_vapor_pressure(t):  # Magnus-Tetens Approximation
    e = 6.11 * np.exp((17.27 * t) / (t + 237.3))  # Vapor Pressure in hPa
    return e


def calc_relative_humidity_era5(t_a, t_d):
    e = _calc_vapor_pressure(t_d)  # vapor pressure from dew point temp
    e_sat = _calc_vapor_pressure(t_a)  # saturation vapor pressure
    rh = 100 * e / e_sat  # Clausius-Clapeyron equation
    return rh


rh = calc_relative_humidity_era5(era5.t2m_C, era5.d2m_C)

wbgt_abm = calc_abm_wbgt(era5.t2m_C, rh)
wbgt_abm

<xarray.DataArray (latitude: 105, longitude: 237)> Size: 199kB
array([[17.87026107, 16.17965608, 14.57478675, ..., 19.19009504,
        18.33482943, 18.17073414],
       [18.5690197 , 17.13604585, 17.36420689, ..., 18.53771728,
        18.47195479, 18.54166586],
       [17.96797569, 18.64702521, 18.87962489, ..., 18.72776049,
        18.70938808, 18.84343994],
       ...,
       [21.67176102, 21.69203306, 21.68911165, ..., 30.19523038,
        29.97646297, 29.82487547],
       [21.93360605, 21.91177247, 21.90737553, ..., 30.34311432,
        30.14710387, 30.0325734 ],
       [22.19253977, 22.14410015, 22.12947582, ..., 30.43313854,
        30.35434565, 30.2302265 ]])
Coordinates:
  * longitude  (longitude) float32 948B -125.0 -124.8 -124.5 ... -66.25 -66.0
  * latitude   (latitude) float32 420B 50.0 49.75 49.5 49.25 ... 24.5 24.25 24.0
    time       datetime64[ns] 8B 1995-07-14T12:00:00

Plotting Chicago ABM WBGT#

fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())

c = plt.contourf(era5.longitude, era5.latitude, wbgt_abm, cmap='inferno')

ax.coastlines()

cbar = plt.colorbar(c, ax=ax, orientation='horizontal')
cbar.set_label('Wet Bulb Global Temperature' + '\N{DEGREE SIGN}' + 'C')

ax.set_title('July 14, 1995 noon - ABM WBGT')

# Annotate location of Chicago
ax.plot(lon_chicago, lat_chicago, 'k*');

/home/runner/micromamba/envs/geocat-applications/lib/python3.12/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
../../_images/3d4430b9bc073183298435050ecd064b357cf9b108702098c311aa6d9256fc4d.png

Bernard#

Bernard’s semi-empirical formula approximates Natural Wet Bulb temperature based on heat exchange of a wetted wick exposed to sun and wind based on measurements of common United States summertime environmental conditions[5].

This is considered an indoor WBGT temperature because it does not include any strong radiative sources in the calculation.

\( WBGT_{id} = \begin{cases} 0.7T_{pwb} + 0.3T_a & \text{if } v > 3 m/s\\ 0.67T_{pwb} + 0.33T_a − 0.048 log_10v (T_a − T_{pwb}) & \text{if } 0.3 \geq v \leq 3 m/s \end{cases} \)

In Bernard’s analysis, wind speeds less than 0.3 m/s are not included since the field of humid heat metrics is primarily concerned with workers, and an outdoor worker is unlikely to be stationary. Apparent wind speeds are assumed to be at least 1 m/s[5].

This formula utilizes thermodynamic Wet Bulb Temperature (\(T_{pwb}\)), which is a wet bulb temperature in the shade and fanned or rotated. This is the wet bulb typically used for dew point calculations, and can be iteratively derived from temperature (\(T_a\)) and dewpoint (\(T_d\)).

# Bernard formula for WBGT
def calc_bernard_wbgt(t_a, t_pwb, v):
    if np.all(v < 0.3):  # m/s
        return np.nan  # Return NaN where velocity is below the threshold
    elif np.all((0.3 <= v) & (v <= 3)):
        wbgt = (0.67 * t_pwb) + (0.33 * t_a) - (0.48 * np.log10(v) * (t_a - t_pwb))
    else:
        wbgt = (0.7 * t_pwb) + (0.3 * t_a)

    return wbgt

$T_{pwb}$ is iteratively solved from McPherson’s formula[6]:

\(1556 e_d - 1.484 e_d * T_{pwb} - 1556 e_w + 1.484 * e_w * T_{pwb} + 1010 * (t_a - t_pwb) = 0\)

where \(e_d = 6.106 * exp(17.27 * T_d / (237.3 + T_d))\)
and \(e_w = 6.106 * exp(17.27 * T_{pwb} / (237.3 + T_{pwb}))\)

Here we use a Newton-Raphson iterative method for the iterative solve for \(t_{pwb}\).

\(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\)

Essentially, Newton-Raphson is a root finding method that plugs your initial guess into the equation in question and the derivative of that equation in order to get a more accurate guess [7]. This is repeated until your new guess is sufficiently close. Thankfully we have scipy.optimize.newton() to handle this solve for us.

# Scalar function to compute t_pwb
def _calc_tpwb_scalar(t_a, t_d):
    def f(t_pwb):
        e_d = 6.106 * np.exp(17.27 * t_d / (237.3 + t_d))  # hPa
        # e_w = 6.106 * np.exp(17.27 * t_pwb / (237.3 + t_pwb))
        func = (
            1556 * e_d
            - 1.484 * e_d * t_pwb
            - 1556 * (6.106 * np.exp(17.27 * t_pwb / (237.3 + t_pwb)))
            + 1.484 * (6.106 * np.exp(17.27 * t_pwb / (237.3 + t_pwb))) * t_pwb
            + 1010 * (t_a - t_pwb)
        )
        return func

    def f_prime(t_pwb, h=1e-5):  # numerical derivative
        return (f(t_pwb + h) - f(t_pwb - h)) / (2 * h)

    # Use the Newton-Raphson method with scipy's newton function
    t_pwb_0 = t_d  # initial guess
    t_pwb = optimize.newton(f, t_pwb_0, fprime=f_prime, tol=1e-6, maxiter=100)
    return t_pwb


# Apply function over grid
def _calc_tpwb(t_a, t_d):
    return xr.apply_ufunc(
        _calc_tpwb_scalar,
        t_a,
        t_d,
        vectorize=True,
        dask="parallelized",
        output_dtypes=[float],
    )

v = np.sqrt(era5.u10**2 + era5.v10**2)  # combine u and v wind components
t_pwb = _calc_tpwb(era5.t2m_C, era5.d2m_C)

wbgt_bernard = calc_bernard_wbgt(era5.t2m_C, t_pwb, v)
wbgt_bernard

<xarray.DataArray (latitude: 105, longitude: 237)> Size: 199kB
array([[13.69517637, 11.8926209 , 10.1376864 , ..., 14.99941759,
        14.13418851, 13.97024568],
       [14.41539187, 12.91478894, 13.15184086, ..., 14.35260232,
        14.28806213, 14.3599888 ],
       [13.78971097, 14.48718552, 14.72253851, ..., 14.55732155,
        14.54146397, 14.67783079],
       ...,
       [17.48845595, 17.50781776, 17.50468557, ..., 24.95087155,
        24.79824126, 24.69033199],
       [17.73964687, 17.71867202, 17.71418707, ..., 25.05317386,
        24.91487184, 24.82648568],
       [17.98647294, 17.94030902, 17.92631711, ..., 25.11072483,
        25.04516283, 24.94564003]])
Coordinates:
  * longitude  (longitude) float32 948B -125.0 -124.8 -124.5 ... -66.25 -66.0
  * latitude   (latitude) float32 420B 50.0 49.75 49.5 49.25 ... 24.5 24.25 24.0
    time       datetime64[ns] 8B 1995-07-14T12:00:00

Plotting Chicago Bernard WBGT#

fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())

c = plt.contourf(era5.longitude, era5.latitude, wbgt_bernard, cmap='inferno')

ax.coastlines()

cbar = plt.colorbar(c, ax=ax, orientation='horizontal')
cbar.set_label('Wet Bulb Global Temperature' + '\N{DEGREE SIGN}' + 'C')

ax.set_title('July 14, 1995 noon - Bernard WBGT method')

# Annotate location of Chicago
ax.plot(lon_chicago, lat_chicago, 'k*');
../../_images/c2d8a1d76f8f1d040caf7afe6394d2e35d1d43932349fe0debd4b65f85cf9b23.png

Comparing methods#

When comparing our output from both ABM and Bernard, ABM tends to estimate at higher WBGT by 4 - 7 degrees Celsius.

fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()

diff = wbgt_abm - wbgt_bernard

c = plt.contourf(era5.longitude, era5.latitude, diff, cmap='Reds')

cbar = plt.colorbar(c, ax=ax, orientation='horizontal')
cbar.set_label(
    '\N{GREEK CAPITAL LETTER DELTA} Wet Bulb Global Temperature'
    + '\N{DEGREE SIGN}'
    + 'C'
)
ax.set_title('July 14, 1995 noon - WBGT Difference (ABM - Bernard)')

# Annotate location of Chicago
ax.plot(lon_chicago, lat_chicago, 'k*');
../../_images/b09998b8b57a21ce3a3b46cf2207407f07f61ea39d9c305f5b1dd04409c62c74.png

References:#

Contents