Course: Dynamics II

Linear stability analysis

Consider the continuous dynamical system described by \[\dot x=f(x,\lambda)\quad \] A bifurcation occurs at \[(x_E,\lambda_0)\] if the Jacobian matrix \[ \textrm{d}f/dx (x_E,\lambda_0)\] has an Eigenvalue with zero real part.

Galerkin approximation: Get a low-order model

Scheme

\[ \mbox{ Saltzman (1962): Expand } \Psi, \Theta \mbox{ in double Fourier series in x and z: } \]

\[\Psi (x,z,t) \, = \, \sum_{k=1}^\infty \sum_{l=1}^\infty \Psi_{k,l} (t) \, \, \sin \left(\frac{k \pi a}{H} x \right) \, \times \, \sin \left(\frac{ l \pi}{H} z \right) \] \[\Theta (x,z,t) \, = \, \sum_{k=1}^\infty \sum_{l=1}^\infty \Theta_{k,l} (t) \cos \left(\frac{k \pi a}{H} x \right) \, \times \, \sin \left( \frac{l \pi}{H} z \right) \]

From partial differential equations to ordinary differential equations for \(\Psi_{k,l}(t)\) and \(\Theta_{k,l}(t)\).

Approximation: Just 3 Modes X(t), Y(t), Z(t)

\[ \frac{a}{1+a^2} \, \kappa \, \Psi = X \sqrt{2} \sin\left(\frac{\pi a}{H} x \right) \sin\left(\frac{\pi}{H} z \right) \]

\[ \pi \frac{R_a}{R_c} \frac{1}{\Delta T} \, \Theta = Y \sqrt{2} \cos\left(\frac{\pi a}{H} x\right) \sin\left(\frac{\pi}{H} z \right) - Z \sin\left(2 \frac{\pi}{H} z \right) \]

Rayleigh number Ra: Buoyancy & Viscosity

\[ \mbox{Motion develops if } \quad R_a = \frac{g \alpha H^3 \Delta T}{\nu \kappa} \quad \mbox{exceeds a critical } \quad R_c \]

As the Rayleigh number increases, the gravitational force becomes more dominant. The critical Rayleigh number can be obtained analytically for a number of different boundary conditions by doing a perturbation analysis on the linearized equations in the stable state.

\(R_c = \pi^4 \frac{\left(1+a^2\right)^3}{a^2} = \pi^4 \frac{27}{4}= 657.51\) with \(a^2 = 1/2\)

\[\mbox{When } R_a < R_c,\mbox{ heat transfer is due to conduction} \]

\[\mbox{When } R_a > R_c, \mbox{ heat transfer is due to convection.} \]

Bernard Cells

Lorenz system

Bifurcation at \[ r = R_a/R_c = 1\]

 
Famous low-order model:

\[ \dot X = -\sigma X + \sigma Y \]

\[ \dot Y = r X - Y - X Z \]

\[ \dot Z = -b Z + X Y \]

North Atlantic Current & Gulfstream

Gulf Stream & North Atlantic Current

Part of deep ocean

brings warm water northward where it cools.

returns southward as a cold, deep, western-boundary current.

Gulf Stream carries 40 Sv of 18°C water northward.

Of this, 15 Sv return southward in the deep western boundary current at a temperature of 2°C.

Ocean Conveyor Belt

Conveyor

Conveyor Belt: Industry

Conveyor

Conveyor belt circulation

The the conveyor is driven by deepwater formation in the northern North Atlantic.

The conveyor belt metaphor necessarily simplifies the ocean system, it is of course not a full description of the deep ocean circulation.

Broecker’s concept provides a successful approach for global ocean circulation, although several features can be wrong like the missing Antarctic bottom water, the upwelling areas etc..

metaphor inspired new ideas of halting or reversing the ocean circulation and put it into a global climate context.

interpretation of Greenland ice core records indicating different climate states with different ocean modes of operation (like on and off states of a mechanical maschine).

Thermohaline ocean circulation

Overturning

Modelled meridional overturning streamfunction in Sv 10^6 = m^3 /s in the Atlantic Ocean. Grey areas represent zonally integrated smoothed bathymetry

