9.3 Extreme Events

Climate extremes -- heat waves, floods, droughts, and tropical cyclones -- cause the most acute societal and ecological impacts. Extreme value theory provides the statistical framework for analyzing rare events, while attribution science quantifies how climate change alters their probability and intensity. A warming climate shifts the entire distribution of weather, making previously rare extremes more frequent and introducing unprecedented events.

Extreme Value Theory (EVT)

The Fisher-Tippett-Gnedenko theorem states that the distribution of block maxima (e.g., annual maximum daily temperature) converges to the Generalized Extreme Value (GEV) distribution regardless of the parent distribution:

$$F(x) = \exp\left\{-\left[1 + \xi\left(\frac{x - \mu}{\sigma}\right)\right]^{-1/\xi}\right\}$$

$\mu$ = location, $\sigma$ = scale, $\xi$ = shape parameter

The shape parameter $\xi$ determines the tail behavior and defines three sub-families:

Gumbel ($\xi = 0$)

Light-tailed (exponential tail decay). Common for temperature maxima in mid-latitudes. CDF simplifies to:

$F(x) = \exp\{-\exp[-(x-\mu)/\sigma]\}$

Frechet ($\xi > 0$)

Heavy-tailed (polynomial tail decay). Common for precipitation and wind speed extremes. No finite upper bound. Higher $\xi$ means heavier tails.

Weibull ($\xi < 0$)

Bounded upper tail. Finite maximum value at $\mu - \sigma/\xi$. Relevant for quantities with physical upper limits (e.g., relative humidity capped at 100%).

The return period T (in years) for an event exceeding level x is the inverse of the exceedance probability:

$$T = \frac{1}{1 - F(x)} \quad \Longleftrightarrow \quad x_T = \mu + \frac{\sigma}{\xi}\left[\left(-\ln\left(1 - \frac{1}{T}\right)\right)^{-\xi} - 1\right]$$

$x_T$ = return level: the value expected to be exceeded once every T years on average

Generalized Pareto Distribution for Threshold Exceedances

An alternative to block maxima, the peaks-over-threshold (POT) approach models exceedances above a high threshold u using the Generalized Pareto Distribution (GPD). If block maxima follow a GEV, then excesses $y = x - u$ above a sufficiently high threshold follow:

$$H(y) = 1 - \left(1 + \frac{\xi\,y}{\tilde{\sigma}}\right)^{-1/\xi} \qquad y > 0, \;\; \tilde{\sigma} > 0$$

$\tilde{\sigma} = \sigma + \xi(u - \mu)$ is the reparameterized scale; $\xi$ is the same shape as the GEV

The GPD is more data-efficient than the block maxima approach because it uses all values above the threshold rather than only annual maxima. The threshold u is chosen by examining mean residual life plots and parameter stability. For $\xi = 0$, the GPD reduces to the exponential distribution. The N-year return level is:

$$x_N = u + \frac{\tilde{\sigma}}{\xi}\left[(N\,n_y\,\zeta_u)^{\xi} - 1\right]$$

where $n_y$ = observations per year, $\zeta_u$ = probability of exceeding threshold u

Event Attribution Science

Probabilistic event attribution quantifies the role of climate change in specific extreme events. The Fraction of Attributable Risk (FAR) compares the probability of an event in the actual (human-influenced) climate versus a counterfactual (natural-only) climate:

$$\text{FAR} = 1 - \frac{P_0}{P_1}$$

$P_0$ = probability in counterfactual world, $P_1$ = probability in factual world. FAR = 1 means entirely attributable.

The probability ratio (risk ratio) RR = $P_1/P_0$ quantifies how much more likely the event has become. The conditional probability approach uses Bayes' theorem:

$$P(\text{anthropogenic} | \text{event}) = \frac{P(\text{event} | \text{anthropogenic}) \cdot P(\text{anthropogenic})}{P(\text{event})}$$

Attribution uses three complementary approaches:

Observational

Fit non-stationary GEV to observed record with time or GMST as covariate. Compare return period then vs now.

Model-based

Compare large ensembles: factual (all forcings) vs counterfactual (natural forcings only). Count exceedances in each.

Physical reasoning

Apply known physical relationships (Clausius-Clapeyron scaling, thermodynamics) to estimate the anthropogenic contribution.

Heatwaves and Heat Stress

Heatwaves are prolonged periods of anomalously high temperature. Their impacts depend on both temperature and humidity, which together determine the body's ability to cool through evaporation.

