Dynamics II
Lecture: April 27, 2026 (Monday), 14:00 Prof. Dr. Gerrit Lohmann
Preparation
trailer Cellules de Bénard (1min)
Rayleigh–Bénard convection: cooking oil and small aluminium particles (5 min)
Simulations:
Rayleigh Benard Thermal Convection with LBM (5 min)
Rayleigh Benard
Thermal Convection 3D Simulation (2 min)
Content for today
April 13, 14:00: Lecture (G. Lohmann)
Content for today: Fluid Dynamics, Non-dimensional parameters, Dynamical similarity, Elimination of the pressure term and vorticity, bifurcations, applications
Non-dimensional parameters: The Reynolds number
For the case of an incompressible flow in the Navier-Stokes equations, assuming the temperature effects are negligible and external forces are neglected.
conservation of mass \[ \nabla \cdot \mathbf{u} = 0 \] conservation of momentum \[ \partial_t \mathbf{u} + ( \mathbf{u} \cdot \nabla) \mathbf{u} = - \frac{1}{\rho_0} \nabla p + \nu \nabla^2 \mathbf{u} \]
The equations can be made dimensionless by a length-scale L, determined by the geometry of the flow, and by a characteristic velocity U.
For analytical solutions, numerical results, and experimental measurements, it is useful to report the results in a dimensionless system (concept of dynamic similarity).
Goal: replace physical parameters with dimensionless numbers, which completely determine the dynamical behavior
representative values for velocity \((U),\) time \((T),\) distances \((L)\)
Using these values, the values in the dimensionless-system (written with subscript d) can be defined: \[ u =U \cdot u_d \] \[ t = T \cdot t_d \] \[ x = L \cdot x_d \] with \(U = L/T\).
From these scalings, we can also derive \[ \partial_t = \frac{\partial}{\partial t } = \frac{1}{T} \cdot \frac{\partial}{\partial t_d } \] \[ \partial_x = \frac{\partial}{\partial x} = \frac{1}{L} \cdot \frac{\partial}{\partial x_d } \] Note furthermore the units of \([\rho_0] = kg/m^3\), \([p] = kg/(m s^2)\), and \([p]/[\rho_0]= m^2/s^2\). Therefore the pressure gradient term has the scaling \(U^2/L\).
\[ \nabla_d \cdot \mathbf{u_d} = 0 \] and conservation of momentum
…
The dimensionless parameter \[ Re=UL/ \nu \]
is the Reynolds number and the only parameter left!
For large Reynolds numbers, the flow is turbulent.
In most practical flows \(Re\) is rather large \((10^4-10^8),\) large enough for the flow to be turbulent. A large Reynolds number allows the flow to develop steep gradients locally.
In the literature, the term “equations have been made dimensionless”, means that this procedure is applied and the subscripts d are dropped.
Characterising flows by dimensionless numbers
These numbers can be interpreted as follows:
- Re: (stationary inertial forces)/(viscous forces)
- Sr: (non-stationary inertial forces)/(stationary inertial forces)
- Fr: (stationary inertial forces)/(gravity)
- Fo: (heat conductance)/(non-stationary change in enthalpy)
- Pe: (convective heat transport)/(heat conductance)
- Ec: (viscous dissipation)/(convective heat transport)
- Ma: (velocity)/(speed of sound): objects moving faster than 0.8 produce shockwaves
- Pr and Nu are related to specific materials.
Rayleigh-Benard convection (and Rayleigh number)
Experiments:
Rayleigh–Bénard convection: cooking oil and small aluminium particles (5 min),
Rayleigh Benard
Thermal Convection 3D Simulation (2 min)
Convection in the Rayleigh-Benard system
Rayleigh (1916) temperature difference between the upper- and lower-surfaces \[ T(x, y, z=H) = \, T_0 \] \[ T(x, y, z=0) \, = \, T_0 + \Delta T \]
Furthermore \[ \rho = \rho_0 = const. \] except in the buoyancy term, where:
\[ \varrho = \varrho_0 (1 - \alpha(T-T_0)) \mbox{ with } \alpha > 0 \quad . \]
common feature of geophysical flows
No Convection Equilibrium: Diffusion
Rayleigh–Bénard System (Streamfunction Formulation)
We consider the 2D Boussinesq system
\[ D_t u = -\frac{1}{\rho_0} \partial_x p + \nu \nabla^2 u \]
\[ D_t w = -\frac{1}{\rho_0} \partial_z p + \nu \nabla^2 w + g \alpha \Theta \]
\[ \partial_x u + \partial_z w = 0 \]
\[ D_t T = \kappa \nabla^2 T \]
Equilibrium (pure diffusion)
\[ u = w = 0 \]
\[ T_{eq}(z) = T_0 + \left(1 - \frac{z}{H}\right)\Delta T \]
We introduce the temperature perturbation
\[ T = T_{eq}(z) + \Theta \]
Streamfunction
Using incompressibility:
\[ u = -\partial_z \Psi, \qquad w = \partial_x \Psi \]
and the vorticity:
\[ \omega = \nabla^2 \Psi \]
Vorticity equation
Eliminating pressure gives
\[ D_t (\nabla^2 \Psi) = \nu \nabla^4 \Psi + g \alpha \frac{\partial \Theta}{\partial x} \]
Temperature equation
From
\[ D_t T = \kappa \nabla^2 T \]
we obtain
\[ D_t \Theta = \frac{\Delta T}{H} \frac{\partial \Psi}{\partial x} + \kappa \nabla^2 \Theta \]
Non-dimensionalization
We introduce
\[ x = H x', \quad z = H z', \quad t = \frac{H^2}{\kappa} t' \]
\[ \Psi = \kappa \Psi', \quad \Theta = \Delta T \, \Theta' \]
Dimensionless equations
This yields
\[ D_t (\nabla^2 \Psi) = \sigma \nabla^4 \Psi + \sigma Ra \frac{\partial \Theta}{\partial x} \]
\[ D_t \Theta = \frac{\partial \Psi}{\partial x} + \nabla^2 \Theta \]
Dimensionless parameters
\[ Ra = \frac{g \alpha \Delta T H^3}{\nu \kappa} \qquad \sigma = \frac{\nu}{\kappa} \] ## Galerkin Approximation: Low-Order Model
To derive a reduced model of convection, we expand the fields in Fourier modes following Saltzman (1962).
Spectral expansion
The streamfunction and temperature perturbation are written as:
\[ \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) \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) \sin\left(\frac{l \pi}{H} z\right) \]
Low-order truncation
Keeping only the dominant modes leads to a 3-variable system:
\[ \Psi(x,z,t) = \frac{1+a^2}{a\,\kappa} \, X(t)\,\sqrt{2} \sin\left(\frac{\pi a}{H} x\right) \sin\left(\frac{\pi}{H} z\right) \]
\[ \Theta(x,z,t) = \frac{\Delta T}{\pi} \frac{R_c}{R_a} \left[ Y(t)\,\sqrt{2} \cos\left(\frac{\pi a}{H} x\right) \sin\left(\frac{\pi}{H} z\right) - Z(t)\,\sin\left(\frac{2\pi}{H} z\right) \right] \]
This yields a reduced dynamical system in terms of three amplitudes: \[ X(t), \quad Y(t), \quad Z(t) \]
These represent: - \(X\):
circulation strength
- \(Y\): horizontal temperature
variation
- \(Z\): vertical temperature
stratification
Rayleigh Number and Onset of Convection
The key control parameter is the Rayleigh number:
\[ R_a = \frac{g \alpha \Delta T H^3}{\nu \kappa} \]
It measures the competition between: - buoyancy forcing (destabilizing) - diffusion and viscosity (stabilizing)
Critical Rayleigh number
Convection sets in when
\[ R_a > R_c \]
with
\[ R_c = \pi^4 \frac{(1+a^2)^3}{a^2} \]
The minimum value occurs for
\[ a^2 = \frac{1}{2} \]
giving
\[ R_c \approx 657.5 \]
Physical interpretation
For \(R_a < R_c\):
heat transport is purely diffusive (conduction)For \(R_a > R_c\):
convective motion develops
Interpretation of the Galerkin model
The Galerkin truncation reduces the full PDE system to a small set of
ODEs.
This is the basis of the classical Lorenz model of convection.
It shows that even simple low-order systems can capture:
instability of the conductive state
onset of convection
nonlinear dynamics and chaos
Lorenz system
Bifurcation at \[ r = R_a/R_c = 1\]
Geometry constant \[b = 4(1+a^2)^{-1}\]
Famous low-order model:
\[ \dot X = -\sigma X + \sigma Y \]
\[ \dot Y = r X - Y - X Z \]
\[ \dot Z = -b Z + X Y \]
\[\mbox{dimensionless time } \quad t_d = \pi^2 H^{-2} (1+a^2) \kappa t,\]
\[ \mbox{ Prandtl number } \quad \sigma = \nu \kappa^{-1}, \]
Lorenz system r=24
left: 52%
r=24
s=10
b=8/3
dt=0.01
x=0.1
y=0.1
z=0.1
vx<-c(0)
vy<-c(0)
vz<-c(0)
for(i in 1:10000){
x1=x+s*(y-x)*dt
y1=y+(r*x-y-x*z)*dt
z1=z+(x*y-b*z)*dt
vx[i]=x1
vy[i]=y1
vz[i]=z1
x=x1
y=y1
z=z1}
plot(vx,vy,type="l",xlab="x",ylab="y")
plot(vy,vz,type="l",xlab="y",ylab="z")
Lorenz system r=0.9
r=0.9
s=10
b=8/3
dt=0.01
x=1.1
y=0.1
z=11.1
vx<-c(0)
vy<-c(0)
vz<-c(0)
for(i in 1:100){
x1=x+s*(y-x)*dt
y1=y+(r*x-y-x*z)*dt
z1=z+(x*y-b*z)*dt
vx[i]=x1
vy[i]=y1
vz[i]=z1
x=x1
y=y1
z=z1}
plot(vx,type="l",xlab="time",ylab="x")
plot(vy,type="l",xlab="time",ylab="y")
Lorenz system r=3.5
r=3.5
s=10
b=8/3
dt=0.01
x=1.1
y=0.1
z=11.1
vx<-c(0)
vy<-c(0)
vz<-c(0)
for(i in 1:1000){
x1=x+s*(y-x)*dt
y1=y+(r*x-y-x*z)*dt
z1=z+(x*y-b*z)*dt
vx[i]=x1
vy[i]=y1
vz[i]=z1
x=x1
y=y1
z=z1}
plot(vx,type="l",xlab="time",ylab="x")
plot(vy,type="l",xlab="time",ylab="y")
North Atlantic Current & Gulfstream
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.
How much heat is transported northward ?
Calculation:
\[ \underbrace{ c_p}_{4.2 \cdot 10^3 Ws/(m^3 kg)} \, \cdot \, \underbrace{ \rho }_{10^3 kg/m^3 } \, \cdot \, \underbrace{\Phi}_{15 \cdot 10^6 m^3/s} \, \cdot \, \underbrace{\Delta T}_{(18-2) K } = 1 \cdot 10^{15} W \]
The flow carried by the conveyor belt loses 1 Petawatts (PW), close to estimates of Rintoul and Wunsch (1991)
The deep bottom water from the North Atlantic is mixed upward in other regions and ocean, and it makes its way back to the Gulf Stream and the North Atlantic. Thus most of the water that sinks in the North Atlantic must be replaced by water from the far South Atlantic and Pacific Ocean.
Ocean Conveyor Belt
Conveyor Belt: Industry
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
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}. \]
# Conveyor belt: climate change
left: 55%
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
Vorticity Dynamics of the Meridional Overturning (y,z)
We consider the meridional–vertical (y,z) plane.
Momentum equations
\[ \frac{\partial v}{\partial t} = - \frac{1}{\rho_0} \frac{\partial p}{\partial y} - f u - \kappa v \]
\[ \frac{\partial w}{\partial t} = - \frac{1}{\rho_0} \frac{\partial p}{\partial z} - \frac{g}{\rho_0} (\rho - \rho_0) - \kappa w \]
where \(\kappa\) represents Rayleigh friction.
Continuity and streamfunction
The continuity equation is
\[ \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 \]
We introduce a streamfunction \(\Phi(y,z,t)\):
\[ v = \partial_z \Phi, \qquad w = -\partial_y \Phi \]
Vorticity equation
Eliminating pressure yields
\[ \frac{\partial}{\partial t} \nabla^2 \Phi = - f \frac{\partial u}{\partial z} + \frac{g}{\rho_0} \frac{\partial \rho}{\partial y} - \kappa \nabla^2 \Phi \]
The Coriolis term is typically absorbed into the damping term for simplicity.
Galerkin Approximation
We expand the streamfunction as
\[ \Phi(y,z,t) = \sum_{k=1}^{\infty}\sum_{l=1}^{\infty} \Phi_{k,l}(t) \sin\left(\frac{\pi k y}{L}\right) \sin\left(\frac{\pi l z}{H}\right) \]
Single-mode truncation
We retain only the dominant mode:
\[ \Phi(y,z,t) = \Phi_{max}(t) \sin\left(\frac{\pi y}{L}\right) \sin\left(\frac{\pi z}{H}\right) \]
This satisfies no-normal-flow boundary conditions.
Projection onto the dominant mode
We project the vorticity equation onto the chosen basis function:
\[ \int_0^L dy \int_0^H dz \, (\cdot) \]
Time tendency term
\[ \int \frac{\partial}{\partial t} \nabla^2 \Phi = 4 L H \left(\frac{1}{L^2} + \frac{1}{H^2}\right) \frac{d \Phi_{max}}{dt} \]
Density forcing term
\[ \int \frac{g}{\rho_0} \frac{\partial \rho}{\partial y} = \frac{g}{\rho_0} H (\rho_{north} - \rho_{south}) \]
Friction term
\[ \int \kappa \nabla^2 \Phi = 4 L H \kappa \left(\frac{1}{L^2} + \frac{1}{H^2}\right) \Phi_{max} \]
Final low-order model
Combining all terms gives
\[ \frac{d}{dt} \Phi_{max} = \frac{a}{\rho_0} (\rho_{north} - \rho_{south}) - \kappa \Phi_{max} \]
with
\[ a = \frac{g L H^2}{4 (L^2 + H^2)} \]
Diagnostic limit
Assuming fast adjustment (e.g. Kelvin waves), we obtain
\[ \Phi_{max} = \frac{a}{\rho_0 \kappa} (\rho_{north} - \rho_{south}) \]
Connection to Box Models
The density difference can be written as
\[ \rho_{north} - \rho_{south} = - \rho_0 (\alpha \Delta T - \beta \Delta S) \]
leading to the classical Stommel-type relation
\[ \Phi = - c (\alpha \Delta T - \beta \Delta S) \]
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
Box model of the thermohaline circulation
hemispheric (Stommel-type) or interhemispheric (Rooth-type)
Schematic picture of the hemispheric two box model (a) and of the interhemispheric box model
- The Atlantic surface density is mainly related to temperature
differences.
- But the pole-to-pole differences are caused by salinity differences. }
Box model of the thermohaline circulation
hemispheric (Stommel-type) or interhemispheric (Rooth-type)
\[ \Phi \, = \, - c \, ( \alpha \Delta T - \beta \Delta S ) \qquad , \] where \(\alpha\) and \(\beta\) are the thermal and haline expansion coefficients,
\(\Delta\) denotes the meridional difference operator applied to temperature T and salinity S, respectively.
single hemispheric view:
The meridional density differences are clearly dominated by temperature
Salinity difference breakes the temperature difference
both hemispheres:
densities at high northern and southern latitudes are close,
the pole-to-pole differences are caused by salinity differences.
Ocean: multiple equilibria in GCM
Ocean: multiple equilibria in GCM
Ocean response: Hosing the North Atlantic
- 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%
halting or reversing the ocean circulation
interpretation of Greenland ice core records
climate states with different ocean modes
Abrupt climate change, termination of ice sheets, Climate System II
Ocean box model
Ocean: Temperature
Euler forward for ocean temperature
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
Box model
3) Boxmodel (Version using R-shiny)
How to setup and run the Boxmodel R-shiny app
password: DynamicsII2020
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
Download link password: DynamicsII2020
go to the folder where the boxmodel is located (ui.R): setwd(’ … ’)
install.packages(“shiny”)
library(shiny)
runApp()
a browser window should open that displays the boxmodel controls
first click on the button Spin-Up to create a new spin up
a new window opens, click on Run SpinUp
afterwards click on Apply, then Cancel
define parameters for the model simulation
run the simulation by clicking on Start Simulation