Part III: Climate & Oceans | Chapter 3

Ocean Circulation

Ekman transport, Sverdrup balance, thermohaline circulation, and AMOC

3.1 Ekman Transport

V. Walfrid Ekman (1905) solved the problem of wind-driven ocean currents, inspired by Fridtjof Nansen's observation that ice drifts at 20-40$°$ to the right of the wind in the Arctic.

Derivation 1: The Ekman Spiral

For steady, horizontally uniform flow driven by surface wind stress $\boldsymbol{\tau}$, the balance between Coriolis and viscous forces gives:

$$-fv = A_z \frac{\partial^2 u}{\partial z^2}, \quad fu = A_z \frac{\partial^2 v}{\partial z^2}$$

where $A_z$ is the vertical eddy viscosity. With boundary conditions of wind stress at the surface and vanishing velocity at depth, the solution is:

$$u(z) = V_0 e^{z/D_E} \cos\left(\frac{\pi}{4} + \frac{z}{D_E}\right)$$

$$v(z) = V_0 e^{z/D_E} \sin\left(\frac{\pi}{4} + \frac{z}{D_E}\right)$$

where $D_E = \sqrt{2A_z/|f|}$ is the Ekman depth and $V_0 = \tau/(\rho\sqrt{A_z|f|})$is the surface current speed. The surface current is directed 45$°$ to the right of the wind (Northern Hemisphere). The current direction rotates clockwise with depth while the speed decreases exponentially.

The depth-integrated Ekman transport is exactly perpendicular to the wind:

$$M_E = \int_{-\infty}^{0} \rho \mathbf{v}_E \, dz = \frac{\hat{k} \times \boldsymbol{\tau}}{f}$$

This means the total mass transport in the Ekman layer is 90$°$ to the right of the wind in the Northern Hemisphere. For typical ocean conditions ($A_z \approx 0.01$-0.1 m$^2$/s, $f \approx 10^{-4}$ s$^{-1}$), the Ekman depth is $D_E \approx$ 15-50 m.

3.2 Ekman Pumping and Sverdrup Balance

Derivation 2: Ekman Pumping Velocity

Spatial variations in wind stress cause convergence or divergence of the Ekman transport, which drives vertical motion at the base of the Ekman layer. The Ekman pumping velocity is:

$$w_E = \nabla \cdot \left(\frac{\boldsymbol{\tau}}{\rho f}\right) = \frac{1}{\rho f}\left(\frac{\partial \tau_y}{\partial x} - \frac{\partial \tau_x}{\partial y}\right) + \frac{\beta \tau_x}{\rho f^2}$$

Simplifying for zonal wind only ($\tau_y = 0$):

$$w_E = -\frac{1}{\rho}\frac{\partial}{\partial y}\left(\frac{\tau_x}{f}\right) = \frac{1}{\rho f}\text{curl}_z(\boldsymbol{\tau})$$

Negative wind stress curl (as in subtropical regions) drives downwelling (Ekman pumping), depressing the thermocline and creating the subtropical gyres. Positive curl in subpolar regions drives upwelling.

Derivation 3: The Sverdrup Balance

Harald Sverdrup (1947) showed that the depth-integrated meridional transport in the ocean interior is determined by the wind stress curl. Starting from the vorticity equation for a barotropic ocean:

$$\beta V = f \frac{\partial w}{\partial z} = \frac{1}{\rho}\text{curl}_z(\boldsymbol{\tau})$$

where $V$ is the depth-integrated meridional transport per unit zonal width and $\beta = df/dy$. This is the Sverdrup relation:

$$V = \frac{1}{\rho \beta}\text{curl}_z(\boldsymbol{\tau})$$

The Sverdrup transport can be integrated from the eastern boundary to give the total meridional flow. In the subtropical North Atlantic, the wind stress curl drives ~30 Sv of equatorward Sverdrup transport. This transport must be compensated by a narrow, intense western boundary current: the Gulf Stream.

3.3 Western Boundary Currents

Derivation 4: The Stommel Model of Western Intensification

Henry Stommel (1948) explained western intensification by including bottom friction in the vorticity balance. The depth-integrated vorticity equation for a rectangular basin is:

$$\beta \frac{\partial \psi}{\partial x} + R \nabla^2 \psi = \frac{1}{\rho}\text{curl}_z(\boldsymbol{\tau})$$

