Boussinesq Approximation

The Boussinesq approximation (often called the Oberbeck–Boussinesq approximation) is a foundational analytical and numerical simplification in continuum fluid mechanics for buoyancy-driven and density-stratified flows in which fractional density variations are small compared to a reference background density (|ρ′| / ρ₀ ≪ 1). Formulated independently by Anton Oberbeck in 1879 for thermal conduction in liquids and by Joseph Boussinesq in 1903 in his Théorie analytique de la chaleur, the approximation replaces the spatially and temporally varying fluid density ρ(x, t) with a constant reference density ρ₀ in all inertial acceleration terms and horizontal pressure gradients, while retaining the density perturbation ρ′ = ρ − ρ₀ solely in the vertical gravitational body-force term (ρ′g).

Within geophysical fluid dynamics, physical oceanography, and mesoscale meteorology, applying the Boussinesq approximation to the compressible Navier-Stokes equations replaces the compressible continuity equation with the divergence-free incompressibility condition (∇ · u = 0). This mathematical filtering eliminates fast-propagating compressible acoustic (sound) waves from the dynamical system while preserving the full physics of thermal convection, internal gravity waves, Kelvin–Helmholtz shear instabilities, Ekman transport, and planetary Rossby waves.

In a stratified fluid such as the global ocean or the lower troposphere, total fluid density ρ(x, t) and pressure p(x, t) can be decomposed into a constant reference density ρ₀, a horizontally uniform hydrostatically balanced background profile ρ̄(z) and p̄(z), and dynamical perturbations ρ′(x, y, z, t) and p′(x, y, z, t):

where the reference hydrostatic state satisfies dp̄ / dz = −[ρ₀ + Δρ̄(z)]g. Throughout the entire world ocean, potential density varies between ~1,024 kg m⁻³ in warm tropical surface waters and ~1,028 kg m⁻³ in abyssal polar bottom waters, so |Δρ̄ + ρ′| / ρ₀ ≈ 4 × 10⁻³ ≪ 1.

Momentum Equation and the Buoyancy Variable

Substituting the density decomposition into the rotating Navier-Stokes equations yields the momentum balance:

Because |Δρ̄ + ρ′| / ρ₀ ∼ 10⁻³, neglecting the density anomaly in the inertial and Coriolis terms on the left-hand side introduces an error of less than 0.4%. However, in the vertical momentum balance, gravitational acceleration g ≈ 9.81 m s⁻² is three to six orders of magnitude larger than typical vertical fluid accelerations (Dw/Dt ∼ 10⁻⁵–10⁻² m s⁻²). Consequently, the product (ρ′/ρ₀)g is of the exact same order of magnitude as Dw/Dt and −(1/ρ₀)∂p′/∂z, and must be retained. Dividing through by ρ₀ gives the Boussinesq momentum equation:

where ν = μ/ρ₀ is the kinematic viscosity and b is the buoyancy (with dimensions of acceleration, m s⁻²):

For a fluid with a linearized thermal and haline equation of state ρ = ρ₀ [1 − α(T − T₀) + β(S − S₀)]—where α = −(1/ρ₀)(∂ρ/∂T)_p,S is the isobaric thermal expansion coefficient (~2 × 10⁻⁴ K⁻¹ for seawater, 1/T₀ ≈ 3.4 × 10⁻³ K⁻¹ for an ideal gas) and β = (1/ρ₀)(∂ρ/∂S)_p,T is the haline contraction coefficient (~7.6 × 10⁻⁴ kg g⁻¹)—buoyancy becomes b = g[α(T − T₀) − β(S − S₀)].

Acoustic Filtering and Solenoidal Mass Continuity

In the exact continuity equation, (1/ρ) Dρ/Dt + ∇ · u = 0, scaling shows that the Lagrangian rate of fractional density change (1/ρ) Dρ/Dt scales as max(Ma², H/H_ρ) times the velocity divergence ∇ · u, where Ma = U/c_s is the Mach number and H_ρ = |(ρ₀⁻¹ dρ̄/dz)⁻¹| is the density scale height (~200 km in the ocean, ~8.5 km in the atmosphere). When Ma² ≪ 1 and vertical motion scales H ≪ H_ρ, the continuity equation reduces to the divergence-free (solenoidal) constraint:

Meanwhile, density anomalies evolve via a prognostic thermodynamic advection–diffusion equation for buoyancy (or potential temperature θ):

