Gent–McWilliams Isopycnal Eddy Parameterization

Gent–McWilliams Isopycnal Eddy Parameterization (commonly abbreviated as the GM parameterization or GM90) is a mathematical closure formulated in 1990 by Peter R. Gent and James C. McWilliams to represent the adiabatic, buoyancy-conserving transport of heat, salinity, and passive tracers by unresolved mesoscale baroclinic eddies in coarse-resolution ocean general circulation models (OGCMs). Because the first baroclinic Rossby radius of deformation in the stratified ocean ranges from roughly 30–50 km in the subtropics to less than 10 km at high latitudes, climate-scale ocean models operating on horizontal grid spacings of 50 to 100 km cannot explicitly resolve the mesoscale eddies generated by baroclinic instability. Prior to Gent and McWilliams's formulation, z-coordinate ocean models parameterized subgrid-scale mixing via horizontal Laplacian diffusion, which artificially mixed heat and salt across tilted constant-density (isopycnal) surfaces and generated spurious diapycnal upwelling known as the Veronis effect. By introducing a non-divergent, eddy-induced advection velocity—termed the bolus velocity—driven by the downgradient diffusion of isopycnal layer thickness, the Gent–McWilliams scheme acts as an adiabatic sink of large-scale available potential energy (APE) that slumps tilted density surfaces toward the horizontal without inducing spurious cross-isopycnal mixing.

Spurious Diapycnal Diffusion in Depth-Coordinate Ocean Models

In early three-dimensional ocean general circulation models based on the Bryan–Cox–Semtner fixed-depth (z-coordinate) architecture, subgrid-scale eddy fluxes of potential temperature (θ) and salinity (S) were represented as downgradient horizontal diffusion with a constant horizontal diffusivity A_H. In the real ocean, however, baroclinic fronts such as the Gulf Stream, the Kuroshio Extension, and the Southern Ocean fronts along the Antarctic Circumpolar Current exhibit steeply sloping isopycnal surfaces maintained in thermal-wind balance with vertical shear and the Coriolis effect.

As demonstrated by George Veronis in 1975, applying purely horizontal diffusion across steeply tilted isopycnals in regions such as western boundary currents forces a strong fictitious flux of buoyancy perpendicular to density surfaces (diapycnal mixing). To balance this spurious diapycnal buoyancy gain in the steady-state thermodynamic budget, coarse-resolution models developed an unphysical vertical circulation cell—the Veronis effect—characterized by exaggerated upwelling on the shoreward flank of western boundary currents and unrealistic warming of the deep abyssal ocean.

Redi Isopycnal Rotation Tensor and Its Dynamical Incompleteness

To eliminate spurious cross-isopycnal diffusion, Martha H. Redi (1982) and Lawrence J. Cox (1987) rotated the subgrid-scale tracer diffusion tensor from horizontal-vertical coordinates into local isopycnal-diapycnal coordinates aligned with the neutral density surface. Under the small-slope approximation—valid throughout the ocean interior where the horizontal-to-vertical aspect ratio of isopycnal slopes is of order 10⁻³ to 10⁻²—define the two-dimensional isopycnal slope vector S in terms of the local horizontal and vertical gradients of potential density ρ:

S = (S_x, S_y) = − ∇_hρ / (∂ρ / ∂z)

where ∇_h = (∂/∂x, ∂/∂y) denotes the horizontal gradient operator and ∂ρ/∂z < 0 in a stably stratified water column. While the Redi tensor successfully restricts passive tracer stirring to lie along constant-density surfaces (so that the symmetric Redi flux of density itself vanishes identically, J_Redi(ρ) = 0), it leaves the mean density field ρ(x, y, z, t) and the large-scale pressure field unaltered. Consequently, Redi isopycnal diffusion alone cannot represent the primary dynamical role of baroclinic instability in two-layer baroclinic turbulence: namely, the adiabatic slumping of tilted isopycnal surfaces that converts mean available potential energy into eddy kinetic energy.

Downgradient Diffusion of Isopycnal Layer Thickness

Gent and McWilliams (1990) reasoned from the continuity equation formulated in isopycnal coordinates—where the mass thickness h = −∂z/∂ρ of an infinitesimal density layer is conserved following adiabatic flow—that mesoscale baroclinic eddies mix layer thickness downgradient along isopycnal surfaces with a positive thickness diffusivity κ_GM (units of m² s⁻¹). When transformed back into standard Cartesian height coordinates (x, y, z), this downgradient isopycnal thickness flux produces an effective three-dimensional eddy-induced transport velocity (or bolus velocity) v = (u, v, w) that advects both active buoyancy (potential temperature and salinity) and all passive biogeochemical tracers.