where $\psi$ is the transport stream function and $R$ is a linear friction coefficient. The $\beta$ term breaks the east-west symmetry: the planetary vorticity gradient means that northward flow gains vorticity ($\beta v > 0$), while friction dissipates it.

The western boundary layer width scales as:

$$\delta_W \sim \frac{R}{\beta}$$

For the Gulf Stream, $\delta_W \approx$ 100 km. The Munk (1950) model uses lateral friction instead, giving $\delta_W \sim (A_H/\beta)^{1/3}$, where $A_H$ is horizontal eddy viscosity. Both models correctly predict western intensification but differ in boundary layer structure.

Major western boundary currents include the Gulf Stream (North Atlantic, ~30 Sv), Kuroshio (North Pacific, ~25 Sv), Agulhas Current (Indian Ocean), Brazil Current, and East Australian Current. These currents transport enormous amounts of heat poleward, influencing regional and global climate.

3.4 Thermohaline Circulation and AMOC

Derivation 5: Stommel's Two-Box Model of THC

Stommel (1961) proposed a simple two-box model to understand the stability of the thermohaline circulation. Two well-mixed boxes (tropical and polar) exchange water at a rate proportional to the density difference:

$$q = k(\rho_2 - \rho_1) = k[\beta_S(S_2 - S_1) - \alpha_T(T_2 - T_1)]$$

where $\alpha_T$ and $\beta_S$ are the thermal expansion and haline contraction coefficients. Temperature is restored to surface values on a fast timescale, while salinity is controlled by net freshwater flux. The steady-state equation for salinity difference $\Delta S$:

$$|q| \Delta S = F_w$$

This model has multiple equilibria:

Thermally-driven mode ($q > 0$): Temperature dominates density contrast; deep water forms at high latitudes (like present AMOC)
Salinity-driven mode ($q < 0$): Freshwater input raises salinity contrast enough to reverse circulation

The system exhibits hysteresis: if freshwater input exceeds a critical threshold, the AMOC collapses and cannot recover until freshwater input drops well below the threshold. This is a potential tipping point in the climate system. The AMOC transports approximately 17 Sv of warm water northward and ~1.3 PW of heat, maintaining northern Europe about 5-10$°$C warmer than equivalent latitudes.

3.5 Deep Water Formation and Abyssal Circulation

Deep water forms at only a few locations: the Greenland-Norwegian Sea (North Atlantic Deep Water, NADW), the Weddell and Ross Seas (Antarctic Bottom Water, AABW), and marginally in the Labrador Sea. The deep ocean is filled by these cold, dense water masses that sink and spread along the ocean floor.

The Stommel-Arons (1960) model of abyssal circulation shows that uniform upwelling through the thermocline, combined with the $\beta$-effect, drives deep western boundary currents. The abyssal turnover time is approximately 1000-2000 years, meaning that changes in deep water formation have long-lasting effects on climate.

Modern observations from the RAPID array at 26.5$°$N show the AMOC strength has been approximately 17 Sv since 2004, with significant variability on seasonal to decadal timescales. Paleoceanographic evidence from $\delta^{13}$C and Pa/Th ratios suggests the AMOC was substantially weakened during the last glacial period and collapsed during Heinrich events, contributing to abrupt climate changes.

The Bipolar Seesaw

When the AMOC weakens, northward heat transport decreases, cooling the North Atlantic while warming the South Atlantic. This "bipolar seesaw" is evident in the anti-phased temperature records from Greenland and Antarctic ice cores during Dansgaard-Oeschger events. The thermal seesaw timescale depends on the AMOC adjustment time (~100-200 years for circulation changes to propagate through the Atlantic basin).

The AMOC is considered a potential tipping element in the climate system. Climate models suggest that continued global warming and freshwater input from Greenland ice sheet melting could weaken the AMOC by 25-50% by 2100, with a small but non-negligible probability of complete collapse. Such an event would cool northwestern Europe by several degrees, shift tropical rainfall belts, and disrupt marine ecosystems globally.

3.6 Oceanic Heat Transport and Climate

The ocean absorbs approximately 93% of the excess heat from anthropogenic forcing. The total oceanic heat transport peaks at ~2 PW (petawatts) near 20$°$N, with roughly equal contributions from the wind-driven (shallow) and thermohaline (deep) circulations. The meridional heat transport can be decomposed:

$$H = \rho C_P \int_0^L \int_{-D}^{0} v T \, dz \, dx = H_{\text{mean}} + H_{\text{eddy}} + H_{\text{overturning}}$$

where the overturning component dominates in the Atlantic (the AMOC carries warm surface water northward and cold deep water southward) and the gyre component dominates in the Pacific. The deep ocean is warming at a rate of ~0.5 W/m$^2$, contributing to sea level rise through thermal expansion (~1.4 mm/yr).

Ocean Acidification

As the ocean absorbs CO$_2$, carbonic acid forms and pH decreases. The key reactions are:

$$\text{CO}_2 + \text{H}_2\text{O} \rightleftharpoons \text{H}_2\text{CO}_3 \rightleftharpoons \text{H}^+ + \text{HCO}_3^- \rightleftharpoons 2\text{H}^+ + \text{CO}_3^{2-}$$

Ocean surface pH has decreased from ~8.18 (pre-industrial) to ~8.07 (present), a 25% increase in hydrogen ion concentration. The aragonite saturation state$\Omega_A = [\text{Ca}^{2+}][\text{CO}_3^{2-}]/K_{sp}$ is declining, threatening coral reefs and shell-forming organisms. If CO$_2$ continues rising, surface waters in polar regions will become undersaturated ($\Omega_A < 1$) by mid-century.

Computational Lab: Ocean Circulation

Ekman Spiral, Sverdrup Transport, and Thermohaline Box Model

Python

script.py378 lines

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

# ============================================================
# OCEAN CIRCULATION SIMULATION
# ============================================================
print("OCEAN CIRCULATION SIMULATION")
print("=" * 70)

Omega = 7.292e-5
rho_w = 1025.0

# ============================================================
# 1. Ekman Spiral
# ============================================================
print()
print("1. Ekman Spiral Profile")
print("-" * 60)

lat = 45.0
f = 2 * Omega * np.sin(np.radians(lat))
Az = 0.01  # eddy viscosity (m^2/s)
DE = np.sqrt(2 * Az / abs(f))
tau_x = 0.1  # wind stress (N/m^2)
V0 = tau_x / (rho_w * np.sqrt(Az * abs(f)))

print(f"  Latitude: {lat} deg, f = {f:.4e} s^-1")
print(f"  Eddy viscosity: {Az} m^2/s")
print(f"  Ekman depth: {DE:.1f} m")
print(f"  Surface current speed: {V0*100:.2f} cm/s")
print()
print(f"  {'Depth (m)':14} {'u (cm/s)':12} {'v (cm/s)':12} {'Speed (cm/s)':14} {'Direction (deg)'}")
print("  " + "-" * 60)

for z in [0, -5, -10, -15, -20, -30, -50, -75, -100]:
    u = V0 * np.exp(z/DE) * np.cos(np.pi/4 + z/DE) * 100
    v = V0 * np.exp(z/DE) * np.sin(np.pi/4 + z/DE) * 100
    speed = np.sqrt(u**2 + v**2)
    direction = np.degrees(np.arctan2(v, u)) % 360
    print(f"  {z:12} {u:10.3f} {v:10.3f} {speed:12.3f} {direction:10.1f}")

# Ekman transport
ME = tau_x / abs(f)  # kg/(m*s)
print(f"  Total Ekman transport: {ME:.2f} kg/(m*s) = {ME/rho_w*1000:.2f} m^2/s per m of coastline")
print(f"  Direction: 90 deg to right of wind (to the east for northward wind)")

# ============================================================
# 2. Sverdrup Transport
# ============================================================
print()
print()
print("2. Sverdrup Transport in Subtropical Gyre")
print("-" * 60)

# Wind stress: tau_x = -tau_max * cos(pi*y/L)
a_earth = 6.371e6
beta = 2 * Omega * np.cos(np.radians(30)) / a_earth
tau_max = 0.1  # N/m^2
L_basin = 60.0  # degrees latitude extent
L_m = L_basin * np.pi / 180 * a_earth  # meters

print(f"  beta = {beta:.2e} 1/(m*s)")
print(f"  Basin: 10N to 50N, width ~ 60 deg longitude")
print()
print(f"  {'Latitude':12} {'tau_x (N/m2)':14} {'curl_tau':14} {'V_sv (Sv/deg)':16} {'V_sv (m/s)'}")
print("  " + "-" * 60)

