2.3 Density & Stratification

Step 1: Taylor Expansion of Density

We begin by expanding $\rho$ about a reference state $(S_0, T_0, p_0)$ in a Taylor series. To first order in the perturbations $\Delta T = T - T_0$, $\Delta S = S - S_0$, $\Delta p = p - p_0$:

$$\rho \approx \rho_0 + \frac{\partial \rho}{\partial T}\bigg|_{S,p}\Delta T + \frac{\partial \rho}{\partial S}\bigg|_{T,p}\Delta S + \frac{\partial \rho}{\partial p}\bigg|_{T,S}\Delta p$$

Step 2: Define Thermodynamic Coefficients

The three partial derivatives define the thermal expansion coefficient $\alpha$, the haline contraction coefficient $\beta$, and the isothermal compressibility $\kappa$:

$$\alpha = -\frac{1}{\rho}\frac{\partial \rho}{\partial T}\bigg|_{S,p}, \quad \beta = \frac{1}{\rho}\frac{\partial \rho}{\partial S}\bigg|_{T,p}, \quad \kappa = \frac{1}{\rho}\frac{\partial \rho}{\partial p}\bigg|_{T,S}$$

Note the sign conventions: $\alpha$ is positive because warming decreases density (negative $\partial\rho/\partial T$), while $\beta$ is positive because increasing salinity increases density.

Step 3: Linearized Equation of State

Substituting the coefficient definitions back into the Taylor expansion, the linearized equation of state becomes:

$$\rho = \rho_0\left[1 - \alpha(T - T_0) + \beta(S - S_0) + \kappa(p - p_0)\right]$$

This is the form used in many analytical ocean models (the Boussinesq approximation). Typical values:$\alpha \approx 2 \times 10^{-4}$ K$^{-1}$,$\beta \approx 7.5 \times 10^{-4}$ (g/kg)$^{-1}$,$\kappa \approx 4.5 \times 10^{-10}$ Pa$^{-1}$.

Step 4: UNESCO Polynomial (EOS-80)

The full nonlinear equation of state is a high-order polynomial fit to laboratory measurements. The UNESCO formula (Millero et al., 1980) expresses density at atmospheric pressure as:

$$\rho(S, T, 0) = \rho_w(T) + A(T)\,S + B(T)\,S^{3/2} + C\,S^2$$

where $\rho_w(T)$ is pure water density (a 5th-order polynomial in $T$), and $A(T)$,$B(T)$ are polynomial functions of temperature. The $S^{3/2}$ term reflects the Debye-Huckel limiting law for electrolyte solutions. Pressure dependence is added via the secant bulk modulus$K(S, T, p)$:

$$\rho(S, T, p) = \frac{\rho(S, T, 0)}{1 - p / K(S, T, p)}$$

Step 5: From In-Situ Density to $\sigma_t$

The density anomaly $\sigma_t$ simply evaluates the equation of state at $p = 0$ and subtracts 1000 to make the numbers more convenient:

$$\sigma_t = \rho(S, T, 0) - 1000 \text{ kg/m}^3$$

This removes the dominant (and dynamically uninteresting) pressure contribution to density, leaving only the T and S effects. Values typically range from 23 to 28 kg/m$^3$.

Step 6: TEOS-10 and the Gibbs Function

The modern TEOS-10 standard (2010) replaces the empirical polynomial with a thermodynamically consistent Gibbs potential $g(S_A, T, p)$. All thermodynamic properties are derived from this single function:

$$\rho = \left(\frac{\partial g}{\partial p}\right)^{-1}, \quad \alpha = \frac{1}{\rho}\frac{\partial^2 g}{\partial T \partial p}\left(\frac{\partial g}{\partial p}\right)^{-1}, \quad c_p = -T\frac{\partial^2 g}{\partial T^2}$$

TEOS-10 also introduces Absolute Salinity $S_A$ (mass fraction in g/kg) and Conservative Temperature $\Theta$ as the recommended state variables, replacing Practical Salinity and potential temperature for improved accuracy.

Step 1: Displace a Fluid Parcel

