Part III: Climate & Oceans | Chapter 2

Atmospheric Physics

Radiative transfer, lapse rates, Hadley cells, and the Coriolis effect

3.1 Radiative Transfer in the Atmosphere

Radiation passing through the atmosphere is absorbed and emitted by greenhouse gases. The fundamental equation governing this process is the Schwarzschild equation of radiative transfer.

Derivation 1: The Schwarzschild Equation

Consider monochromatic radiation of specific intensity $I_\nu$ passing through a slab of atmosphere with absorption coefficient $k_\nu$ and density $\rho$. The change in intensity along path length $ds$ includes absorption and thermal emission:

$$\frac{dI_\nu}{ds} = -k_\nu \rho I_\nu + k_\nu \rho B_\nu(T)$$

where $B_\nu(T)$ is the Planck function (Kirchhoff's law: emissivity equals absorptivity at LTE). Defining the optical depth $d\tau_\nu = k_\nu \rho \, ds$:

$$\frac{dI_\nu}{d\tau_\nu} = -I_\nu + B_\nu(T)$$

The formal solution looking upward through the atmosphere from the surface is:

$$I_\nu(0) = B_\nu(T_s) e^{-\tau_\nu^*} + \int_0^{\tau_\nu^*} B_\nu(T(\tau')) e^{-\tau'} d\tau'$$

where $\tau_\nu^*$ is the total optical depth. The first term is attenuated surface emission; the second is atmospheric emission. The outgoing longwave radiation (OLR) is obtained by integrating over all frequencies and solid angles.

Historical Context

Karl Schwarzschild derived the radiative transfer equation in 1906 for stellar atmospheres. Its application to planetary atmospheres was developed by Chandrasekhar (1950) and Goody (1964). Modern climate models use correlated-k methods or line-by-line calculations with databases like HITRAN containing millions of spectral lines.

3.2 Atmospheric Temperature Structure

Derivation 2: The Dry Adiabatic Lapse Rate

When an air parcel rises adiabatically, it expands and cools. From the first law of thermodynamics for an ideal gas with no heat exchange:

$$C_P dT = \frac{dP}{\rho}$$

Combining with hydrostatic balance $dP = -\rho g \, dz$:

$$\Gamma_d = -\frac{dT}{dz} = \frac{g}{C_P} = \frac{9.81}{1004} = 9.8 \text{ K/km}$$

This is the dry adiabatic lapse rate (DALR). The actual lapse rate in the troposphere (~6.5 K/km) is less than the DALR because of latent heat release from condensation.

Derivation 3: The Moist Adiabatic Lapse Rate

When saturated air rises, condensation releases latent heat $L_v$, reducing the cooling rate. The saturated adiabatic lapse rate is:

$$\Gamma_s = \Gamma_d \frac{1 + \frac{L_v w_s}{R_d T}}{1 + \frac{L_v^2 w_s}{C_P R_v T^2}}$$

where $w_s$ is the saturation mixing ratio, $R_d = 287$ J/(kg$\cdot$K) is the gas constant for dry air, $R_v = 461$ J/(kg$\cdot$K) for water vapor, and$L_v \approx 2.5 \times 10^6$ J/kg. The saturated mixing ratio follows the Clausius-Clapeyron equation:

$$\frac{de_s}{dT} = \frac{L_v e_s}{R_v T^2}$$

which gives approximately 7% increase in saturation vapor pressure per kelvin of warming. At typical tropical surface conditions (T = 300 K), $\Gamma_s \approx 4$-5 K/km, while at cold polar conditions the moist rate approaches the dry rate.

Atmospheric Stability

The atmosphere is statically stable when the environmental lapse rate $\Gamma_e$is less than the adiabatic rate: $\Gamma_e < \Gamma_d$. Convective instability occurs when $\Gamma_e > \Gamma_d$. The Brunt-Vaisala frequency:

$$N^2 = \frac{g}{\theta}\frac{d\theta}{dz} = \frac{g}{T}(\Gamma_d - \Gamma_e)$$

where $\theta = T(P_0/P)^{R_d/C_P}$ is potential temperature. When $N^2 > 0$, displaced parcels oscillate; when $N^2 < 0$, they accelerate away (convection).

3.3 The General Circulation

Derivation 4: The Coriolis Effect

In a rotating reference frame, a moving object experiences an apparent deflection called the Coriolis acceleration. For motion with velocity $\mathbf{v}$ on a rotating Earth with angular velocity $\boldsymbol{\Omega}$:

$$\mathbf{a}_{\text{Cor}} = -2\boldsymbol{\Omega} \times \mathbf{v}$$

The vertical component of Earth's rotation at latitude $\phi$ gives the Coriolis parameter:

$$f = 2\Omega \sin\phi$$

where $\Omega = 7.292 \times 10^{-5}$ rad/s. Horizontal motion is deflected to the right in the Northern Hemisphere and to the left in the Southern Hemisphere. At the equator,$f = 0$ and the Coriolis effect vanishes.

Hadley Cells and the Jet Stream

George Hadley (1735) proposed that solar heating at the equator drives rising motion, poleward flow aloft, sinking at the poles, and return flow at the surface. However, the Coriolis effect limits the meridional extent of direct thermal circulation.

The angular momentum constraint limits the Hadley cell. Air moving poleward conserves angular momentum $M = (u + \Omega a \cos\phi)a\cos\phi$, where $a$ is Earth's radius. Starting from rest at the equator, the zonal wind at latitude $\phi$ would be:

$$u(\phi) = \Omega a \frac{\cos^2\phi_0 - \cos^2\phi}{\cos\phi}$$

This gives unrealistically large winds at high latitudes, indicating that the Hadley cell cannot extend to the pole. The Held-Hou (1980) model predicts the Hadley cell extent by balancing angular momentum and thermal constraints.

Derivation 5: Geostrophic Balance and Thermal Wind

At synoptic scales, the dominant balance in the atmosphere is between pressure gradient and Coriolis forces (geostrophic balance):

$$f u = -\frac{1}{\rho}\frac{\partial P}{\partial y}, \quad f v = \frac{1}{\rho}\frac{\partial P}{\partial x}$$

The thermal wind relation connects vertical wind shear to horizontal temperature gradients:

$$f \frac{\partial u_g}{\partial z} = -\frac{g}{T}\frac{\partial T}{\partial y}$$

This explains the subtropical jet stream: the strong equator-to-pole temperature gradient in the upper troposphere drives fast westerly winds, peaking at the tropopause at ~200 hPa (~12 km altitude) with speeds of 30-70 m/s. The jet stream is strongest in winter when the temperature gradient is largest.

3.4 Atmospheric Windows and Spectroscopy

The atmosphere is opaque at most infrared wavelengths due to absorption by H$_2$O, CO$_2$, O$_3$, CH$_4$, and N$_2$O. The key atmospheric window is 8-12 $\mu$m, through which surface radiation escapes to space. CO$_2$ has a strong absorption band centered at 15 $\mu$m (the $\nu_2$ bending mode), H$_2$O absorbs across the far-infrared, and O$_3$ absorbs at 9.6 $\mu$m.

The logarithmic dependence of forcing on CO$_2$ concentration arises because the center of the 15 $\mu$m band is already saturated. Additional CO$_2$ mainly widens the absorption band into the wings, where the increase in optical depth grows logarithmically with concentration. This was first understood by Arrhenius (1896) and formalized by Manabe and Wetherald (1967).

3.5 Rossby Waves and Midlatitude Weather

Carl-Gustaf Rossby (1939) showed that the variation of the Coriolis parameter with latitude ($\beta = df/dy = 2\Omega\cos\phi/a$) supports a class of planetary-scale waves. The dispersion relation for barotropic Rossby waves is:

$$\omega = -\frac{\beta k}{k^2 + l^2}$$

where $k$ and $l$ are zonal and meridional wavenumbers. Rossby waves always propagate westward relative to the mean flow. The phase speed is:

$$c_x = \bar{u} - \frac{\beta}{k^2 + l^2}$$

Stationary Rossby waves (forced by topography and land-sea contrasts) occur when$c_x = 0$, giving the stationary wavenumber$K_s = \sqrt{\beta/\bar{u}}$. For typical midlatitude winds ($\bar{u} \approx 15$ m/s), this gives wavenumber ~4-6, corresponding to the familiar pattern of ridges and troughs in the jet stream that controls weather patterns.

Baroclinic Instability

The midlatitude cyclones (low-pressure systems) that dominate weather in temperate regions arise from baroclinic instability of the jet stream. Jule Charney (1947) and Eric Eady (1949) independently showed that horizontal temperature gradients in a rotating fluid are unstable to wave-like perturbations. The most unstable wavelength in the Eady model is:

$$\lambda_{\max} = \frac{2\pi \times 0.38 f_0 H}{N}$$

where $H$ is the depth of the fluid, $f_0$ is the Coriolis parameter, and $N$ is the Brunt-Vaisala frequency. For typical midlatitude values,$\lambda_{\max} \approx 4000$ km, matching the observed scale of extratropical cyclones. The growth rate is proportional to the temperature gradient, explaining why storms are strongest in winter and at the polar front.

The Eady growth rate provides a fundamental timescale for midlatitude weather:

$$\sigma_{\max} = 0.31 \frac{f_0}{N}\frac{\partial U}{\partial z} \approx 0.31 \frac{f_0}{N}\frac{U}{H}$$

giving e-folding times of ~2-3 days, consistent with the typical development time of midlatitude cyclones from initial perturbation to mature storm.

3.6 Tropical Circulation and Monsoons

The tropical atmosphere is governed by different dynamics than midlatitudes. The Coriolis parameter is small, so geostrophic balance breaks down and the gradient wind or nonlinear balance becomes important. The Walker circulation (east-west overturning along the equator) interacts with the Hadley cell, and its variability drives the El Nino-Southern Oscillation (ENSO).

The Held-Hou Model of Hadley Cell Width

Held and Hou (1980) derived the poleward extent of the Hadley cell $\phi_H$ by requiring that the angular momentum-conserving wind equals the thermal wind at the cell boundary:

$$\phi_H = \left(\frac{5}{3}\frac{\Delta_h g H}{2\Omega^2 a^2}\right)^{1/2}$$

where $\Delta_h$ is the fractional equator-to-pole temperature difference,$H$ is the tropopause height, and $a$ is Earth's radius. This gives$\phi_H \approx 30°$, consistent with the observed Hadley cell extent. The model also predicts that the Hadley cell widens as the planet warms (decreasing temperature gradient) or rotates more slowly.

The Bjerknes Feedback and ENSO

Jacob Bjerknes (1969) identified the positive feedback between SST gradients and atmospheric circulation that drives ENSO. In the normal state, cold SSTs in the eastern Pacific drive easterly trade winds, which maintain upwelling and cold SSTs (a positive feedback). During El Nino, weakened trades allow warm water to spread eastward, further weakening the trades. The system oscillates with a period of 2-7 years due to delayed negative feedback from oceanic wave dynamics (Kelvin and Rossby waves propagating across the Pacific basin).

Computational Lab: Atmospheric Physics

Atmospheric Structure, Radiative Transfer, and General Circulation

Python

script.py331 lines

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

# ============================================================
# ATMOSPHERIC PHYSICS SIMULATION
# ============================================================
print("ATMOSPHERIC PHYSICS SIMULATION")
print("=" * 70)

# ============================================================
# 1. Atmospheric Temperature Profile
# ============================================================
print()
print("1. Standard Atmosphere Temperature Profile")
print("-" * 60)

g = 9.81
Cp = 1004.0
Rd = 287.0
Rv = 461.0
Lv = 2.5e6
P0 = 101325.0
T0 = 288.15

DALR = g / Cp * 1000  # K/km
print(f"  Dry adiabatic lapse rate: {DALR:.2f} K/km")
print()

print(f"  {'Alt (km)':12} {'T (K)':10} {'P (hPa)':12} {'rho (kg/m3)':14} {'Layer'}")
print("  " + "-" * 60)

std_atm = [
    (0, 288.15, 1013.25, "Troposphere"),
    (2, 275.15, 795, "Troposphere"),
    (4, 262.15, 616, "Troposphere"),
    (6, 249.15, 472, "Troposphere"),
    (8, 236.15, 356, "Troposphere"),
    (10, 223.15, 265, "Troposphere"),
    (11, 216.65, 227, "Tropopause"),
    (15, 216.65, 121, "Stratosphere"),
    (20, 216.65, 55.3, "Stratosphere"),
    (25, 221.65, 25.5, "Stratosphere"),
    (30, 226.65, 12.0, "Stratosphere"),
    (40, 250.35, 2.87, "Stratosphere"),
    (50, 270.65, 0.80, "Stratopause"),
    (60, 247.0, 0.22, "Mesosphere"),
    (80, 198.6, 0.011, "Mesopause"),
]

for alt, T, P, layer in std_atm:
    rho = P * 100 / (Rd * T)
    print(f"  {alt:10} {T:8.1f} {P:10.2f} {rho:12.4f} {layer}")

# ============================================================
# 2. Moist Adiabatic Lapse Rate
# ============================================================
print()
print()
print("2. Moist Adiabatic Lapse Rate vs Temperature")
print("-" * 55)

print(f"  {'T (K)':10} {'T (C)':10} {'es (hPa)':12} {'ws (g/kg)':12} {'MALR (K/km)'}")
print("  " + "-" * 50)

for T in [240, 250, 260, 270, 275, 280, 285, 290, 295, 300, 305, 310]:
    # Clausius-Clapeyron (Bolton 1980)
    T_C = T - 273.15
    es = 6.112 * np.exp(17.67 * T_C / (T_C + 243.5))  # hPa
    ws = 0.622 * es / (1013.25 - es)  # kg/kg

numerator = 1 + Lv * ws / (Rd * T)
    denominator = 1 + Lv**2 * ws / (Cp * Rv * T**2)
    MALR = DALR * numerator / denominator

print(f"  {T:8} {T_C:8.1f} {es:10.2f} {ws*1000:10.2f} {MALR:10.2f}")

# ============================================================
# 3. Coriolis Parameter and Geostrophic Wind
# ============================================================
print()
print()
print("3. Coriolis Parameter and Rossby Number")
print("-" * 60)

Omega = 7.292e-5
a_earth = 6.371e6

print(f"  {'Latitude':12} {'f (s^-1)':14} {'Inertial T (hr)':18} {'Beta':16} {'Rossby #'}")
print("  " + "-" * 65)

v_wind = 10.0  # m/s
L_scale = 1000e3  # 1000 km

for lat in [0, 5, 10, 15, 20, 30, 45, 60, 75, 90]:
    f = 2 * Omega * np.sin(np.radians(lat))
    beta = 2 * Omega * np.cos(np.radians(lat)) / a_earth
    T_inert = 2 * np.pi / f / 3600 if f > 0 else float('inf')
    Ro = v_wind / (f * L_scale) if f > 0 else float('inf')
    f_str = f"{f:.4e}" if f > 0 else "0"
    T_str = f"{T_inert:.1f}" if T_inert < 1000 else "inf"
    Ro_str = f"{Ro:.3f}" if Ro < 100 else "inf"
    print(f"  {lat:10} {f_str:12} {T_str:16} {beta:14.2e} {Ro_str}")

# ============================================================
# 4. Thermal Wind: Jet Stream Profile
# ============================================================
print()
print()
print("4. Thermal Wind: Zonal Wind Profile")
print("-" * 55)

# Simplified temperature gradient: dT/dy ~ -1 K per 100 km near 45N
dTdy = -1.0 / 100e3  # K/m
f_45 = 2 * Omega * np.sin(np.radians(45))
T_mean = 250.0  # mean temperature (K)

print(f"  Latitude = 45N, dT/dy = {dTdy*1000:.1f} K/km, f = {f_45:.4e} s^-1")
print()
print(f"  {'Altitude (km)':16} {'P (hPa)':12} {'u_g (m/s)':14} {'Notes'}")
print("  " + "-" * 48)

# Integrate thermal wind from surface
u_g = 0.0  # surface wind
dz = 1000.0  # 1 km steps

for z_km in range(0, 18):
    z = z_km * 1000
    dug_dz = -(g / (f_45 * T_mean)) * dTdy
    P_approx = 1013 * np.exp(-z / 8500)
    note = ""
    if z_km == 0:
        note = "Surface"
    elif z_km == 11:
        note = "Tropopause"
    elif z_km == 12:
        note = "Jet core"
    print(f"  {z_km:14} {P_approx:10.0f} {u_g:12.1f} {note}")

if z_km < 11:
        u_g += dug_dz * dz
    elif z_km >= 11:
        u_g -= 1.0  # weakening above tropopause

# ============================================================
# 5. Planck Function and Atmospheric Windows
# ============================================================
print()
print()
print("5. Planck Function at Key Temperatures")
print("-" * 65)

h = 6.626e-34
c = 3e8
kB = 1.381e-23

def planck_wavelength(lam_um, T):
    lam = lam_um * 1e-6
    return 2*h*c**2/lam**5 / (np.exp(h*c/(lam*kB*T)) - 1) * 1e-6  # W/(m^2.um)

print(f"  {'Lambda (um)':14}", end="")
for T in [220, 255, 288, 300]:
    print(f"{'B(' + str(T) + 'K)':16}", end="")
print()
print("  " + "-" * 65)

for lam in [3, 5, 8, 9.6, 10, 12, 15, 20, 30, 50]:
    print(f"  {lam:12}", end="")
    for T in [220, 255, 288, 300]:
        B = planck_wavelength(lam, T)
        print(f"  {B:14.4f}", end="")
    print(" W/(m^2.um)")

# ============================================================
# PLOTS
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# Plot 1: Temperature profile with stability
ax = axes[0, 0]
z_arr = np.array([a[0] for a in std_atm])
T_arr = np.array([a[1] for a in std_atm])
ax.plot(T_arr, z_arr, 'b-o', linewidth=2, markersize=4, label='Standard atmosphere')

# DALR from surface
z_dalr = np.linspace(0, 15, 50)
T_dalr = T0 - DALR/1000 * z_dalr * 1000
ax.plot(T_dalr, z_dalr, 'r--', linewidth=1.5, label='Dry adiabat')

# MALR from surface
T_malr = [T0]
z_malr = [0]
T_curr = T0
for dz_step in range(1, 16):
    T_C = T_curr - 273.15
    es = 6.112 * np.exp(17.67 * T_C / (T_C + 243.5))
    ws = 0.622 * es / (1013.25 - es)
    num = 1 + Lv * ws / (Rd * T_curr)
    den = 1 + Lv**2 * ws / (Cp * Rv * T_curr**2)
    malr = DALR * num / den / 1000
    T_curr -= malr * 1000
    T_malr.append(T_curr)
    z_malr.append(dz_step)

ax.plot(T_malr, z_malr, 'g--', linewidth=1.5, label='Moist adiabat')

ax.axhline(y=11, color='orange', linestyle=':', alpha=0.5, label='Tropopause')
ax.axhline(y=50, color='purple', linestyle=':', alpha=0.5, label='Stratopause')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Altitude (km)')
ax.set_title('Atmospheric Temperature Profile')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 85)

# Plot 2: Planck functions and atmospheric absorption
ax = axes[0, 1]
lam_range = np.linspace(3, 60, 500)

for T, color, label in [(288, 'red', '288K (surface)'),
                          (255, 'blue', '255K (effective)'),
                          (220, 'green', '220K (tropopause)')]:
    B_arr = [planck_wavelength(l, T) for l in lam_range]
    ax.plot(lam_range, B_arr, color=color, linewidth=2, label=label)

# Mark absorption bands
bands = [(13, 17, 'CO2 15um'), (5, 8, 'H2O'), (9, 10, 'O3 9.6um')]
for l1, l2, name in bands:
    ax.axvspan(l1, l2, alpha=0.15, color='gray')
    ax.text((l1+l2)/2, max([planck_wavelength(10, 288)])*0.9, name,
            ha='center', fontsize=7, color='gray')

ax.fill_between([8, 12], 0, max([planck_wavelength(10, 288)])*0.5,
                alpha=0.1, color='yellow', label='Atmospheric window')
ax.set_xlabel('Wavelength (um)')
ax.set_ylabel('Spectral Radiance (W/m^2/um)')
ax.set_title('Planck Functions and Atmospheric Bands')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# Plot 3: Hadley cell - angular momentum conserving wind
ax = axes[1, 0]
lats = np.linspace(0.1, 60, 200)
u_AM = Omega * a_earth * (np.cos(np.radians(0))**2 - np.cos(np.radians(lats))**2) / np.cos(np.radians(lats))

ax.plot(lats, u_AM, 'b-', linewidth=2, label='Angular momentum conserving')
ax.axhline(y=0, color='gray', linewidth=0.5)

# Observed approximate jet stream speed
lats_obs = [0, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60]
u_obs = [0, 5, 10, 15, 25, 35, 30, 25, 20, 15, 12, 10]
ax.plot(lats_obs, u_obs, 'r-o', linewidth=2, markersize=4, label='Observed (approx)')

ax.axvline(x=30, color='green', linestyle=':', alpha=0.5, label='Hadley cell edge')
ax.set_xlabel('Latitude (degrees)')
ax.set_ylabel('Zonal Wind Speed (m/s)')
ax.set_title('Upper Tropospheric Winds')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 60)
ax.set_ylim(0, 200)

