Module 0

Radiative Transfer & Energy Balance

From stellar radiation to greenhouse warming: the physics that determines planetary temperature

0.1 Blackbody Radiation & the Planck Function

All matter with temperature above absolute zero emits electromagnetic radiation. A blackbody is an idealized object that absorbs all incident radiation and emits with a spectrum determined solely by its temperature. Max Planck derived the spectral radiance in 1900, launching quantum theory:

\( B_\nu(T) = \frac{2h\nu^3}{c^2} \cdot \frac{1}{e^{h\nu/(k_BT)} - 1} \)

Planck function (spectral radiance per unit frequency)

In terms of wavelength \(\lambda\), using \(\nu = c/\lambda\) and\(|d\nu| = c/\lambda^2 |d\lambda|\):

\( B_\lambda(T) = \frac{2hc^2}{\lambda^5} \cdot \frac{1}{e^{hc/(\lambda k_BT)} - 1} \)

Wien's Displacement Law

The wavelength of maximum emission is found by setting \(\partial B_\lambda / \partial \lambda = 0\). Defining \(x = hc/(\lambda k_B T)\), the condition becomes:

\( (5 - x) e^x = 5 \)

\( x \approx 4.965 \implies \lambda_{\max} T = \frac{hc}{4.965\, k_B} \approx 2898\;\mu\text{m}\cdot\text{K} \)

For the Sun (\(T \approx 5778\) K): \(\lambda_{\max} \approx 0.50\;\mu\text{m}\) (green-yellow visible). For Earth (\(T \approx 255\) K): \(\lambda_{\max} \approx 11.4\;\mu\text{m}\) (thermal infrared).

Emissivity and Kirchhoff's Law

Real objects are not perfect blackbodies. The emissivity \(\varepsilon_\lambda\) describes the ratio of actual emission to blackbody emission. Kirchhoff's law states that at thermal equilibrium, absorptivity equals emissivity at each wavelength:

\( \alpha_\lambda = \varepsilon_\lambda \)

Kirchhoff's law of thermal radiation

Integrating the Planck Function: Stefan-Boltzmann Derivation

The total power radiated per unit area is obtained by integrating \(B_\lambda\) over all wavelengths and the hemisphere of solid angle:

\( F = \int_0^\infty \pi B_\lambda(T)\,d\lambda = \int_0^\infty \frac{2\pi hc^2}{\lambda^5(e^{hc/\lambda k_BT}-1)}\,d\lambda \)

\( \text{Let } x = \frac{hc}{\lambda k_BT} \implies d\lambda = -\frac{hc}{k_BT x^2}\,dx \)

\( F = \frac{2\pi k_B^4 T^4}{h^3 c^2}\int_0^\infty \frac{x^3}{e^x-1}\,dx = \frac{2\pi k_B^4 T^4}{h^3 c^2}\cdot\frac{\pi^4}{15} = \sigma T^4 \)

The integral \(\int_0^\infty x^3/(e^x-1)\,dx = \pi^4/15\) is evaluated using the Riemann zeta function: \(\Gamma(4)\zeta(4) = 6 \cdot \pi^4/90 = \pi^4/15\).

Atmospheric Windows

Not all infrared wavelengths are absorbed by the atmosphere. The atmospheric window (8-12 \(\mu\)m) allows surface radiation to escape directly to space, except for the 9.6 \(\mu\)m ozone band. Key absorption bands:

H\(_2\)O: 6.3 \(\mu\)m (bending), rotation band \(>20\;\mu\)m
CO\(_2\): 15 \(\mu\)m (bending), 4.3 \(\mu\)m (asymmetric stretch)
O\(_3\): 9.6 \(\mu\)m (within the window)
CH\(_4\): 7.7 \(\mu\)m
N\(_2\)O: 7.8, 17 \(\mu\)m

The closing of the atmospheric window by increasing greenhouse gas concentrations is the primary mechanism of enhanced greenhouse warming. Each molecule's radiative forcing depends on where its absorption overlaps with the window and with other absorbers.

0.2 Stefan-Boltzmann Law & Effective Temperature

Integrating the Planck function over all wavelengths and solid angles yields the total power radiated per unit area:

\( F = \sigma T^4 \quad \text{where} \quad \sigma = \frac{2\pi^5 k_B^4}{15 c^2 h^3} \approx 5.67 \times 10^{-8}\;\text{W m}^{-2}\text{K}^{-4} \)

Derivation of Effective Temperature

Earth intercepts solar radiation over its cross-section \(\pi R_E^2\) and reflects a fraction\(\alpha\) (planetary albedo \(\approx 0.30\)). The absorbed power is:

\( P_{\text{abs}} = S \cdot \pi R_E^2 \cdot (1 - \alpha) \)

where \(S \approx 1361\;\text{W/m}^2\) is the solar constant. Earth re-radiates as a sphere of area \(4\pi R_E^2\):

\( P_{\text{emit}} = \sigma T_{\text{eff}}^4 \cdot 4\pi R_E^2 \)

Setting \(P_{\text{abs}} = P_{\text{emit}}\) and solving:

\( S \pi R_E^2 (1-\alpha) = \sigma T_{\text{eff}}^4 \cdot 4\pi R_E^2 \)