Consider a fluid parcel at depth $z_0$ with density $\rho_{\text{parcel}}(z_0)$ equal to the environmental density $\bar{\rho}(z_0)$. Displace it vertically by a small distance$\delta z$ (upward, so $\delta z > 0$ with z pointing up). The parcel arrives at$z_0 + \delta z$.

$$\text{Parcel at } z_0 + \delta z, \quad \text{environment density: } \bar{\rho}(z_0 + \delta z)$$

Step 2: Adiabatic Adjustment of the Parcel

As the parcel moves to a region of lower pressure (upward displacement), it expands adiabatically. Its density changes according to the adiabatic compressibility. For an incompressible treatment (Boussinesq), the parcel retains its original density. More precisely, the parcel density at the new location is:

$$\rho_{\text{parcel}}(z_0 + \delta z) \approx \bar{\rho}(z_0) - \frac{\rho_0 g}{c_s^2}\,\delta z$$

where $c_s$ is the speed of sound. The second term accounts for the adiabatic density change due to the pressure difference. In the Boussinesq limit, we neglect this compressibility correction.

Step 3: Compute the Buoyancy Force

The net buoyancy force on the displaced parcel (per unit volume) is the difference between the gravitational force on the parcel and the pressure gradient force from the environment. Using Archimedes' principle:

$$F = -g\left[\rho_{\text{parcel}}(z_0 + \delta z) - \bar{\rho}(z_0 + \delta z)\right]$$

The environmental density at the new depth can be Taylor-expanded:$\bar{\rho}(z_0 + \delta z) \approx \bar{\rho}(z_0) + \frac{\partial \bar{\rho}}{\partial z}\delta z$.

Step 4: Substitute and Simplify

In the Boussinesq approximation, the parcel retains its original density $\bar{\rho}(z_0)$. Substituting:

$$F = -g\left[\bar{\rho}(z_0) - \bar{\rho}(z_0) - \frac{\partial \bar{\rho}}{\partial z}\delta z\right] = g\frac{\partial \bar{\rho}}{\partial z}\,\delta z$$

In a stably stratified ocean, $\partial\bar{\rho}/\partial z < 0$ (density increases downward, z positive upward), so the force opposes the displacement -- it is a restoring force.

Step 5: Apply Newton's Second Law

The acceleration of the parcel (per unit mass, dividing by $\rho_0$) gives the equation of motion:

$$\frac{d^2(\delta z)}{dt^2} = \frac{F}{\rho_0} = \frac{g}{\rho_0}\frac{\partial \bar{\rho}}{\partial z}\,\delta z = -N^2\,\delta z$$

where we define:

$$\boxed{N^2 \equiv -\frac{g}{\rho_0}\frac{\partial \bar{\rho}}{\partial z}}$$

Step 6: Oscillatory Solution

The equation $\ddot{\delta z} + N^2\,\delta z = 0$ is a simple harmonic oscillator when $N^2 > 0$(stable stratification):

$$\delta z(t) = A\cos(Nt) + B\sin(Nt)$$

The parcel oscillates about its equilibrium depth with angular frequency $N$ and period$\tau = 2\pi/N$. When $N^2 < 0$ (unstable), the solution is exponentially growing -- convective overturning ensues.

Step 7: Decomposition into T and S Contributions

Using the linearized equation of state, the density gradient can be decomposed into thermal and haline contributions:

$$\frac{\partial \rho}{\partial z} = \frac{\partial \rho}{\partial T}\frac{\partial T}{\partial z} + \frac{\partial \rho}{\partial S}\frac{\partial S}{\partial z} = -\rho_0\alpha\frac{\partial T}{\partial z} + \rho_0\beta\frac{\partial S}{\partial z}$$

Substituting into the definition of $N^2$:

$$N^2 = g\left(\alpha\frac{\partial T}{\partial z} - \beta\frac{\partial S}{\partial z}\right)$$

This decomposition is crucial for understanding double-diffusive processes, where the thermal and haline contributions to stability can have opposite signs.

Step 1: Thermal and Haline Density Gradients

Using the linearized equation of state, define the thermal and haline contributions to the vertical density gradient separately. The thermal contribution to buoyancy frequency is:

$$N_T^2 = g\alpha\frac{\partial \theta}{\partial z}, \quad N_S^2 = -g\beta\frac{\partial S}{\partial z}$$