where κ is the thermal or molecular diffusivity. Because ∇ · u = 0 eliminates temporal pressure derivatives from the continuity equation, compressible acoustic waves (c_s ≈ 1,500 m s⁻¹ in seawater; ~340 m s⁻¹ in air) are filtered out (effectively setting c_s → ∞). In numerical general circulation models, this acoustic filtering removes the restrictive Courant–Friedrichs–Lewy (CFL) time-step penalty of sound waves, allowing integration time steps orders of magnitude longer.

Baroclinic Vorticity Generation

Taking the curl of the Boussinesq momentum equation reveals how buoyancy gradients generate vorticity ω = ∇ × u (extending the barotropic vorticity equation): because −∇(1/ρ₀ p′) has zero curl, the baroclinic torque (−1/ρ²) ∇ρ × ∇p simplifies under the Boussinesq approximation to ∇b × k̂ = (∂b/∂y, −∂b/∂x, 0). Thus, any horizontal gradient of temperature or salinity (∂b/∂x ≠ 0 or ∂b/∂y ≠ 0) directly generates horizontal overturning vorticity—the fundamental mechanism driving land–sea breezes, atmospheric fronts, and the oceanic thermal wind balance (∂u_g/∂z = (f⁻¹) k̂ × ∇_hb).

When a fluid is stably stratified, density decreases with height (dρ̄/dz < 0, or equivalently background buoyancy increases with height, db̄/dz > 0). If a fluid parcel at equilibrium height z₀ is displaced vertically upward by a small distance ζ without mixing, its buoyancy anomaly relative to the new surroundings is b′ = −(db̄/dz)ζ. Neglecting pressure perturbations and viscosity for a vertical parcel oscillation yields the harmonic oscillator equation:

Here N (measured in rad s⁻¹) is the Brunt–Väisälä frequency (or buoyancy frequency), discovered independently by David Brunt (1927) and Vilho Väisälä (1925). For a compressible atmosphere, ρ̄ is replaced by potential temperature θ̄ = T(p₀/p)^R/c_p, giving N² = (g/θ̄)(dθ̄/dz). In the oceanic seasonal pycnocline and lower stratosphere, N ≈ 10⁻² rad s⁻¹ (oscillation period 2π/N ≈ 10 minutes); in the abyssal ocean, N ≈ 10⁻³ rad s⁻¹ (period ~1.7 hours). When N² < 0 (heavy fluid overlying light fluid), the stratification is statically unstable and perturbations grow exponentially at rate |N| via convective overturning.

Dispersion Relation of Internal Gravity Waves

Linearizing the inviscid, non-rotating Boussinesq equations around a rest state with uniform stratification N² and seeking plane-wave solutions proportional to exp[i(kx + ly + mz − ωt)]—where k = (k, l, m) is the three-dimensional wavevector and k_h = √(k² + l²) is the horizontal wavenumber—yields the Boussinesq internal gravity wave dispersion relation:

where ϑ is the angle between the three-dimensional wavevector k and the horizontal plane. This relation exhibits three remarkable physical properties:

1. Frequency Cutoff: Internal gravity waves can only propagate at frequencies below the buoyancy frequency (0 < ω ≤ N, or f ≤ ω ≤ N when planetary rotation f is included).
2. Angle Selected solely by Frequency: The propagation angle ϑ = arccos(ω/N) depends purely on the oscillation frequency ω and is independent of the wavelength 2π/|k|.
3. Orthogonality of Phase and Group Velocity: Because ∇ · u = 0 requires ik · u = 0, fluid parcels oscillate parallel to the wave crests (u ⊥ k). Differentiating ω(k) with respect to k shows that the group velocity vector c_g = ∇_kω (which carries wave energy) is strictly perpendicular to the phase velocity vector c_p = (ω/|k|²)k (c_g · c_p = 0) and has an opposite vertical component: when wave crests (phase) propagate downward (m < 0), wave energy (c_g) propagates upward into the middle atmosphere or upper ocean.

The Rayleigh and Prandtl Numbers

When a horizontal fluid layer of thickness H is heated from below (T_bottom = T₀ + ΔT) and cooled from above (T_top = T₀), thermal expansion creates an unstable top-heavy density gradient (N² < 0). For small ΔT, viscous friction (ν) and thermal diffusion (κ) stabilize the motionless conductive state. Lord Rayleigh (1916) performed a linear stability analysis of the Oberbeck–Boussinesq equations and demonstrated that convective overturning cells (Rayleigh–Bénard convection) erupt when the dimensionless Rayleigh number (Ra) exceeds a critical threshold Ra_c:

The Rayleigh number measures the ratio of destabilizing buoyancy forcing to stabilizing viscous and thermal diffusion, while the Prandtl number (Pr = ν/κ) measures the ratio of momentum diffusivity to thermal diffusivity. For stress-free horizontal boundaries, Rayleigh proved analytically that Ra_c = 27π⁴ / 4 ≈ 657.5 (at horizontal wavenumber k_c = π/(√2 H)); for rigid no-slip boundaries (relevant to laboratory experiments), Jeffreys (1928) found Ra_c ≈ 1,708 (with k_c ≈ 3.117/H). In Earth's silicate mantle (H ≈ 2,900 km, Pr ≈ 10²³, Stokes flow regime), Ra ≈ 10⁷ ≈ 10⁴ Ra_c, driving vigorous mantle plumes and plate tectonics.

Deep Atmospheric Limits and the Anelastic Approximation

While the Boussinesq approximation is extraordinarily accurate throughout the entire depth of the ocean (H ≈ 4 km ≪ H_ρ ≈ 200 km), in the atmosphere air density drops by ~70% between the surface and the tropopause (H ≈ 10 km > H_ρ ≈ 8.5 km). Consequently, the Boussinesq approximation is strictly valid in meteorology only for shallow boundary-layer phenomena (H ≲ 1–2 km). For deep atmospheric cumulonimbus convection and tropospheric gravity waves, Yoshimitsu Ogura and Norman Phillips (1962) developed the anelastic approximation, which retains the vertically varying reference density profile ρ̄(z) inside the continuity constraint:

Like the Boussinesq system, the anelastic equations filter out fast acoustic waves while accurately capturing the adiabatic expansion and vertical velocity amplification (w ∝ ρ̄(z)^−1/2) of air parcels ascending through multiple atmospheric scale heights.

How can a fluid satisfy ∇ · u = 0 (incompressibility) while still having density variations that drive convection?

In the Boussinesq regime, the thermal expansion coefficient α is small (~2 × 10⁻⁴ K⁻¹ in water), so a temperature difference of ΔT = 10 K changes fluid density by only 0.2% (|ρ′|/ρ₀ = 0.002). In the mass continuity equation (∇ · u = −(1/ρ) Dρ/Dt), this 0.2% volume change makes a negligible contribution to the kinematic velocity balance compared to spatial velocity gradients ∂u/∂x and ∂w/∂z, allowing ∇ · u = 0 to hold to 99.8% accuracy. Yet when that same 0.2% density anomaly is multiplied by Earth's large gravitational acceleration (g = 9.81 m s⁻²), it produces a buoyancy acceleration b = gαΔT ≈ 0.02 m s⁻², which is large enough to drive powerful convective plumes over minutes to hours.

Why do internal gravity waves propagate at an angle rather than horizontally like ocean surface waves?

Surface water waves are trapped on a two-dimensional interface (the sea surface) and can only propagate horizontally. Internal gravity waves propagate through a continuously stratified three-dimensional fluid interior where restoring buoyancy acts strictly in the vertical direction (k̂), while the Boussinesq continuity condition (∇ · u = 0) requires fluid motions to form closed loops in the plane spanned by the wavevector k and the vertical axis. Balancing the vertical buoyancy restoring force against the total three-dimensional inertia of the tilted parcel oscillation introduces the geometric factor cos ϑ = k_h/|k|, locking the wave propagation angle ϑ directly to the ratio of wave frequency to the Brunt–Väisälä frequency (ω = N cos ϑ).

What is the non-Boussinesq Steric Sea-Level Effect in climate modeling?

Because a Boussinesq ocean model enforces strictly zero volume divergence (∇ · u = 0) everywhere, the total volume of the global ocean in the model remains constant even when the ocean absorbs anthropogenic heat and warms. In reality, warming seawater expands thermally (thermosteric sea-level rise, accounting for ~40% of observed modern global mean sea-level rise). Ocean general circulation models using the Boussinesq approximation compensate for this by diagnosing the spatially integrated steric density change post-hoc and adding a uniform global barotropic sea-surface height correction (the Greatbatch correction).

1. Boussinesq, J. (1903). Théorie analytique de la chaleur, mise en harmonie avec la thermodynamique et avec la théorie mécanique de la lumière (Vol. 2). Gauthier-Villars, Paris. doi:10.3931/e-rara-105749
2. Spiegel, E. A., & Veronis, G. (1960). "On the Boussinesq approximation for a compressible fluid." The Astrophysical Journal, 131, 442–447. doi:10.1086/146849
3. Rayleigh, Lord (1916). "On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side." The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 32(192), 529–546. doi:10.1080/14786441608635602
4. Vallis, G. K. (2017). Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation (2nd ed.). Cambridge University Press. doi:10.1017/9781107588417