# Plot 4: Rossby wave dispersion
ax = axes[1, 1]
beta = 2 * Omega * np.cos(np.radians(45)) / a_earth
k_range = np.linspace(0.1e-6, 5e-6, 200)
u_mean_vals = [5, 10, 15, 20, 30]

for u_mean in u_mean_vals:
    c_phase = u_mean - beta / k_range**2
    ax.plot(k_range * 1e6, c_phase, linewidth=2, label=f'U = {u_mean} m/s')

ax.axhline(y=0, color='red', linewidth=2, linestyle='--', label='Stationary (c=0)')
ax.set_xlabel('Zonal wavenumber k (x10^-6 m^-1)')
ax.set_ylabel('Phase speed (m/s)')
ax.set_title('Rossby Wave Dispersion at 45N')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_ylim(-30, 40)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Held-Hou Hadley Cell Width")
print("-" * 55)

delta_h = 1.0/3.0  # fractional equator-pole T difference
H_trop = 15e3     # tropopause height (m)
a_HH = 6.371e6

print(f"  Parameters: delta_h = {delta_h:.2f}, H = {H_trop/1000:.0f} km")
print()
print(f"  {'Omega (Earth=1)':18} {'phi_H (deg)':14} {'u_jet (m/s)':14} {'Notes'}")
print("  " + "-" * 55)

