Linear Shallow Water Equations: From Core Principles to Global Forecasting

The vast, chaotic movements of the Earth's oceans and atmosphere can seem bewilderingly complex, yet they are governed by a set of core physical principles. The shallow water equations are a cornerstone of geophysical fluid dynamics, providing a powerful lens through which to understand this planetary-scale motion. However, their full nonlinear form can be mathematically challenging. To build a foundational understanding, we can simplify our approach by examining the linear shallow water equations, which describe small disturbances to a fluid at rest. This simplification peels back layers of complexity to reveal the fundamental physics of waves, rotation, and adjustment that orchestrate the behavior of our fluid world.

This article provides a journey into the heart of these elegant equations. In the first chapter, "Principles and Mechanisms," we will explore the fundamental physics, starting with the simple feedback between water height and velocity that creates a wave. We will then see how planetary rotation complicates this picture, introducing new wave types, fundamental length scales, and the powerful concept of geostrophic balance. Finally, we will uncover how a conserved quantity, potential vorticity, provides a unifying framework for understanding this adjustment process. In the following chapter, "Applications and Interdisciplinary Connections," we will bridge theory and practice. We will see how these principles explain the destructive speed of tsunamis, underpin the algorithms that power supercomputer simulations of the ocean, and guide the sophisticated methods used in modern weather and climate forecasting. Our exploration begins with the core principles that form the bedrock of our understanding of planetary fluids.

To truly appreciate the symphony of the oceans and atmosphere, we must first learn to listen to its most fundamental notes. These notes are waves, and their orchestra is governed by a beautifully simple set of rules: the shallow water equations. While the full score can be dauntingly complex and nonlinear, we can understand the essential harmony by studying a linearized version, which describes small disturbances to a vast, calm body of water. This simplified model strips away the complexities to reveal the pure physics at play, offering a surprisingly accurate picture of everything from tsunamis to the large-scale structure of weather systems.

Imagine a perfectly still, infinitely long canal of water with a uniform depth, which we'll call $H$. What happens if you momentarily create a small hump of water? Your intuition tells you the hump won't just sit there; it will spread out. But how?

This is where the physics begins. A hump in the water creates a pressure gradient. The water at the peak of the hump is higher, so it's under slightly more pressure from the weight above it than the water at its flanks. This pressure difference pushes the water outwards, away from the peak. So, a spatial variation in height ($\frac{\partial \eta}{\partial x}$) creates an acceleration, a change in velocity over time ($\frac{\partial u}{\partial t}$).

But the story doesn't end there. As the water begins to flow, it must go somewhere. Where the flow converges ($\frac{\partial u}{\partial x} 0$), water piles up, increasing the height. Where it diverges ($\frac{\partial u}{\partial x} > 0$), the water level drops. So, a spatial variation in velocity creates a change in height over time ($\frac{\partial h}{\partial t}$).

This beautiful feedback loop is the essence of a wave: height differences drive flows, and flows create height differences. It is a perpetual dance between potential energy (stored in the height of the water) and kinetic energy (in the motion of the water). The linearized shallow water equations capture this dance with elegant simplicity:

Here, $\eta$ is the tiny perturbation of the height from its mean value $H$, $u$ is the fluid velocity, and $g$ is the acceleration due to gravity. The first equation is a statement of mass conservation (the divergence of flow changes the height), and the second is Newton's second law (the pressure gradient from the height slope accelerates the fluid).

By combining these two equations, we can ask a simple question: how fast does the disturbance travel? The answer turns out to be astonishingly simple. There are two characteristic speeds, $\lambda = \pm\sqrt{gH}$. This means any disturbance will split into two waves, one traveling right and one traveling left, both with a speed of $c = \sqrt{gH}$. This speed depends only on the depth of the water and gravity. For a typical ocean depth of $4$ km, this speed is about $200$ m/s, or $720$ km/h—the speed of a jet airliner. This is the speed of tsunamis crossing the open ocean. The mathematical term for such a system, where information propagates at real and distinct speeds, is ​​strictly hyperbolic​​.