\( \boxed{T_{\text{eff}} = \left[\frac{S(1-\alpha)}{4\sigma}\right]^{1/4} \approx 255\;\text{K} \approx -18°\text{C}} \)

The observed mean surface temperature is about \(288\) K (\(+15\)°C), implying a \(33\) K greenhouse warming.

One-Layer Greenhouse Model

Consider a single atmospheric layer that is transparent to shortwave solar radiation but absorbs all outgoing longwave (infrared) radiation. Let \(T_s\) = surface temperature and\(T_a\) = atmospheric layer temperature.

Top-of-atmosphere balance: The atmosphere must radiate\(\sigma T_a^4\) upward to space, matching the absorbed solar flux:

\( \sigma T_a^4 = \frac{S(1-\alpha)}{4} = \sigma T_{\text{eff}}^4 \implies T_a = T_{\text{eff}} \)

Atmospheric layer balance: The layer absorbs \(\sigma T_s^4\) from below and emits \(\sigma T_a^4\) both upward and downward:

\( \sigma T_s^4 = 2\sigma T_a^4 \)

Surface temperature:

\( \boxed{T_s = 2^{1/4} T_{\text{eff}} \approx 1.189 \times 255 \approx 303\;\text{K}} \)

The one-layer model predicts \(\sim 303\) K, overshooting by \(\sim 15\) K because real greenhouse absorption is not 100%. Multi-layer and partial-absorption models refine this further.

Surface Energy Balance

At the surface, the full energy balance includes both radiative and non-radiative terms:

\( S_{\text{abs}} + F_{\text{atm}\downarrow} = \sigma T_s^4 + F_{\text{SH}} + F_{\text{LH}} + G \)

\(G\) = ground heat flux, \(F_{\text{SH}}\) = sensible, \(F_{\text{LH}}\) = latent heat

The Bowen ratio \(\beta_B = F_{\text{SH}}/F_{\text{LH}}\) characterizes the surface energy partitioning: over oceans \(\beta_B \sim 0.1\) (evaporation dominates), over deserts \(\beta_B > 5\) (sensible heat dominates). This ratio plays a crucial role in determining local climate and is sensitive to soil moisture and vegetation changes.

The latent heat flux is parameterized using bulk aerodynamic formulas:

\( F_{\text{LH}} = \rho L_v C_E |\mathbf{u}|(q_s - q_a) \)

\(L_v \approx 2.5 \times 10^6\) J/kg, \(C_E \approx 10^{-3}\), \(q_s - q_a\) = specific humidity difference

0.3 The Schwarzschild Equation of Radiative Transfer

The fundamental equation governing how radiation propagates through an absorbing and emitting medium:

\( \frac{dI_\nu}{ds} = -\kappa_\nu \left( I_\nu - B_\nu(T) \right) \)

Schwarzschild equation for monochromatic intensity

where \(I_\nu\) is the specific intensity, \(\kappa_\nu\) is the absorption coefficient (m\(^{-1}\)), \(B_\nu(T)\) is the Planck function at local temperature, and \(s\) is the path length.

Optical Depth

Define the optical depth along the path:

\( d\tau_\nu = \kappa_\nu \, ds = \kappa_\nu \rho \, ds \quad \implies \quad \tau_\nu = \int_0^s \kappa_\nu \rho \, ds' \)

In terms of optical depth, the Schwarzschild equation becomes:

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

Beer-Lambert Law

For a non-emitting medium (\(B_\nu = 0\) or cold background), the equation\(dI_\nu/d\tau = -I_\nu\) integrates to:

\( I_\nu(\tau) = I_\nu(0) \, e^{-\tau_\nu} \)

\( \mathcal{T}_\nu = e^{-\tau_\nu} \)

Beer-Lambert transmission

The general formal solution with emission included:

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

The first term represents attenuated background radiation; the second is the contribution of emission along the path, each element attenuated by the remaining optical depth.

0.4 CO\(_2\) Absorption & Radiative Forcing

Carbon dioxide is a linear triatomic molecule (O=C=O) with three vibrational modes. The bending mode at \(\sim 15\;\mu\text{m}\) (\(667\;\text{cm}^{-1}\)) is the most important for climate because it overlaps with Earth's peak thermal emission:

Symmetric stretch (\(\nu_1\)): \(\sim 7.2\;\mu\text{m}\) — infrared inactive (no dipole change)
Bending mode (\(\nu_2\)): \(\sim 15\;\mu\text{m}\) — doubly degenerate, strong absorber
Asymmetric stretch (\(\nu_3\)): \(\sim 4.3\;\mu\text{m}\) — active but less climatically important

Line-by-Line vs Band Models

Line-by-line (LBL) calculations compute absorption for each spectral line using Voigt profiles (Lorentzian pressure broadening + Doppler broadening). Each line has a strength \(S_j\) and half-width \(\gamma_j\):

\( \kappa(\nu) = \sum_j S_j \cdot f(\nu - \nu_j; \gamma_j) \)

Band models simplify by averaging over spectral intervals. The Malkmus random-band model gives the mean transmission:

\( \overline{\mathcal{T}} = \exp\!\left[-\frac{\pi \alpha}{\delta}\left(\sqrt{1 + \frac{S u}{\pi \alpha}} - 1\right)\right] \)

Radiative Forcing from CO\(_2\)

Because the 15 \(\mu\)m band centre is already saturated, additional CO\(_2\) primarily broadens the absorption wings. This leads to a logarithmic dependence of radiative forcing on concentration:

\( \Delta F = 5.35 \cdot \ln\!\left(\frac{C}{C_0}\right) \quad \text{[W/m}^2\text{]} \)

Myhre et al. (1998) — validated against LBL calculations

Derivation sketch: The absorption cross-section in the Lorentzian wings falls as \(\kappa \propto (\nu - \nu_0)^{-2}\). The optical depth reaches unity (new absorption becomes effective) at a wing distance proportional to \(\sqrt{C}\). The integral of the Planck function over this incrementally absorbed bandwidth gives:

\( \Delta \nu_{\text{abs}} \propto \sqrt{C} \implies \Delta F \propto \sqrt{C} \)

\( \text{But with multiple overlapping lines: } \Delta F \propto \ln(C) \)

For a CO\(_2\) doubling from 280 to 560 ppm: \(\Delta F = 5.35 \ln(2) \approx 3.7\;\text{W/m}^2\). This is the benchmark for equilibrium climate sensitivity calculations.

0.5 Earth's Full Energy Budget

The complete energy budget includes not just radiative fluxes but also convective (sensible) and evaporative (latent) heat transport from the surface:

\( \frac{S}{4}(1-\alpha) = \varepsilon \sigma T_s^4 - F_{\text{atm}\downarrow} + F_{\text{SH}} + F_{\text{LH}} \)

Key flux values (Trenberth, Fasullo & Kiehl 2009):

\(\bullet\) Incoming solar: \(S/4 = 340\;\text{W/m}^2\)
\(\bullet\) Reflected solar: \(\sim 100\;\text{W/m}^2\) (albedo \(\alpha \approx 0.29\))
\(\bullet\) Absorbed solar (surface): \(\sim 185\;\text{W/m}^2\)
\(\bullet\) Absorbed solar (atmosphere): \(\sim 75\;\text{W/m}^2\)
\(\bullet\) Surface longwave up: \(\sim 398\;\text{W/m}^2\)
\(\bullet\) Atmospheric longwave down: \(\sim 342\;\text{W/m}^2\)
\(\bullet\) Sensible heat: \(\sim 17\;\text{W/m}^2\)
\(\bullet\) Latent heat: \(\sim 80\;\text{W/m}^2\)
\(\bullet\) OLR (top of atmosphere): \(\sim 239\;\text{W/m}^2\)

Albedo Components

The planetary albedo \(\alpha \approx 0.29\) arises from multiple contributors:

Clouds: \(\sim 0.15\) (largest contributor, especially low clouds)
Surface ice/snow: \(\sim 0.06\) (ice-albedo feedback)
Land surfaces: \(\sim 0.04\) (deserts high, forests low)
Rayleigh scattering: \(\sim 0.03\)
Aerosols: \(\sim 0.01\)

As discussed in the Climate & Biodiversity course, changes to albedo from ice loss and vegetation shifts create powerful climate feedbacks that amplify initial warming.

Earth Energy Balance Diagram

Earth energy balance based on Trenberth, Fasullo & Kiehl (2009). Values in W/m\(^2\).

Greenhouse Gas Overlap & Global Warming Potential

The radiative impact of each greenhouse gas depends on spectral overlap with other absorbers, particularly water vapor. The per-molecule forcing efficiency varies enormously:

CO\(_2\): \(\sim 1.4 \times 10^{-5}\) W/m\(^2\) per ppb — dominant due to high concentration
CH\(_4\): \(\sim 3.6 \times 10^{-4}\) W/m\(^2\) per ppb — 28x stronger per molecule
N\(_2\)O: \(\sim 3.0 \times 10^{-3}\) W/m\(^2\) per ppb — absorbs in a clear window
CFCs: \(\sim 0.1{-}0.3\) W/m\(^2\) per ppb — absorb in the 8-12 \(\mu\)m window

The Global Warming Potential (GWP) integrates radiative forcing of a pulse emission over time horizon \(H\), relative to CO\(_2\):

\( \text{GWP}_i(H) = \frac{\int_0^H a_i \cdot c_i(t)\,dt}{\int_0^H a_{\text{CO}_2} \cdot c_{\text{CO}_2}(t)\,dt} \)

For CH\(_4\): GWP\(_{20} = 81\), GWP\(_{100} = 27\) (lifetime \(\sim 12\) yr). For N\(_2\)O: GWP\(_{100} = 273\) (lifetime \(\sim 121\) yr).

Simulation: Planck Spectra at Different Temperatures

Compare the Planck spectral radiance for the Sun (\(\sim 5778\) K, scaled down) and Earth (\(\sim 255\) K and \(\sim 288\) K), showing the separation of shortwave and longwave windows:

Python

script.py1 lines

import numpy as np\nimport matplotlib.pyplot as plt\n\n# Constants\nh = 6.626e-34  # Planck constant\nc = 3.0e8      # Speed of light\nk_B = 1.381e-23  # Boltzmann constant\n\ndef planck_lambda(lam, T):\n    """Spectral radiance B_lambda(T) in W/m^2/sr/m"""\n    x = h * c / (lam * k_B * T)\n    return 2 * h * c**2 / lam**5 / (np.exp(x) - 1)\n\n# Wavelength ranges\nlam_sw = np.linspace(0.1e-6, 4e-6, 500)   # shortwave\nlam_lw = np.linspace(3e-6, 50e-6, 500)     # longwave\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))\nfig.patch.set_facecolor('#0a0a1a')\n\n# Solar spectrum (shortwave)\nfor T, color, label in [(5778, '#fbbf24', 'Sun (5778 K)'),\n                        (5000, '#f59e0b', 'Cooler star (5000 K)'),\n                        (6500, '#fef3c7', 'Hotter star (6500 K)')]:\n    B = planck_lambda(lam_sw, T) * 1e-6  # convert to W/m^2/sr/um\n    ax1.plot(lam_sw * 1e6, B, color=color, linewidth=2, label=label)\n    # Wien peak\n    lam_max = 2898e-6 / T\n    B_max = planck_lambda(lam_max, T) * 1e-6\n    ax1.axvline(lam_max * 1e6, color=color, linestyle='--', alpha=0.5)\n\nax1.set_xlabel('Wavelength (\u03bcm)', color='white', fontsize=12)\nax1.set_ylabel('Spectral Radiance (W/m\u00b2/sr/\u03bcm)', color='white', fontsize=12)\nax1.set_title('Solar (Shortwave) Spectra', color='#818cf8', fontsize=14, fontweight='bold')\nax1.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax1.set_facecolor('#0f0f2a')\nax1.tick_params(colors='white')\nax1.set_xlim(0.1, 4)\nax1.grid(True, alpha=0.2, color='#4f46e5')\n\n# Earth spectrum (longwave)\nfor T, color, label in [(288, '#f87171', 'Earth surface (288 K)'),\n                        (255, '#60a5fa', 'Effective T (255 K)'),\n                        (220, '#38bdf8', 'Tropopause (220 K)')]:\n    B = planck_lambda(lam_lw, T) * 1e-6\n    ax2.plot(lam_lw * 1e6, B, color=color, linewidth=2, label=label)\n    lam_max = 2898e-6 / T\n    ax2.axvline(lam_max * 1e6, color=color, linestyle='--', alpha=0.5)\n\n# Mark CO2 absorption band\nax2.axvspan(13, 17, alpha=0.15, color='#ef4444', label='CO\u2082 15 \u03bcm band')\nax2.axvspan(5.5, 7.5, alpha=0.1, color='#3b82f6', label='H\u2082O 6.3 \u03bcm band')\n\nax2.set_xlabel('Wavelength (\u03bcm)', color='white', fontsize=12)\nax2.set_ylabel('Spectral Radiance (W/m\u00b2/sr/\u03bcm)', color='white', fontsize=12)\nax2.set_title('Earth (Longwave) Spectra', color='#818cf8', fontsize=14, fontweight='bold')\nax2.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax2.set_facecolor('#0f0f2a')\nax2.tick_params(colors='white')\nax2.set_xlim(3, 50)\nax2.grid(True, alpha=0.2, color='#4f46e5')\n\nplt.tight_layout()\nplt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')\nprint('Planck spectra generated successfully')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: CO\(_2\) Absorption Spectrum

Modelling the CO\(_2\) bending mode absorption using a simplified Lorentzian line shape centred at 667 cm\(^{-1}\) (15 \(\mu\)m):

Python

script.py1 lines

import numpy as np\nimport matplotlib.pyplot as plt\n\n# Simplified CO2 absorption model\n# Main bending mode at 667 cm^-1 (15 um) with satellite lines\nwavenumber = np.linspace(500, 850, 2000)  # cm^-1\n\n# Lorentzian line shape\ndef lorentzian(nu, nu0, gamma, S):\n    return S * gamma / (np.pi * ((nu - nu0)**2 + gamma**2))\n\n# Main CO2 absorption lines (simplified)\nlines = [\n    (618, 2.0, 0.8), (624, 3.0, 1.2), (632, 4.5, 1.5),\n    (640, 6.0, 2.0), (648, 8.0, 2.5), (655, 12.0, 3.0),\n    (661, 18.0, 3.5), (667, 25.0, 4.0),  # Q branch center\n    (673, 18.0, 3.5), (679, 12.0, 3.0), (686, 8.0, 2.5),\n    (694, 6.0, 2.0), (702, 4.5, 1.5), (710, 3.0, 1.2),\n    (718, 2.0, 0.8), (725, 1.5, 0.6), (735, 1.0, 0.5),\n]\n\nfig, axes = plt.subplots(2, 2, figsize=(14, 10))\nfig.patch.set_facecolor('#0a0a1a')\n\n# Panel 1: Absorption coefficient\nax = axes[0, 0]\nkappa = np.zeros_like(wavenumber)\nfor nu0, S, gamma in lines:\n    kappa += lorentzian(wavenumber, nu0, gamma, S)\nax.plot(wavenumber, kappa, color='#818cf8', linewidth=1.5)\nax.fill_between(wavenumber, 0, kappa, alpha=0.3, color='#6366f1')\nax.set_xlabel('Wavenumber (cm\u207b\u00b9)', color='white', fontsize=11)\nax.set_ylabel('Absorption coefficient (arb.)', color='white', fontsize=11)\nax.set_title('CO\u2082 Absorption Coefficient', color='#818cf8', fontsize=13, fontweight='bold')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 2: Transmission for different path lengths\nax = axes[0, 1]\nfor path_mult, color, label in [(0.3, '#86efac', 'Thin path'),\n                                  (1.0, '#fbbf24', 'Standard'),\n                                  (3.0, '#f87171', 'Thick path')]:\n    tau = kappa * path_mult\n    T_trans = np.exp(-tau)\n    ax.plot(wavenumber, T_trans, color=color, linewidth=1.5, label=label)\nax.set_xlabel('Wavenumber (cm\u207b\u00b9)', color='white', fontsize=11)\nax.set_ylabel('Transmittance', color='white', fontsize=11)\nax.set_title('Beer-Lambert Transmission', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.set_ylim(-0.05, 1.05)\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 3: Emission to space (brightness temperature)\nax = axes[1, 0]\ndef planck_wn(nu_cm, T):\n    nu = nu_cm * 100 * 3e8  # cm^-1 to Hz\n    h, c, kB = 6.626e-34, 3e8, 1.381e-23\n    return 2*h*nu**3/c**2 / (np.exp(h*nu/(kB*T)) - 1)\n\nT_surface = 288\nT_tropo = 220\nB_surface = planck_wn(wavenumber, T_surface)\nB_tropo = planck_wn(wavenumber, T_tropo)\ntau_atm = kappa * 1.5\ntrans = np.exp(-tau_atm)\nI_toa = B_surface * trans + B_tropo * (1 - trans)\n\nax.plot(wavenumber, B_surface / B_surface.max(), color='#f87171', linewidth=1, alpha=0.5, label='Surface (288 K)')\nax.plot(wavenumber, B_tropo / B_surface.max(), color='#60a5fa', linewidth=1, alpha=0.5, label='Tropopause (220 K)')\nax.plot(wavenumber, I_toa / B_surface.max(), color='#fbbf24', linewidth=2, label='TOA emission')\nax.fill_between(wavenumber, I_toa/B_surface.max(), B_surface/B_surface.max(), alpha=0.15, color='#ef4444')\nax.set_xlabel('Wavenumber (cm\u207b\u00b9)', color='white', fontsize=11)\nax.set_ylabel('Normalized Radiance', color='white', fontsize=11)\nax.set_title('Top-of-Atmosphere Emission', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 4: Radiative forcing vs CO2\nax = axes[1, 1]\nC0 = 280  # preindustrial ppm\nC_range = np.linspace(280, 1200, 200)\ndelta_F = 5.35 * np.log(C_range / C0)\nax.plot(C_range, delta_F, color='#f87171', linewidth=2.5)\nax.axhline(y=5.35*np.log(2), color='#fbbf24', linestyle='--', alpha=0.7, label=f'2x CO\u2082: {5.35*np.log(2):.1f} W/m\u00b2')\nax.axhline(y=5.35*np.log(560/280), color='#fbbf24', linestyle='--', alpha=0.3)\nax.axvline(x=420, color='#86efac', linestyle=':', alpha=0.7, label='Current (420 ppm)')\nax.axvline(x=560, color='#fbbf24', linestyle=':', alpha=0.5, label='2x preindustrial')\nax.set_xlabel('CO\u2082 Concentration (ppm)', color='white', fontsize=11)\nax.set_ylabel('Radiative Forcing (W/m\u00b2)', color='white', fontsize=11)\nax.set_title('\u0394F = 5.35 ln(C/C\u2080)', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\nplt.tight_layout()\nplt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')\nprint('CO2 absorption simulation complete')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: Energy Balance Model

A zero-dimensional energy balance model with greenhouse feedback and ice-albedo feedback:

\( C \frac{dT}{dt} = \frac{S}{4}(1-\alpha(T)) - \varepsilon(T)\,\sigma\,T^4 \)

Python

script.py1 lines

import numpy as np\nimport matplotlib.pyplot as plt\n\n# 0D Energy Balance Model with feedbacks\nsigma = 5.67e-8\nS = 1361.0  # Solar constant W/m^2\nC = 4.0e8   # Heat capacity (J/m^2/K) ~ 100m mixed layer\ndt = 86400 * 30  # 30-day timestep\nyears = 200\nsteps = int(years * 365.25 * 86400 / dt)\n\ndef albedo(T):\n    """Ice-albedo feedback: higher albedo when cold"""\n    if T < 250:\n        return 0.6  # Snowball Earth\n    elif T > 290:\n        return 0.25  # Ice-free\n    else:\n        return 0.6 - 0.35 * (T - 250) / 40\n\ndef emissivity(T, CO2_ppm=400):\n    """Effective emissivity including greenhouse effect"""\n    # Base emissivity reduced by greenhouse gases\n    eps0 = 0.62  # base with preindustrial CO2\n    # Logarithmic CO2 effect\n    delta_eps = -0.02 * np.log(CO2_ppm / 280)\n    # Water vapor feedback (Clausius-Clapeyron)\n    delta_wv = -0.005 * (T - 288)  # stronger greenhouse when warmer\n    return max(0.3, min(1.0, eps0 + delta_eps + delta_wv))\n\nfig, axes = plt.subplots(2, 2, figsize=(14, 10))\nfig.patch.set_facecolor('#0a0a1a')\n\n# Panel 1: Temperature evolution for different CO2 levels\nax = axes[0, 0]\nfor co2, color, label in [(280, '#60a5fa', '280 ppm (preindustrial)'),\n                           (420, '#86efac', '420 ppm (current)'),\n                           (560, '#fbbf24', '560 ppm (2x)'),\n                           (840, '#f87171', '840 ppm (3x)')]:\n    T = 260.0  # initial\n    T_history = [T]\n    for i in range(steps):\n        a = albedo(T)\n        eps = emissivity(T, co2)\n        dTdt = (S/4 * (1 - a) - eps * sigma * T**4) / C\n        T += dTdt * dt\n        T_history.append(T)\n    time = np.arange(len(T_history)) * dt / (365.25 * 86400)\n    ax.plot(time, T_history, color=color, linewidth=2, label=label)\n\nax.set_xlabel('Time (years)', color='white', fontsize=11)\nax.set_ylabel('Global Mean Temperature (K)', color='white', fontsize=11)\nax.set_title('EBM: CO\u2082 Scenarios', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 2: Equilibrium T vs solar constant\nax = axes[0, 1]\nS_range = np.linspace(800, 2000, 300)\nT_eq_up = []\nT_eq_down = []\n# Sweep up\nT = 200.0\nfor s in S_range:\n    for _ in range(5000):\n        a = albedo(T)\n        eps = emissivity(T)\n        dTdt = (s/4 * (1 - a) - eps * sigma * T**4) / C\n        T += dTdt * dt\n    T_eq_up.append(T)\n# Sweep down\nT = 350.0\nfor s in reversed(S_range):\n    for _ in range(5000):\n        a = albedo(T)\n        eps = emissivity(T)\n        dTdt = (s/4 * (1 - a) - eps * sigma * T**4) / C\n        T += dTdt * dt\n    T_eq_down.append(T)\nT_eq_down.reverse()\nax.plot(S_range, T_eq_up, color='#f87171', linewidth=2, label='Warming branch')\nax.plot(S_range, T_eq_down, color='#60a5fa', linewidth=2, label='Cooling branch')\nax.axvline(x=1361, color='#fbbf24', linestyle='--', alpha=0.7, label='Current S=1361')\nax.set_xlabel('Solar Constant (W/m\u00b2)', color='white', fontsize=11)\nax.set_ylabel('Equilibrium Temperature (K)', color='white', fontsize=11)\nax.set_title('Hysteresis: Ice-Albedo Feedback', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 3: Climate sensitivity\nax = axes[1, 0]\nCO2_range = np.linspace(200, 1500, 200)\nT_eq_co2 = []\nfor co2 in CO2_range:\n    T = 285.0\n    for _ in range(8000):\n        a = albedo(T)\n        eps = emissivity(T, co2)\n        dTdt = (S/4 * (1 - a) - eps * sigma * T**4) / C\n        T += dTdt * dt\n    T_eq_co2.append(T)\nT_eq_co2 = np.array(T_eq_co2)\nax.plot(CO2_range, T_eq_co2 - 273.15, color='#f87171', linewidth=2.5)\nax.axvline(x=280, color='#60a5fa', linestyle='--', alpha=0.7, label='Preindustrial')\nax.axvline(x=420, color='#86efac', linestyle='--', alpha=0.7, label='Current')\nax.axvline(x=560, color='#fbbf24', linestyle='--', alpha=0.7, label='2x CO\u2082')\n# ECS\nT_280 = np.interp(280, CO2_range, T_eq_co2)\nT_560 = np.interp(560, CO2_range, T_eq_co2)\nECS = T_560 - T_280\nax.annotate(f'ECS = {ECS:.1f} K', xy=(560, T_560-273.15), xytext=(700, T_560-273.15+1),\n            color='#fbbf24', fontsize=12, fontweight='bold',\n            arrowprops=dict(arrowstyle='->', color='#fbbf24'))\nax.set_xlabel('CO\u2082 (ppm)', color='white', fontsize=11)\nax.set_ylabel('Temperature (\u00b0C)', color='white', fontsize=11)\nax.set_title('Equilibrium Climate Sensitivity', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 4: Greenhouse effect decomposition\nax = axes[1, 1]\nT_arr = np.linspace(220, 320, 200)\nOLR_no_ghg = sigma * T_arr**4\nOLR_with_ghg = np.array([emissivity(T, 400) * sigma * T**4 for T in T_arr])\nabs_solar = S/4 * np.array([1 - albedo(T) for T in T_arr])\n\nax.plot(T_arr, OLR_no_ghg, color='#94a3b8', linewidth=1.5, linestyle='--', label='No greenhouse')\nax.plot(T_arr, OLR_with_ghg, color='#f87171', linewidth=2, label='With greenhouse')\nax.plot(T_arr, abs_solar, color='#fbbf24', linewidth=2, label='Absorbed solar')\nax.fill_between(T_arr, OLR_with_ghg, OLR_no_ghg, alpha=0.1, color='#ef4444')\nax.set_xlabel('Surface Temperature (K)', color='white', fontsize=11)\nax.set_ylabel('Flux (W/m\u00b2)', color='white', fontsize=11)\nax.set_title('OLR vs Absorbed Solar', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\nax.set_xlim(220, 320)\n\nplt.tight_layout()\nplt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')\nprint('Energy balance model simulation complete')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

0.6 Climate Sensitivity & Feedback Analysis

The equilibrium climate sensitivity (ECS) is the steady-state global mean temperature change for a doubling of CO\(_2\). From the radiative forcing framework:

\( \Delta T_{\text{eq}} = \frac{\Delta F}{\lambda_0 - \sum_i \lambda_i} \)

where \(\lambda_0 = 4\sigma T_{\text{eff}}^3 \approx 3.3\;\text{W m}^{-2}\text{K}^{-1}\) is the Planck (no-feedback) response parameter, and \(\lambda_i\) are feedback parameters:

Water vapor: \(\lambda_{\text{WV}} \approx +1.8\;\text{W m}^{-2}\text{K}^{-1}\) (strongest positive feedback)
Lapse rate: \(\lambda_{\text{LR}} \approx -0.6\;\text{W m}^{-2}\text{K}^{-1}\) (negative, partially cancels WV)
Ice-albedo: \(\lambda_{\text{ice}} \approx +0.3\;\text{W m}^{-2}\text{K}^{-1}\) (positive)
Cloud: \(\lambda_{\text{cloud}} \approx +0.3{-}0.9\;\text{W m}^{-2}\text{K}^{-1}\) (largest uncertainty)

Without feedbacks: \(\Delta T_0 = \Delta F / \lambda_0 = 3.7/3.3 \approx 1.1\) K. With feedbacks, the amplification (gain) factor:

\( G = \frac{1}{1 - \sum_i \lambda_i / \lambda_0} \approx 2{-}4 \implies \text{ECS} \approx 2.5{-}4.5\;\text{K} \)

Transient Climate Response

The transient climate response (TCR) accounts for ocean heat uptake, which delays the full equilibrium warming:

\( C \frac{dT}{dt} = \Delta F - \lambda T - F_{\text{ocean}} \)

\( \text{TCR} \approx \text{ECS} \times \frac{\lambda}{\lambda + \kappa} \approx 0.5{-}0.8 \times \text{ECS} \)

where \(\kappa\) parameterizes ocean heat uptake efficiency. Current best estimates give TCR \(\approx 1.4{-}2.2\) K (IPCC AR6). These concepts connect directly to the energy balance models explored in the simulation above.

0.7 Multi-Layer Atmospheric Model

Extending the one-layer greenhouse model to \(N\) absorbing layers, each transparent to shortwave but opaque to longwave:

\( T_{\text{surface}} = (N+1)^{1/4} \cdot T_{\text{eff}} \)

Proof by induction: For the \(N\)-layer model, the top layer must radiate \(\sigma T_{\text{eff}}^4\) upward. Each layer emits equally up and down. The\(k\)-th layer from the top satisfies:

\( \sigma T_k^4 = (k+1) \sigma T_{\text{eff}}^4 \implies T_k = (k+1)^{1/4} T_{\text{eff}} \)

\( T_{\text{surface}} = T_{N+1} = (N+1)^{1/4} T_{\text{eff}} \)

For \(N = 1\): \(T_s = 2^{1/4} \times 255 = 303\) K (as derived in Section 0.2). For \(N = 2\): \(T_s = 3^{1/4} \times 255 = 335\) K (too hot, showing that the real atmosphere is not fully opaque at all wavelengths).

A more realistic approach uses partial absorption with emissivity \(\varepsilon < 1\):

\( T_s^4 = T_{\text{eff}}^4 \left(1 + \frac{\varepsilon}{2 - \varepsilon}\right) \)

Single-layer with partial absorption \(\varepsilon\)

Setting \(\varepsilon \approx 0.77\) yields \(T_s \approx 288\) K, matching the observed global mean surface temperature. This effective emissivity encapsulates the combined effect of all greenhouse gases, including CO\(_2\), H\(_2\)O, CH\(_4\), and N\(_2\)O.

Simulation: Feedback Analysis & Multi-Layer Model

Comparing the no-feedback, single-layer, and multi-layer atmospheric models with feedback decomposition:

Python

script.py1 lines

import numpy as np\nimport matplotlib.pyplot as plt\n\nfig, axes = plt.subplots(2, 2, figsize=(14, 10))\nfig.patch.set_facecolor('#0a0a1a')\nsigma = 5.67e-8\nT_eff = 255.0\n\n# Panel 1: Multi-layer model\nax = axes[0, 0]\nN_layers = np.arange(0, 8)\nT_surface = (N_layers + 1)**0.25 * T_eff\nax.bar(N_layers, T_surface - 273.15, color='#818cf8', edgecolor='#4f46e5', alpha=0.8)\nax.axhline(y=15, color='#f87171', linestyle='--', alpha=0.7, label='Observed (15\u00b0C)')\nax.axhline(y=-18, color='#60a5fa', linestyle='--', alpha=0.7, label='No greenhouse (-18\u00b0C)')\nfor i, t in enumerate(T_surface):\n    ax.text(i, t - 273.15 + 2, f'{t:.0f} K', ha='center', color='white', fontsize=9)\nax.set_xlabel('Number of Layers N', color='white', fontsize=11)\nax.set_ylabel('Surface Temperature (\u00b0C)', color='white', fontsize=11)\nax.set_title('Multi-Layer Greenhouse Model', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 2: Partial emissivity model\nax = axes[0, 1]\neps_range = np.linspace(0, 1, 200)\nT_s_partial = T_eff * (1 + eps_range / (2 - eps_range))**0.25\nax.plot(eps_range, T_s_partial, color='#f87171', linewidth=2.5)\nax.axhline(y=288, color='#fbbf24', linestyle='--', alpha=0.7, label='Observed (288 K)')\neps_match = 2 * ((288/T_eff)**4 - 1) / ((288/T_eff)**4)\nax.axvline(x=eps_match, color='#fbbf24', linestyle=':', alpha=0.5, label=f'\u03b5 = {eps_match:.2f}')\nax.set_xlabel('Effective Emissivity \u03b5', color='white', fontsize=11)\nax.set_ylabel('Surface Temperature (K)', color='white', fontsize=11)\nax.set_title('Partial Absorption Model', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\n\n# Panel 3: Feedback decomposition\nax = axes[1, 0]\nfeedbacks = ['Planck\\n(no FB)', 'Water\\nVapor', 'Lapse\\nRate', 'WV+LR\\ncombined', 'Ice-\\nAlbedo', 'Cloud', 'Total\\nNet']\nlambda_vals = [-3.3, 1.8, -0.6, 1.2, 0.3, 0.6, 1.8+(-0.6)+0.3+0.6]\ncolors_fb = ['#60a5fa', '#f87171', '#86efac', '#fbbf24', '#38bdf8', '#c084fc', '#818cf8']\nbars = ax.bar(range(len(feedbacks)), lambda_vals, color=colors_fb, edgecolor='white', linewidth=0.5, alpha=0.85)\nax.axhline(y=0, color='white', linewidth=0.5)\nfor i, v in enumerate(lambda_vals):\n    ax.text(i, v + (0.1 if v >= 0 else -0.2), f'{v:.1f}', ha='center', color='white', fontsize=9)\nax.set_xticks(range(len(feedbacks)))\nax.set_xticklabels(feedbacks, fontsize=9)\nax.set_ylabel('Feedback Parameter (W/m\u00b2/K)', color='white', fontsize=11)\nax.set_title('Climate Feedback Decomposition', color='#818cf8', fontsize=13, fontweight='bold')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5', axis='y')\n\n# Panel 4: ECS probability distribution\nax = axes[1, 1]\nnp.random.seed(42)\nN_samples = 100000\nlam_0 = 3.3\nlam_wv = np.random.normal(1.8, 0.3, N_samples)\nlam_lr = np.random.normal(-0.6, 0.15, N_samples)\nlam_ice = np.random.normal(0.3, 0.1, N_samples)\nlam_cloud = np.random.normal(0.6, 0.4, N_samples)\nlam_total = lam_wv + lam_lr + lam_ice + lam_cloud\nECS = 3.7 / (lam_0 - lam_total)\nECS = ECS[(ECS > 0) & (ECS < 10)]\nax.hist(ECS, bins=80, density=True, color='#818cf8', alpha=0.7, edgecolor='#4f46e5')\nax.axvline(x=np.median(ECS), color='#fbbf24', linestyle='--', linewidth=2,\n           label=f'Median = {np.median(ECS):.1f} K')\nax.axvline(x=np.percentile(ECS, 17), color='#60a5fa', linestyle=':', alpha=0.7,\n           label=f'17th pct = {np.percentile(ECS, 17):.1f} K')\nax.axvline(x=np.percentile(ECS, 83), color='#f87171', linestyle=':', alpha=0.7,\n           label=f'83rd pct = {np.percentile(ECS, 83):.1f} K')\nax.set_xlabel('ECS (K)', color='white', fontsize=11)\nax.set_ylabel('Probability Density', color='white', fontsize=11)\nax.set_title('ECS Distribution (Monte Carlo)', color='#818cf8', fontsize=13, fontweight='bold')\nax.legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#4f46e5', labelcolor='white')\nax.set_facecolor('#0f0f2a')\nax.tick_params(colors='white')\nax.grid(True, alpha=0.2, color='#4f46e5')\nax.set_xlim(0, 8)\n\nplt.tight_layout()\nplt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')\nprint('Feedback analysis simulation complete')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

References

Pierrehumbert, R. T. (2010). Principles of Planetary Climate. Cambridge University Press.
Trenberth, K. E., Fasullo, J. T., & Kiehl, J. (2009). Earth's global energy budget. Bulletin of the American Meteorological Society, 90(3), 311-324.
Myhre, G., Highwood, E. J., Shine, K. P., & Stordal, F. (1998). New estimates of radiative forcing due to well mixed greenhouse gases. Geophysical Research Letters, 25(14), 2715-2718.
Petty, G. W. (2006). A First Course in Atmospheric Radiation (2nd ed.). Sundog Publishing.
Hartmann, D. L. (2016). Global Physical Climatology (2nd ed.). Elsevier.
Goody, R. M. & Yung, Y. L. (1989). Atmospheric Radiation: Theoretical Basis (2nd ed.). Oxford University Press.