The vertical slice: buoyancy, convection and mountain waves
Atmosphere and Ocean Dynamics, Lecture 4: notes
Learning outcomes
After this lecture you can:
- compute the buoyancy frequency \(N\) of a stratified fluid, say whether it is stable, and what \(N^2 < 0\) does;
- write the equations of a vertical slice in vorticity–streamfunction form and explain why a horizontal buoyancy gradient makes vorticity: why a warm bubble rises as a mushroom and a cold one spreads as a density current;
- derive the dispersion relation of internal gravity waves, find where their energy goes, and say when the hydrostatic approximation fails;
- classify the flow of a wind over a mountain with the numbers \(Na/U\) and \(Nh_0/U\), and compute it with linear theory;
- run a vertical-slice model in vorticity–streamfunction form, with the Poisson problem that recovers the flow at every step.
1 Where we are
So far the vertical was either averaged away or balanced away. Lecture 3 averaged it: the shallow water equations describe one layer of fluid of depth \(H\), whose gravity waves travel at \(\sqrt{gH}\). Lectures 2 and 3 balanced it: in the vertical momentum equation the pressure force and gravity are ten thousand times larger than everything else, so models of large-scale flow replace that equation by hydrostatic balance, \(\partial p/\partial z = -\rho g\). Lecture 3 also introduced the Boussinesq approximation (the density is a constant \(\rho_0\) except in the buoyancy \(b = -g(\rho - \rho_0)/\rho_0\)), and lecture 1 two results we use today: Bjerknes’s circulation theorem, by which a coast heated by the sun spins up a sea breeze, and the 1950 ENIAC forecast, which stepped a vorticity equation and solved a Poisson problem at every step to recover the wind.
Air and water also move up and down. Thunderstorms lift air through the whole troposphere in minutes; their cold outflows run along the ground as gust fronts; winds blowing over mountains make waves that can reach the stratosphere, and lenticular clouds downstream (Figure 1). None of this is hydrostatic, and none of it fits in one layer. This lecture builds the simplest model that contains them all, the two-dimensional Boussinesq equations in a vertical \(x\)–\(z\) plane, a vertical slice, and uses it all morning: for buoyancy, internal waves, convection and mountain waves.
2 Layered fluids
Both fluids are stably stratified: lighter fluid lies over denser fluid. In the ocean the density increases downwards, most of it within a thermocline a few hundred metres deep, below which lies cold, nearly homogeneous abyssal water. In the atmosphere the temperature decreases upwards in the troposphere, by about 6.5 K per km, yet the atmosphere is stable, because air lifted without heating cools faster, by 9.8 K per km (lecture 3). The right measure is the potential temperature, \[ \theta = T\left(\frac{p_0}{p}\right)^{R/c_p}, \qquad p_0 = 1000\ \mathrm{hPa}, \] the temperature a parcel would have if brought adiabatically to the reference pressure \(p_0\). It is conserved by a parcel moving without heating, and it increases with height in a stable atmosphere: by a few kelvin per km in the troposphere and much faster in the stratosphere (Figure 2).
04c-vertical-modes.ipynb: buoyancy periods of about 12 minutes near the top and 2 hours in the abyss.
3 The buoyancy frequency
Lift a small parcel of fluid by \(\delta z\) from height \(z_0\), quickly enough that it does not mix with its surroundings, so that it keeps its density \(\bar\rho(z_0)\) (for the ocean we neglect compressibility; for the atmosphere the same argument works with \(\theta\)). Its new surroundings have the density \(\bar\rho(z_0 + \delta z) \approx \bar\rho(z_0) + \bar\rho'(z_0)\,\delta z\). The net upward force per unit mass, the buoyancy, is \[ b = -g\,\frac{\rho_{\mathrm{parcel}} - \rho_{\mathrm{around}}}{\rho_0} = \frac{g}{\rho_0}\frac{d\bar\rho}{dz}\,\delta z , \] with \(\rho_0\) a constant reference density (the Boussinesq approximation: density differences matter only where they are multiplied by \(g\)). Newton’s law for the parcel is then \[ \frac{d^2\,\delta z}{dt^2} = -N^2\,\delta z, \qquad N^2 = -\frac{g}{\rho_0}\frac{d\bar\rho}{dz} \qquad\left(\text{air: } N^2 = \frac{g}{\theta}\frac{d\theta}{dz}\right). \] If \(N^2 > 0\) the parcel oscillates with the buoyancy or Brunt–Väisälä frequency \(N\) (Figure 3); if \(N^2 < 0\) the displacement grows like \(e^{|N|t}\): the fluid overturns, which is convection. Typical values: \(N \approx 10^{-2}\ \mathrm{s^{-1}}\) in the troposphere and in the ocean thermocline (periods of 10–20 minutes), \(N \approx 10^{-3}\ \mathrm{s^{-1}}\) in the deep ocean (periods of hours) (Vallis 2017; Gill 1982). In the stratosphere \(N\) is about twice its tropospheric value: that is the lid under which a thunderstorm’s anvil spreads (Figure 1).
Example: how stiff is the stratification? Each value comes from one gradient, with \(g = 9.81\ \mathrm{m\,s^{-2}}\) (Figure 4). In the troposphere \(\theta\) grows by 3 K per km at \(\theta \approx 290\) K: \(N^2 = (g/\theta)\,d\theta/dz = 1.0\times10^{-4}\ \mathrm{s^{-2}}\), \(N = 0.010\ \mathrm{s^{-1}}\), a period of 10 minutes. In an ocean thermocline the density grows by 2 kg/m³ over 500 m, with \(\rho_0 = 1025\) kg/m³: \(N = 6.2\times10^{-3}\ \mathrm{s^{-1}}\), 17 minutes. A lake in summer, 24 °C at the surface and 12 °C 5 m below, with the linear equation of state \(\rho = \rho_0[1 - \beta_T(T - T_0)]\) and \(\beta_T = 2\times10^{-4}\ \mathrm{K^{-1}}\): \(N^2 = g\beta_T\,dT/dz = 4.7\times10^{-3}\ \mathrm{s^{-2}}\), \(N = 0.069\ \mathrm{s^{-1}}\), 1.5 minutes. The abyss, where the density grows by only 0.2 kg/m³ over 3000 m: \(N = 8\times10^{-4}\ \mathrm{s^{-1}}\), 2.2 hours. Lifting a parcel by \(\delta z\) costs about \(N^2\delta z^2/2\) per unit mass (the work against the restoring force \(N^2 s\) over the distance), so mixing the water vertically takes the most work in the summer lake, whose \(N^2\) is 7000 times that of the abyss: in summer the wind cannot mix the lake’s warm top into its cold bottom.
4 The vertical-slice model
4.1 The equations
Take the Boussinesq equations of lecture 3 and keep two directions only, the horizontal \(x\) and the vertical \(z\): nothing depends on \(y\) and there is no motion along it (Figure 5). Neglect the rotation of the Earth: the motions of today last minutes to hours and span a few to a few tens of kilometres, so their Rossby number is large (a bubble rises in minutes; the Coriolis parameter, \(10^{-4}\ \mathrm{s^{-1}}\), acts over hours). With \(\phi = \delta p/\rho_0\), the departure \(\delta p\) of the pressure from the hydrostatic pressure of the reference density (\(p = \tilde p(z) + \delta p\), \(d\tilde p/dz = -\rho_0 g\)) divided by \(\rho_0\), \[ \frac{Du}{Dt} = -\frac{\partial\phi}{\partial x}, \qquad \frac{Dw}{Dt} = -\frac{\partial\phi}{\partial z} + b, \qquad \frac{Db}{Dt} = 0, \qquad \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0, \] where \(D/Dt = \partial_t + u\,\partial_x + w\,\partial_z\). The buoyancy \(b\) is carried by the parcels (no heating, no diffusion); in a stratified fluid at rest it increases upwards, \(\bar b(z)\) with \(d\bar b/dz = N^2\). The vertical momentum equation keeps its acceleration \(Dw/Dt\): the model is not hydrostatic. The domain of the notebook is periodic in \(x\) and has two flat, free-slip lids, the ground and a top, where \(w = 0\).
The equations have no equation of their own for the pressure: \(\phi\) is whatever keeps the flow divergence-free. Taking the divergence of the momentum equations gives a Poisson problem for \(\phi\), which is one way to compute it (the projection method of incompressible flow). In two dimensions there is a more elegant way: eliminate the pressure altogether.
4.2 Buoyancy makes vorticity
The vorticity of the slice is \[ \xi = \frac{\partial w}{\partial x} - \frac{\partial u}{\partial z}, \] positive when the fluid turns anticlockwise in the usual picture, \(x\) to the right and \(z\) up (it is minus the \(y\) component of \(\nabla\times\mathbf u\), since \(y\) points into the page). Differentiate the \(w\) equation with respect to \(x\), the \(u\) equation with respect to \(z\), and subtract. The pressure drops out, \(-\partial_x\partial_z\phi + \partial_z\partial_x\phi = 0\), and \[ \frac{\partial\xi}{\partial t} + \frac{\partial}{\partial x}\left(u\,\frac{\partial w}{\partial x} + w\,\frac{\partial w}{\partial z}\right) - \frac{\partial}{\partial z}\left(u\,\frac{\partial u}{\partial x} + w\,\frac{\partial u}{\partial z}\right) = \frac{\partial b}{\partial x} . \] Expanding the derivatives of the products, the terms with derivatives of \(\xi\) give \(u\,\partial_x\xi + w\,\partial_z\xi\), and the remaining four combine into \((\partial_x u)(\partial_x w - \partial_z u) + (\partial_z w)(\partial_x w - \partial_z u) = \xi\,(\partial_x u + \partial_z w) = 0\), by continuity. Hence \[ \boxed{\;\frac{D\xi}{Dt} = \frac{\partial b}{\partial x}\;} \] Vorticity is carried by the flow, and it is created only where the buoyancy changes horizontally (Figure 6): where the fluid is warmer (lighter) on the right, it turns anticlockwise, rising on the warm side and sinking on the cold side. This is Bjerknes’s circulation theorem of lecture 1 for a Boussinesq fluid: the sea breeze, which rises over the warm land and sinks over the cool sea, is this equation (Bjerknes 1898; Vallis 2017). There is no stretching term: a two-dimensional flow cannot stretch its vortex lines.
4.3 The streamfunction, and the 1950 equation on its side
A divergence-free flow in a plane has a streamfunction: with \[ u = -\frac{\partial\psi}{\partial z}, \qquad w = \frac{\partial\psi}{\partial x}, \] continuity holds for any \(\psi\) (this is the convention of the horizontal two-dimensional flows of lectures 1 and 7, \(u = -\partial_y\psi\), \(v = \partial_x\psi\), with \((y, v)\) replaced by \((z, w)\)), and the vorticity is \(\xi = \partial_x w - \partial_z u = \nabla^2\psi\). With the Jacobian \(J(a, c) = \partial_x a\,\partial_z c - \partial_z a\,\partial_x c\), the advection of any quantity is \(\mathbf u\cdot\nabla c = J(\psi, c)\), and the slice model becomes \[ \frac{\partial\xi}{\partial t} + J(\psi, \xi) = \frac{\partial b}{\partial x}, \qquad \frac{\partial b}{\partial t} + J(\psi, b) = 0, \qquad \nabla^2\psi = \xi . \] Compare the barotropic vorticity equation of the 1950 forecast (lectures 1 and 7): \(\partial_t\zeta + J(\psi, \zeta + f) = 0\), \(\nabla^2\psi = \zeta\), in a horizontal plane. The vertical slice is the same equation turned on its side, with a source, the horizontal buoyancy gradient, and a second equation that carries the buoyancy. The algorithm is also the same: at each step, \(\xi\) and \(b\) are stepped forward, a Poisson problem gives \(\psi\), and its derivatives give \(u\) and \(w\) for the next step (Charney et al. 1950). Arakawa’s (1966) discretization of the Jacobian, which conserves the discrete energy and enstrophy of the 1950 equation, serves the slice as well (Arakawa 1966). A fluid at rest with \(b = \bar b(z)\) is a solution: \(\partial_x b = 0\), so no vorticity is made.
5 Internal gravity waves
Because each parcel is an oscillator, a stratified fluid carries waves in its interior. Write \(b = N^2 z + b'\) with constant \(N\), and keep only terms linear in the small quantities \(\xi\), \(\psi\), \(b'\) (the Jacobians are products of two small terms and drop out; the advection of the background, \(w\,d\bar b/dz = N^2 w\), stays): \[ \frac{\partial\xi}{\partial t} = \frac{\partial b'}{\partial x}, \qquad \frac{\partial b'}{\partial t} = -N^2 w = -N^2\,\frac{\partial\psi}{\partial x} . \] Differentiate the first in time and use the second, with \(\xi = \nabla^2\psi\): \[ \frac{\partial^2}{\partial t^2}\nabla^2\psi = -N^2\,\frac{\partial^2\psi}{\partial x^2} . \] A plane wave \(\psi \propto e^{i(kx + mz - \omega t)}\) gives \(\omega^2(k^2 + m^2) = N^2 k^2\), that is \[ \omega^2 = N^2\,\frac{k^2}{k^2 + m^2} = N^2\cos^2\alpha , \] where \(\alpha\) is the angle of the wavevector \((k, m)\) to the horizontal (written \(\vartheta\) by Vallis (2017)), and so also the angle of the crests to the vertical. These steps are those of the Rossby waves of lecture 3, where \(\partial_t\nabla^2\psi = -\beta\,\partial_x\psi\) gives \(\omega = -\beta k/(k^2 + l^2)\): there the gradient of planetary vorticity restores the fluid, here the gradient of buoyancy.
The frequency depends on the direction of the wave but not on its length, and it is never larger than \(N\): vertical crests (\(\alpha = 0\), columns of fluid moving up and down) oscillate at \(N\), as the parcel did; nearly horizontal crests oscillate slowly. The fluid moves along the crests (continuity: \(\mathbf u \perp \mathbf k\)), and so does the energy: the group velocity, \[ \mathbf c_g = \left(\frac{\partial\omega}{\partial k}, \frac{\partial\omega}{\partial m}\right) = \frac{N m}{(k^2 + m^2)^{3/2}}\,(m, -k) \quad (k > 0), \] is perpendicular to the wavevector, and its vertical component has the opposite sign of the vertical phase speed: when the crests move up, the energy goes down, and the other way round (Figure 7). A source oscillating at a frequency \(\omega < N\) therefore sends its energy along four beams at the angle \(\alpha = \arccos(\omega/N)\) from the vertical, the “St Andrew’s cross” of laboratory experiments (Mowbray and Rarity 1967); a source oscillating faster than \(N\) makes no waves at all (Figure 8). In the notebook the measured beams lie at 17.3, 30.0 and 53.3 degrees above the horizontal for \(\omega = 0.3N\), \(0.5N\), \(0.8N\), against 17.5, 30.0 and 53.1 predicted. With rotation, \(\omega^2 = N^2\cos^2\alpha + f^2\sin^2\alpha\): the frequencies of inertia–gravity waves lie between \(f\) and \(N\) (Gill 1982; Vallis 2017). Internal waves are seen from space where they change the roughness of the sea surface, as at the Strait of Gibraltar (Figure 9), or where they make clouds, as behind mountains and islands.
04a-vertical-slice.ipynb: a source of buoyancy oscillating at \(\omega = N/2\) in the middle of a fluid with constant \(N\) sends four beams of waves at \(\arccos(1/2) = 60^\circ\) from the vertical (dashed: the prediction).
5.1 Hydrostatic or not
The hydrostatic approximation drops \(Dw/Dt\) from the vertical momentum equation. In the derivation above this removes \(\partial_x w\) from the vorticity, \(\xi \approx -\partial_z u\), and the dispersion relation becomes \(\omega = Nk/m\) (Figure 10). The two agree when \(k \ll m\), that is when the motion is much wider than it is tall: the weather systems of lectures 1–3 (10 km deep, 1000 km wide, \(k/m \approx 0.01\)) and the large-scale ocean. They disagree when the motion is as tall as it is wide: the ratio of the two frequencies is \(\sqrt{1 + k^2/m^2}\), 10% at \(k/m = 0.46\), and the hydrostatic frequency grows without bound while the true one never exceeds \(N\). Convection, the waves over narrow mountains and everything a kilometre-scale model resolves are of this kind, which is why models with grids of a few kilometres solve the non-hydrostatic equations (lecture 9) (Durran 2010).
6 Convection
6.1 A warm bubble
Put a bubble of warm air (\(b > 0\)) in air at rest with \(N = 0\) (a well-mixed layer, as on a summer afternoon). Inside the bubble there is no horizontal gradient of buoyancy at its centre, but on its flanks there is: \(\partial_x b > 0\) on the left flank and \(\partial_x b < 0\) on the right. By \(D\xi/Dt = \partial_x b\), the left flank turns anticlockwise and the right one clockwise: a vortex pair, with rising air between the two vortices and sinking air outside them (Figure 11, Figure 12). The pair lifts the bubble and rolls its edges outwards and down, and the bubble becomes a mushroom; in three dimensions the pair is a vortex ring, as in a smoke ring. A bubble of radius \(R\) and buoyancy \(b\) moves at speeds of order \(\sqrt{bR}\) (its kinetic energy per unit mass, \(w^2/2\), comes from the work \(bR\) of buoyancy over its own size): for \(\theta' = 2\) K at \(\theta_0 = 300\) K, \(b = g\theta'/\theta_0 = 0.065\ \mathrm{m\,s^{-2}}\), and \(R = 2\) km, \(\sqrt{bR} \approx 11\) m/s. In the notebook’s run of this case (Figure 12) the centre of buoyancy rises from 3 to 5 km at 5.6 m/s, about half the scale, while the updraft in the core reaches 16 m/s; without diffusion the rise speed grows exactly as \(\sqrt{\theta'}\), as the scaling says. Warm bubbles of this kind were among the first tests of non-hydrostatic models with semi-implicit time stepping (Robert 1993), and a dry and a moist version are standard benchmarks (Bryan and Fritsch 2002).
04a-vertical-slice.ipynb (Bryan and Fritsch’s dry case, Boussinesq, 100 m grid): the potential temperature anomaly (colour), the vorticity (contours) and the centre of buoyancy (cross) at 0, 5, 10 and 17 minutes.
6.2 A cold bubble: density currents
A bubble of cold air falls, hits the ground, and spreads along it as a density current (or gravity current): heavy fluid running under light fluid (Figure 13). Its front has a raised head; behind it the shear between the current and the air above rolls up into Kelvin–Helmholtz billows; ahead of it the warm air is lifted over the head. Its speed is of order \(\sqrt{g'h}\), where \(g' = g\,|\theta'|/\theta_0\) is the reduced gravity and \(h\) the depth of the current (Benjamin 1968). Density currents are everywhere: the gust fronts of thunderstorms, marked by shelf clouds (Figure 1), where the lifted warm air can trigger new storms; sea-breeze fronts (lecture 1); and, in the ocean, dense water spilling over sills, as the salty Mediterranean water leaving the Strait of Gibraltar along the bottom.
The cold bubble is the benchmark of Straka et al. (1993): in a domain 51.2 km long and 6.4 km high, a bubble 15 K colder than its neutral surroundings, centred 3 km above the ground with radii of 4 km and 2 km, falls and spreads for 900 s. Its front, the \(-1\) K contour at the ground, moves at more than 30 m/s at first, then slows to about 15 m/s as billows form on its top, one every 300 s or so. The billows make the problem hard: without dissipation the solution never converges, because finer grids resolve ever smaller billows. The benchmark therefore fixes a viscosity and diffusivity of \(75\ \mathrm{m^2/s}\), for which a reference solution converged at 25 m spacing exists; at 100 m spacing most of the schemes compared in the paper placed the front within a few per cent of it, at 400 m the errors in the front speed ranged from 2 to 24% (Straka et al. 1993). The notebook runs a Boussinesq version of this test (Figure 14): its front at 900 s is at 16.3, 16.7 and 17.0 km with grids of 400, 200 and 100 m, and converges to about 16.8 km on finer grids (50 and 25 m), 8% ahead of the compressible reference (15.5 km): the difference comes from the equations (the Boussinesq approximation over a 6 km deep atmosphere), not from the grid. Its front moves at about 16 m/s at the end, some \(1.4\sqrt{g'h}\) with the depth and mean temperature deficit of the current measured in the run. Without viscosity the run stays bounded, but its temperature overshoots by more than 10 K: there is no converged solution to compare with.
04a-vertical-slice.ipynb (Boussinesq version, 100 m grid, \(\nu = \kappa = 75\ \mathrm{m^2/s}\)) at 0, 5, 10 and 15 minutes: the bubble falls, reaches the ground and spreads as a density current with billows on its top. Triangles: the front (\(-1\) K at the ground); open triangle: Straka et al.’s compressible reference at 900 s.
6.3 Convection in models
A thunderstorm’s updraft is one to ten kilometres wide. A global model with 9 km cells cannot represent it; instead, at each step, every column is checked, and where it is unstable (\(N^2 < 0\), or unstable for a saturated parcel), a parameterization moves heat and moisture up as the unresolved convection would. The simplest is convective adjustment: mix the unstable part of the column at once. The optional notebook 04d-column-mixing.ipynb (appendix) does this for the ocean’s winter mixed layer, with an implicit vertical step. Between grid spacings of about 10 and 1 km deep convection is partly resolved, and the parameterizations must step back gradually: this range, partly resolved and partly parameterized, is the “grey zone” (lecture 9; a term in common use among modellers, cited here without a single source). The ocean convects too, in winter, in a few places where storms cool and evaporate the surface water enough to make it sink to great depths: the Gulf of Lion in the Mediterranean, the Labrador and Greenland Seas (Marshall and Schott 1999).
7 Flow over mountains
7.1 The lower boundary
Let a uniform wind \(U\), in a fluid of constant \(N\), blow over the ground \(z = h(x)\). The ground is a streamline, so at the ground the flow is tangent to it, \(w = (U + u')\,dh/dx \approx U\,dh/dx\): up on the windward side, down in the lee (Figure 15). The question is what the stratified fluid above does with this push.
7.2 Linear theory
Write the streamfunction as that of the uniform wind plus a small disturbance, \(\psi = -Uz + \psi'\) (so that \(u = U + u'\)), and look for a steady state. Linearized, the advection operator becomes \(U\,\partial_x\), and the slice equations give \(U\partial_x\xi' = \partial_x b'\) and \(U\partial_x b' + N^2\,\partial_x\psi' = 0\). Integrating in \(x\) (the disturbances vanish far upstream), \(b' = -N^2\psi'/U\) and \(\xi' = b'/U\), and with \(\xi' = \nabla^2\psi'\): \[ \nabla^2\psi' + \frac{N^2}{U^2}\,\psi' = 0, \qquad \psi'(x, 0) = U\,h(x), \] the second condition saying that the ground is a streamline, \(\psi = -Uh + \psi' = 0\) at \(z = h \approx 0\). This is a Helmholtz equation. A Fourier transform in \(x\) turns it into an ordinary differential equation for each wavenumber \(k\), with solutions \(\hat\psi(k, z) = U\hat h(k)\,e^{imz}\) and \[ m^2 = \frac{N^2}{U^2} - k^2 . \] The length \(U/N\) (1 km for \(U = 10\) m/s and \(N = 0.01\ \mathrm{s^{-1}}\)) separates two behaviours:
- \(|k| < N/U\), horizontal scales longer than \(2\pi U/N\): \(m\) is real and the disturbance is a wave that propagates upwards. The sign of \(m\) follows from a radiation condition: the energy must go up, away from the mountain, not come down from infinity. In the frame of the air the wave has the frequency \(\omega = -Uk\) (it is steady over the ground), and the group velocity of the previous section points upwards only if \(m\) has the sign of \(k\). The phase lines \(kx + mz = \mathrm{const}\) then tilt upstream with height (Queney 1948; Durran 1990).
- \(|k| > N/U\), narrow features: \(m = i\sqrt{k^2 - N^2/U^2}\) and the disturbance decays with height.
For a bell-shaped ridge, the Witch of Agnesi \(h(x) = h_0 a^2/(x^2 + a^2)\), whose spectrum \(\hat h \propto e^{-|k|a}\) is concentrated at \(k \lesssim 1/a\), one number sorts the flow (Figure 17):
- \(Na/U \ll 1\) (a narrow hill): the disturbance decays with height; no waves aloft.
- \(Na/U \approx 1\): non-hydrostatic waves, which also trail downstream.
- \(Na/U \gg 1\) (a wide ridge): hydrostatic waves stand above the mountain, tilting upstream, with the vertical wavelength \(2\pi U/N\) (about 6 km). For \(Na/U \to \infty\) the solution is known in closed form (Queney 1948): with the Scorer parameter \(l = N/U\) (its value for constant \(U\); here \(l\) is not the \(y\)-wavenumber), \(\psi'/U = h_0 a\,(a\cos lz - x\sin lz)/(x^2 + a^2)\).
The waves carry horizontal momentum upwards, and the mountain feels the opposite force, a drag: for the hydrostatic Witch of Agnesi it is \(\frac{\pi}{4}\rho_0 N U h_0^2\) per unit length of the ridge (Teixeira 2014). Where the waves break, high up, they deposit that momentum and slow the wind; global models cannot resolve most mountains, so they add this drag as a parameterization (orographic gravity-wave drag), another column-by-column scheme. For still wider mountains, \(a \gtrsim U/f\) (about 100 km), the rotation of the Earth enters and the response becomes the inertia–gravity and quasi-geostrophic flows of lecture 3 (Smith 1979).
04b-mountain-waves.ipynb): a hydrostatic case, with the phase lines tilting upstream and the energy going up.
04b-mountain-waves.ipynb.
7.3 Lee waves
When \(N/U\) decreases with height (a stable layer near the ground under less stable air, or a wind that increases with height), waves of intermediate wavelength propagate in the lower layer but decay in the upper one: they are trapped, and reflect between the ground and the upper layer as they are carried downstream. They form a long train of lee waves behind the mountain, with clouds on their crests where the air is lifted enough to condense (Figure 18). Scorer (1949) found the condition: the Scorer parameter \(l\), \(l^2(z) = N^2/U^2 - U''/U\) (\(l = N/U\) for constant \(U\)), must decrease upwards, and for two layers (the lower of depth \(H_1\)) trapped waves exist if \(l_1^2 - l_2^2 > \pi^2/(4H_1^2)\). Their wavelength lies between \(2\pi/l_1\) and \(2\pi/l_2\), typically 5–20 km.
04b-mountain-waves.ipynb.
7.4 High mountains: blocking and breaking
The linear theory needs the disturbance to be small, and the second number of the problem, \(Nh_0/U\), says when it is not. Lifting a parcel by \(h_0\) in a fluid of buoyancy frequency \(N\) costs about \(N^2h_0^2/2\) per unit mass; the parcel brings the kinetic energy \(U^2/2\). When \(Nh_0/U \gtrsim 1\) the low-level air cannot climb the mountain: in two dimensions it is blocked upstream, in three dimensions it goes around the mountain. Above it, the waves steepen until the isentropes overturn and the waves break, and the air crossing the crest can plunge down the lee slope in a downslope windstorm, as the bora on the Adriatic coast or the foehn of the Alps (Figure 19) (Durran 1990; Smith 1979). The Po valley, closed by the Alps and the Apennines, can hold its cold air under a winter inversion for days, with its fog and smog (the example below).
Example: waves or blocking? With \(N = 0.01\ \mathrm{s^{-1}}\) and \(U = 10\) m/s, \(U/N = 1\) km and hydrostatic mountain waves have the vertical wavelength \(2\pi U/N = 6.3\) km (Figure 20). A hill of half-width \(a = 500\) m and height \(h_0 = 50\) m has \(Na/U = 0.5\) and \(Nh_0/U = 0.05\): the air flows over it, with no waves aloft. A ridge 10 km wide and 500 m high has 10 and 0.5: hydrostatic waves stand above it, tilting upstream. The Alps, \(a \approx 100\) km and \(h_0 \approx 3\) km, have 100 and 3: the low air is blocked and goes around them (and rotation matters too, \(U/(fa) \approx 1\)). The Apennines, \(a \approx 20\) km and \(h_0 \approx 1.2\) km, change with the season. In a winter anticyclone the Po valley holds cold air under an inversion, \(N = 0.03\ \mathrm{s^{-1}}\), with a weak wind, \(U = 2\) m/s: \(Na/U = 300\) and \(Nh_0/U = 18\), firmly blocked, and the valley keeps its cold air, its fog and its smog for days. On a summer afternoon, \(N = 0.005\ \mathrm{s^{-1}}\) and \(U = 10\) m/s give 10 and 0.6: the air flows over and makes waves.
8 The notebooks
04a-vertical-slice.ipynb(in class): the slice model in vorticity–streamfunction form (the three equations above; centred differences with Arakawa’s Jacobian, the Poisson problem solved exactly by fast Fourier and sine transforms, the three-stage Runge–Kutta of lectures 5–6), with a check of the dispersion relation (second order), the warm bubble, Straka’s cold bubble, the St Andrew’s cross and, optionally, convection in an unstable layer. Its exercises 4.1 (the beam angles for two other frequencies) and 5.1 (the growth rate of convective cells against their width) are mandatory.04b-mountain-waves.ipynb(demo, or at home): linear mountain waves by FFT, checked against Queney’s solution; the regimes; trapped lee waves; Schär’s mountain; a time-dependent model with and without a sponge layer. Its exercises 2.1 (your own mountain) and 5.1 (how deep and how strong the sponge must be) are mandatory.- Optional, from the appendix:
04c-vertical-modes.ipynband04d-column-mixing.ipynb.
9 What remains open
One model, the vertical slice, carried the whole lecture: the 1950 vorticity equation turned on its side and driven by buoyancy. It explains why a layered fluid carries waves with frequencies below \(N\), why warm bubbles rise as vortex pairs and cold ones spread as density currents, and what a mountain does to the wind: waves aloft, trapped lee waves, or blocking. The numerics remain: how to step such equations stably and accurately (finite differences and finite volumes, lectures 5 and 6); how to solve the elliptic problem of every step fast, on a plane and on the sphere (Fourier series and spherical harmonics, lecture 7); and how operational models treat the vertical: staggered grids, an implicit vertical for sound, a sponge at the top, and mountains (lecture 9).
Carry forward
- Stratification: \(N^2 = -(g/\rho_0)\,d\bar\rho/dz\) (air: \((g/\theta)\,d\theta/dz\)); buoyancy periods of minutes; \(N^2 < 0\): convection.
- The vertical slice: \(D\xi/Dt = \partial b/\partial x\), \(\nabla^2\psi = \xi\): horizontal buoyancy gradients make vorticity; a Poisson problem at every step.
- Internal gravity waves: \(\omega = N\cos\alpha < N\), energy along the crests; the hydrostatic approximation fails when the motion is as tall as it is wide.
- Warm bubbles rise as vortex pairs; cold air spreads as a density current, at a speed of order \(\sqrt{g'h}\).
- Mountains: \(Na/U\) (waves aloft or decay), \(Nh_0/U\) (linear or blocked); waves tilt upstream, with the vertical wavelength \(2\pi U/N\).
Appendix: vertical modes (optional)
This appendix keeps the vertical modes of the earlier draft of the lecture, for self-study and projects, and for lecture 7, whose semi-implicit method solves one shallow water problem per vertical mode.
The linear equations
Consider small motions about a state of rest with stratification \(N^2(z)\), in a layer of depth \(H\) with a flat bottom and a rigid lid, on an \(f\)-plane, under the hydrostatic and Boussinesq approximations. With \(\phi = \delta p/\rho_0\) the pressure perturbation and \(b\) the buoyancy perturbation (the departure from the mean buoyancy \(\bar b(z)\), whose gradient is \(N^2 = d\bar b/dz\); lecture 3), \[ \begin{aligned} &\partial_t u - f v = -\partial_x\phi, \qquad \partial_t v + f u = -\partial_y\phi, \\ &\partial_z\phi = b, \qquad \partial_t b + N^2(z)\,w = 0, \\ &\partial_x u + \partial_y v + \partial_z w = 0, \qquad w = 0 \ \text{at}\ z = -H \ \text{and}\ z = 0 . \end{aligned} \] The first line is the horizontal momentum of the shallow water equations (lecture 3). What is new is the second line: hydrostatic balance with buoyancy, and the buoyancy equation, in which the buoyancy of a parcel changes only because it moves vertically across the mean stratification.
Separation of variables
Look for solutions in which every field is a horizontal pattern times a vertical shape (Figure 21): \[ (u, v, \phi) = (\hat u, \hat v, \hat\phi)(x, y, t)\; C(z), \qquad w = \hat w(x, y, t)\; W(z), \qquad C = \frac{dW}{dz}, \] with the symbols of Vallis (2017) (§3.4.1) for the horizontal shape \(C\) and of Gill (1982) (§6.11, where \(W\) is called \(R\)) for the vertical-velocity shape \(W\). The index \(m\) below labels the vertical modes (it is not the vertical wavenumber of the mountain waves). \[ \] The continuity equation gives \((\partial_x\hat u + \partial_y\hat v)\,C + \hat w\,W' = 0\), and \(W' = C\), that is \(\hat w = -(\partial_x\hat u + \partial_y\hat v)\). Differentiating the hydrostatic equation in time and using the buoyancy equation, \(\partial_t\partial_z\phi = -N^2 w\), gives \(\partial_t\hat\phi\;C' = -N^2\,\hat w\,W\), or \[ \frac{C'}{N^2\,W} = -\frac{\hat w}{\partial_t\hat\phi} . \] The left side depends only on \(z\), the right side only on \((x, y, t)\): both equal a constant, which we call \(-1/c^2\). This gives two problems.
The vertical problem is a Sturm–Liouville eigenvalue problem, \[ W'' + \frac{N^2(z)}{c^2}\,W = 0, \qquad W(-H) = W(0) = 0 . \] Eliminating \(W\) gives the form of Vallis (2017) (§3.4.1), \(\left(C'/N^2\right)' + C/c^2 = 0\) with \(C' = 0\) at the bottom and the lid. Since \(N^2 > 0\), it has an infinite sequence of eigenvalues \(c_1 > c_2 > \dots > 0\), and the eigenfunctions \(W_n\) are orthogonal with weight \(N^2\); mode \(m\) has \(m-1\) zeros inside the layer, and its horizontal velocity, with the shape \(C_m = W_m'\), has \(m\) zeros. For constant \(N\), \(W_m = \sin(m\pi z/H)\) and \(c_m = NH/(m\pi)\).
The horizontal problem is \(\partial_t\hat\phi = c^2\hat w = -c^2(\partial_x\hat u + \partial_y\hat v)\) together with the horizontal momentum equations: exactly the linear shallow water equations of lecture 3, with \(gH\) replaced by \(c_m^2\). Each mode is therefore a shallow water system with an equivalent depth \(H_m = c_m^2/g\), gravity waves of speed \(c_m\), Kelvin and Rossby waves, and a deformation radius \(L_m = c_m/f\). With a free surface instead of a rigid lid there is also the barotropic mode, nearly uniform in the vertical, with \(c_0 \approx \sqrt{gH}\). Figure 22 lists the equivalent depths of an ocean and of the troposphere: the same equations, with depths from 4 km down to a few centimetres.
04c-vertical-modes.ipynb (baroclinic modes computed, barotropic mode \(\sqrt{gH}\) for a free surface); troposphere: constant \(N = 10^{-2}\ \mathrm{s^{-1}}\), \(H = 10\) km, \(c_m = NH/(m\pi)\).
Numbers
For the ocean (Figure 23), the first baroclinic mode travels at 2–3 m/s at mid-latitudes, with an equivalent depth of less than a metre and a deformation radius of 20–30 km at the latitude of Milan; the barotropic mode travels at about 200 m/s, with a deformation radius of about 2000 km. For the troposphere, with \(N = 10^{-2}\ \mathrm{s^{-1}}\) and \(H = 10\) km, \(c_1 \approx 30\) m/s and \(L_1 \approx 300\) km. (The deformation radius of the stratified atmosphere, \(L_d = NH/f \approx 1000\) km as in lecture 3, is \(\pi\) times larger than \(L_1\): it uses the depth \(H\) as the vertical scale, while mode 1 has the vertical wavenumber \(\pi/H\), so \(c_1 = NH/\pi\) and \(L_1 = c_1/f\). Vallis uses both, \(L_d = NH/f\) and the mode speeds \(c_m = NH/(m\pi)\); here \(L_1\) is meant whenever a mode is.) The ocean’s eddies, set by \(L_1\), are therefore about ten times smaller than the atmosphere’s cyclones, and ocean models need grids of a few kilometres to resolve them (lecture 6, \(L_d/\Delta x\)). A useful estimate is the WKB approximation \(c_m \approx \frac{1}{m\pi}\int_{-H}^{0} N\,dz\) (Gill 1982; Vallis 2017).
Numerically
The notebook notebooks/04c-vertical-modes.ipynb discretizes \(W''\) by second-order centred differences on \(J\) interior points, which turns the problem into the generalized symmetric eigenvalue problem \(-D_2\mathbf{W} = \lambda\,\mathrm{diag}(N_j^2)\,\mathbf{W}\), \(\lambda = 1/c^2\), solved by scipy.linalg.eigh. For constant \(N\) the error in \(c_m\) falls by a factor of 4 each time \(J\) doubles (second order, as the difference formula), and grows with \(n\): higher modes, with more zeros, need more points.
04c-vertical-modes.ipynb, with their speeds \(c_m\) and deformation radii \(L_m = c_m/f\) at Milan. Left: constant \(N = 2\times10^{-3}\ \mathrm{s^{-1}}\) (cosines, \(c_m = NH/m\pi\)); right: the ocean with a thermocline of Figure 2, whose modes are concentrated near the surface. Dashed: the barotropic mode, uniform in the vertical, with \(c_0 = \sqrt{gH} = 198\) m/s.
Example: how fine must an ocean model be? For constant \(N\), \(c_1 = NH/\pi\), \(H_1 = c_1^2/g\) and \(L_1 = c_1/f\); at Milan, \(f = 1.04\times10^{-4}\ \mathrm{s^{-1}}\). The ocean, \(N = 2\times10^{-3}\ \mathrm{s^{-1}}\) and \(H = 4\) km: \(c_1 = 2.5\) m/s, \(H_1 = 0.66\) m, \(L_1 = 24\) km. The troposphere, \(N = 10^{-2}\ \mathrm{s^{-1}}\) and \(H = 10\) km: \(c_1 = 32\) m/s, \(H_1 = 103\) m, \(L_1 = 306\) km. For the ocean with a thermocline of Figure 2, \(\int_{-H}^0 N\,dz = 4.0 + 3.2 = 7.2\) m/s and the WKB estimate gives \(c_1 = 2.3\) m/s (the eigenvalue solver of 04c-vertical-modes.ipynb: 2.1 m/s). The radius falls towards the poles as \(f\) grows (Figure 24): \(L_1 = 46\) km at \(20^\circ\)N, 22 km at Milan and 17 km at \(70^\circ\)N. A model sees eddies only if its grid spacing is at most about \(L_1/3\) (a rule of thumb): 15, 7 and 6 km.
Column models
Processes too small to resolve (turbulent mixing, convection, radiation, clouds) are represented in models by parameterizations that act column by column. The notebook notebooks/04d-column-mixing.ipynb builds the simplest one, for the ocean’s mixed layer in winter: vertical diffusion with a turbulent diffusivity that is large near the surface, a surface heat loss, and convective adjustment, which mixes at once any part of the column where denser water lies over lighter water. Explicit time stepping would need steps shorter than \(\Delta z^2/(2K) \approx 40\) s (lecture 5); backward Euler takes steps of an hour, with one tridiagonal solve per step, and, written in finite-volume form, conserves the heat content exactly (lecture 6). After 90 days of losing 200 W/m² the mixed layer deepens to about 160 m (Figure 25).
04d-column-mixing.ipynb: an autumn profile (start) and the profiles after 30, 60 and 90 days of a surface heat loss of 200 W/m², with implicit diffusion and convective adjustment: the mixed layer cools and deepens to about 160 m.
Project ideas
- The density current of the notebook at several resolutions and viscosities: front speed against \(\sqrt{g'h}\); does the solution converge as the grid is refined at fixed viscosity, and as the viscosity goes to zero?
- A sea-breeze front: heat a strip of land in the slice model and follow the front inland.
- Mountain waves over a real Alpine cross-section (from a digital elevation model): linear theory against a nonlinear run of the slice model; foehn or bora conditions.
- A moist bubble: add water vapour and latent heat as a buoyancy source (after Bryan and Fritsch 2002).
- The slice model with a finite-volume, limited advection of buoyancy (lecture 6) against the centred scheme: overshoots, fronts, conservation.
- From the appendix: the vertical modes of a real ocean profile (World Ocean Atlas) against eddy sizes; a column model of the atmospheric boundary layer over a day; the pressure-gradient error of \(\sigma\) coordinates over a seamount.
Further reading
Vallis (2017) for the Boussinesq equations, stratification and internal waves; Gill (1982), chapters 6 and 8, for internal waves and vertical modes; Durran (1990) and Smith (1979) for mountain waves; Straka et al. (1993) for the density-current benchmark; Durran (2010) for the numerics of fast waves and non-reflecting boundaries.