Numerical study of the sharp stratification limit towards bilayer models

This paper numerically and theoretically compares continuously stratified and bilayer models for oceanic flows, proving convergence to bilayer equations for near-piecewise constant profiles while demonstrating via normal mode analysis that Kelvin-Helmholtz instabilities in shear flows fundamentally limit the validity of widely used bilayer Euler and shallow-water models.

The Gist

Imagine the ocean not as a uniform blue soup, but as a giant, layered cake. In the real world, this "cake" is made of water with slightly different densities at different depths (due to temperature and salt). This layering is called stratification.

Scientists use math to predict how waves move through this ocean cake. They have two main ways of doing this:

1. The "High-Definition" Camera (Stratified Models): This approach treats the ocean as a smooth, continuous gradient of density, like a perfect gradient cake. It's incredibly accurate but requires a supercomputer to solve because the math is so complex.
2. The "Blocky" Lego Model (Bilayer Models): This approach simplifies the ocean into just two distinct layers: a top layer of light water and a bottom layer of heavy water, separated by a sharp, invisible line. It's much easier to calculate, like building with Legos instead of sculpting clay.

The Big Question:
Can we safely replace the complex, high-definition ocean with the simple, blocky Lego version? Specifically, what happens when the transition between the layers becomes super thin (almost a sharp line)?

The Discovery:
The author of this paper, Théo Fradin, ran a massive numerical experiment to answer this. He found that the answer depends entirely on whether the water is moving smoothly or if there is a shear flow (where the top layer moves in one direction and the bottom layer moves in another, like a river current sliding over a slower current).

The Two Scenarios

1. The Calm Ocean (No Shear)

Imagine the ocean layers are just sitting there, or moving together in the same direction.

• The Result: The "Lego model" works perfectly! As the transition layer gets thinner, the complex high-definition model slowly morphs into the simple blocky model. The Lego approximation is a valid, accurate shortcut.
• The Analogy: It's like zooming out on a photo of a smooth gradient. At a distance, it looks like two solid blocks of color, and that's a good enough description for most purposes.

2. The Turbulent Ocean (With Shear)

Now, imagine the top layer is rushing east, and the bottom layer is rushing west. They are sliding past each other.

• The Result: Chaos. The simple "Lego model" breaks down completely.
• The Culprit: A phenomenon called Kelvin-Helmholtz Instability.
  • The Metaphor: Think of blowing across the top of a cup of hot coffee. The wind (top layer) moves fast, the coffee (bottom layer) is still. This creates little ripples that grow into bigger waves. In the ocean, if the speed difference is too high, these ripples don't just stay small; they explode into chaotic, violent turbulence.
• The Problem: In the real, high-definition ocean, this instability happens, but it takes a tiny bit of time to grow. However, in the simplified "Lego" model, the math says these instabilities grow infinitely fast.
• The Consequence: Because the Lego model predicts infinite chaos instantly, it becomes mathematically impossible to use it to predict the future of the real ocean in these conditions. The "shortcut" leads to a dead end.

The "Sharp" Limit

The paper investigates what happens when the transition layer (the "pycnocline") gets thinner and thinner, approaching zero thickness.

• Without Shear: As the layer gets thinner, the complex model smoothly converges to the simple model.
• With Shear: As the layer gets thinner, the "explosive" instabilities in the complex model get faster and faster. The simple model predicts they are already infinite. This means you cannot mathematically prove that the simple model is a good approximation of the real one when shear is present.

The Takeaway

This paper is a warning label for oceanographers and mathematicians.

• If you are studying calm, layered water, you can safely use the simple, fast "two-layer" models.
• If you are studying areas with strong currents sliding past each other (like the Gulf Stream or internal waves in turbulent zones), do not trust the simple two-layer models. They miss the crucial physics of how the layers interact, leading to predictions that are mathematically broken.

In short: You can simplify the ocean's "cake" into two layers only if the layers aren't sliding past each other too fast. If they are, the cake doesn't just split; it explodes, and the simple math can't handle the blast.

Technical Summary

1. Problem Statement

The paper investigates the mathematical and numerical validity of approximating continuously stratified oceanic flows using bilayer models (models with two layers of constant density separated by a sharp interface).

• Context: In geophysical fluid dynamics, density stratification is often modeled either by continuous profiles (Stratified Euler Equations) or by simplified bilayer models. Bilayer models are computationally cheaper and analytically tractable but rely on the assumption that the transition layer (pycnocline) is infinitely thin ($\delta \to 0$).
• The Gap: While the convergence of continuous models to bilayer models is well-understood in the absence of shear flow, the situation changes drastically when a shear flow (horizontal velocity varying with depth) is present.
• Core Question: Does the linearized Stratified Euler Equations converge to the Bilayer Euler Equations in the "sharp stratification limit" ($\delta \to 0$) when shear is present? Specifically, how do Kelvin-Helmholtz (KH) instabilities, which render bilayer models ill-posed in Sobolev spaces, affect this convergence?

2. Methodology

The author employs a hybrid approach combining rigorous mathematical analysis with high-fidelity numerical simulations.

A. Mathematical Framework