Estimates of overturning ?

It is observed that water sinks in to the deep ocean in polar regions of the Atlantic basin at a rate of 15 Sv. (Atlantic basin: 80,000,000 km^2 area * 4 km depth.)

– How long would it take to ‘fill up’ the Atlantic basin?

– Supposing that the local sinking is balanced by large-scale upwelling, estimate the strength of this upwelling.

Hint: Upwelling = area * w

– Compare this number with that of the Ekman pumping!

Estimates of overturning: Solution

Timescale T to ‘fill up’ the Atlantic basin: \[ T = \frac{ 80 \cdot 10^{12} \, m^2 \cdot 4000 \, m}{15 \cdot 10^6 \, m^3 s^{-1}} = 2.13 \cdot 10^{10} s = 676 \;years\]

Overturning is balanced by large-scale upwelling:

\[ area \cdot w = 15 \cdot 10^6 \, m^3 s^{-1}\]

\[ w = 0.1875 \cdot 10^{-6} m\;s^{-1} = 5.9 \cdot 10^{-15} m \; y^{-1}. \]

Ekman pumping \[ w_E \simeq 32 \, \, m \; y^{-1}. \]

Stability of the 0-D Energy Balance Model

The zero-dimensional energy balance model (EBM) with ice-albedo feedback is given by: \[ C \frac{dT}{dt} = \frac{S_0}{4} (1 - \alpha(T)) - \epsilon \sigma T^4 \]

Where:

• \(T\): Global mean surface temperature
  
• \(C\): Effective heat capacity of the climate system
  
• \(S_0 \approx 1361 \, \text{W/m}^2\): Solar constant
  
• \(\alpha(T)\): Planetary albedo, dependent on temperature
  
• \(\epsilon \approx 0.61\): Effective emissivity
  
• \(\sigma = 5.67 \times 10^{-8} \, \text{W/m}^2\text{K}^4\): Stefan–Boltzmann constant

Ice-Albedo Feedback

The albedo is parameterized as: \[ \alpha(T) = \alpha_{\text{ice}} + \frac{\alpha_{\text{water}} - \alpha_{\text{ice}}}{1 + \exp\left( \frac{T - T_c}{\Delta T} \right)} \] with typical values:

• \(\alpha_{\text{ice}} = 0.6\): Albedo of ice
  
• \(\alpha_{\text{water}} = 0.3\): Albedo of open water
  
• \(T_c = 273\,\text{K}\): Critical temperature
  
• \(\Delta T = 5\,\text{K}\): Transition width

Equilibrium and Stability

At equilibrium: \[ \frac{S_0}{4} (1 - \alpha(T)) = \epsilon \sigma T^4 \]

The function: \[ f(T) = \frac{S_0}{4}(1 - \alpha(T)) - \epsilon \sigma T^4 \]

Stability is determined by the sign of the derivative: \[ f'(T) = -\frac{S_0}{4} \alpha'(T) - 4\epsilon \sigma T^3 \]

If \(f'(T_e) < 0\), the equilibrium at \(T_e\) is stable; if \(f'(T_e) > 0\), it is unstable.

# Load necessary libraries
library(ggplot2)
library(gridExtra)

# Physical constants
S0 <- 1361                  # Solar constant [W/m^2]
epsilon <- 0.61             # Emissivity
sigma <- 5.67e-8            # Stefan–Boltzmann constant [W/m^2/K^4]

# Albedo model parameters
alpha_ice <- 0.7
alpha_water <- 0.3
T_c <- 260                  # Transition temperature [K]
delta_T <- 5                # Transition width [K]

# Albedo and its derivative
alpha <- function(T) {
  alpha_water + (alpha_ice - alpha_water) / (1 + exp((T - T_c) / delta_T))
}
alpha_prime <- function(T) {
  expT <- exp((T - T_c) / delta_T)
 - (alpha_ice - alpha_water) * expT / ( delta_T * (1 + expT)^2)
}
onemalpha_int <- function(T) {
  expT <- exp((T - T_c) / delta_T)
 (1-alpha_water) * T - (alpha_ice - alpha_water) * delta_T * log( expT/(1 + expT))
}