Under the small-slope approximation, the horizontal and vertical components of the Gent–McWilliams bolus velocity take the form:

u = − ∂(κ_GM S_x) / ∂z,   v = − ∂(κ_GM S_y) / ∂z,   w^* = ∇_h · (κ_GM S) = ∂(κ_GM S_x)/∂x + ∂(κ_GM S_y)/∂y

Because the bolus velocity is derived from the curl of a vector streamfunction Ψ^* = (−κ_GM S_y, κ_GM S_x, 0) = k × (κ_GM S), it is identically non-divergent in three dimensions:

∇ · v = ∂u/∂x + ∂v/∂y + ∂w/∂z = 0

At the sea surface (z = 0) and at the solid ocean floor (z = −H), the kinematic boundary condition requires that the eddy-induced streamfunction vanish (Ψ^* = 0), ensuring zero integrated mass transport across any closed basin boundary:

∫_−H⁰ u dz = 0,   ∫_−H⁰ v dz = 0

Tracer Conservation under Residual-Mean Advection

When the Gent–McWilliams bolus velocity v is combined with the resolved Eulerian-mean velocity v = (u, v, w), any conservative tracer τ (including potential temperature θ, salinity S, and potential density ρ) satisfies the advection-diffusion conservation equation governed by the residual-mean velocity v_res = v + v:

∂τ / ∂t + (v + v^*) · ∇τ = ∇ · (K_Redi ∇τ) + ∂(κ_dia ∂τ/∂z) / ∂z

where K_Redi is the symmetric along-isopycnal Redi stirring tensor (characterized by the isopycnal tracer diffusivity κ_Redi) and κ_dia is the small background diapycnal diffusivity driven by internal wave breaking. Because the GM transport enters the density equation purely as an advection by the incompressible velocity v, it preserves the volume of every density class in the ocean: the volume integral of any arbitrary function F(ρ) over a closed domain with insulating boundaries is invariant under v advection, guaranteeing that all water-mass transformations remain adiabatic.

Available Potential Energy Conversion and Domain-Integrated Sink

A foundational physical requirement of any parameterization of baroclinic instability is that it mimic the baroclinic conversion pathway by extracting available potential energy (APE) from the large-scale stratification. Multiplying the eddy-advected density tendency ∂ρ/∂t = − v^* · ∇ρ by gravitational acceleration g and vertical coordinate z and integrating over the global ocean volume V yields the domain-integrated rate of change of gravitational potential energy E_p = ∫∫∫ ρ g z dV:

(dE_p / dt)_GM = ∫∫∫_V g z (v^* · ∇ρ) dV = − ∫∫∫_V g κ_GM (− ∂ρ/∂z) |S|² dV ≤ 0

Because a stably stratified water column satisfies −∂ρ/∂z > 0 and κ_GM > 0, the integrand is non-negative everywhere and strictly positive wherever isopycnal surfaces are sloped (|S| > 0). Consequently, the Gent–McWilliams parameterization guarantees a monotonic sink of large-scale potential energy at a local volumetric rate proportional to ρ₀ κ_GM N² |S|², where N² = −(g/ρ₀) ∂ρ/∂z is the squared Brunt–Väisälä buoyancy frequency. Geometrically, wherever wind-driven Ekman pumping or surface buoyancy forcing steepens isopycnal surfaces, the bolus circulation generates an overturning cell that lifts dense water on the deep side of the front and depresses light water on the shallow side, slumping the isopycnals toward a flat, minimum-potential-energy state.

The Griffies (1998) Antisymmetric Skew-Flux Tensor

In numerical ocean models, computing the bolus velocity v explicitly by differentiating the noisy isopycnal slope vector S vertically and then advecting tracers can introduce grid-scale computational modes. In 1998, Stephen M. Griffies demonstrated that the advective transport of a tracer τ by the non-divergent bolus velocity v is mathematically equivalent to the divergence of an antisymmetric skew-diffusive flux, because v · ∇τ = ∇ · (Ψ × ∇τ). Under the small-slope approximation, the combined Redi isopycnal stirring and Gent–McWilliams eddy advection can be written as a single three-dimensional flux divergence ∇ · (J_Redi+GM), where J_Redi+GM = − K_total ∇τ and the unified 3 × 3 diffusivity tensor decomposes into symmetric (Redi) and antisymmetric (GM) parts:

When the isopycnal tracer diffusivity equals the thickness diffusivity (κ_Redi = κ_GM), the upper-right off-diagonal entries of the symmetric and antisymmetric matrices cancel, eliminating four derivative terms per tracer step and improving numerical stability near steep fronts and the surface mixed layer.

