Module 4

Drought, Floods & Extreme Precipitation

From thermodynamic constraints on moisture to the statistics of catastrophe

4.1 Clausius-Clapeyron & the 7% Rule

The Clausius-Clapeyron equation is the single most important thermodynamic relation for understanding how precipitation extremes change with warming. It governs the temperature dependence of the saturation vapor pressure:

$ \frac{dP_{\text{sat}}}{dT} = \frac{L_v \cdot P_{\text{sat}}}{R_v \cdot T^2} $

where $L_v \approx 2.5 \times 10^6$ J/kg is the latent heat of vaporisation,$R_v = 461$ J/(kg$\cdot$K) is the gas constant for water vapour, and$T$ is temperature in Kelvin. Integrating from a reference state:

$ P_{\text{sat}}(T) = P_0 \exp\!\left[\frac{L_v}{R_v}\left(\frac{1}{T_0} - \frac{1}{T}\right)\right] $

Near 20°C, evaluating the fractional change:

$ \frac{1}{P_{\text{sat}}}\frac{dP_{\text{sat}}}{dT} = \frac{L_v}{R_v T^2} \approx \frac{2.5 \times 10^6}{461 \times 293^2} \approx 0.063 = 6.3\%/°\text{C} $

Accounting for the temperature dependence of $L_v$, the commonly quoted figure is approximately 7% per °C. This sets a fundamental thermodynamic constraint: the atmosphere can hold ~7% more water vapour for every degree of warming. Since extreme precipitation events are limited by available moisture, they are expected to intensify at roughly the same rate.

Super-Clausius-Clapeyron Scaling

Observations show that sub-daily convective precipitation extremes can scale at up to 14%/°C (roughly 2$\times$ CC), a phenomenon called super-CC scaling. The mechanism involves convective dynamics: in a warmer atmosphere, the additional latent heat release from condensation strengthens updrafts, which intensifies moisture convergence, creating a positive feedback. This is particularly relevant for flash-flood producing thunderstorms in tropical and subtropical regions. For connections to ecosystem impacts, see ourClimate & Biodiversity course.

4.2 Extreme Value Theory & Return Periods

Designing infrastructure (dams, levees, drainage) requires quantifying the probability of extreme events. The Generalized Extreme Value (GEV) distributionis the limiting distribution for block maxima (e.g., annual maximum daily rainfall):

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

where $\mu$ is the location parameter, $\sigma > 0$ is the scale, and $\xi$ is the shape parameter that determines tail behaviour:

$\xi = 0$: Gumbel (Type I) — light tail, exponential decay
$\xi > 0$: Fréchet (Type II) — heavy tail, power-law decay (common for precipitation)
$\xi < 0$: Weibull (Type III) — bounded upper tail

Return Period

The return period $T_R$ of an event with magnitude $x$ is the expected waiting time:

$ T_R = \frac{1}{1 - F(x)} $

A “100-year flood” has $F(x) = 0.99$, meaning a 1% annual exceedance probability. Inverting to find the return level:

$ x_{T_R} = \mu + \frac{\sigma}{\xi}\left[\left(-\ln\!\left(1 - \frac{1}{T_R}\right)\right)^{-\xi} - 1\right] $

IDF Curves & Climate Change

Intensity-Duration-Frequency (IDF) curves plot rainfall intensity against duration for different return periods. Under climate change, the entire IDF family shifts upward: what was a 100-year event in the historical climate becomes a 30$\text{--}$50 year event under 2°C of warming. This is because the location parameter $\mu$ increases with the Clausius-Clapeyron scaling, and the scale parameter $\sigma$ may also increase (wider distribution of extremes).

4.3 Drought Indices & Water Balance

The Water Balance Equation

At any location, the water balance is:

$ P - \text{ET} - R - \Delta S = 0 $

where $P$ is precipitation, ET is evapotranspiration, $R$ is runoff, and $\Delta S$ is the change in soil moisture storage. Drought occurs when the demand side (ET) persistently exceeds the supply side (P), depleting $S$.

Major Drought Indices

Palmer Drought Severity Index (PDSI)

Uses a two-layer soil moisture model with a water balance accounting system. PDSI ranges from $-4$ (extreme drought) to $+4$ (extreme wet). It accounts for both supply (P) and demand (potential ET from the Thornthwaite equation).

Standardized Precipitation Index (SPI)

Fits a gamma distribution to accumulated precipitation over chosen timescales (1, 3, 6, 12 months) and transforms to a standard normal. SPI < $-2$ indicates extreme drought. Advantage: multi-timescale, purely precipitation-based.

Standardized Precipitation-Evapotranspiration Index (SPEI)

Like SPI but uses $P - \text{PET}$ instead of just $P$. Critical under warming because PET increases with temperature through the vapour pressure deficit (VPD).

Vapour Pressure Deficit (VPD) Under Warming

The VPD drives evaporative demand:

$ \text{VPD} = e_s(T) - e_a = e_s(T)\,(1 - \text{RH}) $

