Appendix B: Formulations

1. Heat Index (HI)

The Heat Index (HI) estimates the "apparent temperature" by combining air temperature and relative humidity. We use the Rothfusz regression formula. If the temperature is provided in Celsius, it should first be converted to Fahrenheit before applying the formula, and then converted back if needed.

def heat_index(T_f, RH):
    """
    Calculate the Heat Index (HI) using the Rothfusz regression formula.
    
    Parameters:
        T_f (float): Air temperature in Fahrenheit.
        RH (float): Relative Humidity (0-100).
        
    Returns:
        float: Heat Index in Fahrenheit.
        
    Note:
        To convert from Celsius to Fahrenheit:
            T_f = (T_c * 9/5) + 32
        To convert back to Celsius:
            T_c = (HI - 32) * 5/9
    """
    HI = (-42.379 +
          2.04901523 * T_f +
          10.14333127 * RH -
          0.22475541 * T_f * RH -
          6.83783e-3 * T_f**2 -
          5.481717e-2 * RH**2 +
          1.22874e-3 * T_f**2 * RH +
          8.5282e-4 * T_f * RH**2 -
          1.99e-6 * T_f**2 * RH**2)
    return HI

# Example usage:
T_celsius = 35.0  # example air temperature in Celsius
RH = 70.0         # relative humidity in percent
T_fahrenheit = (T_celsius * 9/5) + 32  # convert to Fahrenheit
hi_value = heat_index(T_fahrenheit, RH)
print(f"Heat Index (F): {hi_value:.2f}")

2. Wet-Bulb Globe Temperature (WBGT)

The WBGT index assesses overall heat stress by combining air temperature, wet-bulb temperature (which accounts for humidity), and globe temperature (which accounts for radiant heat). The following function calculates a simplified WBGT for shaded conditions.

def wet_bulb_globe_temperature(T, T_wb, T_g):
    """
    Calculate the Wet-Bulb Globe Temperature (WBGT) for shaded environments.
    
    Parameters:
        T (float): Air temperature (e.g., in Celsius).
        T_wb (float): Wet-bulb temperature (Celsius).
        T_g (float): Globe temperature (Celsius).
        
    Returns:
        float: Approximated WBGT value.
    """
    return 0.7 * T_wb + 0.2 * T_g + 0.1 * T

# Example usage:
T = 35.0    # air temperature in Celsius
T_wb = 28.0 # wet-bulb temperature in Celsius
T_g = 30.0  # globe temperature in Celsius
wbgt_value = wet_bulb_globe_temperature(T, T_wb, T_g)
print(f"WBGT: {wbgt_value:.2f}")

3. Standardized Precipitation Evapotranspiration Index (SPEI)

The SPEI quantifies drought conditions by comparing precipitation (P) to potential evapotranspiration (PET). The basic approach involves calculating the water balance D=P−PETD = P - PETD=P−PET and then standardizing it. For simplicity, we’ll use z-score normalization as an approximation.

import numpy as np

def standardized_precipitation_evapotranspiration_index(P, PET, scale=1):
    """
    Calculate an approximation of the Standardized Precipitation Evapotranspiration Index (SPEI).
    
    The function calculates the water balance (D = P - PET) and standardizes it using z-score normalization.
    
    Parameters:
        P (array-like): Precipitation values.
        PET (array-like): Potential Evapotranspiration values.
        scale (int): Aggregation scale (default=1; for higher scales, apply a rolling window aggregation).
        
    Returns:
        numpy.ndarray: Standardized SPEI values.
    """
    P = np.array(P)
    PET = np.array(PET)
    D = P - PET  # Water balance
    D_mean = np.mean(D)
    D_std = np.std(D)
    SPEI = (D - D_mean) / D_std  # Z-score normalization
    return SPEI

# Example usage:
precipitation = [5, 0, 10, 7, 2, 0, 4]  # sample data in mm
evapotranspiration = [3, 2, 8, 6, 3, 1, 3]  # sample data in mm
spei_values = standardized_precipitation_evapotranspiration_index(precipitation, evapotranspiration)
print("SPEI values:", spei_values)