Now, let's take our canal and expand it into a vast ocean on a spinning planet. We must now account for the ​​Coriolis effect​​, an apparent force that deflects any moving object—be it an air parcel or a water parcel—to the right in the Northern Hemisphere and to the left in the Southern Hemisphere. This isn't a magical force; it's a consequence of observing motion from within a rotating frame of reference while momentum is conserved in the absolute sense.

Adding the Coriolis effect, characterized by a parameter $f$, modifies our equations. For waves, it introduces a new restoring force. Not only does gravity try to flatten a height anomaly, but the Coriolis force acts on the resulting motion, deflecting it and creating its own dynamic response. The simple gravity waves now become ​​inertia-gravity waves​​, also known as ​​Poincaré waves​​.

Their behavior is captured in a dispersion relation, $\omega^2 = f^2 + gH|\mathbf{k}|^2$, where $\omega$ is the wave's frequency and $|\mathbf{k}|$ is its wavenumber (inversely related to wavelength). This equation tells us something profound. For very short, fast waves (large $|\mathbf{k}|$), the $f^2$ term is negligible, and they behave just like the simple gravity waves from our canal. But for any wave, no matter how long its wavelength, the frequency $\omega$ can never be less than $f$. The rotation of the planet sets a fundamental low-frequency limit on wave propagation.

This rotation also changes the partitioning of energy. For a simple gravity wave, kinetic and potential energies are, on average, equal. For an inertia-gravity wave, the ratio of time-averaged kinetic energy ($E_K$) to potential energy ($E_P$) is given by $\frac{E_K}{E_P} = \frac{\omega^2+f^2}{\omega^2-f^2}$. As the wave's frequency $\omega$ approaches the inertial frequency $f$, the kinetic energy becomes enormous compared to the potential energy. This hints that slow motions, near this frequency limit, are dominated by currents rather than by height fluctuations.

The low-frequency cutoff for waves poses a fascinating puzzle: what happens if we try to create a disturbance that is very slow, or simply create an imbalance and let it evolve on its own? The system cannot simply radiate the energy away with low-frequency waves, because such waves cannot exist.

Instead, the fluid performs a remarkable process called ​​geostrophic adjustment​​. Any initial "unbalanced" state, like a pile of water sitting at rest, rapidly splits into two components. A portion of the initial energy excites high-frequency inertia-gravity waves that radiate away to infinity, like the ripples from a stone dropped in a pond. What is left behind is a steady, balanced state where the Coriolis force on the moving fluid is in perfect opposition to the pressure gradient force from the sea surface slope. This is ​​geostrophic balance​​, the defining principle of nearly all large-scale weather and ocean current systems. It's why wind on a weather map flows along lines of constant pressure (isobars), not directly from high to low pressure.

But what determines the size of the structures left behind? This is set by a crucial length scale: the ​​Rossby radius of deformation​​, $R_d = \frac{\sqrt{gH}}{f}$. This radius can be thought of as the distance a gravity wave (speed $\sqrt{gH}$) can travel in the time it takes for rotation to become important (roughly $1/f$).

• If you create a disturbance much smaller than the Rossby radius, rotation doesn't have time to act before the disturbance disperses as gravity waves. The initial height anomaly simply flattens out.
• If you create a disturbance much larger than the Rossby radius, rotation is dominant. The disturbance cannot easily disperse; instead, it adjusts, with most of its initial structure settling into a geostrophically balanced flow.

Imagine starting with a rectangular block of elevated water at rest. As geostrophic adjustment proceeds, waves shoot out, and the sharp corners of the block are smoothed away. The system settles into a new, permanent state with a gentle slope in sea surface height and a steady current flowing parallel to the slope, a perfect demonstration of geostrophy. The width of this final, balanced structure is governed by the Rossby radius. On Earth, this scale is tens to hundreds of kilometers in the ocean and thousands of kilometers in the atmosphere, setting the characteristic size of ocean eddies and high-pressure weather systems.

The process of geostrophic adjustment seems almost intelligent. How does the fluid "know" which balanced state to relax into? The secret lies in a deeply powerful conserved quantity known as ​​potential vorticity (PV)​​.