# Energy balance and its derivative
f_T <- function(T) {
  (S0 / 4) * (1 - alpha(T)) - epsilon * sigma * T^4
}
f_prime <- function(T) {
  - (S0 / 4) * alpha_prime(T) - 4 * epsilon * sigma * T^3
}
f_int <- function(T) {
   -(S0 / 4) * onemalpha_int(T) +  epsilon/5 * sigma * T^5
}

# Temperature range and values
T_vals <- seq(220, 310, by = 0.1)
f_vals <- sapply(T_vals, f_T)
f_prime_vals <- sapply(T_vals, f_prime)
f_int_vals <- sapply(T_vals, f_int)

# Find approximate equilibria (sign changes of f)
signs <- sign(f_vals)
eq_indices <- which(diff(signs) != 0)
T_eq <- T_vals[eq_indices]
f_eq <- sapply(T_eq, f_T)
f_prime_eq <- sapply(T_eq, f_prime)
stability <- ifelse(f_prime_eq < 0, "Stable", "Unstable")

# Create data frames
df_f <- data.frame(T = T_vals, f = f_vals)
df_fp <- data.frame(T = T_vals, f_prime = f_prime_vals)
df_fint <- data.frame(T = T_vals, f_int = f_int_vals)
df_eq <- data.frame(T = T_eq, f = f_eq, f_prime = f_prime_eq, Stability = stability)

