Non-Hydrostatic Projection Operators

• Non-hydrostatic projection operators are mathematical tools that decompose pressure fields and enforce divergence-free conditions in complex fluid models.
• They extend classical projection methods to address vertical accelerations, dispersive effects, and tensorial anisotropies in geophysical and relativistic flows.
• Their numerical implementation using discrete stencils and spectral decompositions enhances stability, conservation, and accuracy in simulations of complex topographies.

Non-hydrostatic projection operators are fundamental computational and analytical constructs for enforcing constraints within non-hydrostatic fluid models. They arise in geophysical flows, computational fluid dynamics, meteorological modeling, oceanography, and relativistic hydrodynamics, wherever it is necessary to maintain divergence-free velocity fields or to decompose dynamics into component modes orthogonal with respect to conserved invariants. While their historical origin can be traced to classical projection and operator-splitting methods for incompressible flows, non-hydrostatic projection operators are specifically adapted to regimes where vertical accelerations, dispersive effects, or tensorial anisotropies render the hydrostatic approximation invalid.

1. Mathematical Foundation and General Framework

In a typical non-hydrostatic extension of shallow water or Boussinesq-type equations, the pressure field is decomposed as

$$P(x, z, t) = \rho\, g\,(\eta(x, t) - z) + P^{\rm nh}(x, z, t),$$

where $P^{\rm nh}$ is the non-hydrostatic pressure deviation. The non-hydrostatic projection operator acts as a Lagrange multiplier to enforce an elliptic (often Poisson-type) constraint on an auxiliary field—usually the pressure or its vertical moment—such that the divergence condition or its generalization holds at the discrete or continuous level.

For instance, in the depth-averaged non-hydrostatic model with moving bottom topography, a non-hydrostatic pressure projection operator $\mathcal{P}$ is defined via the elliptic problem

$$-\bigl(H(x)\,p^{\rm nh}_x\bigr)_x + \alpha(x)\,p^{\rm nh} = R(x),$$

where $H(x) > 0$, $\alpha(x) > 0$, and $R(x)$ is a right-hand side derived from intermediate variables and constraint residuals. This construction extends to tensorial decompositions in relativistic hydrodynamics and to the separation of wave and geostrophic modes in non-hydrostatic stratified flows (Firdaus et al., 31 Oct 2024, Escalante et al., 2017, Early et al., 22 Nov 2025, Subias, 2015, Aissiouene et al., 2015, Florkowski et al., 2011).

2. Operator Construction and Discretization Techniques

Projection operators are constructed to enforce instantaneous constraints, such as

• divergence-free depth-integrated velocities,
• incompressibility with non-hydrostatic corrections,
• or the elimination of certain variables via integral-differential constraints in vertical columns.

The canonical Chorin–Temam projection/correction approach involves two steps:

1. Predictor step: The hyperbolic part of the system (usually the hydrostatic approximation) is advanced explicitly, ignoring non-hydrostatic terms.
2. Corrector (projection) step: The non-hydrostatic pressure is determined from an elliptic problem such that the updated fields satisfy the divergence or closure constraint.

On discrete grids, operators such as the discrete gradient $G$, divergence $D$, and Laplacian $L=DG$ are assembled (typically as finite volume or finite difference stencils), yielding symmetric positive-definite systems for pressure. For structured vertical discretizations, as in atmospheric and NWP models, vertical finite-element projection operators based on B-splines are constructed, precisely enforcing algebraic versions of analytic vertical constraints like the so-called $C^1$ constraint

$$\mathcal{G}^*\,\mathcal{S}^* - \mathcal{G}^* - \mathcal{S}^* + \mathcal{N}^* = 0,$$

ensuring that semi-implicit iterations remain consistent and that single-variable reductions are feasible (Subias, 2015).

3. Non-Hydrostatic Projection in Shallow Water and Boussinesq Models

Non-hydrostatic projection methods provide dispersive corrections to standard shallow water systems. For moving bottom problems, the pressure profile can be quadratic in the vertical, leading to relations such as

$$P^{\rm nh}|_{z=-d} = \frac{6}{4 + d_x^2}\,p^{\rm nh} + \frac{d_x}{4 + d_x^2}\,(h\,p^{\rm nh})_x + \phi(x, t),$$

with $\phi(x, t)$ ensuring closure without unknown time derivatives. Substitution into the momentum constraint and divergence-free condition yields a symmetric, coercive block operator system. The projection operator $\mathcal{P}$ thus realizes a map from explicit forcing terms (predictor residuals) to the non-hydrostatic pressure, enabling stable, high-order, and well-posed solutions even in the presence of complex topographies and moving boundaries (Firdaus et al., 31 Oct 2024).