Heat Index

The "feels-like" temperature combining air temperature and relative humidity. Uses the Rothfusz regression equation (NWS). Dangerous above 40°C (104°F). The full formula includes nine terms with cross-products of T and RH.

Wet-Bulb Temperature (T_w)

The lowest temperature achievable through evaporative cooling. At T_w = 35°C, the human body cannot shed metabolic heat even at rest in shade -- a lethal threshold:

$$e(T) = e_s(T_w) - \gamma \, p \, (T - T_w)$$

$\gamma$ = psychrometric constant (~6.6 × 10⁻⁴ /°C)

Urban Heat Island Amplification

Cities are typically 2--8°C warmer than surrounding rural areas due to impervious dark surfaces, reduced vegetation, anthropogenic heat release, and altered geometry that traps longwave radiation. During heatwaves, the UHI effect intensifies because rural areas can cool via soil moisture evaporation while urban areas cannot. This disproportionately affects low-income communities with less access to air conditioning and green space.

Extreme Precipitation and Clausius-Clapeyron Scaling

A warmer atmosphere holds more water vapor according to the Clausius-Clapeyron relation. This provides a baseline expectation for how precipitation extremes change with warming:

$$\frac{de_s}{dT} = \frac{L_v \, e_s}{R_v \, T^2} \quad \Rightarrow \quad \frac{1}{e_s}\frac{de_s}{dT} \approx 7\%/\text{K}$$

Clausius-Clapeyron scaling: saturation vapor pressure increases ~7% per kelvin of warming

Some extreme events show "super-CC" scaling (exceeding 7%/K) due to dynamical feedbacks where enhanced latent heat release intensifies updrafts and moisture convergence.

Intensity-Duration-Frequency (IDF) Curves

IDF curves relate precipitation intensity to storm duration and return period, used extensively for engineering design:

$$i = \frac{a \, T^b}{(d + c)^n}$$

i = intensity (mm/hr), T = return period (yr), d = duration (hr), a,b,c,n = fitted parameters

Climate change is making existing IDF curves obsolete. Non-stationary extreme value analysis is needed to account for temporal trends in GEV parameters.

Drought Indices

Drought is a sustained moisture deficit. Key indices include:

• PDSI: Palmer Drought Severity Index -- combines precipitation, temperature, and soil moisture. Range: -10 (extreme drought) to +10.
• SPI: Standardized Precipitation Index -- precipitation fitted to gamma distribution, transformed to standard normal. SPI < -2 = extreme drought.
• SPEI: Standardized Precipitation-Evapotranspiration Index -- like SPI but includes temperature-driven evaporative demand, capturing warming-driven drought intensification.

Compound Events, Economic Losses, and Nonlinear Responses

Compound events involve multiple drivers and/or hazards that interact to produce impacts greater than any single driver alone:

Concurrent Heat + Drought

Hot conditions drive evaporation, depleting soil moisture. Dry soils reduce evaporative cooling, amplifying heat through a positive land-atmosphere feedback. Wildfire risk escalates dramatically under compound hot-dry conditions.

Storm Surge + Rainfall Flooding

Coastal areas face compound flooding when storm surge prevents rainfall drainage. Rising sea levels raise the baseline for surge while warmer atmospheres increase rainfall intensity. The joint probability may increase faster than either alone.

Economic Damage Functions and Social Cost of Carbon

Economic damages from climate change are often modeled using damage functions relating temperature change to GDP loss:

$$D(\Delta T) = 1 - \frac{1}{1 + \alpha_1 \Delta T + \alpha_2 (\Delta T)^2}$$

The social cost of carbon (SCC) integrates future damages discounted to the present. Current estimates range from $50--$200 per tonne CO₂, but uncertainty is dominated by the discount rate, climate sensitivity, and the functional form of damages at high warming levels.

Why 2°C Is Much Worse Than 1.5°C: Nonlinear Tail Risks

A small shift in the mean of a distribution causes a disproportionately large change in the probability of exceeding a high threshold. At 1.5°C warming, ~14% of the world's population is exposed to severe heatwaves at least once every 5 years; at 2°C this rises to ~37%. Coral reef survival probability drops from ~30% at 1.5°C to near 0% at 2°C. These nonlinearities arise because we are already on the steep part of many impact curves.

Fortran: Monte Carlo Simulation of Extreme Value Statistics