for omega_frac in [0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0]:
    Om = Omega * omega_frac
    phi_H = np.degrees(np.sqrt(5/3 * delta_h * g * H_trop / (2 * Om**2 * a_HH**2)))
    phi_H = min(phi_H, 90)
    # Jet speed at cell boundary
    u_jet = Om * a_HH * (np.cos(np.radians(0))**2 - np.cos(np.radians(phi_H))**2) / np.cos(np.radians(phi_H))
    note = "Earth" if abs(omega_frac - 1.0) < 0.01 else ""
    print(f"  {omega_frac:16.2f} {phi_H:12.1f} {u_jet:12.0f} {note}")

# ENSO indices
print()
print()
print("7. ENSO State Classification")
print("-" * 55)
print(f"  {'Nino3.4 SST anomaly':24} {'Classification':18} {'Typical effects'}")
print("  " + "-" * 65)

enso_states = [
    ("< -1.5 C", "Strong La Nina", "Wet Australia, dry S. America"),
    ("-1.5 to -0.5 C", "La Nina", "Enhanced trade winds"),
    ("-0.5 to +0.5 C", "Neutral", "Normal conditions"),
    ("+0.5 to +1.5 C", "El Nino", "Wet Peru, dry Indonesia"),
    ("> +1.5 C", "Strong El Nino", "Global T spike, coral bleaching"),
    ("> +2.5 C", "Super El Nino", "1997-98, 2015-16 type"),
]