Spatially and Temporally Varying κ_GM Parameterizations

In their original 1990 formulation, Gent and McWilliams prescribed a spatially uniform thickness diffusivity κ_GM on the order of 600 to 1,000 m² s⁻¹. However, satellite altimetry and mooring observations reveal that mesoscale eddy diffusivity varies by more than an order of magnitude across the global ocean—exceeding 2,000 m² s⁻¹ in the core of western boundary current extensions and the Antarctic Circumpolar Current while dropping below 200 m² s⁻¹ in quiescent subtropical gyre interiors.

Building on the Held–Larichev baroclinic eddy diffusivity scaling theory and the linear Eady growth rate σ_Eady = f / √Ri = M² / N (where M² = |(g/ρ₀) ∇_hρ| is the horizontal buoyancy gradient and Ri = N² f² / M⁴ is the balanced Richardson number), Martin Visbeck, John Marshall, Thomas Haine, and Mark Spall (1997) proposed a flow-dependent closure in which κ_GM scales with the vertical average of the baroclinic growth rate and a characteristic eddy mixing length ℓ_eddy:

κ_GM(x, y) = α_V ⟨ M² / N ⟩_z ℓ_eddy²

where α_V ≈ 0.015 is a dimensionless proportionality constant and ⟨·⟩_z denotes a depth average across the main thermocline. Subsequent vertical-structure extensions by Ferrari and colleagues (2008, 2010) replaced the algebraic slope-clipping schemes needed near unstratified surface mixed layers (where N² → 0 causes S to diverge) with a boundary-value problem for the vertical structure of Ψ^* based on baroclinic Rossby wave eigenmodes, smoothly tapering the eddy streamfunction to zero at the sea surface and seafloor.

Transformation of Global Climate and Ocean Modeling

The adoption of the Gent–McWilliams parameterization in the mid-1990s (first demonstrated in the NCAR ocean model by Danabasoglu, McWilliams, and Gent in 1994) transformed global ocean and coupled climate modeling:

1. Elimination of the Veronis Effect and Deep Warming Drift: By replacing horizontal diffusion with adiabatic bolus advection, coarse-resolution OGCMs eliminated fictitious diapycnal upwelling in western boundary currents, cooling the simulated abyssal ocean by 2–3 °C toward observed hydrographic climatologies and stabilizing multi-century coupled climate simulations without artificial heat-flux adjustments.
2. Southern Ocean Overturning and the Deacon Cell: In the latitude band of Drake Passage (55°S–62°S), Eulerian-mean diagnostics of z-coordinate models previously showed a deep wind-driven overturning loop—the Deacon Cell—that appeared to sink across density surfaces in the subantarctic zone. Incorporating the Gent–McWilliams bolus velocity revealed that the southward eddy-induced circulation Ψ^* nearly cancels the northward wind-driven Ekman transport cell Ψ_Ekman, forming the theoretical foundation for the Marshall–Radko residual-mean theory of the Antarctic Circumpolar Current and modern representations of the global thermohaline circulation and Atlantic Meridional Overturning Circulation.

1. Gent, P. R., & McWilliams, J. C. (1990). Isopycnal mixing in ocean circulation models. Journal of Physical Oceanography, 20(1), 150–155. https://doi.org/10.1175/1520-0485(1990)020<0150:IMIOCM>2.0.CO;2020<0150:IMIOCM>2.0.CO;2)
2. Gent, P. R., Willebrand, J., McDougall, T. J., & McWilliams, J. C. (1995). Parameterizing eddy-induced tracer transports in ocean circulation models. Journal of Physical Oceanography, 25(4), 463–474. https://doi.org/10.1175/1520-0485(1995)025<0463:PEITTI>2.0.CO;2025<0463:PEITTI>2.0.CO;2)
3. Danabasoglu, G., McWilliams, J. C., & Gent, P. R. (1994). The role of mesoscale tracer transports in the global ocean circulation. Science, 264(5162), 1123–1126. https://doi.org/10.1126/science.264.5162.1123
4. Visbeck, M., Marshall, J., Haine, T., & Spall, M. (1997). Specification of eddy transfer coefficients in coarse-resolution ocean circulation models. Journal of Physical Oceanography, 27(3), 381–402. https://doi.org/10.1175/1520-0485(1997)027<0381:SOETCI>2.0.CO;2027<0381:SOETCI>2.0.CO;2)
5. Griffies, S. M. (1998). The Gent–McWilliams skew flux. Journal of Physical Oceanography, 28(5), 831–841. https://doi.org/10.1175/1520-0485(1998)028<0831:TGMSF>2.0.CO;2028<0831:TGMSF>2.0.CO;2)