This Fortran program performs Monte Carlo simulation to estimate GEV return levels with bootstrap confidence intervals. It generates synthetic samples from a known GEV distribution, fits the GEV to each bootstrap replicate, and computes confidence intervals for return levels at various return periods.

program monte_carlo_gev
  ! Monte Carlo simulation of GEV return levels with bootstrap CI
  ! Generates GEV samples, fits via L-moments, bootstraps CIs
  ! Compile: gfortran -O2 -o mc_gev monte_carlo_gev.f90
  ! Run:     ./mc_gev
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: n_data = 100       ! years of data
  integer, parameter :: n_boot = 2000      ! bootstrap replicates
  integer, parameter :: n_rp = 6           ! return periods to evaluate

  real(dp) :: data(n_data)                 ! annual maxima
  real(dp) :: boot_sample(n_data)          ! bootstrap sample
  real(dp) :: return_periods(n_rp)
  real(dp) :: rl_boot(n_boot, n_rp)        ! return levels
  real(dp) :: rl_estimate(n_rp)
  real(dp) :: rl_lower(n_rp), rl_upper(n_rp)
  real(dp) :: sorted_rl(n_boot)

  ! GEV parameters
  real(dp) :: mu_true, sigma_true, xi_true
  real(dp) :: mu_fit, sigma_fit, xi_fit

  ! L-moment variables
  real(dp) :: l1, l2, l3, t3
  real(dp) :: c_approx

  ! Random number variables
  real(dp) :: u_rand, gev_quantile
  integer  :: seed_arr(8)
  integer  :: i, j, k, idx, i_lo, i_hi

  ! --- True GEV parameters ---
  mu_true    = 38.0_dp     ! location (degC)
  sigma_true = 3.0_dp      ! scale
  xi_true    = 0.10_dp     ! shape (Frechet type)

  return_periods = (/ 5.0_dp, 10.0_dp, 20.0_dp, 50.0_dp, &
                      100.0_dp, 500.0_dp /)

  ! --- Initialize random seed ---
  call random_seed()

  ! --- Generate synthetic annual maxima from GEV ---
  do i = 1, n_data
    call random_number(u_rand)
    data(i) = gev_quantile_func(u_rand, mu_true, sigma_true, xi_true)
  end do

  ! --- Fit GEV to original data using L-moments ---
  call fit_gev_lmoments(data, n_data, mu_fit, sigma_fit, xi_fit)

  write(*,'(A)') '================================================='
  write(*,'(A)') '  Monte Carlo GEV Return Level Analysis'
  write(*,'(A)') '================================================='
  write(*,'(A,F8.3,A,F8.3,A,F8.4)') &
       '  True:    mu=', mu_true, '  sigma=', sigma_true, '  xi=', xi_true
  write(*,'(A,F8.3,A,F8.3,A,F8.4)') &
       '  Fitted:  mu=', mu_fit, '  sigma=', sigma_fit, '  xi=', xi_fit
  write(*,'(A)') ''

  ! --- Compute return levels from fitted parameters ---
  do k = 1, n_rp
    u_rand = 1.0_dp - 1.0_dp / return_periods(k)
    rl_estimate(k) = gev_quantile_func(u_rand, mu_fit, sigma_fit, xi_fit)
  end do

  ! --- Bootstrap: resample with replacement ---
  write(*,'(A,I5,A)') '  Running ', n_boot, ' bootstrap replicates...'

  do j = 1, n_boot
    ! Draw bootstrap sample (resample with replacement)
    do i = 1, n_data
      call random_number(u_rand)
      idx = int(u_rand * real(n_data, dp)) + 1
      if (idx > n_data) idx = n_data
      boot_sample(i) = data(idx)
    end do

    ! Fit GEV to bootstrap sample
    call fit_gev_lmoments(boot_sample, n_data, mu_fit, sigma_fit, xi_fit)

    ! Compute return levels
    do k = 1, n_rp
      u_rand = 1.0_dp - 1.0_dp / return_periods(k)
      rl_boot(j, k) = gev_quantile_func(u_rand, mu_fit, sigma_fit, xi_fit)
    end do
  end do

  ! Re-fit to original data for point estimates
  call fit_gev_lmoments(data, n_data, mu_fit, sigma_fit, xi_fit)

  ! --- Compute 95% confidence intervals ---
  write(*,'(A)') ''
  write(*,'(A)') '  Return    Point     95% CI        95% CI       True'
  write(*,'(A)') '  Period   Estimate    Lower         Upper       Value'
  write(*,'(A)') '  -------  --------  ----------   ----------   --------'

  do k = 1, n_rp
    ! Sort bootstrap return levels for this return period
    sorted_rl = rl_boot(:, k)
    call sort_array(sorted_rl, n_boot)

    ! 2.5th and 97.5th percentiles
    i_lo = max(1, nint(0.025_dp * real(n_boot, dp)))
    i_hi = min(n_boot, nint(0.975_dp * real(n_boot, dp)))
    rl_lower(k) = sorted_rl(i_lo)
    rl_upper(k) = sorted_rl(i_hi)

    ! True return level
    u_rand = 1.0_dp - 1.0_dp / return_periods(k)
    gev_quantile = gev_quantile_func(u_rand, mu_true, sigma_true, xi_true)

    write(*,'(F8.0, 4F13.3)') return_periods(k), rl_estimate(k), &
         rl_lower(k), rl_upper(k), gev_quantile
  end do

  write(*,'(/,A)') '  Bootstrap confidence intervals bracket true values'
  write(*,'(A)')   '  when the model is correctly specified.'