so that $N^2 = N_T^2 + N_S^2$. When $N_T^2 > 0$, temperature is stabilizing (warm above cold); when $N_S^2 > 0$, salinity is stabilizing (fresh above salty).

Step 2: Define the Density Ratio

The density ratio $R_\rho$ is defined as the ratio of the thermal to haline contributions to the density gradient:

$$R_\rho = \frac{\alpha\,\partial\theta/\partial z}{\beta\,\partial S/\partial z} = \frac{N_T^2}{-N_S^2}$$

When both T and S decrease downward (as in the subtropical thermocline), both numerator and denominator are negative, giving $R_\rho > 0$. If $R_\rho > 1$, the thermal stratification dominates and the column is in the salt-finger regime. If $0 < R_\rho < 1$, salinity dominates and the column is in the diffusive-convective regime.

Step 3: Limitation of $R_\rho$ -- Motivation for the Turner Angle

The density ratio has a singularity when $\beta\,\partial S/\partial z = 0$ (no salinity gradient), which makes it poorly behaved in much of the ocean. Furthermore, $R_\rho$ does not distinguish between different types of instability in a symmetric way. Ruddick (1983) introduced the Turner angle to resolve these issues.

Step 4: Construct the Turner Angle

Define the Turner angle using the four-quadrant arctangent of the sum and difference of the density gradient components:

$$Tu = \arctan\left(\frac{\alpha\,\partial\theta/\partial z + \beta\,\partial S/\partial z}{\alpha\,\partial\theta/\partial z - \beta\,\partial S/\partial z}\right) = \arctan\left(\frac{N_T^2 - N_S^2}{N_T^2 + N_S^2}\right)$$

Note the numerator equals $N_T^2 - N_S^2$ (total buoyancy gradient, proportional to the sum of absolute density contributions) and the denominator equals $N^2$ (the net stratification).

Step 5: Relate Turner Angle to $R_\rho$

The Turner angle and density ratio are related through:

$$\tan(Tu) = \frac{R_\rho + 1}{R_\rho - 1}$$

This can be inverted to give $R_\rho = \frac{\tan(Tu) + 1}{\tan(Tu) - 1}$. The advantage of $Tu$is that it maps the entire range of $R_\rho \in (-\infty, +\infty)$ onto the finite interval$Tu \in (-90°, +90°)$, eliminating the singularity.

Step 6: Stability Regimes in Turner Angle Space

Linear stability analysis of the double-diffusive equations (Stern, 1960; Baines and Gill, 1969) shows that instability occurs when the slower-diffusing component is destabilizing. In terms of Tu:

$$\begin{aligned} 45° < Tu < 90°: &\quad \text{Salt finger regime} \quad (R_\rho > 1) \\ -90° < Tu < -45°: &\quad \text{Diffusive convection} \quad (0 < R_\rho < 1) \\ -45° < Tu < 45°: &\quad \text{Doubly stable (no DD instability)} \\ |Tu| = 90°: &\quad \text{Gravitationally unstable} \end{aligned}$$

The critical angles $\pm 45°$ correspond to $R_\rho = 1$ and $R_\rho = \infty$(or equivalently, one of the two density gradient contributions being zero).

Step 7: Growth Rate of Double-Diffusive Instabilities

From the linearized perturbation equations for temperature and salinity with different diffusivities$\kappa_T$ and $\kappa_S$, the growth rate $\sigma$ of salt fingers satisfies:

$$\sigma = \frac{g\alpha\,\partial\theta/\partial z}{\kappa_T k^2}\left(1 - \frac{1}{R_\rho}\right) \quad \text{for tall, thin fingers } (k \to \infty)$$

where $k$ is the horizontal wavenumber. Growth is fastest when $R_\rho$ is close to 1 (marginally stable). The preferred finger width scales as $\sim (\kappa_T / N)^{1/2} \approx 1$ cm, which matches laboratory observations.

Step 1: Why In-Situ Density Fails for Comparison

Consider two water parcels at different depths with identical T and S. Their in-situ densities differ because of pressure:

$$\rho(S, T, p_1) \neq \rho(S, T, p_2) \quad \text{for } p_1 \neq p_2$$