# Plot f(T)
plot_f <- ggplot(df_f, aes(x = T, y = f)) +
  geom_line(color = "green", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(data = df_eq, aes(x = T, y = f, color = Stability), size = 3) +
  scale_color_manual(values = c("Stable" = "blue", "Unstable" = "orange")) +
  labs(title = "Net Energy Flux f(T)", x = "Temperature [K]", y = expression(f(T)~"[W/m"^2*"]")) +
  theme_minimal()

# Plot f'(T)
plot_fp <- ggplot(df_fp, aes(x = T, y = f_prime)) +
  geom_line(color = "red", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(data = df_eq, aes(x = T, y = f_prime, color = Stability), size = 3) +
  scale_color_manual(values = c("Stable" = "blue", "Unstable" = "orange")) +
  labs(title = "Derivative f'(T)", x = "Temperature [K]", y = expression(f*"'(T)~[W/m"^2*"/K]")) +
  theme_minimal()

# Plot integral of f(T)
plot_fint <- ggplot(df_fint, aes(x = T, y = f_int)) +
  geom_line(color = "red", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(data = df_eq, aes(x = T, y = f_prime, color = Stability), size = 3) +
  scale_color_manual(values = c("Stable" = "blue", "Unstable" = "orange")) +
  labs(title = "Integral f(T)", x = "Temperature [K]", y = expression(V(T)~"[W/m"^2*"K]")) +
  theme_minimal()

# Combine the two plots
grid.arrange(plot_f, plot_fp, ncol = 2)

Scaling: Rotating frame of reference

The Coriolis effect is one of the dominating forces for the large-scale dynamics.

\[ \underbrace{\frac{\partial \mathbf{v}}{\partial t}}_{ U/T \sim 10^{-8} } + \underbrace{\mathbf{v} \cdot \nabla \mathbf{v}}_{ U^2/L \sim 10^{-8} }= {\underbrace{- \frac{1}{\rho} \nabla p}_{ \bf \delta P/(\rho L) \sim 10^{-5} } + \underbrace{2 \mathbf{\Omega \times v}}_{ \bf f_0 U \sim 10^{-5} } + \underbrace{fric}_{ \nu U/H^2 \sim 10^{-13}}} \quad \]

where fric is due to eddy stress divergence (\(\sim \nu \nabla^2 \mathbf{v}\)).

Values from the ocean -> exercise.

Because of the continuity equation \[ U/L \sim W/H \] and horizontal scales are orders of magnitude larger than the vertical ones, \[ W << U . \]

The timescales are related to \[ T \sim L/U \sim H/W \]

Exception for small scales (e.g., ocean convection or cumuls clouds): \[ H \sim L \quad \rightarrow \quad W \sim U \]

Rossby number Ro

\[ \underbrace{\frac{\partial \mathbf{v}}{\partial t}}_{ U/T \sim 10^{-8} } + \underbrace{\mathbf{v} \cdot \nabla \mathbf{v}}_{ U^2/L \sim 10^{-8} } = {\underbrace{- \frac{1}{\rho} \nabla p}_{ \bf \delta P/(\rho L) \sim 10^{-5}} + \underbrace{2 \mathbf{\Omega \times v}}_{ \bf f_0 U \sim 10^{-5} } + \underbrace{fric}_{ \nu U/H^2 \sim 10^{-13}}} \]

\[ Ro = \frac{ \mbox{Inertial (the left hand side)} }{ \mbox{Coriolis term } } \]

\[ Ro = \frac{(U^2/L)}{(f U)} = \frac{U}{f L} \quad \]

characterizes the importance of Coriolis acceleration

Vorticity is the rotation of the fluid

\[ \zeta \equiv \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \]

or in 3D:

\[ \equiv \nabla \times \boldsymbol{ u} \]

Example: Vorticity from shear

Cape

Tomczak & Godfrey: Regional Oceanography

\[u=0, v=v\left(x\right)\]

\[\zeta= \partial v\left(x\right)/\partial x \]

Planetary and relative vorticity

\[ \mbox{Absolute Vorticity }\equiv\left(\zeta+f\right) \]

\[ \frac{Du}{Dt}-f\;v = -\frac{1}{\rho}\frac{\partial p}{\partial x} \] \[ \frac{Dv}{Dt}+f\;u = -\frac{1}{\rho}\frac{\partial p}{\partial y} \]

\[ \mbox{subtract } \partial/\partial y \mbox{ of (u-equation) from } \partial /\partial x \mbox{ of (v-equation) } \]

Use \[ \frac{D}{Dt} f = v \, \partial_y f: \]

to obtain

\[ \underline{ \frac{D}{Dt}\left(\zeta+f\right) + \left(\zeta + f\right)\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0 } \quad \]

Examples for Vorticity: Ocean with depth h(x,y)

\[ \mbox{Because of the continuity equation } \quad \partial_x \left( u h \right) + \partial_y \left( v h \right) \quad = \quad 0 \]

\[ \quad \frac{D}{Dt} h + h \left( \partial_x u + \partial_y v \right) = \quad 0 \]

Therefore, \[ \underline{ \frac{D}{Dt}\left(\zeta+f\right) + \left(\zeta + f\right)\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0 } \quad \]

\[ \frac{D}{Dt}\left(\zeta +f\right)-\frac{\left(\zeta+f\right)}{h}\frac{Dh}{Dt}=0 \]

\[ \frac{1}{h} \frac{D}{Dt}\left(\zeta+f\right) - \left(\zeta + f\right) \frac{D_t h}{h^2} =0 \]

\[ \underline{ \frac{D}{Dt}\left( \frac{\zeta+f}{h}\right) = 0 } \quad \]

Potential vorticity is conserved along a fluid trajectory.

Potential vorticity: Examples

\[ \frac{D}{Dt}\left(\zeta+f\right) + \left(\zeta + f\right)\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0 \quad \]

Ocean/Atmosphere with depth h(x,y)

\[ \frac{D}{Dt}\left( \frac{\zeta+f}{h}\right) = 0 \quad \]

 

Couples depth, vorticity, latitude

– Changes in the depth results in change in \(\zeta\).

– Changes in latitude require a corresponding change in \(\zeta\).

Cape

Dietrich et al. (1980)

---

Cape

Steward, Oceanography

 

Angular Momentum and Hadley Cell

Hadely Cell

Tropical air rises to tropopause & moves poleward

Deflected eastward by the Coriolis force

Subtropical jet: forms at poleward limit of Hadley Cell

It tends to conserve angular momentum, friction small

equatorward moving air: westward component

What drives the ocean currents?

Friction: transfer of momentum from atmosphere to oceanic Ekman layer

Vorticity dynamics for the ocean and include the wind stress term

\[ D_t u - f v = - \frac{1}{\rho} \frac{\partial p}{\partial x} + \frac{1}{\rho} \partial_z \tau_{xz} \] \[ D_t v + f u = - \frac{1}{\rho} \frac{\partial p}{\partial y} + \frac{1}{\rho} \partial_z \tau_{yz} \]

\[ \frac{D}{Dt} \left( {\zeta+f}\right) - \frac{\left(\zeta+f \right)}{h} \frac{D}{Dt} h \, = \, \frac{1}{\rho} \underbrace{\left( \frac{\partial}{\partial x} \, \partial_z \tau_{yz} - \frac{\partial}{\partial y}\, \partial_z \tau_{xz} \right)}_{curl \, \partial_z \tau} \quad . \]

\[ \frac{D}{Dt} \left( \frac{\zeta+f}{h}\right) = \frac{1}{\rho \, h} \, \mbox{curl} \, \partial_z \tau \, \]

Sverdrup transport

\[ \beta v = \frac{1}{\rho } \, \mbox{curl} \, \partial_z \tau \, \]

\[ \int_{-H}^0 dz \, \beta v = \frac{1}{\rho } \, \int_{-H}^0 dz \, \mbox{curl} \, \partial_z \tau \, = \frac{1}{\rho } \, \mbox{curl} \, \tau \, \]

\[ V = \frac{1}{\rho \beta} \, \left( \frac{\partial \tau_{yz} }{\partial x} \, - \frac{\partial \tau_{xz}}{\partial y}\, \right) = \frac{1}{\rho \beta} \, \, \operatorname{curl} \, \tau \]

Sverdrup

applied globally using the wind stress from Hellerman and Rosenstein (1983). Contour interval is \(10\) Sverdrups (Tomczak and Godfrey, 1994).

Ekman Pumping & Sverdrup Transport

Ekman

The center of a subtropical gyre is a high pressure zone: clockwise on the Northern Hemisphere

Ekman surface currents towards the center of the gyre

The Ekman vertical velocity balanced by \[ w_E=w_g \] vertical geostrophic current in the interior

geostrophic flow towards the equator

returned flow towards the pole in western boundary currents

Ekman Pumping: vertical velocity at the bottom of the Ekman layer E

\(w_E\) as the Ekman vertical velocity the bottom of the Ekman layer \[ w_E = - \int_{-E}^0 \frac{\partial w}{\partial z} dz = \frac{\partial}{\partial x} U_E + \frac{\partial}{\partial y} V_E \]

\(\operatorname{curl} \mathbf{\tau}\) produces a divergence of the Ekman transports leading to \(w_E\) at the bottom E

\[ w_E = \, \frac{\partial }{\partial x} \left( \frac{ \tau_{y}}{\rho \;f }\, \right) - \frac{\partial }{\partial y}\, \left( \frac{ \tau_{x}}{\rho \;f }\, \right) =\operatorname{curl}\left(\frac{\mathbf{\tau}}{\rho\;f}\right) \simeq \frac{1}{\rho\;f} \, \operatorname{curl} \mathbf{\tau} \]

The order of magnitude of the Ekman vertical velocity:

typical wind stress variation of \(0.2 N m^{-2}\) per 2000 km in y-direction:

\[ w_E \simeq - \frac{ \Delta \tau_{x}}{\rho \;f_0 \Delta y}\, \simeq \frac{1 }{10^3 kg m^{-3}} \frac{0.2 N m^{-2} }{10^{-4} s^{-1}\, \, 2 \cdot 10^6 m} \simeq 32 \, \, \frac{m}{yr} \]

Conveyor belt: climate change

left: 55%

THC

halting or reversing the ocean circulation

interpretation of Greenland ice core records

climate states with different ocean modes

Overturning Circulation (model)

Modelled meridional overturning streamfunction in Sv 10^6 = m^3 /s in the Atlantic Ocean. Grey areas represent zonally integrated smoothed bathymetry

Atlantic water masses

THC

North

Estimates of overturning

Timescale T to ‘fill up’ the Atlantic basin:

\[ T = \frac{ 80 \cdot 10^{12} \, m^2 \cdot 4000 \, m}{15 \cdot 10^6 \, m^3 s^{-1}} = 2.13 \cdot 10^{10} s = 676 \;years\]

Overturning is balanced by large-scale upwelling:

\[ area \cdot w = 15 \cdot 10^6 \, m^3 s^{-1}\]

\[ w = 0.1875 \cdot 10^{-6} m\;s^{-1} = 5.9 \cdot 10^{-15} m \; y^{-1}. \]

compare to Ekman pumping \[ w_E \simeq 32 \, \, m \; y^{-1}. \]

Vorticity dynamics of meridional overturning (y,z)

\[ \frac{\partial}{\partial t} v \quad = \quad - \frac{1}{\rho_0} \frac{\partial p}{\partial y} \quad - \quad f u \quad - \quad \kappa v \]

\[ \frac{\partial}{\partial t} w \quad = \quad - \frac{1}{\rho_0} \frac{\partial p}{\partial z} \quad - \quad \frac{g}{\rho_0} (\rho -\rho_0) \quad - \quad \kappa w \] \(\kappa\) as parameter for Rayleigh friction.

Amplitude of overturning

\[ \frac{d}{dt} \Phi_{max} \, = \, \frac{a}{\rho_0} (\rho_{north} - \rho_{south}) \, \, - \, \, \kappa \Phi_{max} \] with \[ a = g L H^2/4(L^2 + H^2) \, \]

This shows that the overturning circulation depends on the density differences on the right and left boxes.

It is simplified to a diagnostic relation

\[ \Phi_{max} = \frac{a}{\rho_0 \, \kappa} \, \, (\rho_{north} - \rho_{south}) \quad \]

because the adjustment of \(\Phi_{max}\) is quasi-instantaneous due to adjustment processes, e.g. Kelvin waves.

North-south density gradient in an ocean basin

Primary north-south gradient in balance with an eastward geostrophic current: generates a secondary high & low pressure system, northward current

THC1

=======================================================

THC2

Schematic picture of the hemispheric two box model (a) and of the interhemispheric box model

=======================================================

THC3

1. The Atlantic surface density is mainly related to temperature differences.
   2. But the pole-to-pole differences are caused by salinity differences. }

Box model of the thermohaline circulation

Box model of the thermohaline circulation

Ocean: multiple equilibria in GCM

THC4

Ocean: multiple equilibria in GCM

THC5

Ocean response: Hosing the North Atlantic

THC6

1. AMOC indices of North Atlantic hosing for different hosing areas. Units are Sv. Black line represents the unperturbed LGM experiment. Hosing is for the period 840–990. (b) Annual mean sea surface salinity anomaly between LGM and the perturbation experiment LGM with 0.2 Sv for the model years 900–950.

-> Multi-scale Ocean GCM

Conveyor belt: climate change

left: 55%

THC7

halting or reversing the ocean circulation

interpretation of Greenland ice core records

climate states with different ocean modes

Deglaciation

Abrupt climate change, termination of ice sheets, Climate System II

Ocean box model

THC

Ocean: Temperature

THC

      Tnln = Tnl + dts * ((Hfnl)/(rcz2)-(Tnl-Tml)*phi/Vnl);
      
      Tmln = Tml + dts * ((Hfml)/(rcz1)-(Tml-Tsl)*phi/Vml);
      
      Tsln = Tsl + dts * ((Hfsl)/(rcz2)-(Tsl-Td)*phi/Vsl);
      
      Tdn = Td + dts * (-(Td-Tnl)*(phi/Vd));

Atmosphere

THC

Box model

Box model with jupyter

for the box model, you can download sevenbox_jupyter.ipynb

conda create -n jupyter-R

conda activate jupyter-R

conda install -y -c conda-forge pip notebook nb_conda_kernels jupyter_contrib_nbextensions

conda install -y -c conda-forge r r-irkernel r-ggplot2 r-dplyr

jupyter-notebook sevenbox_jupyter.ipynb

you have to download conda (or miniconda). Here a good web site: https://conda.io/projects/conda/en/latest/user-guide/install/index.html

Boxmodel with R-shiny