contains

  function gev_quantile_func(p, mu, sigma, xi) result(x)
    ! GEV quantile function (inverse CDF)
    real(dp), intent(in) :: p, mu, sigma, xi
    real(dp) :: x, yp
    yp = -log(p)
    if (abs(xi) < 1.0e-8_dp) then
      ! Gumbel limit
      x = mu - sigma * log(yp)
    else
      x = mu + (sigma / xi) * (yp**(-xi) - 1.0_dp)
    end if
  end function

  subroutine fit_gev_lmoments(x, n, mu, sigma, xi)
    ! Fit GEV distribution using L-moments method
    integer, intent(in) :: n
    real(dp), intent(in) :: x(n)
    real(dp), intent(out) :: mu, sigma, xi
    real(dp) :: sorted(n), b0, b1, b2, l1, l2, l3, t3, c_val, gam_val
    integer :: i

    sorted = x
    call sort_array(sorted, n)

    ! Probability-weighted moments
    b0 = 0.0_dp; b1 = 0.0_dp; b2 = 0.0_dp
    do i = 1, n
      b0 = b0 + sorted(i)
      b1 = b1 + sorted(i) * real(i - 1, dp) / real(n - 1, dp)
      b2 = b2 + sorted(i) * real(i-1, dp) * real(i-2, dp) / &
           (real(n-1, dp) * real(n-2, dp))
    end do
    b0 = b0 / real(n, dp)
    b1 = b1 / real(n, dp)
    b2 = b2 / real(n, dp)

    ! L-moments
    l1 = b0
    l2 = 2.0_dp * b1 - b0
    l3 = 6.0_dp * b2 - 6.0_dp * b1 + b0
    t3 = l3 / l2  ! L-skewness

    ! Approximate GEV shape parameter from L-skewness
    c_val = 2.0_dp / (3.0_dp + t3) - log(2.0_dp) / log(3.0_dp)
    xi = 7.859_dp * c_val + 2.9554_dp * c_val**2

    ! GEV scale and location from L-moments
    gam_val = exp(log_gamma(1.0_dp + xi))
    sigma = l2 * xi / (gam_val * (1.0_dp - 2.0_dp**(-xi)))
    mu = l1 - sigma * (gam_val - 1.0_dp) / xi
  end subroutine

  function log_gamma(x) result(lg)
    ! Stirling approximation for ln(Gamma(x))
    real(dp), intent(in) :: x
    real(dp) :: lg
    lg = 0.5_dp*log(6.2831853_dp/x) + x*(log(x + 1.0_dp/(12.0_dp*x &
         - 1.0_dp/(10.0_dp*x))) - 1.0_dp)
  end function

  subroutine sort_array(arr, n)
    ! Simple insertion sort
    integer, intent(in) :: n
    real(dp), intent(inout) :: arr(n)
    real(dp) :: temp
    integer :: i, j
    do i = 2, n
      temp = arr(i); j = i - 1
      do while (j >= 1 .and. arr(j) > temp)
        arr(j+1) = arr(j); j = j - 1
      end do
      arr(j+1) = temp
    end do
  end subroutine

end program monte_carlo_gev

Interactive Simulation: Return Period Analysis (GEV)