For our simple system, the linearized potential vorticity anomaly is given by $q' = \frac{\zeta}{H} - \frac{f\eta}{H^2}$, where $\zeta$ is the relative vorticity (the local "spin" of the fluid) and $\eta$ is the height perturbation. In the absence of friction, every parcel of fluid conserves its value of $q'$ as it moves.

When an initial disturbance is created, it has a certain PV field. As the fast inertia-gravity waves radiate away, they carry energy and momentum, but they do not alter the PV of the fluid parcels left behind. Therefore, the final, steady, geostrophically balanced state must have the same potential vorticity field as the initial state.

This gives us an incredible tool. Knowing the initial state (e.g., a height anomaly with zero velocity) allows us to calculate the initial PV field. Since we know this field is conserved, we can find the final balanced state by solving for the unique geostrophically balanced flow that possesses this PV field. This procedure is called ​​PV inversion​​. Mathematically, it takes the form of a Helmholtz equation, $\nabla^2 \psi - \frac{1}{R_d^2} \psi = H q'$, where $\psi$ is the streamfunction describing the flow. This single equation beautifully encapsulates the entire principle: the final balanced flow ($\psi$) is determined entirely by the conserved PV ($q'$) and the fundamental length scale of the system ($R_d$).

Our story so far has assumed a constant Coriolis parameter $f$. But what happens near the equator, where $f$ is zero and changes sign? The variation of $f$ with latitude, denoted by the parameter $\beta$, turns the equator into a remarkable ​​waveguide​​, trapping energy and creating new types of waves that can travel for thousands of kilometers across an ocean basin.

Two of the main actors on this equatorial stage are Kelvin waves and Rossby waves.

The ​​Equatorial Kelvin Wave​​ is a marvel of simplicity and importance. It is a wave that has no north-south motion whatsoever ($v=0$). North of the equator, the eastward-flowing water is deflected rightward (northward) by the Coriolis force. This deflection is perfectly balanced by a pressure gradient force from a sea surface that slopes down to the north. South of the equator, the same eastward flow is deflected leftward (southward), balanced by a sea surface that slopes down to the south. The result is a hump of water, trapped at the equator, that can only propagate eastward. Amazingly, it does so without changing its shape and at the pure gravity wave speed $c=\sqrt{gH}$. These waves are the primary way that large-scale disturbances, such as those associated with El Niño, communicate across the Pacific Ocean.

The other leading characters are the ​​Equatorial Rossby Waves​​. Unlike Kelvin waves, their very existence depends on the variation of the Coriolis parameter, $\beta$. They are slow, westward-propagating waves that are crucial for the long-term adjustment of the ocean to changes in wind forcing. Together with Kelvin waves and other trapped modes like the mixed Rossby-gravity wave, they form a complete set of tools that the equatorial ocean and atmosphere use to move energy and information around, orchestrating climate variability on a global scale.

From a simple dance of height and velocity in a canal, we have journeyed to a rotating world of intricate adjustments, powerful conserved quantities, and a special zoo of equatorial waves. The linear shallow water equations, in their simplicity, provide the fundamental key to understanding the majestic and complex motions of our planet's fluid envelope.

Having acquainted ourselves with the principles and mechanisms of the linear shallow water equations, we might be tempted to file them away as a neat, but perhaps overly simplified, academic exercise. Nothing could be further from the truth. These elegant equations are not a final destination but a gateway. They are the "hydrogen atom" of geophysical fluid dynamics—a simple, solvable system that unlocks a profound understanding of the complex machinery of our planet's oceans and atmosphere. Embarking on a journey with these equations, we will see how they explain the awesome power of tsunamis, dictate the very rules by which we build our supercomputer simulations, and even guide the sophisticated art of weather and ocean forecasting.

The most immediate and striking application of our equations is their ability to predict the speed of long waves. In an idealized channel, if we were to suddenly raise the water level on one side, creating a "dam break," the disturbance would not spread instantaneously. Instead, it would propagate outwards as two distinct wave fronts. The speed of these fronts, as derived from the very heart of the shallow water equations, is a beautifully simple expression:

where $g$ is the acceleration due to gravity and $H_{0}$ is the undisturbed water depth. This isn't just a formula; it is a fundamental truth about our world, revealing itself in some of nature's most dramatic events.

Consider a tsunami, born from a submarine earthquake. In the vast, open Pacific, the ocean depth $H_{0}$ is on the order of 4000 meters. A quick calculation reveals a wave speed $c$ of about 200 meters per second, or over 700 kilometers per hour—the cruising speed of a modern jetliner. This explains the tsunami's terrifying nature: it traverses entire ocean basins in a matter of hours, its long wavelength and small amplitude making it nearly undetectable in the deep ocean, only to rear up into a destructive wall of water as it enters the shallow coastal regions where $H_{0}$ decreases.

The same physics governs the more gentle, daily rhythm of the tides. As the gravitational pull of the Moon and Sun creates a vast bulge of water, this tidal wave propagates across the oceans. Its wavelength is so immense that even the deep ocean is "shallow" relative to it. Its movement into bays, estuaries, and rivers is a magnificent, real-world demonstration of shallow water wave dynamics, with the speed given by our faithful formula.

Of course, real waves in a real harbor or lake do not slosh back and forth forever. Anyone who has watched the ripples from a thrown stone die out knows that nature has its frictions. By adding a simple linear damping term to our momentum equation, we can account for the drag exerted by the seabed. This allows us to model phenomena like a seiche, a standing wave in an enclosed or semi-enclosed body of water, and correctly predict that its energy will slowly dissipate, its amplitude decaying exponentially over time. This simple addition of a friction term is a crucial step from an idealized physical model to one with greater real-world fidelity.

To move beyond simple channels and basins and model the intricate dance of water in real oceans, we turn to the power of computers. But a computer cannot grasp the elegant, continuous form of our equations. It understands only numbers. We must therefore translate our continuous world into a discrete one, a process called discretization. We lay a grid over our ocean and solve the equations at a finite number of points, stepping forward in finite increments of time. It is here that the shallow water equations reveal their second, deeper layer of influence: they dictate the very rules of their own simulation.

The most fundamental rule is the Courant-Friedrichs-Lewy (CFL) condition. It tells us that for an explicit numerical simulation to remain stable—to not "explode" into a chaos of meaningless numbers—the time step $\Delta t$ we choose is limited by the grid spacing $\Delta x$ and the fastest-moving wave in the system. The relationship is, once again, beautifully simple:

where $C_{\max}$ is the Courant number, a value typically less than or equal to one. The physical wave speed $c = \sqrt{gH}$ that governs tsunamis now reappears as a strict speed limit for our simulation! To capture the physics accurately, the numerical domain of influence (what a grid point "sees" in one time step) must be at least as large as the physical domain of influence. The physics dictates the algorithm. This connection between the physical properties of the system and the stability of its numerical representation is a cornerstone of computational science.

The art of building a good numerical model, however, goes deeper. A naive placement of all variables at the same grid points can lead to spurious, checkerboard-like patterns in the solution—a kind of numerical noise. A far more elegant solution, used in the vast majority of modern ocean and atmosphere models, is the staggered grid, such as the Arakawa C-grid. In this clever arrangement, water level ($\eta$) is stored at the center of a grid cell, while velocities ($u,v$) are stored on the cell faces. This may seem like a mere bookkeeping trick, but its effect is profound. It creates the most natural, robust coupling between the pressure gradient (which depends on differences in $\eta$) and the velocity divergence (which depends on differences in $u$ and $v$), effectively preventing the numerical decoupling that can plague simpler grids.

Even with these clever tricks, our digital ocean is not a perfect replica of reality. Every discretization introduces errors. One of the most insidious is numerical dispersion. In the real world, all shallow water waves, regardless of their wavelength, travel at the same speed, $c = \sqrt{gH}$. In a numerical model, this is no longer true: the wave speed becomes a function of the wavelength. A careful analysis shows that the ratio of the numerical phase speed to the true phase speed is always less than one, meaning the model makes waves travel too slowly. The error is most pronounced for the shortest waves that the grid can resolve. This is a humbling and essential lesson for any modeler: the digital world has its own peculiar physics, and understanding these artifacts is as important as understanding the underlying physical laws themselves.

