🔢 mathematics

A consistent-splitting generalized scalar auxiliary variable scheme for the perturbed Boussinesq system

This paper proposes and analyzes a second-order consistent-splitting scheme based on the generalized scalar auxiliary variable approach for the perturbed Boussinesq system, proving its unconditional weak stability and optimal error estimates while demonstrating its ability to accurately capture internal-wave dynamics and hydrostatic relaxation through numerical experiments.

Original authors: M Nader Alhomsi, Jiahong Wu, Xiaoming Zheng

Published 2026-07-01

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: M Nader Alhomsi, Jiahong Wu, Xiaoming Zheng

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the Earth's atmosphere and oceans as a giant, swirling pot of soup. Sometimes, this soup is calm and layered (like a perfectly stratified drink where the heavy stuff is at the bottom and light stuff is at the top). Sometimes, you stir it, creating waves, swirls, and turbulence. Scientists use a set of complex math rules called the Boussinesq system to predict how this "soup" moves, especially when heat and gravity are involved.

This paper introduces a new, smarter way to solve these math rules on a computer. Here is the breakdown in everyday terms:

1. The Problem: Stirring the Pot

The authors are studying a specific version of this "soup" called the perturbed Boussinesq system.

The Setup: Imagine the soup is already settled into a stable, layered state (like a calm ocean).
The Disturbance: Now, you drop a hot rock or a cold ice cube into it. This creates a "perturbation"—a ripple or a wave that tries to mix the layers.
The Challenge: The math describing how these ripples move, bounce, and eventually settle back down is incredibly difficult. It involves two main things fighting each other:
1. The Flow: The fluid moving around (like traffic).
2. The Heat: Temperature differences causing the fluid to rise or sink (like a hot air balloon).

2. The Solution: A "Splitting" Trick

Computers are bad at solving giant, tangled knots of math all at once. If you try to calculate the flow and the heat simultaneously at every tiny moment, the computer gets overwhelmed or takes forever.

The authors propose a Consistent-Splitting Scheme. Think of this like a chef who doesn't try to chop vegetables, boil water, and fry meat all in the same second. Instead, they do it in a specific, quick sequence:

Step A: Calculate where the water wants to go next based on where it is now.
Step B: Calculate how the heat changes based on that movement.
Step C: Adjust the pressure to make sure the water doesn't magically appear or disappear (keeping the "soup" volume constant).

By breaking the problem into these small, separate steps, the computer can solve simple, linear puzzles instead of one giant, impossible one.

3. The Secret Sauce: The "GSAV" Variable

To make sure the computer doesn't crash or give nonsense answers (like infinite energy), they use a tool called the Generalized Scalar Auxiliary Variable (GSAV).

The Metaphor: Imagine a "safety monitor" or a "thermostat" attached to the simulation. This monitor constantly checks the total energy of the system.
How it works: If the math starts to drift and the energy looks like it's going to explode, this "thermostat" gently nudges the numbers back into a safe zone. It ensures that even if the computer takes a big step, the simulation stays physically realistic and doesn't blow up.

4. What They Found (The Results)

The authors proved two main things about their new method:

It's Accurate: They showed that if you make the time steps smaller, the answer gets closer to the truth very quickly (specifically, "second-order" accuracy). It's like zooming in on a blurry photo; the more you zoom, the sharper the picture gets, and this method sharpens it faster than older methods.
It's Stable (Mostly): They proved the method won't crash for any reasonable time step size, provided the "safety monitor" (the constant C) is set high enough.

5. The Catch: The "Viscosity" Trap

There is one major limitation they discovered, which is like a "speed limit" on how fast the computer can run when the fluid gets very slippery.

