$$ \def\ab{\boldsymbol{a}} \def\bb{\boldsymbol{b}} \def\eb{\boldsymbol{e}} \def\fb{\boldsymbol{f}} \def\gb{\boldsymbol{g}} \def\kb{\boldsymbol{k}} \def\nb{\boldsymbol{n}} \def\ub{\boldsymbol{u}} \def\vb{\boldsymbol{v}} \def\xb{\boldsymbol{x}} \def\Xb{\boldsymbol{X}} \def\Ab{\boldsymbol{A}} \def\Bb{\boldsymbol{B}} \newcommand{\Cb}{\boldsymbol{C}} \newcommand{\Db}{\boldsymbol{D}} \def\Eb{\boldsymbol{E}} \def\Fb{\boldsymbol{F}} \def\Jb{\boldsymbol{J}} \def\Ub{\boldsymbol{U}} \def\xib{\boldsymbol{\xi}} \def\evx{\boldsymbol{e}_x} \def\evy{\boldsymbol{e}_y} \def\evz{\boldsymbol{e}_z} \def\evr{\boldsymbol{e}_r} \def\evt{\boldsymbol{e}_\theta} \def\evp{\boldsymbol{e}_r} \def\evf{\boldsymbol{e}_\phi} \def\evb{\boldsymbol{e}_\parallel} \def\omb{\boldsymbol{\omega}} \def\dA{\;\mathrm{d}\Ab} \def\dS{\;\mathrm{d}\boldsymbol{S}} \def\dV{\;\mathrm{d}V} \def\dl{\mathrm{d}\boldsymbol{l}} \def\bfzero{\boldsymbol{0}} \def\Rey{\mathrm{Re}} \newcommand{\dds}[1]{\frac{\mathrm{d}{#1}}{\mathrm{d}s}} \newcommand{\ddy}[2]{\frac{\partial{#1}}{\partial{#2}}} \newcommand{\ddt}[1]{\frac{\mathrm{d}{#1}}{\mathrm{d}t}} \newcommand{\DDt}[1]{\frac{\mathrm{D}{#1}}{\mathrm{D}t}} \newcommand{\gradH}{\nabla_\mathrm{H}} \newcommand{\grad}{\nabla} \newcommand{\Dl}{\mathrm{\Delta}} \newcommand{\Omb}{\boldsymbol{\mathrm{\Omega}}} \newcommand{\Omrm}{\mathrm{\Omega}} \newcommand{\Ld}{\mathcal{L}_d } \newcommand{\nd}[1]{{#1}^*} \newcommand{\vh}{\nd{\boldsymbol{v}}} \newcommand{\bh}{\nd{b}} \newcommand{\ph}{\nd{p}} \newcommand{\fh}{\nd{\fb}} \newcommand{\tha}{\nd{t}} \newcommand{\wh}{\nd{w}} \newcommand{\zh}{\nd{z}} \newcommand{\gradHh}{\nd{\nabla}_H} \newcommand{\psih}{\nd{\psi}} \newcommand{\Ro}{\mathrm{Ro}} \newcommand{\ddz}[1]{\frac{\mathrm{d}{#1}}{\mathrm{d}z}} $$
2 Stratification
2.1 Buoyancy
A heavy object sinks down in a fluid because of gravity, but lighter objects float because of an upward force that the fluid exerts on them opposite to the direction of gravity. This force – known as buoyancy – is due to the weight of displaced fluid.
Typically, fluids in the atmosphere and ocean are stratified. That is, the density of air in the atmosphere and of water in the ocean varies with height. In a stably stratified medium (where lower density layers are above higher density ones), when a fluid parcel is pushed down vertically, we expect that it will try to come back up due to buoyancy force, and vice-versa. By considering the motion of such a fluid parcel, we can infer its buoyancy which is an important property of stratified fluids.
Consider a fluid at rest with density \(\tilde\rho(z)\) and hydrostatic pressure \(\tilde p\), i.e., \(d\tilde p/dz = -g \tilde\rho\). Then consider raising a parcel of fluid from a height \(z^*\) to \(z^* + \Delta z\) (see Figure 2.1). At \(z^*\), the environment and the parcel are in hydrostatic balance and so \[ \frac{d\tilde p}{dz}\bigg |_{z^*} = -g \tilde\rho(z^*) = -g \rho_p(z^*), \] where we denote by a subscript \(p\) a variable following the parcel and a tilde denotes an environment variable.
The environment at \(z^*+\Delta z\) is also in hydrostatic balance: \[ \frac{d\tilde p}{dz}\bigg |_{z^*+\Delta z} = -g \tilde\rho(z^* +\Delta z). \]
However, the fluid parcel at \(z^*+\Delta z\) is now being acted on by a force given by \[ -\frac{d\tilde p}{dz}\bigg |_{z^*+\Delta z} -g \rho_p(z^* +\Delta z) \] since the pressure of the parcel is the same as the pressure of the environment. Note the pressure gradient acts upwards and gravity downwards. Then, since \(-g \rho_p(z^* +\Delta z) \neq -g\tilde\rho(z^* +\Delta z) = {d\tilde p}/{dz}\bigg |_{z^*+\Delta z}\), there is a net force acting on the fluid parcel (per unit volume) which can be written as \[ F_b = g \tilde\rho(z^* +\Delta z) - g \rho_p( z^* +\Delta z) \tag{2.1}\]
We can then apply Newton’s Second Law to the parcel at \(z^*+\Delta z\) to give \[ \rho_p(z^* +\Delta z) \frac{d^2 \Delta z}{dt^2} = F_b \] and so \[ \frac{d^2 \Delta z}{dt^2} = \frac{g\left[ \tilde \rho(z^* +\Delta z) - \rho_p( z^* +\Delta z) \right]}{\rho_p(z^* +\Delta z) }. \tag{2.2}\] Next, if we assume the density of the parcel is conserved then we have that \(\rho_p(z^* +\Delta z)=\rho_p(z^*)=\tilde\rho(z^*)\). Then we can write Equation 2.2 as \[ \frac{d^2 \Delta z}{dt^2} = \frac{g\left[ \tilde \rho(z^* +\Delta z) - \tilde\rho( z^*) \right]}{\tilde\rho(z^*) }\] \[ \approx\frac{g}{\tilde\rho(z^*)}\frac{d\tilde\rho}{dz}\bigg |_{z^*}\Delta z \tag{2.3}\] where we have used a Taylor series for \(\tilde\rho(z^* +\Delta z)\) to arrive at the final expression. We see that if \({d\tilde\rho}/{dz}<0\) then solutions to Equation 2.3 are oscillatory with frequency \(N(z^*)\) where in general \[ N^2(z) = -\frac{g}{\tilde\rho}\frac{d\tilde\rho}{dz}. \tag{2.4}\] The quantity \(N(z)\) is known as the buoyancy frequency or Brunt–Väisälä frequency. It gives the theoretical frequency of vertical oscillations of a small parcel of fluid in a stratified medium.
If \({d\tilde\rho}/{dz}>0\) then \(N^2<0\) and the stratification is unstable meaning that if we raise a fluid parcel it will continue to rise.
2.2 Boussinesq equations
2.2.1 Density and temperature variation in the ocean and atmosphere
In the lower atmosphere and most parts of the ocean, the density (or temperature) does not change significantly at different horizontal locations \(\xb=(x,y)\) or in time. This motivates the Boussinesq approximation which involves writing the density (and pressure) as a sum of a base profile plus small perturbation, where the base profile only depends on \(z\). We will formulate this approximation mathematically in the next section, but before that it will be insightful to look at few observations of density and temperature variations in the ocean and atmosphere.
Figure 2.2 shows a combination of two types of measurements in a region in the arctic ocean. In panel (a), there are four red dots that are mooring measurements obtained by a collection of devices connected to a wire and anchored on the sea floor. The other dots are measurements taken by a cruise that sends diving devices down into the ocean. The details of these observations are not important for us, but they show that several locations are sampled by different measuring devices and at different times. Panel (b) shows the density profiles of the cruise measurements at different locations. It is clear that the density profile has a (horizontal) mean that is a function of \(z\). The perturbations to this mean by changing the horizontal location are much smaller than the mean value itself. The variation of density from the base profile is better quantified in the panel (c). The solid black curve shows the mean (base) profile of the density next to a single cast with dashed black curve, which marginally deviates from the mean. Panel (d) shows the density profiles for a fixed location (location D in panel (a)), but at different times. It can be inferred from this plot that changing time does not affect the mean density profile very much. Looking at different part of the ocean, we get a more general base profile (a smoother mean). Nevertheless, the underlying message of this example still holds: the density is well approximated with a base profile that is a function of depth only plus small perturbations.
Figure 2.3 displays the (daily) temperature profiles over a region in India. Atmospheric temperature is directly related to the density at each altitude (which has a fixed pressure). Hence, the relative magnitude of the mean temperature to perturbations can translate to the density as well. It is evident from this figure that changes in temperature at each altitude are relatively small for lower atmosphere (below 15km). There is a jump in temperature that is known as “tropopause”. Near the tropopause, the deviations from the mean temperature profile are larger, and the Boussinesq approximation is not justified.
Measurements similar to the examples just discussed are abundant in the literature thanks to several decades of observational campaigns. Considering them all together, we can conclude that in most parts of the ocean, and in the lower atmosphere, the variations in density at a fixed \(z\) are small relative to their means.
2.2.2 Boussinesq approximation
As we have seen, the density in the ocean and in the lower layer of atmosphere varies by only a few percent across the depth or altitude. We can therefore exploit the smallness of these variations and write \[ \rho(\xb,t) = \bar\rho + \delta\rho(\xb,t) \] where \(\bar\rho\) is the average density (a constant) and \(\delta\rho\) are (small) perturbations away from this mean. The perturbations are assumed small so that \(|\delta\rho|\ll\bar\rho\). Substituting this into the inviscid, incompressible momentum equation (i.e., Equation 1.7 with \(\nu=0\)) and again assuming gravity is the only body force acting so that \(\fb_b=-g\eb_z\), we find \[ (\bar\rho + \delta\rho(\xb,t)) \frac{D\ub}{Dt} = -\nabla p - (\bar\rho + \delta\rho(\xb,t))g\eb_z. \] Since \(|\delta\rho|\ll\bar\rho\) we may be tempted to ignore all the terms involving \(\delta\rho\). Instead, we neglect all the terms involving \(\delta\rho\) except in the term involving gravity (here there is only one but in more complex cases there may be others). This is because while \(|\delta\rho|\ll\bar\rho\), gravity is strong and so the buoyancy term (the one involving gravity) is retained. Making this approximation, known as the Boussinesq approximation, leads to \[ \bar\rho \frac{D\ub}{Dt} = -\nabla p - (\bar\rho + \delta\rho(\xb,t))g\eb_z. \] Then, dividing by \(\bar\rho\) and combining with equation (1.4) and (1.6), we obtain the so called Boussinesq equations: \[\frac{\partial\ub}{\partial t} + (\ub\cdot\nabla)\ub = -\frac{1}{\bar\rho}\nabla p -\frac{\rho}{\bar\rho}g\eb_z. \tag{2.6}\] \[ \nabla\cdot\ub = 0, \tag{2.7}\] \[\frac{D\rho}{Dt} = \frac{\partial \rho}{\partial t} + (\ub\cdot\nabla)\rho= 0, \tag{2.8}\]
Although we have derived this set of Boussinesq equations starting from the assumption of an incompressible fluid, it is possible to show that the Boussinesq approximation applies to a more general compressible fluid.
The Boussineq approximation is valid everywhere in the ocean since density variations are small. However, it is also valid for any motions in the atmosphere that are confined to a layer that is thin compared with the scale height of the atmosphere. For the Earth’s atmosphere, we have seen that the scale height is approximately \(8\)km and so the Boussinesq approximation should apply to motions that are confined to layers of about \(1\)km or less.
2.2.3 Buoyancy as a variable
Next, let’s consider splitting our variables into a reference state and perturbations to this state. We will consider a static, time-independent reference state with a vertical density stratification \(\rho=\tilde\rho(z)\). Equation 2.7 and Equation 2.8 are automatically satisfied by this reference state and Equation 2.6 demands \[ p=\tilde p(z) \quad \text{with} \quad \frac{d\tilde p}{dz} = -\tilde\rho g. \tag{2.9}\] Decomposing our variables in this way gives \[ \ub = \ub'(\xb,t), \quad \rho = \tilde\rho(z) + \rho'(\xb,t), \quad p = \tilde p(z) + p'(\xb,t) \tag{2.10}\] where a tilde denotes a reference state variable and the primed quantities are perturbations to these states (for now no assumptions are made about the size of the perturbations).
Now, substituting (2.10) into Equation 2.6 gives \[ \frac{\partial\ub'}{\partial t} + (\ub'\cdot\nabla)\ub' = -\frac{1}{\bar\rho}\nabla p' -\frac{\rho'}{\bar\rho}g\eb_z = -\frac{1}{\bar\rho}\nabla p' +b\eb_z, \] where we have used Equation 2.9 to cancel reference state terms and introduced the buoyancy variable \[ b = -\frac{\rho'}{\bar\rho} g . \tag{2.11}\] Note that this expression is analogous to the buoyancy force per unit mass \(F_b/\bar\rho\) introduced in Equation 2.1.
Substituting (2.10) into Equation 2.7 simply gives \(\nabla\cdot\ub'=0\). Finally, substituting (2.10) into Equation 2.8 gives \[ \frac{\partial(\tilde \rho(z)+\rho')}{\partial t} + (\ub'\cdot\nabla)(\tilde\rho(z)+\rho')=0\] \[\Rightarrow \frac{\partial\rho'}{\partial t} + (\ub'\cdot\nabla)\rho' + w'\frac{d\tilde\rho}{dz}=0. \tag{2.12}\] Multiplying Equation 2.12 by \(-\frac{g}{\bar\rho}\) allows us to form \[\begin{align*} \frac{Db}{Dt}+\left(\frac{-g}{\bar\rho}\frac{d\tilde\rho}{dz}\right)w'&=0,\\ \Rightarrow\frac{Db}{Dt}+N^2(z)w'&=0 \end{align*}\] where \[ N^2(z) = -\frac{g}{\bar\rho}\frac{d\tilde\rho}{dz}. \] Here \(N(z)\) is the buoyancy frequency in the Boussinesq approximation. Note it is subtly different to the definition in Equation 2.4.
Gathering together the Boussinesq equations for the perturbation quantities we have the following equation set:
where we have dropped the primes from perturbation quantities as is customary.
Note here that \(p\) is a perturbation to the background pressure state. In general, when \(p\) appears along side buoyancy in the governing equations then this indicates it is a perturbation pressure. This is in contrast to the Boussinesq equations where \(p\) in Equation 2.6 is taken to mean the full pressure (likewise in Equation 2.6, \(\rho\) is the full density).
The buoyancy (Brunt–Väisälä) frequency that appears in Equation 2.15 is an important quantity in different regimes of geophysical flows. Figure 2.4 shows the values of \(N^2\) for a region in the atmosphere/ocean as functions of altitude/depth. In the left plot, \(N^2\) over the region of Gadanki in India is shown. There is a noticeable jump at around 18km altitude (similar to what we observed in figure Figure 2.3, marking the “tropopause”, which is the border between two layers of atmosphere (troposphere and stratosphere). Note that the height of tropopause varies to some degrees in different locations (that is why the jump in Figure 2.3 and Figure 2.4 appears at slightly different altitudes). Around the tropopause, constant \(N\) or even the Boussinesq approximation are not valid, but in the lower atmosphere (below the tropopause) \(N\) can be assumed constant. Horizontal bars show standard deviations obtained while considering 43 different measurements during April 2006 to March 2007. The low variability in the lower atmosphere (relatively small bars) justifies the Boussinesq approximation, whereas the high variability near tropopause shows this approximation is no longer valid around this altitude.
The right plot is obtained from a mooring observation (“NISKINe” campaign). For the deep ocean (lower than 800 meters), it is safe to use constant \(N\). However, in the upper ocean the variation of \(N\) with depth cannot be neglected. A comparison between the two figures shows that \(N\) is an order of magnitude smaller in the ocean.
2.2.4 Hydrostatic approximation
A common approximation in geophysical fluid dynamics is to assume that the horizontal scale is large compared to the vertical scale. In this case, the vertical component of the momentum equation is simply replaced by our hydrostatic balance equation. This is known as the hydrostatic approximation and leads to the following hydrostatic Boussinesq equations
where \(\vb=(u,v)\) is the horizontal velocity and \(\nabla_H=\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}\right)\).
Note that the hydrostatic approximation has resulted in the vertical momentum equation being replaced by hydrostatic balance.
The Boussinesq equations (2.13) - (2.15) (and their hydrostatic counterparts (2.16) - (2.19)) are the fundamental equations that are solved numerically in weather forecast and climate models. Because we observe and measure everything from the rotating Earth, these equations are solved in the rotating frame of reference (which we will come to in Chapter 4). Most climate models use the Hydrostatic equations as they require far fewer computational resources. Hence, because climate models (unlike weather models) are run for very long times (several years or decades), it is worth reducing the computational cost at the expense of making the hydrostatic assumption.
In the geophysical fluids literature, equations (2.16) - (2.19) are often known as the primitive equations. This is somewhat confusing as the non-Hydrostatic Boussinesq equation set is more general. However, the Hydrostatic Boussinesq equations are still more accurate than some other reduced equations such as the shallow water equations and quasigeostrophic equations that we will study in Chapters 3 and 5, respectively.