for anom, cls, effects in enso_states:
    print(f"  {anom:22} {cls:16} {effects}")

print()
print()
print("[Chart saved: T profile, Planck functions, Hadley winds, Rossby dispersion]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Atmospheric Stability and Baroclinic Instability Analysis

Python

script.py74 lines

import numpy as np

print("ATMOSPHERIC STABILITY AND DYNAMICS")
print("=" * 70)

g = 9.81
Omega = 7.292e-5
a = 6.371e6
Cp = 1004.0
Rd = 287.0

# Brunt-Vaisala frequency
print()
print("1. Atmospheric Stability: Brunt-Vaisala Frequency")
print("-" * 60)

DALR = g / Cp * 1000  # K/km

print(f"  {'Env. lapse (K/km)':20} {'N^2 (s^-2)':14} {'N (s^-1)':12} {'Period (min)':14} {'Stability'}")
print("  " + "-" * 65)

T_ref = 250.0  # K

for gamma in [4, 5, 6, 6.5, 7, 8, 9, 9.8, 10, 12]:
    N2 = g / T_ref * (DALR - gamma) / 1000 * 1000
    if N2 > 0:
        N = np.sqrt(N2)
        period = 2 * np.pi / N / 60
        stab = "Stable"
    elif N2 == 0:
        N = 0
        period = float('inf')
        stab = "Neutral"
    else:
        N = np.sqrt(-N2)
        period = 0
        stab = "UNSTABLE"
    N_str = f"{N:.4f}" if N2 > 0 else ("0" if N2 == 0 else f"i*{N:.4f}")
    per_str = f"{period:.1f}" if period < 1000 else "inf"
    print(f"  {gamma:18.1f} {N2:12.2e} {N_str:10} {per_str:12} {stab}")

# Baroclinic instability
print()
print()
print("2. Eady Baroclinic Instability")
print("-" * 55)

H_trop = 10000  # m
N_bv = 0.01     # s^-1 (typical troposphere)

print(f"  H = {H_trop/1000:.0f} km, N = {N_bv} s^-1")
print()
print(f"  {'Latitude':12} {'f0':12} {'U (m/s)':12} {'lambda (km)':14} {'Growth rate (1/day)':20} {'e-fold (days)'}")
print("  " + "-" * 75)

for lat in [30, 35, 40, 45, 50, 55, 60]:
    f0 = 2 * Omega * np.sin(np.radians(lat))
    # Typical thermal wind: U ~ 20 m/s at 45N
    U = 20 * np.sin(np.radians(lat)) / np.sin(np.radians(45))
    # Eady most unstable wavelength
    lam = 2 * np.pi * 0.38 * f0 * H_trop / N_bv / 1000  # km
    # Eady growth rate
    sigma = 0.31 * f0 * U / (N_bv * H_trop)
    sigma_day = sigma * 86400
    e_fold = 1.0 / sigma / 86400
    print(f"  {lat:10} {f0:10.4e} {U:10.1f} {lam:12.0f} {sigma_day:18.3f} {e_fold:10.1f}")

print()
print("  Key insight: most unstable wavelength ~4000 km matches")
print("  observed extratropical cyclone scale.")

print()
print("Complete: atmospheric dynamics analysis finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Climate Science Next: Ocean Circulation →