Free-surface Stokes problem: stability estimates and time-step improvements

Technical Summary: Free-surface Stokes problem: stability estimates and time-step improvements

Problem Statement
Accurate simulation of highly viscous free-surface flows, such as ice sheet dynamics, mantle convection, and lava flow, requires solving the coupled Stokes/free-surface equations. These problems are characterized by a moving domain where the fluid velocity and pressure are coupled with the evolution of the free-surface height. While standard discretizations using finite elements in space and Euler methods in time are common, they often suffer from severe stability constraints. Specifically, explicit and semi-implicit Euler schemes frequently require prohibitively small time steps to remain stable, leading to high computational costs. Existing stabilization techniques, such as the Free-Surface Stabilization Algorithm (FSSA), improve time step sizes but lack a priori theoretical bounds on stability and do not guarantee unconditional stability or volume conservation.

Methodology
The authors present a rigorous theoretical analysis of the stability and conservation properties of the coupled Stokes/free-surface system for both Newtonian and non-Newtonian (power-law viscosity) fluids. The analysis is conducted at two levels:

1. Continuous Level: The authors derive a total energy balance that couples the energy of the Stokes momentum equations with the energy of the free-surface height equation. This derivation relies on establishing specific relations between the free-surface advection term and surface integrals, and subsequently linking these to the gravity term in the Stokes equations.
2. Discrete Level: The study examines fully discrete schemes using finite elements in space combined with explicit, semi-implicit, and implicit Euler time-stepping.
   • Regularization: To address the lack of inherent regularity in the free-surface height function (which is required for the well-posedness of the Stokes problem), the authors introduce a high-order edge stabilization term, $J(s^{n+1}, w)$, treated implicitly in time.
   • Stabilization Proposal: Motivated by the identification of "spurious" terms in the energy estimates of the explicit Euler method, the authors propose a new stabilization term added to the Stokes weak form. This term penalizes the normal velocity at the free surface and is designed to precisely cancel the destabilizing terms arising from the explicit treatment of the free-surface equation.

Key Contributions

• Total Energy Balance: The paper derives a novel total energy balance for the coupled Stokes/free-surface system, coupling the momentum energy with the surface height energy. This is identified as a fundamental building block for understanding stability, which was previously missing in the literature.
• Theoretical Stability Proofs:
  • Implicit Euler: Proven to be unconditionally stable with respect to the time step size (assuming Lipschitz continuous surface height) and to conserve domain volume.
  • Explicit Euler (Unstabilized): Proven to be unstable in the total energy norm for any time step size due to spurious terms related to the normal surface velocity.
  • Semi-implicit Euler (FSSA): Shown to be only conditionally stable and, critically, proven not to conserve domain volume.
  • Proposed Stabilized Explicit Euler: The authors prove that their new stabilization term renders the explicit Euler method unconditionally stable with respect to the time step size (under the Lipschitz continuity assumption) while preserving domain volume.
• Volume Conservation: The paper provides theoretical proofs and numerical verification that the proposed stabilized explicit scheme conserves the domain volume, whereas the FSSA-stabilized semi-implicit scheme does not.

Results
Numerical experiments were conducted using the open-source library FEniCS for both Newtonian and non-Newtonian fluid cases, including a complex geometry based on Greenland ice sheet data.

• Stability: The unstabilized explicit Euler scheme was unstable for all tested time steps. The FSSA-stabilized schemes (both semi-implicit and explicit) were conditionally stable, failing at larger time steps. In contrast, the proposed EE-Stabilized scheme remained stable for arbitrarily large time steps in both Newtonian and non-Newtonian cases.
• Volume Conservation: The EE-Stabilized scheme conserved domain volume up to machine precision. The SIE-FSSA scheme failed to conserve volume, confirming the theoretical findings.
• Accuracy: Convergence tests showed that the proposed stabilization terms are consistent with the discretized problem. The EE-Stabilized scheme maintained first-order accuracy in surface height and velocity under simultaneous mesh and time-step refinement, with errors comparable to or smaller than other stabilized schemes.
• Oscillations: The implicit edge stabilization term effectively suppressed high-frequency oscillations in the free-surface height without introducing additional time-step restrictions.

Significance and Claims
The paper claims that the proposed stabilized explicit Euler scheme offers a robust and efficient alternative to computationally expensive implicit solvers. By theoretically guaranteeing unconditional stability with respect to the time step size and ensuring volume conservation, the method enables stable simulations of free-surface flows even with complex geometries and source terms. The authors emphasize that their stabilization is the first to theoretically guarantee arbitrary time step sizes for the explicit Euler discretization of this coupled problem. The approach is noted for being symmetric (beneficial for iterative solvers), simple to implement, and free of tunable parameters. The authors conclude that while their analysis assumes Lipschitz continuity of the surface height, the derived $L^2$-norm stability estimates are a necessary foundation for well-posedness once appropriate regularization is applied. Future work is suggested to explore regularization terms for well-posedness and to extend this stability analysis to Navier-Stokes/free-surface problems and higher-order time integrators.