• Equations: The study focuses on the linearized Stratified Euler Equations in a 2D strip ($T_L \times [-H, 0]$).
• Coordinates: The equations are reformulated in isopycnal coordinates (coordinates following density surfaces). This transformation is crucial as it maps the moving interface in bilayer models to a fixed coordinate in the continuous model, facilitating comparison.
• Modal Decomposition: The author derives a normal mode decomposition based on the Taylor-Goldstein equation. The solution is expanded into a series of vertical modes (eigenfunctions of a Sturm-Liouville problem) coupled via linear operators.
  • This reduces the PDE system to a system of coupled ODEs for the mode coefficients.
  • The coupling is driven by the non-constant Brunt-Väisälä frequency ($N^2$) and the shear profile.

B. Numerical Strategy

• Discretization:
  • Vertical: A spectral-like truncation of the modal decomposition (keeping $N$ modes). The eigenvalues and eigenfunctions of the Sturm-Liouville problem are computed using finite differences.
  • Horizontal: Truncated Fourier series (spectral method).
  • Time: Fourth-order Runge-Kutta method.
• Dispersion Relation Analysis: The primary numerical tool is the computation of the dispersion relation (phase velocities $c$ vs. wavenumber $k$). By solving the eigenvalue problem of the discretized system matrix, the author identifies stable and unstable modes.
• Profiles Tested:
  • Density: A smooth profile converging to a step function as $\delta \to 0$ (pycnocline thickness $\delta$).
  • Shear: A smooth shear profile $V^\delta$ converging to a piecewise constant shear (discontinuity at the interface) as $\delta \to 0$.

3. Key Contributions

1. Rigorous Convergence without Shear ($V=0$)

The author provides a full justification of the bilayer limit in the absence of shear.

• Result: Proved that solutions to the linear Stratified Euler equations converge to the linear Bilayer Euler equations in Sobolev spaces as $\delta \to 0$.
• Rate: The error is bounded by $O(\delta^{1/2} |\log \delta|)$.
• Significance: This confirms that bilayer models are valid approximations for stable, non-sheared stratified flows.

2. Numerical Evidence of KH Instabilities in Continuous Models

In the presence of shear, the bilayer Euler equations are known to be ill-posed in Sobolev spaces due to KH instabilities (high-frequency modes grow exponentially).

• Finding: The numerical simulations of the continuous Stratified Euler equations (with a thin but non-zero $\delta$) reveal that KH instabilities are indeed present.
• Mechanism: As $\delta \to 0$, the continuous model exhibits a range of unstable wavenumbers $[k_{min, \delta}, k_{max, \delta}]$.
  • Low frequencies: Stable (asymptotically stable as $\delta \to 0$).
  • Intermediate frequencies: Unstable (KH instabilities).
  • High frequencies: Surprisingly, the continuous model shows stability for very high frequencies ($k > k_{max, \delta}$), unlike the bilayer model where all high frequencies are unstable.

3. Quantitative Scaling Laws

The paper provides conjectured scaling laws for the instability features as $\delta \to 0$:

• Upper Threshold: $k_{max, \delta} \approx \delta^{-1}$. The range of unstable frequencies extends to infinity as the interface sharpens.
• Growth Rate: The maximum growth rate of the unstable modes scales as $\text{Im}(\omega^*_\delta) \approx \delta^{-1}$.
• Implication: The growth rate becomes arbitrarily fast as the pycnocline thins.

4. Results and Implications

A. Failure of Full Justification in Sobolev Spaces

The most significant result is the negative answer to the justification of bilayer models in the presence of shear within Sobolev spaces.

• Reasoning: Because the growth rate of the unstable modes scales as $\delta^{-1}$, any solution starting with generic initial data in a Sobolev space $H^s$ will blow up in a time interval that shrinks to zero as $\delta \to 0$.
• Conclusion: One cannot prove that the continuous solution stays close to the bilayer solution on a time interval independent of $\delta$. The "full justification" fails.

B. The Shallow-Water Paradox

The paper addresses a counter-intuitive scenario involving the Bilayer Shallow-Water Equations (BSWE).

• Context: BSWE are well-posed in Sobolev spaces (unlike Bilayer Euler) and are often derived as the limit of Stratified Euler when both $\delta \to 0$ (sharp stratification) and $\mu \to 0$ (shallow water parameter).
• Finding: Even though the reduced BSWE model is well-posed, the continuous Stratified Euler equations leading to it are not. The KH instabilities in the continuous model grow so fast ($\sim \delta^{-1}$) that they prevent the convergence to the BSWE in Sobolev spaces, regardless of how small $\mu$ is.
• Nuance: The paper suggests that justification might still be possible in analytic function spaces (where Fourier coefficients decay exponentially), as the instability growth might be tamed by the rapid decay of initial data, but this remains a conjecture.

5. Significance

• Theoretical Impact: The work challenges the standard practice of using bilayer models for sheared flows in geophysical contexts without acknowledging the underlying ill-posedness of the limit. It highlights a fundamental disconnect between continuous and discontinuous models in the presence of shear.
• Numerical Insight: It demonstrates that numerical simulations of continuous models with thin interfaces can accurately capture the onset of KH instabilities, serving as a proxy for the ill-posed bilayer limit.
• Future Directions: The paper calls for rigorous proofs of the conjectured dispersion relation features (Conjecture 6.5) and suggests extending the modal decomposition methods to non-linear settings, potentially requiring fast Sturm-Liouville transforms.

In summary, Fradin's study rigorously establishes that while bilayer models are excellent approximations for stable stratified flows, they fail as rigorous limits for sheared flows in standard Sobolev spaces due to the explosive nature of Kelvin-Helmholtz instabilities, a phenomenon that persists even in the continuous, smoothed-out version of the problem.