SciencePediaThe Earth's oceans are not uniform bodies of water but are carefully layered, or stratified, with lighter water on top of denser water. The slow, energy-intensive process of mixing these layers, known as diapycnal mixing, is a fundamental driver of global ocean circulation and long-term climate. However, when scientists attempt to replicate this process in numerical climate models, they encounter a significant challenge: the emergence of a "phantom" mixing that does not exist in the real world. This artifact, known as spurious diapycnal mixing, arises from the inherent conflict between the ocean's natural physics and the artificial structure of a model's grid. It represents a critical knowledge gap that can compromise the accuracy of climate predictions.
This article provides a comprehensive exploration of this pivotal issue in ocean modeling. By examining its causes, consequences, and the ingenious solutions developed to overcome it, you will gain a deeper understanding of the complexities behind modern climate simulation. The following chapters will first delve into the "Principles and Mechanisms" of how spurious mixing is born from numerical choices and fundamental ocean physics. We will then explore the "Applications and Interdisciplinary Connections," revealing how modelers combat this artifact and why doing so is essential for accurately forecasting phenomena from deep ocean currents to El Niño.
To understand the challenge of modeling our oceans, we must first appreciate their fundamental structure. Imagine a giant, planetary-scale layer cake. The ocean is not a uniform tub of water; it is meticulously stratified, with lighter, warmer water sitting stably atop denser, colder, saltier water. This stratification is the ocean's natural state of rest. Mixing these layers, much like trying to stir oil and water, requires a tremendous amount of energy. The surfaces separating these layers, where the density is constant, are called isopycnals. In a perfectly quiet ocean, water would flow almost exclusively along these surfaces, finding it nearly impossible to cross them. This slow, difficult process of crossing the layers is called diapycnal mixing, and it plays a crucial role in setting the long-term climate of our planet.
Now, imagine you are tasked with building a miniature, digital version of this ocean inside a computer. This is the world of the climate modeler.
A numerical model must divide the smooth, continuous ocean into a finite number of discrete blocks, or grid cells, much like a digital photograph is composed of pixels. The most intuitive way to do this is to use a simple Cartesian grid, a scaffold of horizontal and vertical lines—what oceanographers call a z-level model. Herein lies a profound difficulty. In the real ocean, especially in regions with active currents and eddies, isopycnal surfaces are not flat; they slope and undulate. But our model's grid lines are rigidly horizontal and vertical.
The consequence is a fundamental mismatch: the natural pathways of ocean flow are not aligned with the artificial highways of the numerical grid. This geometric conflict is the birthplace of a phantom that haunts ocean modeling: spurious diapycnal mixing. It is a form of mixing that doesn't happen in the real world but appears in the model as an artifact of our numerical choices.
This phantom mixing arises from several sources, all rooted in the same fundamental misalignment between physics and the grid.
First, let's consider how models represent the small-scale turbulence that we cannot explicitly resolve. A common approach is to add a diffusion term to the equations, which acts to smooth out sharp gradients, much like a drop of ink spreads out in a glass of water. Modelers often apply a large "horizontal diffusion" to represent the vigorous stirring by eddies along isopycnals, and a very small "vertical diffusion" for the weak diapycnal mixing. But what does "horizontal" mean to a computer model built on a z-level grid? It means along the model's horizontal grid lines.
When isopycnals are tilted, these horizontal grid lines cut directly through them. Consequently, an instruction to "mix horizontally" becomes an unintended instruction to mix water of different densities. This is not a physical process, but a numerical error—a ghost in the machine. We can even write down the magnitude of this error. If a model has a physical vertical diffusivity and a horizontal diffusivity , the total effective diapycnal diffusivity, , on an isopycnal with slope turns out to be approximately:
The first term, , is the real, physical mixing we intended. The second term, , is the spurious contribution. Since the horizontal diffusivity is typically thousands or even millions of times larger than the physical diapycnal diffusivity , even a minuscule slope ( on the order of ) can generate spurious mixing that completely overwhelms the real signal.
The same problem plagues the way models handle advection—the simple act of moving a tracer with the flow. The mathematical schemes used to compute this movement are not perfect. Many simple schemes, like the "first-order upwind" method, have an unfortunate side effect: they implicitly introduce a small amount of numerical diffusion along the grid axes to maintain stability. Just like the explicit diffusion operator, this numerical diffusion acts along grid lines and, when they cross tilted isopycnals, creates spurious diapycnal mixing.
Even the fundamental laws of motion can be corrupted. When models try to represent sloped ocean floors with a series of "stair-steps," they can introduce errors in the calculation of the Pressure Gradient Force. This can create phantom currents that don't exist in reality but flow up and down the grid's staircases, artificially pushing water across density layers.
The problems we've discussed so far seem to suggest a simple solution: if the grid is the problem, why not build a model whose grid layers are themselves aligned with the ocean's density surfaces? This is the idea behind isopycnal coordinate models. But this approach leads us to a deeper, more beautiful question: what, precisely, is a density surface?
The density of seawater is a tricky thing. It depends not only on temperature and salinity, but also on pressure. As a parcel of water sinks, it is compressed, and its density increases, even if no heat or salt is exchanged. To get around this, oceanographers invented potential density (), which is the density a parcel would have if it were moved adiabatically to a reference pressure (usually the sea surface, ). For decades, it was believed that ocean mixing happens primarily along these surfaces of constant potential density.
But nature has a subtle and wonderful surprise in store, an effect known as thermobaricity. In simple terms, the way temperature affects density changes with pressure. The coefficient of thermal expansion, , which tells us how much water expands when heated, is itself a function of pressure—it gets larger in the deep ocean. The consequence is startling: a surface of constant potential density (referenced to the surface) is not a truly "neutral" surface at depth. If you were to take a parcel of water in the deep ocean and slide it along a surface of constant potential density, it could become buoyant or heavy relative to its new surroundings. This means that mixing along potential density surfaces—our supposedly "correct" physical pathway—can itself introduce an artificial buoyancy flux, another form of spurious mixing.
This discovery forces us to ask an even more fundamental question: If potential density surfaces are flawed, what is the true surface of neutral buoyancy along which water mixes? Such a surface is called a neutral surface. It is defined at any point as the plane where a small displacement of a water parcel results in no change in its buoyancy.
And here, we arrive at one of the most elegant and frustrating truths in physical oceanography: in general, it is mathematically impossible to connect these local neutral planes to form a single, globally consistent set of surfaces. The reason lies in the complex, nonlinear equation of state for seawater. The vector field that defines the neutral direction has a mathematical property known as non-zero curl. This means that if you start on what you think is a "level" neutral path and follow it on a long journey around an ocean basin, you can arrive back at your starting longitude and latitude, but at a different "level"—you will have spiraled up or down.
There is no perfect, global map of neutral surfaces. They are inherently "slippery". This isn't a failure of our models; it's a fundamental property of the real ocean.
Faced with these challenges—misaligned grids, imperfect numerics, and the non-existence of perfect mixing surfaces—ocean modelers have developed a suite of remarkably clever strategies. They have learned to tame the phantom, even if they cannot completely exorcise it.
Rather than diffusing along rigid horizontal lines, models can be taught to calculate the local slope of the isopycnal and rotate the diffusion tensor to act along that slope.
Parameterizations like the Gent-McWilliams (GM) scheme don't try to mix along isopycnals directly. Instead, they simulate the primary effect of eddies, which is to release potential energy by causing tilted isopycnals to slump back toward being flat. This is an adiabatic process that reduces the very slopes that are the source of so much spurious mixing.
To deal with the problem of thermobaricity, scientists have developed approximate neutral density variables (like ) that provide a much better, though still imperfect, guide for mixing in the deep ocean than potential density does.
And finally, we can even diagnose the magnitude of the phantom's influence. By tracking the variance of a passive tracer on a density surface, we can precisely calculate how much of its decay is due to explicit, physical diffusion. Any residual decay beyond that must be the work of the phantom—the spurious diapycnal mixing generated by the model's numerics.
The story of spurious diapycnal mixing is a perfect example of the scientific process. It is a journey that begins with a practical engineering problem—how to build a simulation—and leads to deep and beautiful insights into the fundamental physics of the ocean itself. It teaches us that even our most basic concepts, like a "level surface," are filled with a surprising and wonderful complexity.
Having grappled with the fundamental principles of spurious diapycnal mixing, we now arrive at a fascinating question: what do we do about it? If our numerical models, by their very design, are prone to creating fictitious mixing that corrupts the physics, how can we ever hope to simulate the Earth’s oceans with any fidelity? This is not merely a technical problem for computer scientists; it is a central challenge in physical oceanography and climate science. The solutions reveal a beautiful interplay between physical intuition, mathematical ingenuity, and a deep understanding of the phenomena we seek to capture.
This journey from a numerical artifact to its real-world consequences is a story in three parts. First, we will explore the clever schemes modelers have invented to "teach" a model to respect the physics of stratification, even when its grid does not. Second, we will see how these ideas have led to entirely new ways of building model grids. And finally, we will witness why this entire endeavor is so critical, by connecting these seemingly esoteric details to climate phenomena that shape our world.
Imagine the task of describing the vast, complex, three-dimensional ocean on a finite computer grid. You must make a choice about how to slice it up in the vertical. Do you use simple, horizontal layers of fixed depth, like floors in a building? This is the essence of a -coordinate model. Or perhaps you use flexible layers that stretch and squeeze to follow the contours of the seafloor? This is a terrain-following or -coordinate model. Or maybe, in a stroke of physical intuition, you decide your layers should follow the natural stratification of the ocean itself, surfaces of constant density? This is an isopycnal coordinate model.
Each choice represents a different "worldview," and each comes with profound trade-offs. The simple -coordinate grid is computationally straightforward, and the all-important pressure gradient force is easy to calculate. But as we've seen, its rigid horizontal levels cut across sloping density surfaces, leading to severe spurious diapycnal mixing. The terrain-following -grid brilliantly represents the complex shape of the seafloor, which is vital for understanding bottom currents. However, over steep undersea mountains and continental slopes, the mathematical transformation required to do this can introduce large errors in the pressure gradient calculation, creating powerful phantom currents where none should exist. The isopycnal grid seems like the perfect solution—by definition, motion along its layers is adiabatic, and spurious mixing is dramatically reduced. But what happens in the turbulent surface mixed layer, where density is nearly uniform? The coordinate system breaks down. And how does it handle the ocean bottom, which is certainly not a surface of constant density?
This dilemma—that no single coordinate system is perfect everywhere—has spurred decades of innovation. If you cannot build a perfect grid, perhaps you can invent a way to correct for its flaws.
The most common ocean models today use -coordinates for their simplicity and robustness. The challenge, then, is to counteract their tendency to generate spurious mixing. This has led to the development of "subgrid-scale parameterizations," which are essentially rules added to the model's code to represent physical processes that the grid cannot resolve. Two of the most important are the Redi and Gent-McWilliams (GM) schemes.
Rotating the Compass: Isoneutral Diffusion
In the real ocean, small-scale turbulence mixes tracers. This mixing is highly anisotropic: it happens thousands or millions of times more easily along isopycnal surfaces than across them. A simple -coordinate model applying diffusion horizontally and vertically gets this spectacularly wrong, as the "horizontal" direction on the grid is not the "isopycnal" direction in the ocean.
The Redi isoneutral diffusion scheme is a clever fix. It tells the model: before you calculate the diffusive mixing, first figure out the local orientation of the true isopycnal surface. Then, project the tracer gradients onto that surface and apply a large diffusion coefficient along it, and a tiny one across it. In essence, you build an anisotropic diffusion tensor whose principal axes are aligned with the physics, not the grid. This ensures that the diffusive flux of buoyancy, , has a near-zero component across isopycnals, suppressing the artificial mixing that would otherwise plague the model.
The Ghost in the Machine: Eddy-Induced Advection
But the story doesn't end with diffusion. In the ocean, a significant portion of tracer transport is done by mesoscale eddies—swirling vortices of water tens to hundreds of kilometers across. These eddies are too small to be resolved by typical global climate models. So, their effect must also be parameterized.
Crucially, these eddies are born from instabilities that feed on the potential energy stored in sloping isopycnals. As they churn the ocean, they tend to flatten these slopes, releasing energy. Like a ball rolling downhill, this is a fundamentally advective process—it moves water parcels from one place to another. And, because these eddies are large-scale motions, they are largely adiabatic; they stir water along isopycnals.
The Gent-McWilliams (GM) parameterization is designed to mimic this effect. It introduces a fictitious "eddy-induced velocity," often called a "bolus velocity" , into the tracer advection equation. This velocity is mathematically constructed to be non-divergent (so it conserves mass) and to be directed in a way that flattens isopycnal slopes. Its magnitude is proportional to the slope itself—the steeper the slope, the stronger the eddy-induced transport.
It is critical to understand that GM and Redi represent different physics. Redi is a symmetric, diffusive operator; it mixes tracers and reduces their variance. GM is an antisymmetric, advective operator; it rearranges tracers without destroying their variance. Redi represents the stirring and filamentation of tracers along isopycnals, while GM represents the large-scale slumping of the density field that releases available potential energy. Together, they provide a powerful one-two punch to represent the physics of unresolved eddies in a way that respects the ocean's stratification.
While parameterizations are powerful, an alternative philosophy is to improve the grid itself. This has led to hybrid coordinate models, which are masterpieces of pragmatism. These models use -coordinates near the surface, where they are best suited to handle the complex interactions with the atmosphere and the well-mixed layer. But deeper down, in the vast, stably stratified ocean interior, the model transitions to an isopycnal coordinate system. This design plays to the strengths of each approach, creating a more physically faithful representation of the full water column.
Even in these advanced models, the quest for perfection continues. The coordinate surfaces might not align perfectly with the true neutral surfaces of the ocean. How much misalignment is tolerable? A revealing calculation shows that for a typical large horizontal mixing coefficient of and a target physical diapycnal mixing of , the maximum allowable angle between the model's grid surface and the true isopycnal is given by:
This angle is minuscule—about 0.006 degrees! This demonstrates the extraordinary precision required in modern ocean modeling. If the alignment is off by any more than this, the numerical error will overwhelm the subtle but climatically crucial physical mixing processes.
The story gets even deeper. The very definition of an "isopycnal" or "neutral" surface depends on the equation of state for seawater. For decades, models used simplified equations. The adoption of the new thermodynamic standard, TEOS-10, provides a more accurate definition of neutral surfaces. Using these more accurate definitions can further reduce the subtle misalignments between the model's assumptions and reality, leading to a measurable reduction in spurious mixing and a purer simulation of ocean physics.
At this point, you might be wondering if this obsessive focus on numerical minutiae is truly necessary. The answer is an emphatic yes. Getting the mixing right is not an academic exercise; it is fundamental to simulating some of the most important features of our climate system.
Case Study 1: The Fate of Deep Waters
Consider the dense, salty water that spills over the sill of Gibraltar into the Atlantic, or the cold, dense water that cascades off the Antarctic continental shelf. These dense overflows are like vast undersea rivers that carry water into the abyss, forming the deep limbs of the global ocean circulation. To model them correctly, a model must preserve their high density as they descend.
Now imagine trying to simulate this with a model prone to spurious diapycnal mixing. As the dense water plume flows down the steep continental slope, a -coordinate model with simple horizontal diffusion would violently mix the plume with the surrounding lighter water. A -coordinate model would be even worse; the coordinate surfaces are aligned with the steep topography, creating a massive misalignment with the gently sloping plume, and the resulting spurious mixing would be catastrophic. In such a model, the overflow would effectively "vanish," its defining density signature erased by numerical error before it ever reaches the deep ocean. Only a model that explicitly aligns its mixing with the isopycnal surfaces—either an isopycnal-coordinate model or a -model with rotated isoneutral diffusion—can hope to simulate this crucial process correctly.
Case Study 2: Predicting El Niño
On a shorter timescale, these numerical choices have a direct impact on our ability to predict phenomena like the El Niño–Southern Oscillation (ENSO). ENSO's behavior is governed by a delicate feedback loop in the tropical Pacific between the atmosphere and the ocean, known as the Bjerknes feedback. A key player in this is the thermocline, the sharp density gradient separating the warm surface waters from the cold deep ocean.
In a -coordinate model, spurious diapycnal mixing tends to weaken and deepen this thermocline. A weaker thermocline has two major consequences. First, it reduces the effectiveness of upwelling in the eastern Pacific; the water brought to the surface is not as cold as it should be, damping the temperature anomalies that drive the atmospheric response. Second, the speed of oceanic Kelvin waves, which communicate signals across the equator and set the timing of the ENSO cycle, is dependent on the strength of the stratification. A weaker thermocline leads to slower waves. The combined effect is a change in the fundamental character of ENSO in the model—its amplitude, its frequency, and its predictability are all compromised. This illustrates a direct, traceable path from a numerical artifact to the simulation of a global climate pattern that affects weather, agriculture, and economies worldwide.
The challenge of spurious diapycnal mixing forces us to think like scientists and engineers. There is no single "best" ocean model, only the best model for a specific scientific question. The choice of which tool to use is a strategic decision based on a clear understanding of these trade-offs.
If the goal is to study rapid SST variability in response to a hurricane, and the deep ocean's memory is irrelevant, a simple slab model that captures only the mixed layer's heat budget is the perfect tool—efficient and fit for purpose.
If the goal is to trace the ventilation pathways of water masses from the surface into the ocean interior over decades, following intricate pathways along surfaces of constant density, then an isopycnal model is the undisputed champion. Its coordinate system is intrinsically aligned with the adiabatic nature of the flow.
And if the goal is to simulate the full complexity of century-scale climate change, including the slow, diabatic uptake of heat into the deep ocean, then a robust, general-purpose -coordinate model—fortified with sophisticated parameterizations like Redi and GM to correct for its inherent flaws—is often the workhorse of choice.
The struggle against spurious diapycnal mixing is a microcosm of the entire enterprise of climate modeling. It is a field where fidelity to physics must be balanced with computational reality, and where deep physical insight is required to see through the fog of numerical artifacts. It is a testament to the ingenuity of the scientific community that we can take a collection of equations, discretize them on a grid, and produce a simulation that not only looks like our world but can teach us how it works.