Python

Generates synthetic annual maximum precipitation data and fits Generalized Extreme Value distributions. Computes return period curves and shows how climate change shifts the distribution, increasing 10-yr, 50-yr, and 100-yr return levels.

gev_return_period_analysis.py80 lines

import numpy as np
from scipy.stats import genextreme
import matplotlib.pyplot as plt

np.random.seed(42)

# Generate synthetic annual maximum precipitation data (mm)
n_years = 60
mu_hist, sigma_hist, xi_hist = 85.0, 15.0, 0.15
data_hist = genextreme.rvs(-xi_hist, loc=mu_hist, scale=sigma_hist, size=n_years)

# Climate-shifted distribution (+10% intensity, heavier tail)
mu_shift = mu_hist * 1.10
sigma_shift = sigma_hist * 1.15
data_future = genextreme.rvs(-xi_hist, loc=mu_shift, scale=sigma_shift, size=n_years)

# Fit GEV to both datasets
c_h, loc_h, sc_h = genextreme.fit(data_hist)
c_f, loc_f, sc_f = genextreme.fit(data_future)

# Return period curves
rp = np.logspace(0.3, 3, 300)
prob = 1.0 - 1.0 / rp
rl_hist = genextreme.ppf(prob, c_h, loc=loc_h, scale=sc_h)
rl_future = genextreme.ppf(prob, c_f, loc=loc_f, scale=sc_f)

# Specific return levels
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')
for ax in [ax1, ax2]:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    ax.xaxis.label.set_color('#cbd5e1')
    ax.yaxis.label.set_color('#cbd5e1')
    ax.title.set_color('#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)

# Return period plot
ax1.semilogx(rp, rl_hist, '#38bdf8', lw=2.5, label='Historical')
ax1.semilogx(rp, rl_future, '#f87171', lw=2.5, label='Climate-shifted')
ax1.fill_between(rp, rl_hist, rl_future, alpha=0.15, color='#f87171')
for T in [10, 50, 100]:
    p_val = 1 - 1/T
    rl_h = genextreme.ppf(p_val, c_h, loc=loc_h, scale=sc_h)
    rl_f = genextreme.ppf(p_val, c_f, loc=loc_f, scale=sc_f)
    ax1.plot(T, rl_h, 'o', color='#38bdf8', ms=8, zorder=5)
    ax1.plot(T, rl_f, 's', color='#f87171', ms=8, zorder=5)
    ax1.annotate(f'{T}-yr', xy=(T, rl_f), xytext=(T*1.3, rl_f+5),
                 color='#cbd5e1', fontsize=9)
ax1.set_xlabel('Return Period (years)')
ax1.set_ylabel('Precipitation (mm)')
ax1.set_title('Return Period Curves (GEV)')
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# PDF comparison
x = np.linspace(40, 180, 300)
pdf_h = genextreme.pdf(x, c_h, loc=loc_h, scale=sc_h)
pdf_f = genextreme.pdf(x, c_f, loc=loc_f, scale=sc_f)
ax2.plot(x, pdf_h, '#38bdf8', lw=2.5, label='Historical')
ax2.fill_between(x, pdf_h, alpha=0.2, color='#38bdf8')
ax2.plot(x, pdf_f, '#f87171', lw=2.5, label='Climate-shifted')
ax2.fill_between(x, pdf_f, alpha=0.2, color='#f87171')
ax2.set_xlabel('Annual Max Precipitation (mm)')
ax2.set_ylabel('Probability Density')
ax2.set_title('GEV Distribution Shift')
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor=fig.get_facecolor())

print("=== Return Period Analysis (GEV) ===")
print(f"Historical GEV: mu={loc_h:.1f}, sigma={sc_h:.1f}, xi={-c_h:.3f}")
print(f"Shifted GEV:    mu={loc_f:.1f}, sigma={sc_f:.1f}, xi={-c_f:.3f}")
for T in [10, 50, 100]:
    p_val = 1 - 1/T
    rl_h = genextreme.ppf(p_val, c_h, loc=loc_h, scale=sc_h)
    rl_f = genextreme.ppf(p_val, c_f, loc=loc_f, scale=sc_f)
    print(f"{T}-yr return level: Historical={rl_h:.1f} mm, Shifted={rl_f:.1f} mm (+{(rl_f/rl_h-1)*100:.1f}%)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← 9.2 Climate Change 9.4 Mitigation →