4. Atmospheric Stability Metrics: CAPE and CIN

4.1 Convective Available Potential Energy (CAPE)

CAPE quantifies the energy available for convection by integrating the buoyancy of an air parcel from the level of free convection (LFC) to the equilibrium level (EL).

def compute_CAPE(heights, T_parcel, T_env, g=9.81):
    """
    Compute Convective Available Potential Energy (CAPE) by numerical integration.
    
    CAPE is calculated as:
        CAPE = ∫ (from LFC to EL) g * ((T_parcel - T_env) / T_env) dz
    Only layers where T_parcel > T_env are considered.
    
    Parameters:
        heights (array-like): Altitudes (in meters) from LFC to EL.
        T_parcel (array-like): Temperature profile of the air parcel (in Kelvin).
        T_env (array-like): Environmental temperature profile (in Kelvin).
        g (float): Acceleration due to gravity (default is 9.81 m/s²).
        
    Returns:
        float: CAPE in J/kg.
    """
    heights = np.array(heights)
    T_parcel = np.array(T_parcel)
    T_env = np.array(T_env)
    
    # Consider only layers where the parcel is warmer than the environment
    mask = T_parcel > T_env
    integrand = g * ((T_parcel - T_env) / T_env) * mask
    CAPE = np.trapz(integrand, heights)
    return CAPE

# Example usage:
altitudes = np.linspace(1000, 3000, 21)  # Altitudes from 1000m to 3000m
T_parcel_profile = 300 - 0.005 * (altitudes - 1000)  # Hypothetical profile (K)
T_env_profile = 300 - 0.006 * (altitudes - 1000)     # Hypothetical profile (K)
cape_value = compute_CAPE(altitudes, T_parcel_profile, T_env_profile)
print(f"CAPE: {cape_value:.2f} J/kg")

4.2 Convective Inhibition (CIN)

CIN quantifies the energy barrier that must be overcome for convection to initiate. It is computed over layers where the air parcel is cooler than the environment.

def compute_CIN(heights, T_parcel, T_env, g=9.81):
    """
    Compute Convective Inhibition (CIN) by numerical integration.
    
    CIN quantifies the energy barrier for convection, calculated over layers where T_parcel < T_env.
    
    Parameters:
        heights (array-like): Altitudes (in meters) where T_parcel < T_env.
        T_parcel (array-like): Temperature profile of the air parcel (in Kelvin).
        T_env (array-like): Environmental temperature profile (in Kelvin).
        g (float): Acceleration due to gravity (default is 9.81 m/s²).
        
    Returns:
        float: CIN in J/kg.
    """
    heights = np.array(heights)
    T_parcel = np.array(T_parcel)
    T_env = np.array(T_env)
    
    # Only consider layers where the air parcel is cooler than the environment
    mask = T_parcel < T_env
    if np.sum(mask) == 0:
        return 0.0
    integrand = g * ((T_env - T_parcel) / T_env)
    CIN = np.trapz(integrand[mask], heights[mask])
    return CIN

# Example usage:
cin_value = compute_CIN(altitudes, T_parcel_profile, T_env_profile)
print(f"CIN: {cin_value:.2f} J/kg")

5. Summary

This appendix provides the mathematical formulations and corresponding Python code snippets for key derived indices as examples for heatwave prediction system:

  • Heat Index (HI): Computes the apparent temperature based on air temperature and relative humidity using the Rothfusz regression formula.

  • Wet-Bulb Globe Temperature (WBGT): Estimates overall heat stress by combining air, wet-bulb, and globe temperatures.

  • Standardized Precipitation Evapotranspiration Index (SPEI): Approximates drought conditions by standardizing the water balance (precipitation minus potential evapotranspiration).

  • Atmospheric Stability Metrics (CAPE and CIN): Quantifies the potential for convection (CAPE) and the energy barrier (CIN) using numerical integration over vertical temperature profiles.

Last updated

Was this helpful?