total_sv = 0
for lat_i in range(10, 52, 5):
    y_frac = (lat_i - 10) / 40.0
    tau_x_local = -tau_max * np.cos(np.pi * y_frac)
    # d(tau_x)/dy approximation
    dy = 5 * np.pi / 180 * a_earth
    tau_x_next = -tau_max * np.cos(np.pi * (y_frac + 5/40))
    curl_tau = -(tau_x_next - tau_x_local) / dy
    f_local = 2 * Omega * np.sin(np.radians(lat_i))
    beta_local = 2 * Omega * np.cos(np.radians(lat_i)) / a_earth

# Sverdrup meridional transport per unit width
    V_sv = curl_tau / (rho_w * beta_local)  # m^2/s
    # Multiply by basin width to get Sv
    basin_width = 60 * np.pi / 180 * a_earth * np.cos(np.radians(lat_i))
    V_total_Sv = V_sv * basin_width / 1e6

print(f"  {lat_i:10} {tau_x_local:12.4f} {curl_tau:12.2e} {V_total_Sv:14.1f} {V_sv:12.2e}")
    total_sv += abs(V_total_Sv)

print(f"  Approximate Gulf Stream transport: ~30 Sv (observed)")

# ============================================================
# 3. Stommel Two-Box Model
# ============================================================
print()
print()
print("3. Stommel Two-Box Model: THC Stability")
print("-" * 60)

alpha_T = 2e-4  # thermal expansion (1/K)
beta_S = 7.6e-4  # haline contraction (1/psu)
k_flow = 1.5e-6  # flow constant (1/s)
dT = 20.0  # temperature difference (K)

print(f"  alpha_T = {alpha_T}, beta_S = {beta_S}")
print(f"  Temperature diff = {dT} K")
print()
print(f"  {'Fw (Sv)':12} {'dS (psu)':12} {'q (Sv)':12} {'Regime':20} {'Stable?'}")
print("  " + "-" * 60)

# Fw parameterized as fraction of reference value
for Fw_sv in [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4]:
    # Find steady states: q * dS = Fw, q = k*(alpha_T*dT - beta_S*dS)
    # k*(alpha_T*dT - beta_S*dS)*dS = Fw
    # Let x = beta_S*dS, then k*(alpha_T*dT - x)*x/beta_S = Fw
    # Quadratic: -k*x^2/beta_S + k*alpha_T*dT*x/beta_S = Fw
    a_coeff = -k_flow / beta_S
    b_coeff = k_flow * alpha_T * dT / beta_S
    c_coeff = -Fw_sv * 1e6  # convert to m^3/s
    discriminant = b_coeff**2 - 4*a_coeff*c_coeff
    if discriminant >= 0:
        x1 = (-b_coeff + np.sqrt(discriminant)) / (2*a_coeff)
        x2 = (-b_coeff - np.sqrt(discriminant)) / (2*a_coeff)
        for x in [x1, x2]:
            dS = x / beta_S
            q = k_flow * (alpha_T * dT - beta_S * dS) * 1e6 / 1e6  # rough Sv scaling
            q_val = alpha_T * dT - x
            regime = "Thermal" if q_val > 0 else "Haline"
            stable = "stable" if (regime == "Thermal" and dS > 0) or (regime == "Haline") else "unstable"
            print(f"  {Fw_sv:10.2f} {dS:10.2f} {q_val*1e6/1e6:10.4f} {regime:18} {stable}")
    else:
        print(f"  {Fw_sv:10.2f} {'no real solution'}")

# ============================================================
# 4. Western Boundary Current Width
# ============================================================
print()
print()
print("4. Western Boundary Current Properties")
print("-" * 55)

currents = [
    ("Gulf Stream", 30, 25, 100, 1000),
    ("Kuroshio", 25, 30, 80, 800),
    ("Agulhas", -35, 70, 60, 2200),
    ("Brazil", -20, 15, 50, 500),
    ("E. Australian", -30, 10, 40, 300),
]

print(f"  {'Current':18} {'Lat':8} {'Transport (Sv)':16} {'Width (km)':12} {'Depth (m)'}")
print("  " + "-" * 58)