This operator is critical for dispersive wave modeling (e.g., tsunamis, landslides) since lower-order closures (linear or simplified quadratic profiles) compromise dispersive accuracy and may lead to ill-posedness on variable bathymetry.

4. Vertical Operator Theory and Spectral Decomposition

In models requiring modal or spectral decomposition (e.g., for analyzing wave–geostrophic energy exchanges in rotating, stratified, non-hydrostatic flows), non-hydrostatic projection operators serve a basis-separating function. The approach projects the state vector onto orthogonal eigenmodes determined by energy and enstrophy inner products, such as

$$\langle \psi_1, \psi_2 \rangle_E = \int (u_1 u_2 + v_1 v_2 + w_1 w_2 + N^2 \eta_1 \eta_2) \, dV,$$

separating wave-like and balanced (geostrophic) motions via projections $A_0^{k\ell j}$ and $A_\pm^{k\ell j}$, constructed by inner-product integration against modal vertical structures (Early et al., 22 Nov 2025).

In NWP vertical columns, B-spline-based finite element spaces are built to preserve invertibility and exact algebraic constraints on the vertical operators. Specialized bases (e.g., $\Xi_{i,k}$ and $\Sigma_{i,k}$) enable machine-precision enforcement of algebraic constraints, leading to banded or Krylov-iterative solvers with robust convergence in operational codes (Subias, 2015).

5. Algorithmic Implementation and Computational Aspects

Implementation follows a predictor–corrector splitting:

1. Compute predictor by integrating the nonlinear (typically hydrostatic) equations.
2. Assemble right-hand side for the elliptic system (involving explicit velocity, water depth, and bottom geometry).
3. Solve elliptic projection problem for the non-hydrostatic pressure (matrix-free on GPUs (Escalante et al., 2017), direct or iterative solution for banded matrices in vertical projection (Subias, 2015)).
4. Update velocity and other prognostic fields using the projected pressure to enforce divergence or associated constraint.

Boundary conditions for projection operators are chosen to ensure well-posedness (Dirichlet for open boundaries, Neumann for impermeable walls), and compatibility is ensured by checking that the right-hand-side has zero integral when necessary.

For efficiency in large-scale systems (GPU or multi-core), sparse, matrix-free, and preconditioned iterative methods (CG, BiCGStab) dominate. Formal second-order accuracy in both space and time is routinely achieved, with well-balancing and depth-positivity in shallow water models (Escalante et al., 2017, Aissiouene et al., 2015).

6. Extensions: Tensorial and Relativistic Projections

Beyond incompressible and stratified flow, projection operator methodology generalizes to relativistic and anisotropic fluid dynamics. In boost-invariant, cylindrically symmetric flows, the projection operator formalism decomposes tensors using a locally orthonormal tetrad $(u^\mu, x^\mu, y^\mu, z^\mu)$, with one-, two-, and four-index projectors orthogonalizing expansion, shear, and stress tensors: $$\Delta^{\mu\nu} = g^{\mu\nu} - u^\mu u^\nu,\quad X^{\mu\nu} = x^\mu x^\nu, \ldots$$ enabling direct reduction of the governing PDEs to closed systems of scalar equations, as in the Israel–Stewart and ADHYDRO frameworks. Such operators ensure orthogonality and closure even in strong dissipation or high anisotropy regimes, with explicit matching procedures between models of differing kinetic structure (Florkowski et al., 2011).

7. Applications, Limitations, and Comparative Accuracy

Non-hydrostatic projection operators are deployed in tsunami generation, landslide-induced waves, operational weather forecast models (ALADIN-HIRLAM), global and regional ocean circulation, and relativistic heavy-ion collision models. Their success depends on the consistent enforcement of algebraic constraints (for accuracy and stability), compatibility with energy and enstrophy conservation, and ability to resolve dispersive effects beyond the hydrostatic limit.

Limitations are set by approximations in vertical structure (e.g., linear or quadratic pressure/velocity profiles), discretization errors (notably in the presence of severe topographic gradients or stratification jumps), and assumptions required for modal completeness or boundary closure. Loss of $C^1$ constraint in vertical discretization, neglect of nonlinear or small-scale effects in mode-projected energy budgets, and incompatibility with ill-posed pressure relations on non-flat topography are key challenges (Firdaus et al., 31 Oct 2024, Subias, 2015, Early et al., 22 Nov 2025).

A plausible implication is that, for physically and numerically consistent non-hydrostatic models, exact algebraic enforcement via projection operators is essential, and the choice of operator form directly governs the stability, accuracy, and conservation properties of the overall scheme.