Since $e_s(T)$ increases at ~7%/°C (Clausius-Clapeyron) and relative humidity has remained approximately constant over land, VPD increases exponentially with warming. Each 1°C of warming increases VPD by ~7%, increasing evaporative stress on vegetation, exacerbating agricultural and ecological droughts even without changes in precipitation.

4.4 Atmospheric Rivers

Atmospheric rivers (ARs) are narrow corridors of concentrated moisture transport in the lower troposphere, typically 300$\text{--}$500 km wide and 2000+ km long. They are responsible for $\sim$90% of poleward moisture transport in the midlatitudes and can deliver precipitation equivalent to major river discharge.

Integrated Vapour Transport (IVT)

The key diagnostic for ARs is the vertically integrated vapour transport:

$ \text{IVT} = \frac{1}{g}\int_{300}^{p_{\text{sfc}}} q \cdot |\mathbf{v}|\, dp $

where $q$ is specific humidity, $\mathbf{v}$ is the horizontal wind vector, and the integral runs from 300 hPa to the surface. The magnitude has units of kg/(m$\cdot$s). An AR is typically defined as a region where IVT exceeds 250 kg/(m$\cdot$s) over a length $\geq$ 2000 km.

The Atmospheric River Ranking Scale

AR Category

IVT (kg/m/s)

Impact

AR 1

250$\text{--}$500

Beneficial

AR 2

500$\text{--}$750

Mostly beneficial

AR 3

750$\text{--}$1000

Balance of beneficial/hazardous

AR 4

1000$\text{--}$1250

Mostly hazardous

AR 5

$\geq$ 1250

Exceptionally hazardous

The Pineapple Express is a well-known AR pattern that taps tropical moisture near Hawaii and delivers it to the US West Coast. Under warming, AR IVT increases at CC rates (~7%/°C), and AR frequency may also increase as the jet stream shifts.

4.5 Flood Physics & Urbanisation

Manning's Equation

Open-channel flow in rivers and drainage systems is described by Manning's equation:

$ Q = \frac{1}{n} A R_h^{2/3} S^{1/2} $

where $Q$ is discharge (m$^3$/s), $n$ is Manning's roughness coefficient, $A$ is cross-sectional flow area, $R_h = A/P_w$is the hydraulic radius (area divided by wetted perimeter), and $S$ is the energy slope. Typical $n$ values: concrete channel 0.013, natural stream 0.035, floodplain with vegetation 0.10.

Unit Hydrograph

The unit hydrograph $u(t)$ gives the runoff response to a unit pulse of effective rainfall. The total hydrograph for arbitrary rainfall $P(t)$ is obtained by convolution:

$ Q(t) = \int_0^t P(\tau)\, u(t - \tau)\, d\tau $

Urbanisation & Flood Risk

Urbanisation transforms the hydrologic response in several critical ways:

Reduced infiltration: impervious surfaces (roads, roofs, parking) increase runoff coefficient from $\sim 0.2$ (forest) to $\sim 0.9$ (urban)
Faster response time: smooth surfaces and storm drains accelerate flow, reducing the time to peak discharge
Higher peak discharge: combined effect can increase peak flows by 2$\text{--}$5$\times$
Urban heat island: warmer urban surfaces enhance convective initiation, potentially increasing local rainfall intensity

4.6 Compound Extreme Events

The most devastating impacts arise from compound events: concurrent or sequential extremes whose combined effect exceeds the sum of individual impacts. Common compound events include drought + heatwave + wildfire, coastal flooding + river flooding (from ARs hitting during high tide), and tropical cyclone + extreme precipitation + storm surge (see Module 3).

Joint Probability & Dependence

For independent events: $P(A \cap B) = P(A) \times P(B)$. But under climate change, extremes become positively dependent:

$ P(A \cap B) > P(A) \times P(B) \quad \text{(positive dependence)} $

For example, heatwaves and droughts share common drivers (persistent high-pressure blocking, soil moisture-temperature feedback), making them far more likely to co-occur than independence would suggest.

Copula Models

Copulas model the dependence structure between marginal distributions. By Sklar's theorem, any joint CDF $F(x,y)$ can be decomposed:

$ F(x, y) = C\!\left(F_X(x),\, F_Y(y)\right) $

where $C(u, v)$ is the copula function. Common choices include the Gaussian copula (symmetric tail dependence), Clayton copula (lower tail dependence, useful for drought-heat joint modelling), and Gumbel copula (upper tail dependence, useful for precipitation-flood modelling). Climate change can alter the copula parameters, increasing tail dependence and therefore compound event probability. As discussed in ourClimate & Biodiversity course, compound events are particularly devastating for ecosystems.

4.7 The Wildfire-Drought Connection

Drought and wildfire are intimately coupled. The fire weather index (FWI) system quantifies fire danger based on meteorological variables. The fine fuel moisture code (FFMC) decays exponentially during dry periods:

$ m(t) = m_{\text{eq}} + (m_0 - m_{\text{eq}})\exp(-t/\tau) $

where $m_{\text{eq}}$ is the equilibrium moisture content set by ambient conditions (RH, temperature), $m_0$ is the initial moisture, and$\tau$ is the drying timescale ($\sim$hours for fine fuels, weeks for duff). Under warming, $m_{\text{eq}}$ decreases because VPD increases, making fuels drier and extending the fire season by weeks to months.

The area burned in western North American forests scales exponentially with VPD:

$ A_{\text{burned}} \propto \exp(\beta \cdot \text{VPD}) $

with $\beta \approx 1.0\text{--}1.5$ (kPa$^{-1}$). Given$\sim$7%/°C VPD increase, each degree of warming approximately doubles the area burned in fire-prone ecosystems. Post-fire landscapes are then more vulnerable to flooding (hydrophobic soils, reduced interception), creating compound hazard cascades as discussed in Section 4.6.

IDF Curves & Atmospheric River Schematic

Intensity-Duration-Frequency Curves with Climate Change Shift

Atmospheric River Schematic

Simulation: Clausius-Clapeyron Moisture Scaling

Clausius-Clapeyron and Extreme Precipitation

Python

Saturation vapor pressure, precipitable water, and precipitation scaling with temperature

script.py111 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Clausius-Clapeyron
Lv = 2.5e6     # J/kg
Rv = 461.0     # J/(kg*K)

def e_sat(T_C):
    """Saturation vapor pressure (Pa) from Bolton (1980)"""
    return 611.2 * np.exp(17.67 * T_C / (T_C + 243.5))

def q_sat(T_C, P=101325):
    """Saturation specific humidity (kg/kg)"""
    es = e_sat(T_C)
    return 0.622 * es / (P - 0.378 * es)

T = np.linspace(-10, 45, 300)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0a0a1a')

# Panel 1: Saturation vapor pressure
ax = axes[0, 0]
ax.set_facecolor('#0a0a2e')
ax.plot(T, e_sat(T) / 100, color='#6688ff', linewidth=2.5)  # convert to hPa
ax.set_xlabel('Temperature (°C)', color='white', fontsize=11)
ax.set_ylabel('Saturation Vapor Pressure (hPa)', color='white', fontsize=11)
ax.set_title('Clausius-Clapeyron Relation', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2, color='#4444aa')

# Annotate 7% rule
for T0 in [10, 20, 30]:
    e1 = e_sat(T0) / 100
    e2 = e_sat(T0 + 1) / 100
    pct = (e2 - e1) / e1 * 100
    ax.annotate(f'{pct:.1f}%/°C', xy=(T0, e1), xytext=(T0 - 8, e1 + 15),
                arrowprops=dict(arrowstyle='->', color='#ffaa44'), color='#ffaa44', fontsize=10)

# Panel 2: Fractional rate of change
ax = axes[0, 1]
ax.set_facecolor('#0a0a2e')
frac_rate = np.gradient(e_sat(T), T[1] - T[0]) / e_sat(T) * 100
ax.plot(T, frac_rate, color='#ff6644', linewidth=2.5, label='Exact CC rate')
ax.axhline(y=7, color='#44ccaa', linewidth=1.5, linestyle='--', label='7%/°C rule of thumb')

# Super-CC region
ax.fill_between(T, 7, 14, alpha=0.15, color='#ff4444', label='Super-CC range')
ax.axhline(y=14, color='#ff4444', linewidth=1, linestyle=':', alpha=0.5)
ax.text(35, 14.5, '2x CC (convective)', color='#ff6666', fontsize=9)