for name, lat_c, transport, width, depth in currents:
    print(f"  {name:16} {lat_c:6} {transport:14} {width:10} {depth}")

# ============================================================
# 5. Seawater Equation of State
# ============================================================
print()
print()
print("5. Seawater Density at Key Locations")
print("-" * 65)

print(f"  {'Location':20} {'T (C)':10} {'S (psu)':10} {'rho (kg/m3)':14} {'sigma-t'}")
print("  " + "-" * 58)

locations = [
    ("Tropical surface", 28, 35.0),
    ("Subtropical", 20, 36.0),
    ("N. Atlantic deep", 2, 34.9),
    ("AABW", -0.5, 34.66),
    ("Mediterranean", 13, 38.4),
    ("Arctic surface", -1.8, 33.0),
    ("Labrador Sea", 3, 34.85),
]

for name, T, S in locations:
    # UNESCO equation of state (simplified)
    rho = 999.842594 + 6.793952e-2*T - 9.095290e-3*T**2 + 1.001685e-4*T**3
    rho += (0.824493 - 4.0899e-3*T + 7.6438e-5*T**2) * S
    rho += (-5.72466e-3 + 1.0227e-4*T) * S**1.5 + 4.8314e-4 * S**2
    sigma_t = rho - 1000
    print(f"  {name:18} {T:8.1f} {S:8.2f} {rho:12.2f} {sigma_t:8.2f}")

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

# Plot 1: Ekman Spiral (plan view)
ax = axes[0, 0]
z_arr = np.linspace(0, -4*DE, 200)
u_arr = V0 * np.exp(z_arr/DE) * np.cos(np.pi/4 + z_arr/DE) * 100
v_arr = V0 * np.exp(z_arr/DE) * np.sin(np.pi/4 + z_arr/DE) * 100

ax.plot(u_arr, v_arr, 'b-', linewidth=2)
ax.plot(u_arr[0], v_arr[0], 'ro', markersize=10, label='Surface')
# Mark depth levels
for depth_mark in [DE, 2*DE, 3*DE]:
    idx = np.argmin(abs(z_arr + depth_mark))
    ax.plot(u_arr[idx], v_arr[idx], 'ko', markersize=6)
    ax.annotate(f'{depth_mark:.0f}m', (u_arr[idx], v_arr[idx]),
                textcoords="offset points", xytext=(5, 5), fontsize=8)

