A potential energy conserving finite element method for turbulent variable density flow: application to glacier-fjord circulation

TL;DR

This work develops a continuous Galerkin finite element method for variable-density, gravity-driven flows that conserves total energy and angular momentum without enforcing strict divergence-free conditions. It introduces a shift-invariant SI-MEEDMAC formulation with a new symmetric tensor viscosity, constructed via a residual-viscosity approach and augmented with high-order dissipation to act as an implicit LES. The method is demonstrated on glacier–fjord circulation problems, including a 2D Sherard Osborn fjord test, where it yields improved energy balance, reduced artificial diffusion, and resolution of turbulent features compared with a high-resolution FV model (MITgcm). The results highlight the potential of unstructured-CG approaches for ice–ocean interactions and motivate future extensions to Coriolis effects, more realistic melt physics, and challenging mesh configurations.

Abstract

We introduce a continuous Galerkin finite element discretization of the non-hydrostatic Boussinesq approximation of the Navier-Stokes equations, suitable for various applications such as coastal ocean dynamics and ice-ocean interactions, among others. In particular, we introduce a consistent modification of the gravity force term which enhances conservation properties for Galerkin methods without strictly enforcing the divergence-free condition. We show that this modification results in a sharp energy estimate, including both kinetic and potential energy. Additionally, we propose a new, symmetric, tensor-based viscosity operator that is especially suitable for modeling turbulence in stratified flow. The viscosity coefficients are constructed using a residual-based shock-capturing method and the method conserves angular momentum and dissipates kinetic energy. We validate our proposed method through numerical tests and use it to model the ocean circulation and basal melting beneath the ice tongue of the Ryder Glacier and the adjacent Sherard Osborn fjord in two dimensions on a fully unstructured mesh. Our results compare favorably with a standard numerical ocean model, showing better resolved turbulent flow features and reduced artificial diffusion.

Figures (6)

• Figure 1: Time evolution of (a) kinetic + potential energy and (b) $L^2$-error for $T_h$ using the SI-MEDMAC formulation (\ref{['eq:ns_EMAC']}) (red) and the SI-MEEDMAC formulation (\ref{['eq:ns_EMAC_potential']}) (blue). No stabilization is used and $\nu = \kappa = 0$. The unstructured mesh consists of $4909$${\mathbb P}_3$ and $2209$${\mathbb P}_2$ nodes.
• Figure 2: Unstructured mesh used for the Ryder simulation at $x \in [18.8km, 21.7 km]$.
• Figure 3: Time-averaged melt rate on the ice-shelf (red line) and time-averaged streamlines (black) superimposed on the density profile. The averages are taken from day 10 until the simulation is ended ($\approx 35$ days).
• Figure 4: From top to bottom, snapshots of the discontinuity indicator for temperature in logarithmic scale (first), temperature (second), velocity magnitude (third) at $t =$ 10 days at the place where most melting occurs ($x \approx 15$ km). The unstructured mesh (bottom) is locally refined close to the ice shelf with $h \approx 3.5 m$. An animation is available in the supplementary material online.
• Figure 5: Time-averaged horizontal velocity (a), temperature (b) and salinity (c) at $x =$ 21 km as a function of depth. The average is taken from day 10 until the simulation is ended ($\approx 35$ days). Initial and open-ocean boundary condition profiles are also displayed in (b) and (c).

Theorems & Definitions (8)

• Remark 2.1
• Remark 3.1
• Theorem 4.1
• proof
• Theorem 4.2
• proof
• Proposition 4.1
• proof