ax.set_xlabel('Temperature (°C)', color='white', fontsize=11)
ax.set_ylabel('Fractional Rate (%/°C)', color='white', fontsize=11)
ax.set_title('Rate of Moisture Increase', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.set_ylim(4, 16)
ax.grid(True, alpha=0.2, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

# Panel 3: Column precipitable water
ax = axes[1, 0]
ax.set_facecolor('#0a0a2e')
# Approximate column PW: PW ~ q_sat(T_sfc) * H_scale / g
H = 2500  # scale height for moisture (m)
g = 9.81
PW = 1000 * q_sat(T) * H * 1.2  # kg/m^2 = mm
ax.plot(T, PW, color='#aa66ff', linewidth=2.5)
ax.fill_between(T, 0, PW, alpha=0.15, color='#6644aa')

# Mark current and future climate
for T0, label, color in [(15, 'Current\nmid-lat', '#44ccaa'), (20, '+2°C', '#ffaa44'), (28, 'Tropics', '#ff6644')]:
    pw0 = 1000 * q_sat(np.array([T0])) * H * 1.2
    ax.axvline(x=T0, color=color, linewidth=1, linestyle=':', alpha=0.7)
    ax.text(T0 + 0.5, pw0[0] + 3, label, color=color, fontsize=9)

ax.set_xlabel('Surface Temperature (°C)', color='white', fontsize=11)
ax.set_ylabel('Precipitable Water (mm)', color='white', fontsize=11)
ax.set_title('Column Precipitable Water', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2, color='#4444aa')

# Panel 4: Warming scenarios - extreme precip scaling
ax = axes[1, 1]
ax.set_facecolor('#0a0a2e')
dT = np.linspace(0, 4, 100)
cc_scaling = (1 + 0.07) ** dT  # 7%/°C
supercc = (1 + 0.14) ** dT     # 14%/°C (sub-daily convective)
obs_scaling = (1 + 0.10) ** dT  # ~10% observed average

ax.plot(dT, (cc_scaling - 1) * 100, color='#44ccaa', linewidth=2.5, label='CC scaling (7%/°C)')
ax.plot(dT, (obs_scaling - 1) * 100, color='#ffaa44', linewidth=2.5, linestyle='--', label='Observed (~10%/°C)')
ax.plot(dT, (supercc - 1) * 100, color='#ff4444', linewidth=2.5, label='Super-CC (14%/°C)')
ax.fill_between(dT, (cc_scaling - 1) * 100, (supercc - 1) * 100, alpha=0.1, color='#ff8844')

ax.set_xlabel('Warming (°C)', color='white', fontsize=11)
ax.set_ylabel('Precipitation Intensity Change (%)', color='white', fontsize=11)
ax.set_title('Extreme Precip Scaling with Warming', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

print("At 20°C: e_sat = {:.1f} hPa, q_sat = {:.1f} g/kg".format(e_sat(20)/100, q_sat(20)*1000))
print("At 30°C: e_sat = {:.1f} hPa, q_sat = {:.1f} g/kg".format(e_sat(30)/100, q_sat(30)*1000))
print("10°C warming increases moisture by {:.0f}%".format((e_sat(30)/e_sat(20) - 1)*100))
print("At +2°C warming: CC predicts +{:.0f}% extreme precip".format(((1.07)**2 - 1)*100))
print("At +2°C warming: super-CC predicts +{:.0f}%".format(((1.14)**2 - 1)*100))

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: GEV Return Period Analysis

Generalized Extreme Value Distribution & Return Periods

Python

Fitting GEV to precipitation extremes and computing return levels under climate change

script.py120 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# Generate synthetic annual maximum daily precipitation
# Historical: GEV with mu=60, sigma=15, xi=0.15 (Frechet-type, heavy tail)
mu_hist, sigma_hist, xi_hist = 60, 15, 0.15
n_years = 70

# GEV random samples
def gev_sample(mu, sigma, xi, n):
    u = np.random.uniform(0, 1, n)
    if abs(xi) < 1e-10:
        return mu - sigma * np.log(-np.log(u))
    return mu + (sigma / xi) * ((-np.log(u))**(-xi) - 1)

def gev_cdf(x, mu, sigma, xi):
    t = 1 + xi * (x - mu) / sigma
    t = np.maximum(t, 1e-10)
    return np.exp(-t**(-1/xi))

def gev_return_level(T_return, mu, sigma, xi):
    """Return level for return period T"""
    p = 1 - 1/T_return
    if abs(xi) < 1e-10:
        return mu - sigma * np.log(-np.log(p))
    return mu + (sigma / xi) * ((-np.log(p))**(-xi) - 1)

hist_data = gev_sample(mu_hist, sigma_hist, xi_hist, n_years)

# Future: shifted by CC scaling (+14% for +2°C)
mu_future = mu_hist * 1.14
sigma_future = sigma_hist * 1.10
xi_future = xi_hist  # shape unchanged
future_data = gev_sample(mu_future, sigma_future, xi_future, n_years)

fig, axes = plt.subplots(1, 3, figsize=(16, 5.5))
fig.patch.set_facecolor('#0a0a1a')

# Panel 1: Histograms
ax = axes[0]
ax.set_facecolor('#0a0a2e')
bins = np.linspace(20, 160, 35)
ax.hist(hist_data, bins=bins, density=True, alpha=0.5, color='#4488ff', label='Historical', edgecolor='#3366cc')
ax.hist(future_data, bins=bins, density=True, alpha=0.5, color='#ff6644', label='+2°C', edgecolor='#cc4422')

# Fitted PDFs
x = np.linspace(20, 160, 300)
pdf_hist = np.gradient(gev_cdf(x, mu_hist, sigma_hist, xi_hist), x[1]-x[0])
pdf_future = np.gradient(gev_cdf(x, mu_future, sigma_future, xi_future), x[1]-x[0])
ax.plot(x, pdf_hist, color='#6699ff', linewidth=2)
ax.plot(x, pdf_future, color='#ff8866', linewidth=2)

ax.set_xlabel('Annual Max Daily Precip (mm)', color='white', fontsize=11)
ax.set_ylabel('Probability Density', color='white', fontsize=11)
ax.set_title('GEV Distribution', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=10)

# Panel 2: Return period plot
ax = axes[1]
ax.set_facecolor('#0a0a2e')
T_range = np.logspace(0.3, 3, 200)  # 2 to 1000 years
rl_hist = [gev_return_level(T, mu_hist, sigma_hist, xi_hist) for T in T_range]
rl_future = [gev_return_level(T, mu_future, sigma_future, xi_future) for T in T_range]

ax.semilogx(T_range, rl_hist, color='#4488ff', linewidth=2.5, label='Historical')
ax.semilogx(T_range, rl_future, color='#ff6644', linewidth=2.5, label='+2°C warming')

# Show how return period changes for a fixed level
rl100 = gev_return_level(100, mu_hist, sigma_hist, xi_hist)
# Find new return period for this level under future climate
new_T = 1 / (1 - gev_cdf(rl100, mu_future, sigma_future, xi_future))
ax.axhline(y=rl100, color='#ffaa44', linewidth=1, linestyle=':', alpha=0.7)
ax.plot(100, rl100, 'o', color='#4488ff', markersize=8)
ax.plot(new_T, rl100, 'o', color='#ff6644', markersize=8)
ax.annotate(f'100-yr event\nbecomes {new_T:.0f}-yr', xy=(new_T, rl100),
            xytext=(10, rl100 + 15), color='#ffaa44', fontsize=10, fontweight='bold',
            arrowprops=dict(arrowstyle='->', color='#ffaa44'))

ax.set_xlabel('Return Period (years)', color='white', fontsize=11)
ax.set_ylabel('Return Level (mm)', color='white', fontsize=11)
ax.set_title('Return Period Analysis', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=10)

# Panel 3: Risk multiplier
ax = axes[2]
ax.set_facecolor('#0a0a2e')
T_orig = np.array([5, 10, 20, 50, 100, 200, 500])
for dT_warm, color, label in [(1, '#44ccaa', '+1°C'), (2, '#ffaa44', '+2°C'), (3, '#ff4444', '+3°C')]:
    mu_w = mu_hist * (1 + 0.07 * dT_warm)
    sigma_w = sigma_hist * (1 + 0.05 * dT_warm)
    multipliers = []
    for T in T_orig:
        rl = gev_return_level(T, mu_hist, sigma_hist, xi_hist)
        new_T_val = 1 / (1 - gev_cdf(rl, mu_w, sigma_w, xi_future))
        multipliers.append(T / new_T_val)  # frequency multiplier
    ax.semilogx(T_orig, multipliers, 'o-', color=color, linewidth=2, markersize=6, label=label)

ax.axhline(y=1, color='#666688', linewidth=1, linestyle='--')
ax.set_xlabel('Historical Return Period (years)', color='white', fontsize=11)
ax.set_ylabel('Frequency Multiplier', color='white', fontsize=11)
ax.set_title('How Often Historical Events Recur', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=10)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

print(f"Historical 100-yr return level: {rl100:.1f} mm")
print(f"Under +2°C warming, this becomes a {new_T:.0f}-yr event")
print(f"Frequency multiplier: {100/new_T:.1f}x more frequent")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: Drought Index Time Series

SPI & SPEI Drought Indices Under Warming

Python

Comparing precipitation-only (SPI) and evaporation-adjusted (SPEI) drought indices

script.py125 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(123)

# Simulate 100 years of monthly data
n_years = 100
n_months = n_years * 12
t = np.arange(n_months)
years = 1925 + t / 12

# Precipitation: seasonal cycle + trend + noise
P_mean = 80  # mm/month baseline
P_seasonal = 30 * np.sin(2 * np.pi * t / 12 - np.pi/3)  # seasonal
P_trend = -0.01 * (t - n_months/2)  # slight drying trend
P_noise = np.random.gamma(4, 5, n_months)  # gamma-distributed noise
P = np.maximum(P_mean + P_seasonal + P_trend + P_noise - 20, 0)

# Temperature: warming trend
T_mean = 15 + 10 * np.sin(2 * np.pi * t / 12)  # seasonal
T_trend = 0.02 * t / 12  # 0.02°C/year warming
T = T_mean + T_trend + np.random.normal(0, 1.5, n_months)

# PET using simplified Thornthwaite
def thornthwaite_PET(T_C, month_idx):
    T_pos = np.maximum(T_C, 0)
    I = np.sum([(np.maximum(T_pos[i::12].mean(), 0) / 5)**1.514 for i in range(12)])
    a = 6.75e-7 * I**3 - 7.71e-5 * I**2 + 1.79e-2 * I + 0.49
    PET = 16 * (10 * T_pos / I)**a
    return PET

PET = thornthwaite_PET(T, t)

# Water balance: D = P - PET
D = P - PET

# SPI: standardized precipitation (12-month accumulation)
window = 12
P_accum = np.convolve(P, np.ones(window), mode='same')
D_accum = np.convolve(D, np.ones(window), mode='same')

# Standardize
def standardize(x, ref_end=600):
    mu = x[:ref_end].mean()
    sd = x[:ref_end].std()
    return (x - mu) / sd

SPI = standardize(P_accum)
SPEI = standardize(D_accum)

# Trim edges
mask = (t >= 12) & (t < n_months - 12)

fig, axes = plt.subplots(3, 1, figsize=(14, 12))
fig.patch.set_facecolor('#0a0a1a')

# Panel 1: Temperature and precipitation
ax = axes[0]
ax.set_facecolor('#0a0a2e')
ax2 = ax.twinx()
# Annual means
yearly_P = [P[i*12:(i+1)*12].sum() for i in range(n_years)]
yearly_T = [T[i*12:(i+1)*12].mean() for i in range(n_years)]
year_labels = np.arange(1925, 1925 + n_years)

ax.bar(year_labels, yearly_P, color='#4488ff', alpha=0.6, width=0.8, label='Annual P')
ax2.plot(year_labels, yearly_T, color='#ff6644', linewidth=2, label='Annual T')
ax.set_ylabel('Annual Precipitation (mm)', color='#4488ff', fontsize=11)
ax2.set_ylabel('Mean Temperature (°C)', color='#ff6644', fontsize=11)
ax.set_title('Annual Precipitation & Temperature', color='white', fontsize=13, fontweight='bold')
ax.tick_params(axis='y', colors='#4488ff')
ax2.tick_params(axis='y', colors='#ff6644')
ax.tick_params(axis='x', colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')

# Panel 2: SPI
ax = axes[1]
ax.set_facecolor('#0a0a2e')
yr = years[mask]
spi = SPI[mask]
ax.fill_between(yr, 0, spi, where=spi >= 0, color='#4488ff', alpha=0.6)
ax.fill_between(yr, 0, spi, where=spi < 0, color='#ff4444', alpha=0.6)
ax.axhline(y=0, color='#666688', linewidth=0.5)
ax.axhline(y=-2, color='#ff6644', linewidth=1, linestyle='--', alpha=0.5)
ax.text(1927, -2.2, 'Extreme drought', color='#ff6644', fontsize=9)
ax.set_ylabel('SPI-12', color='white', fontsize=11)
ax.set_title('Standardized Precipitation Index (12-month)', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.set_ylim(-3.5, 3.5)

# Panel 3: SPEI
ax = axes[2]
ax.set_facecolor('#0a0a2e')
spei = SPEI[mask]
ax.fill_between(yr, 0, spei, where=spei >= 0, color='#44cc88', alpha=0.6)
ax.fill_between(yr, 0, spei, where=spei < 0, color='#ff4444', alpha=0.6)
ax.axhline(y=0, color='#666688', linewidth=0.5)
ax.axhline(y=-2, color='#ff6644', linewidth=1, linestyle='--', alpha=0.5)

# 10-year running mean to show trend
kernel = np.ones(120) / 120
spei_smooth = np.convolve(SPEI, kernel, mode='same')
ax.plot(years[mask], spei_smooth[mask], color='#ffffff', linewidth=2, label='10-yr mean')

ax.set_xlabel('Year', color='white', fontsize=11)
ax.set_ylabel('SPEI-12', color='white', fontsize=11)
ax.set_title('SPEI Shows Drying Trend from Increased Evaporative Demand', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.set_ylim(-3.5, 3.5)
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=10)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

# Drought statistics
spi_drought = np.sum(SPI[mask] < -1.5) / np.sum(mask) * 100
spei_drought = np.sum(SPEI[mask] < -1.5) / np.sum(mask) * 100
print(f"Fraction of time in severe drought (SPI < -1.5): {spi_drought:.1f}%")
print(f"Fraction of time in severe drought (SPEI < -1.5): {spei_drought:.1f}%")
print(f"SPEI captures {spei_drought/spi_drought:.1f}x more drought due to warming-driven ET increase")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: Atmospheric River IVT

Integrated Vapour Transport Calculation

Python

Computing IVT for an idealized atmospheric river profile and its response to warming

script.py124 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

g = 9.81  # m/s^2

def compute_IVT(T_sfc, RH_profile, wind_profile, p_levels):
    """Compute IVT from vertical profiles.
    T_sfc: surface temperature (°C)
    RH_profile: relative humidity at each level
    wind_profile: wind speed at each level (m/s)
    p_levels: pressure levels (Pa)
    """
    # Temperature profile (6.5°C/km lapse rate, scale height 8km)
    T_profile = T_sfc + 273.15 - 6.5 * (-8000 * np.log(p_levels / p_levels[0]))
    T_C = T_profile - 273.15
    
    # Saturation specific humidity
    e_sat = 611.2 * np.exp(17.67 * T_C / (T_C + 243.5))
    q_sat = 0.622 * e_sat / (p_levels - 0.378 * e_sat)
    q = RH_profile * q_sat
    
    # IVT = (1/g) * integral(q * |v| * dp)
    integrand = q * wind_profile
    IVT = np.trapezoid(integrand, p_levels) / g
    
    return abs(IVT), q, T_profile

# Pressure levels (surface to 300 hPa)
p_levels = np.linspace(101325, 30000, 50)  # Pa

# Typical AR profile
RH_ar = 0.95 * np.exp(-((p_levels - 85000)**2) / (2 * 20000**2))  # peaked in lower tropo
RH_ar = np.clip(RH_ar + 0.3, 0, 1)
wind_ar = 25 * np.exp(-((p_levels - 85000)**2) / (2 * 25000**2)) + 10  # low-level jet

# Non-AR profile
RH_nonar = np.full_like(p_levels, 0.5)
wind_nonar = np.full_like(p_levels, 8.0)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0a0a1a')

# Panel 1: Vertical profiles
ax = axes[0, 0]
ax.set_facecolor('#0a0a2e')
IVT_ar, q_ar, T_ar = compute_IVT(18, RH_ar, wind_ar, p_levels)
IVT_nonar, q_nonar, T_nonar = compute_IVT(18, RH_nonar, wind_nonar, p_levels)

ax.plot(q_ar * 1000, p_levels / 100, color='#66aaff', linewidth=2.5, label=f'AR (IVT={IVT_ar:.0f})')
ax.plot(q_nonar * 1000, p_levels / 100, color='#888888', linewidth=2, linestyle='--', label=f'Non-AR (IVT={IVT_nonar:.0f})')
ax.set_xlabel('Specific Humidity (g/kg)', color='white', fontsize=11)
ax.set_ylabel('Pressure (hPa)', color='white', fontsize=11)
ax.set_title('Moisture Profile', color='white', fontsize=13, fontweight='bold')
ax.invert_yaxis()
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=10)

# Panel 2: Wind profile
ax = axes[0, 1]
ax.set_facecolor('#0a0a2e')
ax.plot(wind_ar, p_levels / 100, color='#ff8844', linewidth=2.5, label='AR (low-level jet)')
ax.plot(wind_nonar, p_levels / 100, color='#888888', linewidth=2, linestyle='--', label='Non-AR')
ax.set_xlabel('Wind Speed (m/s)', color='white', fontsize=11)
ax.set_ylabel('Pressure (hPa)', color='white', fontsize=11)
ax.set_title('Wind Profile', color='white', fontsize=13, fontweight='bold')
ax.invert_yaxis()
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=10)

# Panel 3: IVT vs surface temperature
ax = axes[1, 0]
ax.set_facecolor('#0a0a2e')
T_range = np.linspace(5, 35, 60)
IVT_vs_T = [compute_IVT(T, RH_ar, wind_ar, p_levels)[0] for T in T_range]
IVT_vs_T_nonar = [compute_IVT(T, RH_nonar, wind_nonar, p_levels)[0] for T in T_range]

ax.plot(T_range, IVT_vs_T, color='#66aaff', linewidth=2.5, label='AR conditions')
ax.plot(T_range, IVT_vs_T_nonar, color='#888888', linewidth=2, linestyle='--', label='Non-AR')
ax.axhline(y=250, color='#ffaa44', linewidth=1.5, linestyle=':', label='AR threshold (250)')

# AR categories
for ivt_thresh, cat, color in [(500, 'AR2', '#44cc88'), (750, 'AR3', '#cccc44'), (1000, 'AR4', '#ff8844'), (1250, 'AR5', '#ff4444')]:
    ax.axhline(y=ivt_thresh, color=color, linewidth=0.8, linestyle=':', alpha=0.5)
    ax.text(6, ivt_thresh + 20, cat, color=color, fontsize=9)

ax.set_xlabel('Surface Temperature (°C)', color='white', fontsize=11)
ax.set_ylabel('IVT (kg/m/s)', color='white', fontsize=11)
ax.set_title('IVT Scales with CC', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

# Panel 4: IVT cross-section of an AR
ax = axes[1, 1]
ax.set_facecolor('#0a0a2e')
y_cross = np.linspace(-500, 500, 200)  # km from AR centre
IVT_cross = 800 * np.exp(-y_cross**2 / (2 * 120**2))  # Gaussian cross-section
IVT_warm = 800 * 1.14 * np.exp(-y_cross**2 / (2 * 130**2))  # warmer: wider + stronger

ax.fill_between(y_cross, 0, IVT_cross, alpha=0.3, color='#4488ff')
ax.plot(y_cross, IVT_cross, color='#4488ff', linewidth=2.5, label='Current climate')
ax.fill_between(y_cross, 0, IVT_warm, alpha=0.2, color='#ff6644')
ax.plot(y_cross, IVT_warm, color='#ff6644', linewidth=2.5, label='+2°C warming')
ax.axhline(y=250, color='#ffaa44', linewidth=1.5, linestyle=':', label='AR threshold')

ax.set_xlabel('Cross-AR Distance (km)', color='white', fontsize=11)
ax.set_ylabel('IVT (kg/m/s)', color='white', fontsize=11)
ax.set_title('AR Cross-Section', color='white', fontsize=13, fontweight='bold')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.15, color='#4444aa')
ax.legend(facecolor='#1a1a3e', edgecolor='#4444aa', labelcolor='white', fontsize=9)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')

print(f"AR IVT at 18°C: {IVT_ar:.0f} kg/m/s")
IVT_28 = compute_IVT(28, RH_ar, wind_ar, p_levels)[0]
print(f"AR IVT at 28°C: {IVT_28:.0f} kg/m/s")
print(f"IVT increase: {(IVT_28/IVT_ar - 1)*100:.0f}%")
print(f"Effective CC scaling: {(IVT_28/IVT_ar - 1)*100/10:.1f}%/°C")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4.7 Precipitation Recycling & Moisture Feedback

Not all precipitation originates from oceanic moisture advection. A significant fraction isrecycled: evapotranspired from the land surface and re-precipitated locally or downstream. The recycling ratio $\rho$ is defined as:

$ \rho = \frac{P_{\text{recycled}}}{P_{\text{total}}} = \frac{E \cdot L}{E \cdot L + F_{\text{in}}} $

where $E$ is the local evapotranspiration rate, $L$ is the characteristic length scale of the domain, and $F_{\text{in}}$ is the moisture influx from external sources. In the Amazon basin, recycling ratios reach$\sim 25\text{--}35\%$, meaning over a quarter of rainfall originates from the forest itself. Deforestation breaks this cycle, reducing downwind rainfall.

The soil moisture-precipitation feedback can be positive or negative depending on the evaporative regime. In energy-limited regions (humid climates), soil moisture has little effect on ET. In moisture-limited regions (semi-arid), soil moisture strongly controls ET and therefore local convective triggering. The Budyko framework classifies catchments by the aridity index $\Phi = \text{PET}/P$:

$ \frac{\text{ET}}{P} = \frac{\Phi \tanh(1/\Phi)(1 - \exp(-\Phi))}{\Phi \tanh(1/\Phi) + (1 - \exp(-\Phi))} $

Under warming, the transition from energy-limited to moisture-limited regimes shifts poleward, expanding the area where land-atmosphere coupling amplifies drought. This has direct implications for the drought indices discussed in Section 4.3 and for the compound events in Section 4.6.

4.8 Flash Drought & Rapid Onset Events

Flash droughts develop rapidly (weeks rather than months), driven by anomalous atmospheric demand rather than prolonged precipitation deficits. They are characterised by:

Rapid soil moisture depletion: persistent high VPD and above-normal temperatures drive anomalous ET, depleting root-zone moisture within 2$\text{--}$4 weeks
Heat wave amplification: as soil dries, sensible heat flux increases (latent flux decreases), raising surface temperature further — a positive feedback loop
Crop failure: onset during critical growth stages (flowering, grain fill) causes disproportionate yield loss
Forecasting difficulty: standard drought indices (SPI, PDSI) respond too slowly; specialised metrics using ET anomalies are needed

The 2012 US Great Plains flash drought caused $\sim$$30 billion in agricultural losses. Flash drought frequency is projected to increase under warming because the VPD increases exponentially with temperature, accelerating the ET-driven soil moisture drawdown. The rate of soil moisture decline during a flash drought can be modelled as:

$ \frac{d\theta}{dt} = \frac{P - \text{ET}(\theta, \text{VPD}) - R(\theta)}{z_r} $

where $\theta$ is volumetric soil moisture, $z_r$ is root zone depth, and ET depends on both soil moisture availability and atmospheric demand (VPD). When VPD is anomalously high and precipitation ceases, $d\theta/dt$ becomes strongly negative, with exponential decay toward the wilting point.

References

Clausius, R. (1850). Über die bewegende Kraft der Wärme. Annalen der Physik, 79, 368-397, 500-524.

Trenberth, K. E. (2011). Changes in precipitation with climate change. Climate Research, 47(1-2), 123-138.

Lenderink, G. & van Meijgaard, E. (2008). Increase in hourly precipitation extremes beyond expectations from temperature changes. Nature Geoscience, 1, 511-514.

Coles, S. (2001). An Introduction to Statistical Modeling of Extreme Values. Springer.

Palmer, W. C. (1965). Meteorological Drought. U.S. Weather Bureau Research Paper No. 45.

Vicente-Serrano, S. M., Beguería, S., & López-Moreno, J. I. (2010). A multiscalar drought index sensitive to global warming: The Standardized Precipitation Evapotranspiration Index. Journal of Climate, 23(7), 1696-1718.

Ralph, F. M., et al. (2019). A scale to characterize the strength and impacts of atmospheric rivers. Bulletin of the American Meteorological Society, 100(2), 269-289.

Zscheischler, J., et al. (2018). Future climate risk from compound events. Nature Climate Change, 8, 469-477.