The Ocean Model That Learns to See in Multiple Dimensions
The Problem with Looking at the Ocean from Only One Angle
Imagine trying to understand a city's traffic patterns by watching only from directly above, or only from street level, or only from a specific altitude. Each view would reveal something different. Some patterns would be invisible from certain angles. Some congestion would blur together. And trying to switch between perspectives mid-observation would mean starting over from scratch.
This is essentially the problem oceanographers have faced with their mathematical models of the ocean. The tools they use to simulate ocean dynamics have long been forced to commit to a single perspective—a single "vertical coordinate system" that determines how the ocean is divided into layers for calculation. And each of these perspectives has blind spots.
A team of researchers at Stanford University and the U.S. Geological Survey has now built something different: an ocean model that can fluidly switch between perspectives mid-simulation, using the best lens for whatever phenomenon it's observing. Their model, published in July 2026, represents a significant step forward in our ability to simulate complex ocean processes—particularly a type of wave that shapes coastlines, moves sediment, and stirs nutrients from the deep.
The waves in question are called internal solitary waves. They're not the waves you see crashing on the beach. They're invisible to the casual observer, rolling through the ocean's interior like slow-motion tsunamis, hundreds of meters below the surface. They can be hundreds of kilometers long and carry enough energy to rattle offshore oil platforms. They reshape the seafloor. They drive mixing that brings nutrients up from the deep, feeding marine ecosystems from the bottom of the food chain upward.
Simulating these waves has historically required a painful tradeoff: you could either track the waves themselves efficiently, or you could capture what happens when they slam into the bottom and create chaos. Doing both simultaneously has been beyond the reach of existing models—until now.
The Science
The new model builds on decades of ocean modeling work, incorporating lessons learned from some of the field's most widely used tools. The researchers started with SUNTANS, a model developed at Stanford that employs an "unstructured grid"—meaning it can accommodate the irregular coastlines and complicated seafloor topography of real ocean basins rather than forcing the geometry into neat rectangles.
SUNTANS uses what modelers call a "z-coordinate" system: it divides the ocean into horizontal layers at fixed heights above a reference point, like stacked pancakes. This approach is intuitive and handles some processes well, but it struggles with stratified fluids—layers of water of different densities stacked on top of each other, like the ocean's interior. When density layers don't align with the fixed-height grid, the model's calculations become inaccurate.
Other popular models use different coordinate systems. ROMS and FVCOM, widely used in coastal modeling, employ "sigma coordinates" or terrain-following systems that hug the shape of the seafloor. These work well near the coast where the bottom is complex, but they introduce their own distortions when modeling deep water. A third family of models, including the isopycnal-coordinate approach developed previously by two of the paper's authors, tracks water masses by their density rather than their depth. This works beautifully for internal waves propagating through the open ocean, but it breaks down near boundaries where density layers can overturn and disappear.
"The vertical coordinate system is critical to the design of an ocean model," the researchers write, citing the consensus among modelers. "These vertical coordinates have been employed in different ocean models... z-level coordinates, terrain-following or sigma-coordinates, or density-following (isopycnal) or rho-coordinates."
Each of these coordinate types is a lens that reveals some aspects of ocean behavior while obscuring others. The Stanford team asked: what if you didn't have to choose?
Their solution is a "generalized vertical coordinate" system, or GVC. Rather than fixing the grid to a single perspective, the GVC framework allows the model to use different coordinate types in different regions of space and time, transitioning smoothly between them as conditions demand.
The mathematical transformation that enables this flexibility is a classic technique adapted from weather modeling. The researchers map the physical ocean space—which has irregular topography, stratified layers, and moving water—onto a simplified computational space where calculations are easier. The mapping is defined by a quantity called the Jacobian, which in this context represents the local thickness of each model layer. By allowing this Jacobian to vary in space and time, and by allowing the vertical coordinate to move, the model can represent any of the traditional coordinate systems—or combinations of them.
The core equations being solved are the Reynolds-Averaged Navier-Stokes equations with the Boussinesq approximation—essentially Newton's laws of motion applied to fluid flow, averaged over small-scale turbulence. The model tracks three velocity components, density (parameterized through salinity in these simulations), and the free surface elevation. Crucially, it includes nonhydrostatic terms, meaning it doesn't assume that pressure at any depth is determined simply by the weight of water above. This matters enormously for internal solitary waves, which are governed by a balance between nonlinear steepening and "nonhydrostatic frequency dispersion"—the very phenomenon the hydrostatic approximation ignores.
In practical terms, the model solves a three-dimensional elliptic equation for nonhydrostatic pressure at each timestep. This is computationally expensive—roughly an order of magnitude more costly than hydrostatic simulations—but it's necessary to capture the full physics of internal waves near their limits of existence.
The unstructured grid in the horizontal plane uses a C-grid arrangement, where velocities are defined at cell edges while pressure and density live at cell centers. This staggered layout helps avoid certain numerical artifacts, though the researchers note that triangular grids can produce "checkerboard noise" that they mitigate through careful implementation. For the vertical direction, the grid is divided into layers that can have variable thickness, with each layer mapped uniformly to computational space even when the physical layer heights differ.
The test case the researchers focus on involves internal solitary waves in a stratified fluid. They set up an idealized ocean scenario with a two-layer density structure—lighter water on top, denser water below—over a flat bottom. They then introduce a perturbation that generates internal waves, and track how these waves evolve, interact with the bottom boundary, and drive instability in the layer immediately adjacent to the seafloor.
What They Found
The researchers first validated their numerical methods with a canonical test case: turbulent flow in a straight channel. This "channel flow" scenario is a workhorse of computational fluid dynamics, with well-established theoretical predictions for velocity profiles near the boundaries. The model was run with three different advection schemes—methods for calculating how momentum moves with the flow—and all three reproduced the expected log-law profile in the boundary layer.
This validation step established that the model's momentum advection scheme, which was specifically designed for unstructured grids with variable layer heights, was working correctly. The scheme ensures that momentum is conserved locally and globally, even when model layers thin to near-zero thickness, which can happen near complex coastlines or in regions of strong vertical shear.
For the internal solitary wave simulations, the researchers compared four different configurations: a hybrid coordinate case (H3), a pure isopycnal case (I3), a pure z-coordinate case (Z3), and a terrain-following sigma-coordinate case (S3). Each used 60 vertical layers. The hybrid configuration used isopycnal coordinates in the upper water column to efficiently resolve the internal wave, transitioned through terrain-following layers in the middle, and switched to fixed z-coordinates near the bottom to capture the boundary layer dynamics.
The results demonstrated the power of the hybrid approach.
When the internal solitary wave passed through, it compressed the upper pycnocline—the sharp density gradient separating the two layers—and depressed the lower one. In the wake of this passage, the bottom boundary layer underwent instability. Near-bottom velocity measurements showed oscillations driven by the wave, with peak amplitudes exceeding the nominal boundary layer thickness.
The key diagnostic the researchers used was background potential energy: the energy required to redistribute the density field back to a reference state. An increase in background potential energy represents irreversible mixing, energy that's lost from the organized wave motion and transferred to turbulence and heat.
Mixing Efficiency by Coordinate Type (Nk=60)
| Label | Value |
|---|---|
| Hybrid Coordinate (H3) | 0.82 |
| Isopycnal Coordinate (I3) | 0.75 |
| Z-Level Coordinate (Z3) | 0.68 |
| Sigma Coordinate (S3) | 0.55 |
Background Potential Energy Evolution
| Label | Value |
|---|---|
| 0 | 0 |
| 0.2 | 0.2 |
| 0.4 | 0.4 |
| 0.6 | 0.6 |
| 0.8 | 0.8 |
| 1.0 | 1 |
The simulations showed that the hybrid coordinate configuration captured the wave evolution and boundary layer instability with the same fidelity as pure isopycnal coordinates for the wave tracking, while simultaneously resolving the bottom boundary layer processes that pure isopycnal coordinates cannot handle without extreme vertical refinement.
In contrast, the pure z-coordinate simulation required far more vertical layers to achieve comparable resolution of the internal wave itself—because z-coordinates waste resolution on unstratified regions and lack the Lagrangian following of isopycnal coordinates that concentrates layers where density gradients concentrate.
The pure sigma-coordinate simulation suffered from "sigma-coordinate error": when the terrain-following layers are steeply sloped, pressure-gradient calculations become inaccurate, leading to spurious velocities. The researchers observed that even with 60 vertical layers, the sigma simulation showed numerical artifacts near the bottom.
The hybrid approach thus achieved something that none of the individual coordinate systems could: efficient resolution of the internal wave while simultaneously capturing the boundary layer response, without numerical artifacts from coordinate distortion.
The model reproduced the theoretical expectation that internal solitary waves in a two-layer system are governed by the KdV (Korteweg-de Vries) equation for small amplitudes, with the nonlinearity and dispersion balancing to produce waves of constant form. The wave amplitude, propagation speed, and structure matched theoretical predictions.
What distinguished the Stanford model was its ability to capture the wave-bottom interaction that previous models struggled with: the oscillation of the bottom boundary layer driven by the wave's passing, and the secondary instability that developed as the wave decayed.
Why This Changes Things
Ocean modeling has long been constrained by the tyranny of the coordinate system. When you're building a telescope, you have to choose your lens before you start observing. But the ocean doesn't cooperate with our choices—it throws density layers, internal waves, and boundary turbulence at us all at once, in the same domain, sometimes in the same location.
Internal solitary waves are ubiquitous in the ocean. They form wherever strong stratification meets tides flowing over topography—off the coasts of California and Peru, in the South China Sea (where some of the largest examples on Earth have been documented), in the Andaman Sea, in the Sulu Sea. They propagate across hundreds of kilometers, interacting with each other and with the boundaries of ocean basins. They drive mixing that influences everything from climate to fisheries.
The South China Sea alone has been the subject of intense research interest because of its spectacular internal wave field. Waves with amplitudes exceeding 100 meters have been observed, generated by tidal flow over the Luzon Strait and propagating westward toward the Chinese mainland. These waves have been implicated in everything from submarine navigation challenges to sediment resuspension on the continental shelf.
Simulating these waves accurately requires nonhydrostatic models, which capture the dispersive character that distinguishes solitary waves from simpler wave theories. But nonhydrostatic models are expensive. The Stanford team cites prior work showing that isopycnal coordinates can reduce computational cost by decreasing the required number of vertical layers, because density-following layers naturally concentrate where the action is.
The hybrid approach takes this insight further. By using isopycnal coordinates where the wave lives, the model can run efficiently through the open ocean. By switching to z-coordinates near the bottom, it captures the boundary layer processes that drive sediment transport, nutrient flux, and marine ecosystem dynamics. By using terrain-following coordinates as transition layers, it avoids the numerical problems that occur when coordinate systems collide.
This is conceptually similar to adaptive mesh refinement in other areas of computational physics, where the computational grid automatically concentrates where it's needed. But here the adaptation is in the coordinate system itself, not just the resolution—making it compatible with existing model frameworks while extending their flexibility.
The practical implications extend beyond internal waves. The same coordinate framework could be applied to other phenomena where different vertical processes dominate in different regions: estuarine circulation, overflows, convective mixing, internal tides. Any problem where you want the model to "zoom in" automatically where interesting things are happening, without requiring the modeler to specify this behavior a priori.
The researchers note that similar grid formulations combining different coordinate types have been used successfully in weather modeling, particularly to reduce noise near topography. Their work extends this approach to nonhydrostatic ocean modeling with unstructured grids, opening the door to realistic domain simulations with complex coastlines and internal wave fields.
The computational efficiency gains could be significant for real-world applications. A model that can use fewer vertical layers where the physics permits—because those layers are placed intelligently—will run faster and consume less memory. For operational ocean forecasting, where model runs must complete within tight time windows, this matters.
For scientific discovery, the ability to simulate internal wave-bottom interactions at the fidelity the hybrid approach enables could unlock new understanding. These interactions drive boundary layer turbulence that matters for sediment transport, biogeochemistry, and the exchanges between the deep ocean and the continental margins.
What's Next
The paper is explicitly labeled as the first in a two-part series. The researchers' immediate next step is to apply the method to internal solitary wave shoaling and breaking using "horizontally hybrid vertical coordinates." This refers to a different kind of hybrid: not mixing coordinate types vertically, but allowing different coordinate systems to dominate in different horizontal regions of the domain.
For wave shoaling—the process by which waves grow in amplitude as they move from deep water toward the coast—the horizontal domain structure matters. In deep water, isopycnal coordinates might be optimal. Near the shelf break where waves steepen and eventually break, terrain-following or z-coordinates might be preferred. The horizontal hybrid approach would allow seamless transitions across these regions.
Several open questions remain. The researchers have demonstrated their method on idealized test cases; applying it to real ocean basins with realistic topography, stratification, and forcing will require further validation. The transition between coordinate systems must be handled carefully to avoid introducing numerical artifacts or violating conservation principles.
The computational cost of the nonhydrostatic pressure solver remains a challenge. While the hybrid approach reduces the number of vertical layers needed, the fundamental expense of solving a 3D elliptic equation at each timestep doesn't disappear. Further optimization, potentially using graphics processing units (GPUs) or other parallel computing approaches, will be important for making these simulations practical at large scales.
The researchers also note that their current implementation uses a linear equation of state relating salinity to density. Real seawater obeys a nonlinear equation of state that depends on both temperature and salinity, which introduces additional complications like double-diffusive convection. Extending the model to handle full equation-of-state effects will be necessary for some applications.
More broadly, the work points toward a future where ocean models are more adaptive and intelligent about where they concentrate their resources. The ocean is not a uniform fluid in a rectangular basin; it's a complex, stratified, evolving system that throws different phenomena at us at different scales and in different places. A modeling framework that can adapt its perspective to match the physics at hand—rather than forcing all physics through a single predetermined lens—is better equipped to capture that complexity.
The Stanford team's contribution is both technical and conceptual. The technical advances—a conservative momentum advection scheme for unstructured grids, a positivity-preserving layer height scheme, locally and globally conservative scalar transport—solve specific numerical challenges that have limited previous approaches. The conceptual advance is the demonstration that hybrid coordinate systems can offer advantages that none of their components achieve alone.
As climate change intensifies the need to understand and predict ocean behavior, from sea level rise to biogeochemical cycling, models that can simulate the ocean's complexity more efficiently and accurately become increasingly valuable. Internal waves are just one piece of that puzzle—but they're a piece that shapes the ocean's vertical structure, drives mixing between layers, and connects the deep circulation to coastal processes.
The waves are already out there, rolling through the ocean's interior, reshaping the seafloor, and stirring the water column. Models that can capture their full complexity—without forcing us to sacrifice the bottom boundary for the wave interior—bring us closer to understanding a system that remains, despite decades of study, remarkably underobserved and surprisingly mysterious.
The code is built on the open-source SUNTANS framework, meaning other researchers can examine, use, and extend the implementation. As the method matures and is applied to more realistic scenarios, it could become a standard tool in the ocean modeler's toolkit—another instrument in the expanding orchestra of computational approaches that together constitute our best efforts to understand the ocean that covers most of our planet's surface and absorbs most of the heat our civilization releases into the atmosphere.