The Issue: In the real world, fluids have "stickiness" (viscosity) and "heat spreading" (diffusivity). If these values are very low (meaning the fluid is very slippery and heat doesn't spread easily, like in the upper atmosphere), the math gets extremely sensitive.
The Analogy: Imagine trying to balance a pencil on its tip. If the table is sticky, it's easy. If the table is perfectly smooth ice, the pencil falls over instantly.
The Result: Their method works great for normal fluids. But as the fluid becomes "slipperier" (viscosity approaches zero), the computer needs to take tinier and tinier steps to stay stable. The math shows that the error constant grows quadruply exponentially (a huge, terrifying number) as the fluid gets less sticky. This means the method isn't "robust" for extremely slippery fluids without using a massive amount of computing power.

6. The Simulation: Watching the Waves

To test their theory, they ran a long-term simulation (like watching a movie of the fluid for a long time).

Phase 1 (The Splash): They dropped a hot and a cold "blob" into the fluid. The fluid surged, converting heat energy into movement.
Phase 2 (The Ripples): The fluid started to oscillate, creating internal waves (like ripples inside a glass of water). These waves bounced back and forth with a predictable rhythm.
Phase 3 (The Calm): Eventually, friction (viscosity) slowed everything down, and the fluid settled back into its calm, layered state.

Summary

The paper presents a clever, fast, and accurate way to simulate how heat and fluid interact in a layered environment. It uses a "divide and conquer" strategy with a safety monitor to keep things stable. While it works beautifully for normal conditions, it struggles when the fluid becomes extremely slippery, requiring the computer to work much harder to avoid errors.

Technical Summary: A Consistent-Splitting Generalized Scalar Auxiliary Variable Scheme for the Perturbed Boussinesq System

Problem Statement
This work addresses the numerical simulation of the two-dimensional perturbed Boussinesq system, which models buoyancy-driven flows in stably stratified environments (e.g., atmospheric and oceanic dynamics). The system is derived by subtracting a stable, linearly-stratified hydrostatic equilibrium from the standard Boussinesq equations. The resulting perturbed system (Equation 6) features a specific coupling where buoyancy drives vertical motion in the momentum equation, while vertical motion advects the linear background temperature gradient in the temperature equation. This coupling generates internal gravity waves and drives the relaxation of disturbances toward hydrostatic balance. The primary mathematical challenge lies in the incompressibility constraint and the saddle-point structure of the velocity-pressure system, compounded by the need for high-order accuracy and stability in the presence of nonlinear convection and linear stratification terms.

Methodology
The authors propose a second-order consistent-splitting scheme based on the Generalized Scalar Auxiliary Variable (GSAV) approach. The methodology integrates several key components:

GSAV Framework: A scalar auxiliary variable $r(t)$ is introduced to mimic the continuous energy law of the system. This allows for the construction of a modified discrete energy that remains unconditionally stable.
Consistent-Splitting Time Discretization: The scheme extends the time-shifted second-order BDF2 (Backward Differentiation Formula) framework developed by Huang and Shen [17] for Navier–Stokes equations.
- Explicit Treatment: Nonlinear convection ( $u \cdot \nabla u$ ), advection ( $u \cdot \nabla \theta$ ), and linear couplings (buoyancy $\theta e_2$ and stratification $\alpha v$ ) are treated explicitly using second-order extrapolation at the shifted time level $t_{n+k}$ .
- Decoupling: This explicit treatment reduces each time step to solving a small number of decoupled linear systems with constant coefficients, avoiding the need for nonlinear solvers or monolithic coupled systems.
Pressure Correction: A pressure Poisson equation is solved using a consistent splitting approach that retains optimal second-order accuracy, avoiding the boundary condition inconsistencies often found in standard projection methods.
Stability and Error Analysis: The authors prove an unconditional weak stability theorem and derive optimal second-order error estimates for velocity, pressure, and temperature. The analysis rigorously tracks the dependence of error constants on the kinematic viscosity ( $\nu$ ) and thermal diffusivity ( $\gamma$ ).

Key Contributions
The paper makes three distinct contributions compared to existing literature on SAV-type schemes for Boussinesq systems:

Perturbed System Analysis: Unlike previous works that analyze the standard Boussinesq system, this paper specifically addresses the perturbed system (6) containing the linear stratification coupling term $\alpha v$ in the temperature equation.
Framework Integration: While Wagner et al. [36] applied the GSAV approach to the standard Boussinesq system, this work adopts the specific time-shifted second-order BDF consistent-splitting framework of Huang and Shen [17] but retains the standard GSAV update mechanism, ensuring full second-order accuracy in both velocity and temperature equations through consistent extrapolation of all coupling terms.
Explicit Error Constant Dependence: A critical contribution is the explicit tracing of the error constant's dependence on $\nu^{-1}$ and $\gamma^{-1}$ . The analysis reveals that the error constant contains negative powers of viscosity and diffusivity, leading to a quadruply-nested exponential dependence on these parameters.

Results

Theoretical Stability: The scheme is proven to be unconditionally weakly stable. The sequences of the auxiliary variable and the scaling factor remain non-negative, and the numerical solutions satisfy uniform bounds independent of the time step size, provided the stabilization constant $C$ exceeds a threshold determined by the forcing terms.
Convergence: Optimal second-order convergence rates are established for the velocity, pressure, and temperature in appropriate norms.
Numerical Experiments:
- Convergence Test: A manufactured solution on a bounded domain confirms the theoretical second-order convergence in both $L^2$ and $H^1$ norms.
- Stratified Flow Simulation: A long-time simulation ( $T=10^4$ $T = 1 0^{4}$ ) of a stratified flow with initial Gaussian temperature perturbations reproduces the expected physical dynamics:
  - A buoyancy transient regime where potential energy converts to kinetic energy.
  - An internal-wave oscillation regime with a period consistent with theoretical predictions ( $T \approx \pi\sqrt{2/\alpha}$ ).
  - An exponential decay regime where perturbations relax toward hydrostatic balance, with decay rates consistent with the analytical bounds derived in Lemma 2.

Significance and Limitations
The paper claims that the proposed scheme offers a computationally efficient method for simulating stratified flows, reducing each time step to decoupled linear solves while maintaining second-order accuracy. However, the authors are modest regarding the scheme's robustness in the limit of vanishing viscosity or diffusivity.

The analysis demonstrates that the error constant grows as a quadruply-nested exponential of $\nu^{-1}$ and $\gamma^{-1}$ . Consequently, the scheme is not robust as $\nu, \gamma \to 0$ . The numerical experiments confirm that for small viscosity and thermal diffusivity, the scheme requires increasingly fine spatial and temporal resolutions to remain stable, and it may blow up at coarser resolutions. The authors conclude that while the scheme is effective for fully dissipative regimes, a fully implicit or stabilized treatment of convection terms would be necessary to obtain error bounds uniform in $\nu$ and $\gamma$ . This limitation highlights a trade-off between the computational efficiency of the explicit consistent-splitting approach and its stability in low-diffusion regimes.