# Wind arrow
ax.annotate('', xy=(0, V0*100*0.8), xytext=(0, 0),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax.text(0.5, V0*100*0.4, 'Wind', color='red', fontsize=10)

# Transport arrow
ax.annotate('', xy=(V0*100*0.5*DE/20, 0), xytext=(0, 0),
            arrowprops=dict(arrowstyle='->', color='green', lw=2))
ax.text(V0*100*0.3, -0.3, 'Transport', color='green', fontsize=10)

ax.set_xlabel('u velocity (cm/s)')
ax.set_ylabel('v velocity (cm/s)')
ax.set_title(f'Ekman Spiral (plan view, {lat}N)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')

# Plot 2: Sverdrup stream function
ax = axes[0, 1]
nx_grid, ny_grid = 50, 40
x_grid = np.linspace(0, 60, nx_grid)  # longitude degrees
y_grid = np.linspace(10, 50, ny_grid)  # latitude degrees
X_grid, Y_grid = np.meshgrid(x_grid, y_grid)

# Wind stress: tau_x = -tau_max * cos(pi*(Y-10)/40)
tau_grid = -tau_max * np.cos(np.pi * (Y_grid - 10) / 40)
# Curl of tau (d/dy only)
dy_grid = (y_grid[1] - y_grid[0]) * np.pi / 180 * a_earth
curl_tau_grid = np.gradient(-tau_grid, dy_grid, axis=0)

# Sverdrup stream function (integrated from east)
psi_sv = np.zeros_like(X_grid)
dx_grid = (x_grid[1] - x_grid[0]) * np.pi / 180 * a_earth * np.cos(np.radians(Y_grid))
beta_grid = 2 * Omega * np.cos(np.radians(Y_grid)) / a_earth

for j in range(ny_grid):
    for i in range(nx_grid-2, -1, -1):
        psi_sv[j, i] = psi_sv[j, i+1] - curl_tau_grid[j, i] / (rho_w * beta_grid[j, i]) * dx_grid[j, i]

psi_sv_Sv = psi_sv / 1e6

cf = ax.contourf(X_grid, Y_grid, psi_sv_Sv, levels=20, cmap='RdBu_r')
plt.colorbar(cf, ax=ax, label='Transport (Sv)')
ax.contour(X_grid, Y_grid, psi_sv_Sv, levels=10, colors='black', linewidths=0.5)
ax.set_xlabel('Longitude (deg from west)')
ax.set_ylabel('Latitude (N)')
ax.set_title('Sverdrup Stream Function')

# Plot 3: THC hysteresis
ax = axes[1, 0]
Fw_range = np.linspace(0, 0.5, 500)
# Find equilibria for each Fw
thermal_branch = []
haline_branch = []

for Fw in Fw_range:
    a_c = -k_flow / beta_S
    b_c = k_flow * alpha_T * dT / beta_S
    c_c = -Fw * 1e6
    disc = b_c**2 - 4*a_c*c_c
    if disc >= 0:
        x1 = (-b_c + np.sqrt(disc)) / (2*a_c)
        x2 = (-b_c - np.sqrt(disc)) / (2*a_c)
        for x_val in [x1, x2]:
            q_val = alpha_T * dT - x_val
            if q_val > 0:
                thermal_branch.append((Fw, q_val*1e7))
            else:
                haline_branch.append((Fw, q_val*1e7))

if thermal_branch:
    tb = np.array(thermal_branch)
    ax.plot(tb[:, 0], tb[:, 1], 'r-', linewidth=2, label='Thermal mode (AMOC on)')
if haline_branch:
    hb = np.array(haline_branch)
    ax.plot(hb[:, 0], hb[:, 1], 'b-', linewidth=2, label='Haline mode (AMOC off)')

ax.axhline(y=0, color='gray', linewidth=0.5)
ax.set_xlabel('Freshwater forcing (Sv)')
ax.set_ylabel('Overturning strength (arb. units)')
ax.set_title('Stommel Box Model: THC Hysteresis')
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 4: T-S diagram
ax = axes[1, 1]
T_plot = np.linspace(-2, 30, 100)
S_plot = np.linspace(33, 38, 100)
T_mesh, S_mesh = np.meshgrid(T_plot, S_plot)

# Density calculation
rho_mesh = 999.842594 + 6.793952e-2*T_mesh - 9.095290e-3*T_mesh**2
rho_mesh += (0.824493 - 4.0899e-3*T_mesh) * S_mesh
sigma_mesh = rho_mesh - 1000

cs = ax.contour(S_mesh, T_mesh, sigma_mesh, levels=np.arange(20, 32, 0.5),
                colors='gray', linewidths=0.5)
ax.clabel(cs, fontsize=7, fmt='%.1f')

# Plot water masses
water_masses = {
    'NADW': (34.9, 2.0),
    'AABW': (34.66, -0.5),
    'Tropical': (35.0, 28.0),
    'Mediterranean': (38.4, 13.0),
    'Labrador': (34.85, 3.0),
    'Arctic': (33.0, -1.8),
}

colors_wm = ['red', 'blue', 'orange', 'purple', 'cyan', 'green']
for (name, (s, t)), color in zip(water_masses.items(), colors_wm):
    ax.plot(s, t, 'o', markersize=10, color=color, label=name)

ax.set_xlabel('Salinity (psu)')
ax.set_ylabel('Temperature (C)')
ax.set_title('T-S Diagram with Water Masses')
ax.legend(fontsize=8, loc='upper left')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Ocean Heat Content and Sea Level Rise")
print("-" * 60)

print(f"  {'Depth range':20} {'OHC change (ZJ)':18} {'SLR contribution (mm/yr)'}")
print("  " + "-" * 55)

ohc_data = [
    ("0-300 m", 150, 0.9),
    ("300-700 m", 80, 0.4),
    ("700-2000 m", 60, 0.3),
    ("2000 m-bottom", 25, 0.1),
    ("Total ocean", 315, 1.7),
]

for depth, ohc, slr in ohc_data:
    print(f"  {depth:18} {ohc:16} {slr:10.1f}")

print()
print("  Ocean pH Change and Acidification:")
print("  " + "-" * 45)

for pco2 in [280, 350, 400, 450, 500, 560, 700, 1000]:
    # Simplified pH calculation
    pH = 8.18 - 0.7 * np.log10(pco2 / 280)
    omega_arag = 4.0 * (280/pco2)**0.5  # approximate
    state = "supersaturated" if omega_arag > 1 else "UNDERSATURATED"
    print(f"    pCO2 = {pco2} ppm: pH = {pH:.2f}, Omega_arag = {omega_arag:.2f} ({state})")

# Ekman depth for different conditions
print()
print()
print("7. Ekman Depth and Transport at Various Latitudes")
print("-" * 60)

print(f"  {'Latitude':12} {'f (s^-1)':14} {'D_E (m)':12} {'M_E (m^2/s)':14}")
print("  " + "-" * 50)

tau_ex = 0.1  # N/m^2
Az_ex = 0.01  # m^2/s

for lat_ek in [5, 10, 15, 20, 30, 45, 60, 75]:
    f_ek = 2 * Omega * np.sin(np.radians(lat_ek))
    DE_ek = np.sqrt(2 * Az_ex / abs(f_ek))
    ME_ek = tau_ex / (rho_w * abs(f_ek))
    print(f"  {lat_ek:10} {f_ek:12.4e} {DE_ek:10.1f} {ME_ek:12.4f}")

print()
print()
print("[Chart saved: Ekman spiral, Sverdrup transport, THC hysteresis, T-S diagram]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Wind-Driven Gyre Dynamics and Deep Water Properties

Python

script.py71 lines

import numpy as np

print("OCEAN GYRE DYNAMICS AND WATER MASS ANALYSIS")
print("=" * 70)

Omega = 7.292e-5
a = 6.371e6
rho = 1025.0

# Stommel boundary layer
print()
print("1. Western Boundary Current Width (Stommel vs Munk)")
print("-" * 60)

print(f"  {'Latitude':12} {'beta (1/ms)':14} {'delta_S (km)':14} {'delta_M (km)':14}")
print("  " + "-" * 50)

R_fric = 1e-6  # bottom friction (1/s)
A_H = 1e4      # lateral viscosity (m^2/s)

for lat in [15, 20, 25, 30, 35, 40, 45]:
    beta = 2 * Omega * np.cos(np.radians(lat)) / a
    delta_S = R_fric / beta / 1000  # Stommel (km)
    delta_M = (A_H / beta)**(1.0/3.0) / 1000  # Munk (km)
    print(f"  {lat:10} {beta:12.2e} {delta_S:12.1f} {delta_M:12.1f}")

# Deep water properties
print()
print()
print("2. Deep Water Mass Properties")
print("-" * 65)

water_masses = [
    ("NADW", 2.0, 34.90, 1.0, 2000, "N. Atlantic"),
    ("AABW", -0.5, 34.66, 5.0, 4000, "Weddell Sea"),
    ("AAIW", 3.5, 34.20, 0.8, 900, "Southern Ocean"),
    ("CDW", 1.5, 34.72, 2.0, 2500, "Circumpolar"),
    ("NPDW", 1.5, 34.68, 1.5, 3000, "N. Pacific"),
    ("Med Water", 11.0, 36.50, 0.5, 1000, "Mediterranean"),
]

print(f"  {'Water Mass':14} {'T (C)':10} {'S (psu)':10} {'O2 (mL/L)':12} {'Depth (m)':12} {'Source'}")
print("  " + "-" * 65)

for name, T, S, O2, depth, source in water_masses:
    print(f"  {name:12} {T:8.1f} {S:8.2f} {O2:10.1f} {depth:10} {source}")

# Ventilation ages
print()
print()
print("3. Deep Ocean Ventilation Ages (from 14C)")
print("-" * 55)

basins = [
    ("N. Atlantic deep", 200),
    ("S. Atlantic deep", 400),
    ("Indian deep", 800),
    ("N. Pacific deep", 1500),
    ("S. Pacific deep", 1200),
]

print(f"  {'Basin':25} {'Ventilation age (yr)':22} {'Turnover'}")
print("  " + "-" * 55)

for name, age in basins:
    turn = "Fast" if age < 300 else ("Moderate" if age < 800 else "Slow")
    print(f"  {name:23} {age:20} {turn}")

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

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Atmospheric Physics Next: Paleoclimatology →