Finally, what about the edges of our model? A regional ocean model, say of the Gulf of Mexico, is not an isolated tub of water; it is connected to the wider Atlantic Ocean. How do we create an "open" boundary that allows waves generated inside our model to pass out freely, without reflecting back, while simultaneously allowing waves from the outside world to enter? The answer lies in the method of characteristics. By rewriting the shallow water equations, we can identify two types of information: one that propagates out of the model, and one that propagates in. A sophisticated radiation boundary condition, such as the Flather condition, uses the model's own interior solution to define the outgoing characteristic, while specifying the incoming one from an external source (like a larger global model). This elegant piece of applied mathematics ensures that our artificial boundaries are as transparent as possible, a critical component for any realistic regional simulation of oceans or atmospheres.

When we scale up to global ocean models that run for decades or centuries, or operational weather forecasts that must be delivered on a strict deadline, two new challenges emerge: computational efficiency and the integration of real-world data. Here, too, the shallow water equations provide the conceptual framework.

The ocean is a multi-layered system. It has the fast, surface-height-driven motions we have been discussing, known as the external or barotropic mode. But it also has much, much slower internal waves that move along density surfaces deep within the ocean, known as internal or baroclinic modes. The speed of the external waves is our familiar $c_0 = \sqrt{gH} \approx 200 \text{ m/s}$. The speed of the fastest internal wave, $c_1$, is typically only a few meters per second. Because the CFL condition is governed by the fastest wave, a straightforward simulation would be forced to take tiny time steps, on the order of minutes, dictated by the barotropic waves. Yet the evolution of the ocean's heat content and internal structure happens on timescales of days to years. This is a computational nightmare.

The solution is a technique called mode splitting. Since the barotropic speed is about 100 times faster than the baroclinic speed, the time step required is about 100 times smaller. Mode splitting brilliantly exploits this separation of timescales. The model is algorithmically split: the fast barotropic motions are calculated separately with the necessary small time step, while the slow, computationally expensive baroclinic motions and tracer transport are updated with a much larger time step. This simple idea, directly motivated by the disparity in wave speeds, has made long-term, high-resolution ocean modeling computationally feasible.

A model running in isolation, no matter how sophisticated, will inevitably drift from reality. To produce a useful forecast, it must be constantly nudged back on track by real-world observations. This process is called data assimilation. But one cannot simply "paste" an observation—say, from a satellite measuring sea surface height—into the model. Doing so would "shock" the system, creating a storm of spurious, high-frequency inertia-gravity waves that would radiate through the model, corrupting the forecast.

The key to avoiding this is the concept of dynamical balance. On the large scales relevant to climate and weather, the rotating ocean and atmosphere are, to a good approximation, in a state of geostrophic balance, where the pressure gradient force is balanced by the Coriolis force. This is a state of slow, graceful evolution, dominated by vorticity (swirling motions) rather than divergence (spreading-out motions). The goal of modern data assimilation is to ensure that the correction applied to the model is itself dynamically balanced. This is achieved by building the physics of geostrophic balance directly into the statistical heart of the assimilation system—the background error covariance matrix, $B$. This matrix tells the system how errors in one variable (like sea surface height) are related to errors in other variables (like velocity). By using our physical understanding to construct these relationships, we guide the assimilation to produce analysis updates that lie on the "slow manifold" of balanced flow, minimizing the generation of noisy gravity waves. Modelers can even diagnose the health of their system by measuring the ratio of divergence to vorticity in the analysis increments, or by projecting the updates onto the system's normal modes to see what fraction of the energy is going into the unwanted gravity wave modes. This is a beautiful synthesis of dynamics, statistics, and control theory, and it all rests on the fundamental distinction between balanced flows and gravity waves, a distinction first made clear by the shallow water equations.

From the speed of a tsunami to the stability of a supercomputer code, from the design of a numerical grid to the statistical fusion of a model with satellite data, the linear shallow water equations are a constant, unifying thread. They remind us that within the most complex systems, there often lies a simple, elegant core whose principles, once understood, reverberate through every layer of our scientific and technological endeavors.