At 4000 m depth, pressure adds roughly 20 kg/m$^3$ to the density, while the T-S signals we care about are only $\sim 0.01$ to $1$ kg/m$^3$. We need a way to compare densities at a common reference pressure.

Step 2: Adiabatic Lapse Rate from the First Law

For a reversible adiabatic (isentropic) process, the first law of thermodynamics gives:

$$dq = 0 = c_p\,dT - \frac{\alpha T}{\rho}\,dp$$

where $\alpha$ is the thermal expansion coefficient. Rearranging gives the adiabatic lapse rate:

$$\Gamma = \left(\frac{\partial T}{\partial p}\right)_\eta = \frac{\alpha T}{\rho c_p}$$

where $\eta$ denotes constant entropy. For seawater, $\Gamma \approx 1.5 \times 10^{-4}$ K/dbar, meaning temperature increases by about 0.15 degrees C per 1000 m of depth due to adiabatic compression alone.

Step 3: Define Potential Temperature

The potential temperature $\theta$ is the temperature a parcel would have if moved adiabatically to a reference pressure $p_r$ (usually 0 dbar, the surface). It is obtained by integrating the lapse rate:

$$\theta(S, T, p) = T - \int_{p_r}^{p} \Gamma(S, T', p')\,dp'$$

For small pressure changes, $\theta \approx T - \Gamma \cdot (p - p_r)$. The integral must be computed numerically for accuracy because $\Gamma$ depends on T, S, and p. A 4th-order Runge-Kutta scheme is typically used.

Step 4: Define Potential Density

Potential density is the density evaluated at the reference pressure using the potential temperature:

$$\rho_\theta = \rho(S, \theta, p_r)$$

The corresponding density anomaly is:

$$\sigma_\theta = \rho(S, \theta, 0) - 1000 \text{ kg/m}^3$$

By using $\theta$ instead of $T$ and evaluating at $p = 0$, we remove both the direct pressure effect on density and the adiabatic temperature change, isolating the true T-S signature of the water.

Step 5: The Thermobaric Problem

Potential density referenced to the surface ($\sigma_0$) fails for deep water because the thermal expansion coefficient $\alpha$ depends on pressure. Consider two parcels with the same $\sigma_0$but different T-S properties:

$$\sigma_0(\text{parcel A}) = \sigma_0(\text{parcel B}) \quad \text{but} \quad \rho(S_A, T_A, 4000\text{ dbar}) \neq \rho(S_B, T_B, 4000\text{ dbar})$$

This occurs because $\alpha$ increases with pressure: cold water is more compressible than warm water. The thermobaric coefficient $\partial\alpha/\partial p \approx 1.7 \times 10^{-12}$ K$^{-1}$ Pa$^{-1}$quantifies this effect. At 4000 dbar, the density error from using $\sigma_0$ can exceed 0.05 kg/m$^3$ -- larger than many deep-ocean density differences.

Step 6: Patched Potential Density

To mitigate the thermobaric problem, oceanographers use different reference pressures for different depth ranges:

$$\sigma_0 \text{ (ref. 0 dbar)}, \quad \sigma_1 \text{ (ref. 1000 dbar)}, \quad \sigma_2 \text{ (ref. 2000 dbar)}, \quad \sigma_4 \text{ (ref. 4000 dbar)}$$

Each $\sigma_n$ is defined as $\sigma_n = \rho(S, \theta_n, n \times 1000 \text{ dbar}) - 1000$, where $\theta_n$ is the potential temperature referenced to $n \times 1000$ dbar. This "patched" approach reduces the thermobaric error but introduces discontinuities at the boundaries between reference levels.

Step 7: Neutral Density as the Solution

The fundamental issue is that no single reference-pressure potential density can perfectly represent the density ordering at all depths. Neutral density $\gamma^n$ solves this by defining surfaces along which parcels can move without doing work against buoyancy:

$$\nabla \gamma^n \propto -\rho_0\left(\alpha\nabla\theta - \beta\nabla S\right) + \frac{\rho_0 g}{c_s^2}\hat{k}$$

The neutral density gradient is aligned with the locally referenced potential density gradient at every point, avoiding the thermobaric errors of any single reference pressure. It is computed as a nonlinear function fitted to a global hydrographic dataset (Jackett and McDougall, 1997).