Obsidian/Zotero/001_artiklid/tominecstabilityanalysisfreesurface2025.md

---

category: literaturenote citekey: tominecstabilityanalysisfreesurface2025 title: Stability analysis of the free-surface Stokes problem and an unconditionally stable explicit scheme authors: "Tominec, Igor; Lundgren, Lukas; Löfgren, André; Ahlkrona, Josefin" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.10447v1" zotero_key: 6VV3UEPS zotero_storage: ME7TJ7D6 collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Tominec, Igor"

status: converted

Stability analysis of the free-surface Stokes problem and an unconditionally stable explicit scheme

IGOR TOMINEC∗,1,² , LUKAS LUNDGREN1,2,³ , ANDRE´ LOFGREN ¨ 1,2 , JOSEFIN AHLKRONA1,2,³

¹ Department of Mathematics, Stockholm University, Sweden ² Bolin Centre for Climate Research, Stockholm University, Sweden ³ Swedish e-Science Research Centre, Sweden

*Corresponding author: igor.tominec@math.su.se

Accurate simulations of ice sheet dynamics, mantle convection, lava flow, and other highly viscous freesurface flows involve solving the coupled Stokes/free-surface equations. In this paper, we theoretically analyze the stability and conservation properties of the weak form of this system for Newtonian fluids and non-Newtonian fluids, at both the continuous and discrete levels. We perform the fully discrete stability analysis for finite element methods used in space with explicit and implicit Euler time-stepping methods used in time. Motivated by the theory, we propose a stabilization term designed for the explicit Euler discretization, which ensures unconditional time stability and permits conservation of the domain volume. Numerical experiments validate and support our theoretical findings.

Keywords: stability, free-surface flow, Stokes problem, time-stepping, finite element method

1. Introduction

Free-surface flows of Newtonian and non-Newtonian fluids in the low Reynolds number regime are commonly modelled by the Stokes equations coupled with a free-surface evolution equation, resulting in a moving domain problem. Real-world applications include ice sheet dynamics [14], mantle convection [18], lava flow [15], interactions between viscous drops [23], and others. These applications demand numerically robust simulations that remain stable across a range of mesh resolutions and time step sizes. However, standard discretizations of the coupled Stokes/free-surface problem may suffer from instabilities that constrain the simulations to small time step sizes, which then induce large computational costs. In this paper, we rigorously analyze the stability of the coupled Stokes/free-surface system and, motivated by our theoretical results, propose new numerical stabilization approaches. These enable stable and efficient simulations of free-surface flows, even in the presence of complex geometries and source terms.

Related work. Finite element methods in space combined with explicit, semi-implicit, and implicit Euler discretizations in time, are a common choice for numerically solving the coupled Stokes/freesurface problem. The implicit Euler method allows for the largest time step sizes, but is computationally expensive as it requires nonlinear iterations coupling the Stokes/free-surface problem. On the other hand, explicit Euler and semi-implicit Euler discretizations in time [9, 12, 18, 31] do not require additional nonlinear iterations. However, these approaches lead to severe restrictions on the maximum stable time step size. Free-Surface Stabilization Algorithm (FSSA) introduced in [18] is an approach for improving the maximum time step size of the semi-implicit Euler discretization in time. The FSSA adds an unsymmetric surface integral term to the weak form Stokes problem to bring the full discretization closer to an implicit method. The approach was first introduced in the context of the mantle convection applications [18] and was later extended to ice sheet dynamics applications in [20, 22, 28]. While the authors report significant time step size improvements, the time step restriction is not a priori known and the overall theoretical understanding of FSSA is still missing. In [4], the FSSA stabilization terms were modified to preserve the symmetry of the linear systems. In [26], the authors give a procedure for computing the largest stable time-step for the Stokes/free-surface system simplified to an associated linear homogeneous eigenvalue problem, by modifying the FSSA stabilization terms [18]. In [5], the

FIG. 1. Plot (a) shows an example of domain Ω ⊂ R ^d with boundary split into Γ^s (free-surface), Γ^b (bottom topography/bedrock), and Γ^l (the lateral boundaries). Plot (b) shows an explicit form free-surface height function h = h(x ^⊥,t), where x ^⊥ is evaluated over free-surface domain Ω^⊥ ⊂ R *d*−1

authors were the first to examine the numerical stability properties of the coupled Navier-Stokes/freesurface problem, discretized using the finite element method and the explicit Euler method. To improve the maximum time step size, the authors supplemented the free-surface equation and the Navier-Stokes problem with additional first-order consistent stabilization terms. Their approach improves the time step sizes up to a CFL-type restriction that can be evaluated up to an unknown constant. Furthermore, the authors consider Newtonian fluids without a source term in the free-surface equation. In [7], the author gives surface height error estimates for the implicit Euler discretized coupled Stokes/free-surface problem for non-Newtonian fluids, however, the stability properties of computationally less expensive approaches were not investigated. Moreover, the estimate coupling the momentum energy balance and the surface height energy through a total energy norm is missing.

Contributions of the present paper. In this paper, we present a novel theoretical study of the stability of the coupled Stokes/free-surface problem for Newtonian and non-Newtonian fluids with a source term, in both the continuous and fully discrete settings. On the continuous level, we derive a total energy balance that couples the energy of the Stokes problem with that of the free-surface equation, which is key to understanding the stability. From this, we establish a stability estimate and the conservation of domain volume. On the discrete level, we show that the implicit Euler method is unconditionally stable and preserves domain volume. By contrast, we provide numerical counterexamples showing that the semi-implicit Euler method does not conserve domain volume. Moreover, the FSSA from [18] only partially mitigates the stability issues of both the explicit and the semi-implicit Euler schemes, resulting in a conditional time step constraint. Next, we prove that the explicit Euler method conserves domain volume, but is not stable in the total energy norm. To address this limitation, we use our stability estimates to derive an optimally scaled stabilization term that penalizes the normal velocity at the free surface in the Stokes problem. This stabilization ensures that the explicit Euler method used for advancing the free-surface equation is unconditionally stable with respect to the time step size. Our formulation also extends to cases with a source term in the free-surface equation. In addition, we prove that the scheme preserves the volume of the domain and numerically verify that it maintains first-order accuracy under simultaneous refinement of the mesh and the time step. Our analysis does not rely on any simplifications of the geometry or the original Stokes/free-surface formulation and is valid in both two and three spatial dimensions.

Paper organization. This paper is structured as follows. In Section 2 we state the Stokes problem coupled to the free-surface equation in strong and weak form and define the domains and their boundary representations. In Section 3 and Section 4 we respectively write the continuous and the discrete theoretical stability analysis. In Section 5 we provide numerical experiments verifying our theory for the Newtonian fluid case and non-Newtonian fluid case (simplified ice sheet dynamics). In Section 6 we give final remarks.

2. Mathematical model

2.1. Governing equations

The fluid is represented as an open and bounded domain Ω = Ω(t) ⊂ R d that evolves in time t ∈ T = [0,*t*ˆ]. The evolving free boundary (one part of ∂Ω) is represented as the surface height h(x ^⊥,t) with x ⊥ evaluated over a horizontal (flat) domain Ω^⊥ ⊂ R *d*−1 . The remaining part of ∂Ω is a topography/bedrock b(x ^⊥). The area between h and b is the fluid thickness *h*−b > 0. The relation between Ω, h, and Ω^⊥ is illustrated in Figure 1. The coupled Stokes/free-surface problem reads. Find the fluid surface height h : Ω^⊥ ×T → R, velocity u : Ω → R d , and pressure π : Ω → R, such that:

$$\begin{cases} \partial_t h + \mathbf{u}^{\perp} \cdot \nabla^{\perp} h - u_z = a, \ -\nabla \cdot (2\mu(\mathbf{u})\mathbf{D}\mathbf{u}) + \nabla \pi = -\rho g \hat{\mathbf{z}}, \ \nabla \cdot \mathbf{u} = 0, \end{cases}$$ (2.1)

where u ^⊥ : Ω^⊥ ×T → R *d*−1 is a vector of the horizontal surface velocities, ∇ ^⊥ is the horizontal gradient operator over Ω⊥, u^z is the vertical surface velocity, a : Ω^⊥ ×T → R is the source term, ∇ is the gradient operator over Ω, Du = ∇u + (∇u) T is the strain rate tensor, g > 0 is the gravitational acceleration, ρ > 0 is the fluid density, and µ : Ω → R ⁺ is a viscosity function. We consider a family of viscosity functions in the power-law form:

$$\mu(\mathbf{u}) = \mu_0 |\mathbf{D}\mathbf{u}|^{p-2}, \quad 1 (2.2)$$

where µ⁰ ∈ R ⁺ is an effective viscosity constant, |Du| = (∑ d i, j=1 (Du)i j(Du)i j) 1/2 is the Frobenius norm of Du. Glaciologists use Glen's flow law exponent n in place of p, where the relation is p = (1/n) +1 [13, 16]. When p = 2, the fluid is Newtonian. When 1 < p < 2, the fluid is non-Newtonian. We supplement (2.1) with an initial condition for the surface h x ^⊥,0 = h⁰ x ⊥ , and the boundary conditions:

$$\mathbf{u} = \mathbf{0} \text{ on } \Gamma_b, \quad \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \Gamma_l, \quad \boldsymbol{\sigma} \cdot \mathbf{n} = \mathbf{0} \text{ on } \Gamma_s,$$ (2.3)

where Γ^b is the bedrock, Γ^s is the free-surface, Γ^l are the lateral boundaries, and σ = 2µ(u)Du − *p*I is the deviatoric stress tensor.

2.2. Weak formulation and function spaces

We analyze the weak formulation of the coupled problem (2.1), which is defined as follows. Find u ∈ V (Ω), π ∈ Q(Ω), and h ∈ Z(Ω⊥)×C 1 (T) such that:

$$\begin{cases} (\partial{t}h, w){\Omega^{\perp}} + \left(\mathbf{u}^{\perp} \cdot \nabla^{\perp}h - u{z}, w\right){\Omega^{\perp}} = (a, w){\Omega^{\perp}}, \ 2\left(\sqrt{\mu} \, \mathbf{D} \mathbf{u}, \sqrt{\mu} \, \mathbf{D} \mathbf{v}\right){\Omega(t)} - (\pi, \nabla \cdot \mathbf{v}){\Omega(t)} = -\rho g(\hat{\mathbf{z}}, \mathbf{v}){\Omega(t)}, \ (\nabla \cdot \mathbf{u}, q)_{\Omega(t)} = 0, \end{cases} (2.4)$$

for all v ∈V (Ω), q ∈ Q(Ω), and w ∈ Z(Ω⊥). We now define the natural solution spacesV , Q, andZ. The Lebesgue space L α (Ω) is defined as [11]: L α (Ω) = {q : Ω → R|q measurable and ∥*q*∥ α L ^α (Ω) < ∞}, α = [1,∞), where an associated L α -norm is ∥*q∥^L* ^α (Ω) = (R Ω |q| ^α dx) 1/α . The Sobolev space W ¹,^α (Ω) and its restriction to vanishing functions over Γ^b ⊂ ∂Ω are respectively defined through:

$$\begin{split} \pmb{W}^{1,\alpha}(\Omega) &= \left{ \pmb{\nu} \in [L^{\alpha}(\Omega)]^d \, \middle| \, | \pmb{\nu} |{L^{\alpha}(\Omega)}^{\alpha} + | \nabla \pmb{\nu} |{L^{\alpha}(\Omega)}^{\alpha} < \infty \right} \ \pmb{W}^{1,\alpha}_{\Gammab}(\Omega) &= \left{ \pmb{\nu} \in \pmb{W}^{1,\alpha}(\Omega) \, \middle| \, \pmb{\nu} |{\Gamma_b} = 0 \right}. \end{split}$$

The solution to the nonlinear Stokes problem in weak form belongs to the velocity space V (Ω) and the pressure space Q(Ω), defined as [6]:

$$V(\Omega) = W_{\Gamma_b}^{1,p}(\Omega), \quad Q = L^{p'}(\Omega), \quad p' = \frac{p}{p-1}, \quad 1 (2.5)$$

where p is a parameter in (2.2).

By the Sobolev trace theorem, if u ∈ *W1,^p* (Ω), then u|∂^Ω ∈ L p (∂Ω). We assume that h belongs to L 2 Ω^⊥ . Accordingly, we define:

$$Z(\Omega^{\perp}) = L^2(\Omega^{\perp}),$$

where an L inner product is (q,w)^Ω = R Ω qwdx.

We approximate the solutions to the Stokes problem and the free-surface equation using the finite element method. The discrete domains Ω^h and Ω^⊥ h are represented by meshes T^h and T ^⊥ h . Each mesh is composed of a finite number of disjoint elements K and K ^⊥ so that Ω^h = ∪*K*∈ThK and Ω^{⊥ h} = ∪*K*∈^{T ⊥} K ⊥.

We denote E^h and E ⊥ h a set of interior d − 1 dimensional faces in T^h and T ^⊥ h respectively. The associated finite element spaces are the lowest-order Taylor-Hood velocity and pressure spaces, and a piecewise linear space for free-surface height:

$$\begin{split} \boldsymbol{V}_{\boldsymbol{h}}(\Omega_h) &= \Big{ \boldsymbol{v}_h \in \mathscr{C}^0(\Omega_h) \, \Big| \, \boldsymbol{v}_h \big|_K \in \mathbb{P}_2(K), \forall K \in \mathscr{T}_h \Big}, \ Q_h(\Omega_h) &= \Big{ q_h \in \mathscr{C}^0(\Omega_h) \, \Big| \, q_h \big|_K \in \mathbb{P}_1(K), \forall K \in \mathscr{T}_h \Big}, \ Z_h(\Omega_h^\perp) &= \Big{ w_h \in \mathscr{C}^0(\Omega_h^\perp) \, \Big| \, wh \big|{K^\perp} \in \mathbb{P}_1(K^\perp), \forall K^\perp \in \mathscr{T}_h^\perp \Big}, \end{split}$$

where P1, P² are multivariate polynomial spaces of order 1 and 2. The analysis naturally extends to higher-order Taylor-Hood elements and higher-order mesh representations.

2.3. Surface representations and integration domains

In this section, we describe how the domain, surface integrals, and normal vectors depend on the surface height h, as this will be required in the subsequent analysis. For d = 2 we define x = (x,z), x ^⊥ = x, whereas for d = 3 we define x = (x, y,z), x ^⊥ = (x, y). Then, the domain and its boundary are defined as:

$$\Omega = \left{ \boldsymbol{x} \in \mathbb{R}^d \mid \boldsymbol{x}^{\perp} \in \Omega^{\perp} \text{ and } b(\boldsymbol{x}^{\perp}) < z < h(\boldsymbol{x}^{\perp}) \right}, \quad \partial \Omega = \Gamma_s \cup \Gamma_b \cup \Gamma_l, \tag{2.6}$$

where:

$$\Gamma{s} = \left{ \left( \boldsymbol{x}^{\perp}, h(\boldsymbol{x}^{\perp}) \right) \in \mathbb{R}^{d} \, | \, \boldsymbol{x}^{\perp} \in \Omega^{\perp} \right}, \quad \Gamma{b} = \left{ \left( \boldsymbol{x}^{\perp}, b(\boldsymbol{x}^{\perp}) \right) \in \mathbb{R}^{d} \, | \, \boldsymbol{x}^{\perp} \in \Omega^{\perp} \right}, \quad \Gamma{l} = \partial \Omega \setminus (\Gamma{s} \cup \Gamma_{b}). \tag{2.7}$$

The volume of the domain is then defined as |Ω(t)| = R Ω(t) dx = R ^Ω⊥ (*h*−b)*dx*⊥.

When computing integrals over Γ^s , it is convenient to change the integration domain to the timeindependent evaluation domain Ω⊥. Let r^s = (x ^⊥,h(x ^⊥,t)) = (x, y,h(x, y,t)) be a parametrization of Γ^s . Then we have:

$$\int_{\Gammas} f(\mathbf{x}) ds = \int{\Omega^{\perp}} f(\mathbf{r}_s) | \partial_x \mathbf{r}_s \times \partial_y \mathbf{r}_s |2 d\mathbf{x}^{\perp} = \int{\Omega^{\perp}} f(\mathbf{r}_s) | (-\nabla^{\perp} h, 1) |_2 d\mathbf{x}^{\perp}, \tag{2.8}$$

where ∥ · ∥² denotes the Euclidean norm. This follows from the fact that ∂xr^s = (1,0,∂xh), ∂yr^s = (0,1,∂yh), and thus ∂xr^s × ∂yr^s = (−∂xh,−∂yh,1) = (−∇ ^⊥h,1). Using that relation, the outward unit normal vector n to Γ^s is:

$$\boldsymbol{n}(\boldsymbol{x}) = \frac{\partial{\boldsymbol{x}} \boldsymbol{r}{\boldsymbol{s}} \times \partial{\boldsymbol{y}} \boldsymbol{r}{\boldsymbol{s}}}{|\partial{\boldsymbol{x}} \boldsymbol{r}{\boldsymbol{s}} \times \partial{\boldsymbol{y}} \boldsymbol{r}{\boldsymbol{s}}|{2}} = \frac{(-\nabla^{\perp} h, 1)}{|(-\nabla^{\perp} h, 1)|{2}}, \quad \boldsymbol{x} \in \Gamma_{\boldsymbol{s}}.$$ (2.9)

Finally, the surface velocity u ^⊥ : Ω^⊥ → R *d*−1 is denoted as u ^⊥ = ux x ^⊥,h ,u^y x ^⊥,h when d = 3 and u ^⊥ = ux x ^⊥,h when d = 2. Note that the relations R Γs f(x)ds = R ^Ω⊥ f(rs)∥(−∇ ^⊥h,1)∥² dx^⊥ and n = (−∇ ^⊥h,1) ∥(−∇⊥h,1)∥² derived above also applies for d = 2.

3. Continuous problem analysis

In this section, we analyze the stability of the problem as stated in (2.4), posed over the infinitedimensional spaces V (Ω), Q(Ω), and Z(Ω⊥). Our goal is to derive the total energy balance, which couples the energy balance of the momentum equations to the free-surface height energy. The framework developed here is also used in later sections to analyze stability and conservation in the discrete setting, where explicit, semi-implicit, and implicit Euler methods are used for time discretization, and finite elements are used in space.

3.1. Relation between the Stokes problem on Ω and the free-surface equation on Ω^⊥

We establish two relations used for the stability analysis of the coupled Stokes/free-surface problem in the sections that follow. The first relates one term in the free-surface equation on Ω^⊥ to an equivalent formulation on Γ^s . The second relates the term formulated on Γ^s , to the right-hand side of the Stokes equation.

Lemma 1 (Free-surface equation to surface integral) Let h : Ω^⊥ × T → R be the surface height function associated with Γ^s (2.7). Then the following relation holds:

$$(-\boldsymbol{u}^{\perp} \cdot \nabla^{\perp} h + u{z}, w){\Omega^{\perp}} = \int{\Gamma{s}(t)} (\boldsymbol{u} \cdot \boldsymbol{n}) w \, ds, \quad \forall w \in Z \left(\Omega^{\perp}\right), \tag{3.1}$$

and when w = h, the relation becomes (−u ^⊥ ·∇ ^⊥h+u^z ,h)Ω^⊥ = R Γs(t) (u · n)zds with (x, y,z) ∈ Ω(t).

Proof We rewrite the integration domain from Ω^⊥ to Γs(t) of the term (−u ^⊥ · ∇ ^⊥h + u^z ,w)Ω^⊥ . We have (−u ^⊥ · ∇ ^⊥h + u^z ,w)Ω^⊥ = R ^Ω⊥ (−u ^⊥ · ∇ ^⊥h + uz)wdx^⊥ = R ^Ω⊥ u · (−∇ ^⊥h,1)wdx^⊥ = R ^Ω⊥ u ⊥ · (−∇ ^⊥h,1)w ∥(−∇ ^⊥h,1)∥² ∥(−∇⊥h,1)∥² *dx*⊥. It follows that: (i) the vector (−^{∇ ⊥}h,1) ∥(−∇⊥h,1)∥² corresponds to the unit normal n as defined in (2.9), and (ii) the remaining factor ∥(−∇ ^⊥h,1)∥² corresponds to the arclength ds of the surface Γs(t), as defined in (2.8). Using that, we rewrite the original integral as:

$$(-\mathbf{u}^{\perp} \cdot \nabla^{\perp} h + u{z}, w){\Omega^{\perp}} = \int{\Omega^{\perp}} \mathbf{u} \cdot \frac{(-\nabla^{\perp} h, 1)}{|(-\nabla^{\perp} h, 1)|{2}} w |(-\nabla^{\perp} h, 1)|{2} d\mathbf{x}^{\perp} = \int{\Gamma_{s}(t)} (\mathbf{u} \cdot \mathbf{n}) w ds, \qquad (3.2)$$

which proves the first statement (3.1) of the lemma. In (3.2) we set w = h. Since Γs(t) is in (2.7) parameterized as z = h(x ^⊥) and since we are integrating over Γ^s , we substitute h = z, which proves the second statement of this proposition. □

Lemma 2 (Stokes problem to surface integral) Let u ∈ V (Ω) solve the Stokes problem in (2.4). Then, the relation between the volume integral −ρg(*z*ˆ,u)Ω(t) and the surface integral over Γs(t) is:

$$-\rho g(\hat{\boldsymbol{z}}, \boldsymbol{u}){\Omega(t)} = -\rho g \int{\Gamma_s(t)} (\boldsymbol{u} \cdot \boldsymbol{n}) z ds, \quad \boldsymbol{x} \in \Omega(t).$$ (3.3)

Proof We write the gravity term as −ρ*gz*ˆ = −ρ*g*∇z. It follows that:

$$-\rho g(\hat{\boldsymbol{z}}, \boldsymbol{u}){\Omega(t)} = -\rho g(\nabla z, \boldsymbol{u}){\Omega(t)} = \rho g(\nabla \cdot \boldsymbol{u}, z){\Omega(t)} - \rho g \int{\partial \Omega(t)} (\boldsymbol{u} \cdot \boldsymbol{n}) z ds = -\rho g \int{\Gamma{s}(t)} (\boldsymbol{u} \cdot \boldsymbol{n}) z ds.$$

In the first step, we integrated (∇z,u)Ω(t) by parts. In the second step, we used the incompressibility condition (∇·u,z)Ω(t) = 0 from (2.4), noting thatz ∈ L 2 (Ω). Finally, the boundary integral was restricted to Γ^s using the boundary conditions. □

3.2. Stability analysis

In this section, we show the stability of the coupled Stokes/free-surface problem (2.4) in the continuous setting, through an energy balance and a corresponding total energy estimate, which is the key to understanding the stability also in discretized settings. In the later sections, we use this total energy estimate to argue about the stability of the numerical schemes.

Theorem 3 (Continuous stability) Let the surface height h and the velocity u solve the coupled system (2.4) on domain Ω = Ω(t), with source term a = a x ^⊥,t , viscosity µ = µ(u), fluid density ρ and the gravitational acceleration g. Then, the following total energy balance holds:

$$\frac{1}{2}\partial{t}|h|{L^{2}(\Omega^{\perp})}^{2} = -\frac{2}{\rho g}|\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^{2}(\Omega(t))}^{2} + (a,h){\Omega^{\perp}}. \tag{3.4}$$

Moreover, the following strong stability estimate in the total energy norm holds:

$$||h(t)||{L^{2}(\Omega^{\perp})}^{2} + \frac{4}{\rho g} \int{0}^{t} ||\sqrt{\mu} \mathbf{D} \mathbf{u}||{L^{2}(\Omega)}^{2} d\tau \leq ||h{0}||{L^{2}(\Omega^{\perp})}^{2} + 2 \int{0}^{t} ||a||{L^{2}(\Omega^{\perp})} \int{0}^{t} ||a||{L^{2}(\Omega^{\perp})} d\tau d\tau + 2 \int{0}^{t} ||a||{L^{2}(\Omega^{\perp})} ||h{0}||_{L^{2}(\Omega^{\perp})} d\tau \qquad \forall t \in T.$$

$$(3.5)$$

Proof

(i) Energy balance. Setting v = u and q = π in (2.4), yields 2∥ √ µ *Du*∥ 2 L ²(Ω(t)) ^{= −}ρg(*z*ˆ,u)^{L 2}(Ω(t)). Applying Lemma 2 to the right-hand side gives 2∥ µ *Du*∥ 2 L ²(Ω(t)) ^{= −}ρ^g R Γs(t) (u · n)zds. We now: (i) apply Lemma 1 to rewrite the surface integral over Γs(t) as an integral over Ω⊥, and (ii) use free-surface equation in (2.4) with w = h to obtain

$$2|\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^{2}(\Omega(t))}^{2} = -\rho g \left(-\mathbf{u}^{\perp} \cdot \nabla^{\perp} h + u{z}, h\right){\Omega^{\perp}} = -\rho g (\partial{t} h, h){\Omega^{\perp}} + \rho g(a, h){\Omega^{\perp}}$$

$$= -\frac{\rho g}{2} \partial{t} |h|{L^{2}(\Omega^{\perp})}^{2} + \rho g(a, h)_{\Omega^{\perp}}.$$ (3.6)

Next, we divide both sides of (3.6) by ρg. Rearranging the terms then proves the first statement (3.4) of this proposition.

(ii) Stability estimate. Using (3.4) we have:

$$\frac{1}{2} \partialt |h|{L^2(\Omega^{\perp})}^2 = |h|_{L^2(\Omega^{\perp})} \, \partialt |h|{L^2(\Omega^{\perp})} \le -\frac{2}{\rho{\mathcal{B}}} |\sqrt{\mu} \, \mathbf{Du}|{L^2(\Omega(t))}^2 + |a|{L^2(\Omega^{\perp})} |h|{L^2(\Omega^{\perp})}. \tag{3.7}$$

In the above relation we use − 2 ρg √ µ *Du*∥ 2 L ²(Ω(t)) ^≤ 0 and then divide the resulting inequality by ∥*h∥^L* ²(Ω⊥) , which gives: ∂*t*∥h(t)∥^{L 2}(Ω(t)) ≤ ∥*a∥^L* ²(Ω⊥) . We integrate the latter relation in time to arrive at: ∥h(t)∥^{L 2}(Ω(t)) ≤ ∥*h0∥^L* ²(Ω(t)) + R t 0 ∥a(τ)∥^{L 2}(Ω⊥) *d*τ. Inserting this on the right-hand-side of (3.7), and then time-integrating and multiplying by 2 on both sides of the inequality, finalizes the proof of the second statement (3.5) of this proposition. □

Remark 1 We note that the coupled problem (2.4) depends on the parameter p through the viscosity function (2.2). Accordingly, the appropriate norm for the strain tensor Du is the one associated with the L^p space. The scaled L^p -norm of the strain tensor is related to its viscosity-weighted L² -norm as follows:

$$\mu0 | \mathbf{D} \mathbf{u} |{L^p(\Omega)}^p = \int_{\Omega} \mu_0 | \mathbf{D} \mathbf{u} |^p d\mathbf{x}$$

$$= \int_{\Omega} \mu0 | \mathbf{D} \mathbf{u} |^{p-2} | \mathbf{D} \mathbf{u} |^2 d\mathbf{x} = \int{\Omega} \mu | \mathbf{D} \mathbf{u} |^2 d\mathbf{x} = \int{\Omega} |\sqrt{\mu} \mathbf{D} \mathbf{u} |^2 d\mathbf{x} = |\sqrt{\mu} \mathbf{D} \mathbf{u} |{L^2(\Omega)}^2,$$ (3.8)

where, in the third step, we used the viscosity definition (2.2). Therefore, the estimates in terms of √ µ *Du*∥ L ²(Ω) remain valid also for p ̸= 2.

3.3. Volume conservation

In the continuous setting, the domain volume |Ω| is conserved over time. For completeness, we briefly revisit the volume conservation analysis. Using the free-surface equation (2.4) with w = 1 we obtain:

$$\partial_t |\Omega(t)| = \partialt \int{\Omega^{\perp}} (h - b) \, d\mathbf{x}^{\perp} = \int_{\Omega^{\perp}} \partial_t h \, d\mathbf{x}^{\perp} = (\partial_t h, 1) = (-\mathbf{u}^{\perp} \cdot \nabla^{\perp} h + uz + a, 1){\Omega^{\perp}}. \tag{3.9}$$

Next, applying Lemma 1 to relate the integration over Ω^⊥ to integration over the surface Γs(t), and using integration by parts together with the incompressibility condition from (2.4) with q = 1 ∈ L (Ω) we obtain:

$$\partial{t}|\Omega(t)| = (-\boldsymbol{u}^{\perp} \cdot \nabla^{\perp} h + u{z} + a, 1){\Omega^{\perp}} = \int{\partial \Omega} \boldsymbol{u} \cdot \boldsymbol{n} \, ds + \int{\Omega^{\perp}} a \, d\boldsymbol{x}^{\perp} = \int{\Omega(t)} \nabla \cdot \boldsymbol{u} \, d\boldsymbol{x} + \int_{\Omega^{\perp}} a \, d\boldsymbol{x}^{\perp}$$

$$= \int_{\Omega^{\perp}} a \, d\boldsymbol{x}^{\perp}.$$ (3.10)

This shows that the domain volume only changes due to the source term a.

4. Discretized Stokes/free-surface problems and their analysis

In this section, we analyze stability and volume conservation when the finite element method is used in space and explicit Euler, semi-implicit Euler, and implicit Euler are used to discretize time. We aim to derive discrete analogues of the continuous total energy estimate (3.5) and volume conservation balance (3.10).

Deriving a discrete counterpart to the continuous total energy estimate (3.5) is a crucial step for having the numerical stability in time. However, this is not sufficient to suppress high-frequency oscillations when the high-order derivatives either entail discontinuities or are under-resolved on a computational mesh, as the total energy of the couple Stokes/free-surface problem does not control the higher-order derivatives on its own. To suppress the high-order oscillations, we use the edge stabilization approach [8, 10]. This is a third-order accurate term, preferred to use on P¹ spaces over the hyperviscosity-based approach that requires higher-order polynomial degrees. To all the considered numerical schemes in this paper we add the edge stabilization term defined as:

$$\sum_{e \in \mathcal{E}_h^{\perp}} \int_e \gamma_e [![ \nabla^{\perp} h ]!] \cdot [![ \nabla^{\perp} w ]!] ds, \quad \gammae = \frac{1}{2} h{K^{\perp}}^2 |\mathbf{u}|, \tag{4.1}$$

where the jump of a vector field over an interior mesh edge e = ∂K⁺ ∩∂K⁻ ∈ E is defined as [[w]] := w⁺ · n⁺ −w⁻ · n⁻ with w^± = w|K^± and n^± denoting the outward-pointing unit normal vector on ∂*K*±. The mesh size hK^⊥ is defined as the cell diameter. We note that this stabilization approach does not stabilize the scheme in time, but is only provided to suppress the high-frequency oscillations that, in some cases, appear during numerical simulations. Some alternative linear high-order stabilization approaches (combined with nonlinear artificial viscosity) are given in [10, 17, 19, 21, 25, 27, 29, 30].

4.1. Implicit Euler discretized problem analysis

The implicit Euler discretization of (2.4), including edge stabilization (4.1), is defined as follows. Given h n and a n+1 , find u ⁿ+¹ ∈ V ^h(Ω n+1 h ), π ⁿ+¹ ∈ Qh(Ω n+1 h ), and h ⁿ+¹ ∈ Zh(Ω^⊥ h ) such that:

$$\begin{cases} \left(\frac{h^{n+1}-h^n}{\Delta t},w\right)_{\Omega_h^\perp} + (\boldsymbol{u}^{\perp,n+1}\cdot\nabla^\perp h^{n+1}-uz^{n+1},w){\Omegah^\perp} + \sum{e\in\mathcal{E}_h^\perp} \int_e \gammae [![\nabla^\perp h^{n+1}]!]\cdot [![\nabla^\perp w]!] ds = (a^{n+1},w){\Omegah^\perp}, \ 2(\sqrt{\mu}\,\boldsymbol{D}\boldsymbol{u}^{n+1},\sqrt{\mu}\,\boldsymbol{D}\boldsymbol{v}){\Omegah^{n+1}} - (\pi^{n+1},\nabla\cdot\boldsymbol{v}){\Omegah^{n+1}} = -\rho\,g(\boldsymbol{\hat{z}},\boldsymbol{v}){\Omegah^{n+1}}, \ (\nabla\cdot\boldsymbol{u}^{n+1},q){\Omega_h^{n+1}} = 0, \end{cases}$$

for all v ∈V ^h(Ω n+1 h ), q ∈ Qh(Ω n+1 h ), and w ∈ Zh(Ω^⊥ h ). In the two propositions below, we show that this scheme:

• 1. is unconditionally stable with respect to ∆t in the total energy norm, mimicing the continuous counterpart (3.5),
• 2. enables the domain volume conservation.

However, the scheme requires nonlinear iterations accross the Stokes/free-surface problem due to the dependence of Ωn+¹ in the Stokes problem to h n+1 in the free-surface equation. In that sense, the scheme is less computationally efficient.

Proposition 4 (Implicit Euler stability) Let (u n+1 ,π n+1 ,h n+1 ) be a solution to (4.2) with source term a ⁿ+¹ = a(x ^⊥,t n+1 ). Then the following stability estimate in the total energy norm holds:

$$|h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} + \frac{4}{\rho g} \Delta t |\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n+1})}^{2} \leq |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2\Delta t |a^{n+1}|{L^{2}(\Omega{h}^{\perp})} |h^{n}|{L^{2}(\Omega{h}^{\perp})} + 2(\Delta t)^{2} |a^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2}.$$

$$+ 2(\Delta t)^{2} |a^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2}.$$

$$(4.3)$$

Proof

(i) Energy balance. We set w = 2h n+1 in the free-surface equation in (4.2) and then apply the identity 2a ² −2ab = a ² −b ² + (*a*−b) ² with a = h n+1 and b = h n , to obtain:

$$\begin{split} |h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} &= |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} - |h^{n+1} - h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2\Delta t (-\boldsymbol{u}^{\perp,n+1} \cdot \nabla h^{n+1} + u{3}^{\perp} + a^{n+1}, h^{n+1}){\Omega{h}^{\perp}} \ &- 2\Delta t \sum{e \in \mathcal{E}{h}^{\perp}} \int{e} \gamma_{e} [![\nabla^{\perp} h^{n+1}]!] \cdot [![\nabla^{\perp} h^{n+1}]!] ds. \end{split} \tag{4.4}$$

We now focus on the term (−u ⊥,n+1 ·∇h ⁿ+¹ +u ⊥,n+1 3 ,h n+1 )Ω^⊥ h . Since all variables in this bilinear form are evaluated at the same time t = t n+1 , the analysis framework for the continuous case as outlined in Lemma 1, Lemma 2, and Theorem 3 applies. In particular, using (3.6) from the proof of Theorem 3 we obtain: 2∥ √ µ *Du*∥ 2 L ²(Ω n+1 h = −ρg(−u ⊥,n+1 ·∇h ⁿ+¹ +u ⊥,n+1 3 ,h n+1 )Ω^⊥ h . Dividing the equation above by ρg, substituting this result into (4.4), and rearranging the terms, gives:

$$\begin{split} |h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} + \frac{4}{\rho g} \Delta t |\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n+1})}^{2} &= |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2\Delta t \, (a^{n+1}, h^{n+1}){\Omega{h}^{\perp}} \ &- 2\Delta t \sum{e \in \mathcal{E}{h}^{\perp}} \int{\mathcal{E}} \gamma{e} [![\nabla^{\perp} h^{n+1}]!] \cdot [![\nabla^{\perp} h^{n+1}]!] ds - |h^{n+1} - h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2}. \end{split}$$

$$(4.5)$$

(ii) An intermediate L 2 -estimate of h. We set w = h n+1 , v = u n+1 , and q = π n+1 in (4.2), and repeat similar steps as in the first part of this proof, to obtain:

$$|h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} = (h^{n}, h^{n+1}){\Omega{h}^{\perp}} - \frac{2}{\rho{g}} \Delta t |\sqrt{\mu} \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n+1})}^{2} + \Delta t (a^{n+1}, h^{n+1}){\Omega{h}^{\perp}} - \Delta t \sum{e \in \mathcal{E}{h}^{\perp}} \int{e} \gamma_{e} [![\nabla^{\perp} h^{n+1}]!] \cdot [![\nabla^{\perp} h^{n+1}]!] ds.$$

$$(4.6)$$

Neglecting the negative terms and applying the Cauchy-Schwarz inequality, we estimate:

$$||h^{n+1}||{L^{2}(\Omega{h}^{\perp})}^{2} \leq (h^{n}, h^{n+1}){\Omega{h}^{\perp}} + \Delta t (a^{n+1}, h^{n+1}){\Omega{h}^{\perp}}$$

$$= (h^{n} + \Delta t a^{n+1}, h^{n+1}){\Omega{h}^{\perp}}$$

$$\leq ||h^{n} + \Delta t a^{n+1}||{L^{2}(\Omega{h}^{\perp})} ||h^{n+1}||{L^{2}(\Omega{h}^{\perp})}.$$ (4.7)

Dividing both sides by ∥h ⁿ+1∥^{L 2}(Ω^⊥ h , and applying the triangle inequality then gives:

$$|h^{n+1}|_{L^2(\Omegah^{\perp})} \le |h^n|{L^2(\Omegah^{\perp})} + \Delta t |a^{n+1}|{L^2(\Omega_h^{\perp})}. \tag{4.8}$$

(iii) Total energy estimate. We start at (4.5) and use −2∆t ∑*e*∈^E h R e γe[[∇ ⊥h n+1 ]] · [[∇ ⊥h n+1 ]]ds ≤ 0 as well as 2∆t(a n+1 ,h n+1 )Ω^⊥ h ≤ 2∆*t*∥a ⁿ+1∥^{L 2}(Ω^⊥ h ∥h ⁿ+1∥^{L 2}(Ω^⊥ h . Then we have:

$$\begin{split} |h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} + \frac{4}{\rho{g}} \Delta t |\sqrt{\mu} \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n+1})}^{2} &\leq |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2\Delta t |a^{n+1}|{L^{2}(\Omega{h}^{\perp})} |h^{n+1}|{L^{2}(\Omega{h}^{\perp})} \ &\leq |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2\Delta t |a^{n+1}|{L^{2}(\Omega{h}^{\perp})} |h^{n}|{L^{2}(\Omega{h}^{\perp})} \ &+ 2(\Delta t)^{2} |a^{n+1}|{L^{2}(\Omega_{h}^{\perp})}^{2}, \end{split} \tag{4.9}$$

where we used (4.8) in the final step. □

Proposition 5 (Implicit Euler conservation) The solution to (4.2) with source term an+¹ = a(x ^⊥,t n+1 conserves domain volume: |Ω n+1 h | = |Ωⁿ h |+∆t R Ω^⊥ a ⁿ+¹ *dx*⊥.

Proof We set w = 1 in the free-surface equation part of (4.2) to obtain:

$$(h^{n+1},1){\Omega{h}^{\perp}} = (h^{n},1){\Omega{h}^{\perp}} + \Delta t \left( -\boldsymbol{u}^{\perp,n+1} \cdot \nabla h^{n+1} + u{3}^{\perp,n+1} + a^{n+1}, 1 \right){\Omega_{h}^{\perp}}. \tag{4.10}$$

Since all terms in the last integrand are evaluated at time t = t n+1 , the analysis in Section 3.3 applies. Following the derivation of the continuous volume conservation in (3.10), we have that:

$$(-\boldsymbol{u}^{\perp,n+1}\cdot\nabla h^{n+1}+u3^{\perp,n+1},1){\Omegah^\perp}=\int{\partial\Omegah^{n+1}}(\boldsymbol{u}\cdot\boldsymbol{n})\,ds=\int{\Omega_h^{n+1}}\nabla\cdot\boldsymbol{u}^{n+1}\,ds=0, \tag{4.11}$$

where the final equality follows from the incompressibility condition (4.2) with q = 1 ∈ Qh(Ω n+1 h ). Combining (4.10) and (4.11), and using the definition of domain volume via the height function, we conclude |Ω n+1 h | − |Ωⁿ h | = (h ⁿ+¹ − b,1) − (h ⁿ − b,1) = (h n+1 ,1) − (h n ,1) = ∆t R Ω^⊥ h a ⁿ+¹ *dx*⊥, which completes the proof. □

4.2. Explicit Euler discretized problem analysis

We now consider the explicit Euler discretization of (2.4) with edge stabilization. Given h n and a n , find u ⁿ ∈ V ^h(Ωⁿ h ), π ⁿ ∈ Qh(Ωⁿ h ), and h ⁿ+¹ ∈ Zh(Ω^⊥ h ) such that:

$$\begin{cases} \left(\frac{h^{n+1}-h^n}{\Delta t},w\right)_{\Omega_h^{\perp}} + \left(\boldsymbol{u}^{\perp,n}\cdot\nabla^{\perp}h^n - uz^n,w\right){\Omegah^{\perp}} + \sum{e\in\mathcal{E}h^{\perp}}\int{e}\gammae^n[![\nabla^{\perp}h^{n+1}]!]\cdot[![\nabla^{\perp}w]!]ds = (a^n,w){\Omegah^{\perp}},\ 2\left(\sqrt{\mu}\,\boldsymbol{D}\boldsymbol{u}^n,\sqrt{\mu}\,\boldsymbol{D}\boldsymbol{v}\right){\Omegah^n} - (\pi^n,\nabla\cdot\boldsymbol{v}){\Omegah^n} = -\rho\,g(\boldsymbol{z},\boldsymbol{v}){\Omegah^n},\ (\nabla\cdot\boldsymbol{u}^n,q){\Omega_h^n} = 0, \end{cases}$$

$$(4.12)$$

for all $$\mathbf{v} \in \mathbf{V}_h(\Omega_h^n)$$ , $q \in Q_h(\Omega_h^n)$ , and $w \in Z_h(\Omega_h^\perp)$ .

This is a numerical scheme central to this paper. The scheme does not require nonlinear iterations across the Stokes/free-surface problems and is, in that sense, computationally efficient. In the two propositions below, we show that this scheme:

• 1. is unstable in the total energy norm for any ∆t, with spurious terms (∆t) R Γn ω (u · n) ² ds and 2(∆t) 2 R Γn (u · n)a ⁿ ds that prevent mimicing the continuous total energy estimate (3.5),
• 2. enables the domain volume conservation.

Proposition 6 (Explicit Euler stability estimate) Let (u n ,π n ,h n+1 ) be a solution to (4.12) with source term aⁿ = a(x ^⊥,t n ). Then the following stability estimate in the total energy norm holds:

$$|h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} + \frac{4}{\rho g} \Delta t |\sqrt{\mu} \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n})}^{2} \leq |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2 \Delta t (a^{n}, h^{n}){\Omega{h}^{\perp}} + (\Delta t)^{2} |a^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + (\Delta t)^{2} \int{\Gamma{s}^{n}} \omega (\mathbf{u} \cdot \mathbf{n})^{2} ds + 2 (\Delta t)^{2} \int{\Gamma{s}^{n}} (\mathbf{u} \cdot \mathbf{n}) a^{n} ds.$$

$$(4.13)$$

Proof

(i) Decomposition of the free-surface equation into explicit and implicit parts. In the first equation

of (4.12), we use h ⁿ+¹ = h n+1 ^E +h n+1 I to then decompose the equation into an explicit part and an implicit part:

$$\left(\frac{hE^{n+1} - h^n}{\Delta t}, w\right){\Omega_h^{\perp}} + \left(\boldsymbol{u}^{\perp, n} \cdot \nabla^{\perp} h^n - uz^n, w\right){\Omegah^{\perp}} = (a^n, w){\Omega_h^{\perp}},\tag{4.14}$$

$$\left(\frac{hI^{n+1}}{\Delta t}, w\right){\Omegah^{\perp}} = -\sum{e \in \mathscr{E}_h^{\perp}} \int_e \gamma_e^n [![\nabla^{\perp} h^{n+1}]!] \cdot [![\nabla^{\perp} w]!] ds. \tag{4.15}$$

In (4.14), we now set w = 2h n and use the identity 2*ab*−2b ² = a ² −b ² −(*a*−b) ² with a = h n+1 E , b = h n , to obtain:

$$|hE^{n+1}|{L^2(\Omegah^{\perp})}^2 = |h^n|{L^2(\Omega_h^{\perp})}^2 + 2\Delta t(-\boldsymbol{u}^{\perp,n} \cdot \nabla^{\perp} h^n + uz^n, h^n){\Omegah^{\perp}} + 2\Delta t(a^n, h^n){\Omega_h^{\perp}} + |hE^{n+1} - h^n|{L^2(\Omega_h^{\perp})}^2. \tag{4.16}$$

(ii) Energy balance. We use Lemma 1 to write term (−u ⊥,n · ∇ ⊥h ⁿ + u n ,h n )Ω^⊥ h as a surface integral, followed by using Lemma 2 to write the resulting surface integral to the scaled Stokes gravity term. This gives: 2∆t(−u ⊥,n · ∇ ⊥h n ,h n )Ω^⊥ h = 2∆t R Γn s (u · n)zds = 2∆t(*z*ˆ,u)Ωⁿ h . Next, we use the Stokes problem (4.12) formulation with v = u n and q = π n to relate the scaled gravity term to the strain tensor norm as:

$$2\Delta t(-\boldsymbol{u}^{\perp,n}\cdot\nabla^{\perp}h^{n},h^{n}){\Omega{h}^{\perp}}=2\Delta t\left(\hat{\boldsymbol{z}},\boldsymbol{u}\right){\Omega{h}^{n}}=-\frac{4}{\rho g}\Delta t\left|\sqrt{\mu}\,\boldsymbol{D}\boldsymbol{u}\right|{L^{2}(\Omega{h}^{n})}^{2}.$$ (4.17)

Inserting this relation into (4.16) and then rearranging, gives the energy balance:

$$|hE^{n+1}|{L^2(\Omega_h^{\perp})}^2 + \frac{4}{\rhog} \Delta t |\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^2(\Omegah^n)}^2 = |h^n|{L^2(\Omega_h^{\perp})}^2 + |hE^{n+1} - h^n|{L^2(\Omegah^{\perp})}^2 + 2\Delta t (a^n, h^n){\Omega_h^{\perp}}. \tag{4.18}$$

(iii) Spurious term estimate. We estimate the spurious term ∥h n+1 ^E − *h n*∥ 2 L ²(Ω^⊥ h in (4.18) by setting w = h n+1 ^E − h n and then estimating: ∥h n+1 ^E − *h n*∥ 2 L ²(Ω^⊥ = ∆t(−u ⊥,n · ∇ ⊥h ⁿ + u n ^z + a n ,h n+1 ^E − h n )Ω^⊥ h ≤ ∆*t*∥ − u ⊥,n · ∇ ⊥h ⁿ + u n ^z + *a n*∥L ²(Ω^⊥ ∥h n+1 ^E − *h n*∥L ²(Ω^⊥ . Dividing both sides by ∥h n+1 ^E − *h n*∥L ²(Ω^⊥ and squaring the resulting inequality gives:

$$|h{E}^{n+1} - h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} \leq (\Delta t)^{2} | - \mathbf{u}^{\perp,n} \cdot \nabla^{\perp} h^{n} + u{z}^{n} + a^{n}|{L^{2}(\Omega{h}^{\perp})}^{2}$$

$$= (\Delta t)^{2} \underbrace{| - \mathbf{u}^{\perp,n} \cdot \nabla^{\perp} h^{n} + u{z}^{n}|{L^{2}(\Omega{h}^{\perp})}^{2}}{\mathbf{Term A}}$$

$$+ 2(\Delta t)^{2} \underbrace{(a^{n}, -\mathbf{u}^{\perp,n} \cdot \nabla^{\perp} h^{n} + u{z}^{n})}{\mathbf{Term B}} + (\Delta t)^{2} |a^{n}|{L^{2}(\Omega{h}^{\perp})}^{2}.$$

$$(4.19)$$

We now rewrite Term A in (4.19) by defining an integral scaling: ω = ∥(−1,∇ ⊥h n )∥² = 1 nz > 0, where n^z denotes the vertical component of an outward unit normal vector n on Γ n . Inserting this relation into (4.19), we have:

$$|-\mathbf{u}^{\perp,n}\cdot\nabla^{\perp}h^{n}+u{z}^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} = \int{\Omega{h}^{\perp}}\frac{\omega}{\omega}\left(-\mathbf{u}^{\perp,n}\cdot\nabla^{\perp}h^{n}+u{z}^{n}\right)\left(-\mathbf{u}^{\perp,n}\cdot\nabla^{\perp}h^{n}+u_{z}^{n}\right)d\mathbf{x}^{\perp}$$

$$= \int{\Gamma{n}^{n}}(\mathbf{u}\cdot\mathbf{n})\left(-\mathbf{u}^{\perp,n}\cdot\nabla^{\perp}h^{n}+u{z}^{n}\right)ds = \int{\Gamma_{n}^{n}}\omega\left(\mathbf{u}\cdot\mathbf{n}\right)^{2}ds,$$

$$(4.20)$$

where we in the second step used: ¹ ω (−u ⊥,n · ∇ ⊥h ⁿ + u n ) = u n · n and ω dx^⊥ = ds (as in the proof of Lemma 1). In the final step, we used −u ⊥,n ·∇ ⊥h ⁿ +u n ^z = ω (u · n). We proceed by rewriting Term B in (4.19) as follows:

$$(a, -\boldsymbol{u}^{\perp,n} \cdot \nabla^{\perp} h^n + uz^n) = \int{\Omega_h^{\perp}} \frac{\omega}{\omega} a \left( -\boldsymbol{u}^{\perp,n} \cdot \nabla^{\perp} h^n + uz^n \right) d\boldsymbol{x}^{\perp} = \int{\Gamma_s^n} (\boldsymbol{u} \cdot \boldsymbol{n}) \, a^n \, ds. \tag{4.21}$$

Alltogether, inserting (4.20) and (4.21) into (4.19), gives the final bound of the spurious term:

$$|hE^{n+1} - h^n|{L^2(\Omegah^{\perp})}^2 \le (\Delta t)^2 \int{\Gammah^n} \omega (\mathbf{u} \cdot \mathbf{n})^2 ds + 2 (\Delta t)^2 \int{\Gammah^n} (\mathbf{u} \cdot \mathbf{n}) a^n ds + (\Delta t)^2 |a^n|{L^2(\Omega_h^{\perp})}^2. \tag{4.22}$$

(iv) Finalization. We now make steps towards estimating ∥h ⁿ+1∥^{L 2}(Ω^⊥ h by writing ∥h ⁿ+1∥^{L 2}(Ω^⊥ h ) = ∥h n+1 ^E + h n+1 I L ²(Ω^⊥ = ∥h n+1 E 2 L ²(Ω^⊥ + ∥h n+1 I L ²(Ω^⊥ + 2(h n+1 E ,h n+1 I ). We continue by testing (4.15) against w = h n+1 ^E +h n+1 I to first expand the term 2(h n+1 E ,h n+1 I )Ω^⊥ h , and then bound out all the resulting negative terms, as follows:

$$|h^{n+1}|_{L^2(\Omega_h^\perp)} = |hE^{n+1}|{L^2(\Omega_h^\perp)}^2 - |hI^{n+1}|{L^2(\Omegah^\perp)}^2 - 2\Delta t \sum{e \in \mathscr{E}h^\perp} \int{e} \gamma_e^n \llbracket \nabla^\perp h^{n+1} \rrbracket \cdot \llbracket \nabla^\perp h^{n+1} \rrbracket \, \mathrm{d}s \leq |hE^{n+1}|{L^2(\Omega_h^\perp)}^2. \tag{4.23}$$

Finally, inserting (4.18) into (4.23), followed by inserting (4.22) into the resulting estimate, proves this proposition. □

Proposition 7 (Explicit Euler domain volume conservation) Solution to the explicit Euler scheme (4.12) with source term aⁿ = a(x ^⊥,t n ), conserves domain volume: |Ω n+1 h | −|Ωⁿ h | = ∆t R Ω^⊥ a ⁿ dx.

Proof Applying all steps of the proof of Proposition 5 to (4.12), proves this proposition. □

4.2.1. Weakly stable time step size for the unstabilized explicit Euler discretized problem The stability estimate (4.13) implies that the corresponding scheme (4.12) is unstable in the total energy norm for any ∆t. In the steps that follow, we determine ∆t that implies stability in a weaker sense, that is, stability without having a full control over the total energy norm. This is disadvantageous as it implies a weaker coupling between the Stokes/free-surface equations and gives rise to increased energy errors.

We first write ⁴ ρg ∆t ∥ √ µ *Du*∥ L ²(Ωⁿ = (1−ε) ρg ∆t ∥ √ µ *Du*∥ L ²(Ωⁿ +ε ρg ∆t ∥ √ µ *Du*∥ L ²(Ωⁿ , where 0 < ε < 1. We insert that relation into (4.13) and rearrange the terms to arrive at a weaker stability estimate:

$$\begin{split} |h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} + (1-\varepsilon) \frac{4}{\rho g} \Delta t \, |\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n})}^{2} &\leq |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2 \Delta t (a^{n}, h^{n}){\Omega{h}^{\perp}} + (\Delta t)^{2} |a^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} \ &+ (\Delta t)^{2} \int{\Gamma{s}^{n}} \omega (\mathbf{u} \cdot \mathbf{n})^{2} \, ds + 2 (\Delta t)^{2} \int{\Gamma{s}^{n}} (\mathbf{u} \cdot \mathbf{n}) \, a^{n} \, ds \ &- \varepsilon \, \frac{4}{\rho g} \, \Delta t \, |\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n})}^{2}. \end{split} \tag{4.24}$$

For weaker stability we require that:

$$(\Delta t)^2 \int_{\Gamma3^n} \omega \left( \boldsymbol{u} \cdot \boldsymbol{n} \right)^2 ds + 2 \left( \Delta t \right)^2 \int{\Gamma3^n} (\boldsymbol{u} \cdot \boldsymbol{n}) \, d^n \, ds - \varepsilon \, \frac{4}{\rho g} \, \Delta t \, | \sqrt{\mu} \, \boldsymbol{D} \boldsymbol{u} |{L^2(\Omega_h^n)}^2 \leq 0.$$

Solving that inequality for ∆t implies a time step restriction that gives stability in a weaker sense:

$$\Delta t^n \le \varepsilon \frac{\frac{4}{\rhoS} |\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^2(\Omegah^n)}^2}{\int{\GammaS^n} \omega \left(\mathbf{u} \cdot \mathbf{n}\right)^2 ds + 2 \int{\Gamma_S^n} (\mathbf{u} \cdot \mathbf{n}) \, a^n \, ds}, \quad 0 < \varepsilon < 1, \quad n = 0, 1, 2, \dots$$

$$(4.25)$$

Using such ∆t, the norm (1 − ε) ρg ∆t ∥ √ µ *Du*∥ 2 L ²(Ωⁿ in (4.24) is now controlled only up to the ε parameter. We use ∆t computed according to (4.25) only for drawing comparisons across numerical schemes in the numerical experiments section. We do not recommend using this approach in practice as it leads to a weak coupling between the Stokes/free-surface equations and increased energy errors. Instead, in Section 4.3, we give a stabilization term that (i) allows a full total energy norm control without ε-dependence, and (ii) enables stability for any choice of ∆t.

4.3. Stabilized explicit Euler discretized coupled problem

Motivated by the stability estimate (4.13) for the explicit Euler method (4.12), we introduce the following stabilized method. Given h n and a n , find h ⁿ+¹ ∈ Zh(Ω^⊥ h ), u ⁿ ∈ V ^h(Ωⁿ h ), π ⁿ ∈ Qh(Ωⁿ h ) such that:

$$\begin{cases} \left(\frac{h^{n+1}-h^n}{\Delta t},w\right)_{\Omega_h^{\perp}}+\left(\mathbf{u}^{\perp,n}\cdot\nabla^{\perp}h^n-uz^n,w\right){\Omegah^{\perp}}+\sum{e\in\mathcal{E}_h^{\perp}}\int_e\gammae[![\nabla^{\perp}h^{n+1}]!]\cdot[![\nabla^{\perp}w]!]ds=(a^n,w){\Omegah^{\perp}},\ 2(\sqrt{\mu}\,\mathbf{D}\mathbf{u}^n,\sqrt{\mu}\,\mathbf{D}\mathbf{v}){\Omegah^n}+\frac{\rho g\Delta t}{2}\int{\Gamma_s^n}!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!$$

for all w ∈ Zh(Ω^⊥ h ), v ∈ V ^h(Ωⁿ h ), and q ∈ Qh(Ωⁿ h ).

This scheme does not require nonlinear iterations accross the Stokes/free-surface problem and is in that sense computationally efficient. The only difference between the unstabilized and stabilized systems (4.12) and (4.26) are the two added stabilization terms ^ρ*g∆^t* R Γn ω (u · n) (v · n)ds, and ρ*g*∆t R Γn a(v ·n)ds. The first term is, up to the scaling, the same as introduced in [26], where it was used for time step size estimation. The scaling that we propose is motivated by our analysis in Proposition 6 without simplifications of the original Stokes/free-surface problem formulation. The added stabilization terms can be seen as of FSSA-type [18], however, compared to the FSSA stabilization terms in (4.29), our stabilization terms are derived such that they precisely cancel the spurious term in the total energy estimate outlined in Proposition 6. Furthermore, these stabilization terms do not break the symmetry of the Stokes problem, which is advantageous when using the iterative linear solvers.

In the two propositions below, we show that the scheme (4.26):

• 1. is unconditionally stable with respect to ∆t in the total energy norm, mimicing the continuous counterpart (3.5),
• 2. enables the domain volume conservation.

Proposition 8 (Stabilized explicit Euler stability estimate) Let (u n ,π n ,h n+1 ) be a solution to the stabilized explicit Euler scheme (4.26) with source term aⁿ = a(x ^⊥,t n ). Then the following stability estimate in the total energy norm holds:

$$|h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} + \frac{4\Delta t}{\rho{g}} |\sqrt{\mu} \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n})}^{2} \leq |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2\Delta t(a^{n}, h^{n}){\Omega{h}^{\perp}} + (\Delta t)^{2} |a^{n}|{L^{2}(\Omega_{h}^{\perp})}^{2}. \tag{4.27}$$

Proof We reuse the proof of Proposition 6. Taking into account the stabilization terms in (4.26) we rewrite (4.17) to: 2∆t(u ⊥,n · ∇ ⊥h ⁿ −u n ,h n )Ω^⊥ = 2∆t(*z*ˆ,u)Ωⁿ h = − 4∆t ρg √ µ *Du*∥ 2 L ²(Ωⁿ h −(∆t) 2 R Γn ω (u · n) ² *ds*−2(∆t 2 ) R Γn (u · n)a ⁿ ds. We now use this relation to reiterate the remaining steps in the proof of Proposition 6 (from (4.17) on) and then obtain the following estimate:

$$\begin{split} &|h^{n+1}|{L^{2}(\Omega{h}^{\perp})}^{2} + \frac{4\Delta t}{\rho g} |\sqrt{\mu} \, \mathbf{D} \mathbf{u}|{L^{2}(\Omega{h}^{n})}^{2} + (\Delta t)^{2} \int{\Gamma{s}^{n}} \omega \, (\mathbf{u} \cdot \mathbf{n})^{2} \, ds + 2 \, (\Delta t)^{2} \int{\Gamma{s}^{n}} (\mathbf{u} \cdot \mathbf{n}) \, a^{n} \, ds \ & \leq |h^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + 2\Delta t (a^{n}, h^{n}){\Omega{h}^{\perp}} + (\Delta t)^{2} |a^{n}|{L^{2}(\Omega{h}^{\perp})}^{2} + (\Delta t)^{2} \int{\Gamma{s}^{n}} \omega \, (\mathbf{u} \cdot \mathbf{n})^{2} \, ds + 2 \, (\Delta t)^{2} \int{\Gamma{s}^{n}} (\mathbf{u} \cdot \mathbf{n}) \, a^{n} \, ds. \end{split}$$

$$(4.28)$$

We subtract (∆t) R Γn ω (u ·n) ² ds+2(∆t) R Γn a n (u ·n)ds on both sides of the inequality to finalize the proof of this proposition. □

Proposition 9 (Stabilized explicit Euler domain volume conservation) Solution to the explicit Euler scheme (4.26) with source term aⁿ = a(x ^⊥,t n ), conserves domain volume: |Ω n+1 h |−|Ωⁿ h | = ∆t R Ω^⊥ a ⁿ dx.

4.4. Semi-implicit Euler discretized coupled problem stabilized by FSSA

The semi-implicit Euler method stabilized using the FSSA approach [18] is formulated as follows. Given h n and a n , find h ⁿ+¹ ∈ Zh(Ω^⊥ h ), u ⁿ ∈ V ^h(Ωⁿ h ), and π ⁿ ∈ Qh(Ωⁿ h ), such that:

$$\begin{cases} \left(\frac{h^{n+1}-h^n}{\Delta t},w\right)_{\Omega_h^{\perp}}+(\boldsymbol{u}^{\perp,n}\cdot\nabla^{\perp}h^{n+1}-uz^n,w){\Omegah^{\perp}}+\sum{e\in\mathcal{E}h^{\perp}}\int{e}\gammae[![\nabla^{\perp}h^{n+1}]!]\cdot[![\nabla^{\perp}w]!]ds=(a^n,w){\Omegah^{\perp}},\ 2(\sqrt{\mu}\boldsymbol{D}\boldsymbol{u}^n,\sqrt{\mu}\boldsymbol{D}\boldsymbol{v}){\Omega_h^n}+\rhog\Delta t\int{\Gamma_s^n}(\boldsymbol{u}\cdot\boldsymbol{n})(\hat{\boldsymbol{z}}\cdot\boldsymbol{v})ds+\rhog\Delta t\int{\Gammas^n}(\boldsymbol{v}\cdot\hat{\boldsymbol{z}}\boldsymbol{n})ads-(\pi^n,\nabla\cdot\boldsymbol{v}){\Omega_h^n}=-\rhog(\hat{\boldsymbol{z}},\boldsymbol{v}){\Omegah^n},\ (\nabla\cdot\boldsymbol{u}^n,q){\Omega_h^n}=0, \end{cases}$$

for all v ∈ V ^h(Ωⁿ h ), q ∈ Qh(Ωⁿ h ), and w ∈ Zh(Ω^⊥ h ). This semi-implicit formulation is a common discretization choice in for instance numerical models of ice sheet dynamics [20, 22, 28]. The formulation can also be written in the explicit Euler setting, but numerical evidence suggests that the semi-implicit formulation allows taking larger stable time step sizes [20, 22]. In both settings, the added non-symmetric term ρ*g*∆t R Γn (u · n) (*z*ˆ · v)ds breaks the symmetry of the Stokes problem, which could be disadvantageous when using the iterative linear solvers.

The stability analysis of (4.29) requires extending the theoretical framework in the proof of Theorem 3 beyond the scope of the present paper. Instead, in Section 5.1 and Section 5.3 we experimentally show that the solution to this scheme:

• 1. is conditionally stable with respect to ∆t in the total energy norm,
• 2. does not enable the domain volume conservation.

5. Numerical experiments

In the sections that follow, we outline numerical experiments verifying our theoretical results on the discrete stability of the Stokes free-surface flow. We use the open-source finite element library FEniCS, version 2019.1 [1, 2, 3] to perform the computations on a laptop with AMD Ryzen 7 PRO 6850U processor and 16 GB RAM.

Boundary conditions. In all considered cases, the boundary conditions are: u = 0 on Γb, u · n = 0 on Γ^l , and σ · n = 0 on Γ^s

Mesh resolution. We denote ∆x and ∆x ^⊥ as mesh sizes on Ω^h and Ω^⊥ h respectively. In some places, we report the results in terms of the number of mesh elements in horizontal direction N^x and vertical direction Ny.

Stability criterion. The stability criterion that mimics the continuous stability estimate in Theorem 3 is defined on the discrete level using the stability estimate stated in Proposition 8. We write the left and right-hand sides of the estimate as E n ^L = ∥h ⁿ+1∥ 2 L ²(Ω^⊥ h + 4∆t ρg √ µ *Du∥ 2 ²(Ωⁿ h* and E n ^R = ∥*h n*∥ L ²(Ω^⊥ h + 2∆t(a n ,h n )Ω^⊥ + (∆t) ²∥*a n*∥ 2 L ²(Ω^⊥ h respectively. Then we define the normalized energy at t n as:

$$\bar{E}^n = \frac{E_L^n - E_R^n}{\max_n |E_R^n|}, \qquad n = 0, 1, 2, \dots$$ (5.1)

where *E*¯ ⁿ ≤ 0 implies stability per time step. The above criterion is stated in the total energy norm. This norm reflects the coupling between u and h and provides a natural measure of stability. In contrast, separate estimates for u and h that ignore their interaction are not sufficient for analyzing the stability of the coupled system.

Abbreviations. To distinguish between the different numerical schemes, we use the following abbreviations: EE-Unstabilized for unstabilized explicit Euler scheme (4.12), EE-Unstabilized-W for EE-Unstabilized with ∆t computed using (4.25) (weakly stable), EE-Stabilized for the stabilized explicit Euler scheme (4.26), SIE-FSSA for the FSSA stabilized semi-implicit Euler scheme (4.29), EE-FSSA for the FSSA stabilized explicit Euler scheme.

5.1. Newtonian fluid dynamics in a tank

In this section, we consider a Newtonian fluid case. This implies a linear Stokes problem. The domain Ω^h is defined as in (2.6), where we set Ω^⊥ = [−1,1], b(x) = 0 and h(x) = 0.5 tanh(2x − 1) + 0.2. In (2.4) we set µ = 1 and −ρ*gz*ˆ = (0,−9.82). In Figure 2 we display the domain shape and a reference solution with a = 0 when the final simulation time is *t*ˆ = 4.

FIG. 2. (a) Initial domain for Newtonian fluid in a tank case. (b) free-surface height evolving in time from t = 0 to t = 4.

5.1.1. Stability comparison of the time stabilization approaches.

We use N^x = 120 and N^y = 120 mesh elements in each direction and use the normalized energy (5.1) as the stability criterion to verify the stability at different time steps when *t*ˆ = 4. Results are given in Figure 3. First, we observe that the EE-Unstabilized scheme is unstable for all the considered time steps, confirming our theoretical observation provided in Proposition 6. The edge stabilization term (4.1) on itself does not correct the stability in time. Second, we observe that the EE-Stabilized scheme, with our new stabilization approach, is always stable, which confirms our theoretical result in Proposition 8. Third, the EE-FSSA scheme and the SIE-FSSA scheme are only stable for some time steps.

We also test the numerical stability when the source term takes values a = a(x ^⊥,t) = 0.2(x ⊥) 2 (0.3+ sin(x ^⊥))sin(2t) . The result is given in Figure 4. We find the interpretation of these results analogous to those when a = 0.

Overall, we find the newly proposed EE-Stabilized scheme (4.26) best in terms of stability as it allows arbitrarily large time step sizes.

Newtonian fluid in a tank Normalized time-step energy difference 1.0 1.0

FIG. 3. Normalized consecutive time step energy difference (stability criterion) for a Newtonian fluid case with source term a = 0 compared across three different choices of time discretization. Other parameters are: viscosity µ = 0.3, uniform mesh N^x = N^y = 120, final time *t*ˆ = 4.

Newtonian fluid in a tank with source term Normalized time-step energy difference

FIG. 4. Normalized consecutive time step energy difference (stability criterion) for a Newtonian fluid case with source term a = 0.2(x (0.3+sin(x ^⊥))sin(2t) compared across three different choices of time discretization. Other parameters are: viscosity µ = 0.3, uniform mesh N^x = N^y = 120, final time *t*ˆ = 4.

5.1.2. Numerical errors under mesh and time step size refinement for different stabilization approaches.

We numerically test the convergence of the numerical errors in u, u ^⊥, and h towards 0. This is only possible when simultaneously refining ∆x → 0, ∆x ^⊥ → 0, and ∆t → 0. We examine how adding different stabilization terms affects the numerical error, i.e., to what extent the stabilizations are consistent with the discretized problem. In all the tested cases we simultaneously refine the mesh size and the time step size according to Table 5. We chose the time step sizes so that they satisfy the time step restriction (4.25), required to use EE-Unstabilized-W. The final simulation time is *t*ˆ= 1. We computed the relative numerical errors against a numerical reference using the EE-Unstabilized-W scheme with mesh sizes ∆x = 0.015, ∆x ^⊥ = 0.0125, and time step ∆t = 0.005.

Newtonian fluid in a tank Table of time steps and the corresponding mesh sizes for the convergence test

∆t	0.5	0.25	0.125	0.0625
⊥ ∆x	0.2	0.1	0.05	0.025
∆x	0.24	0.11	0.06	0.03

FIG. 5. Time step ∆t and the corresponding mesh sizes ∆x ^⊥ on Ω^⊥ h and ∆x on Ω^h for the error convergence test in the Newtonian fluid case.

Newtonian fluid in a tank Convergence of error under simultaneous mesh and time refinement

FIG. 6. Convergence of the relative numerical error (Newtonian fluid case) under simultaneous mesh and time refinement in surface height on Ω⊥, surface velocity on Ω⊥, and velocity on Ω, for different choices of discretization in time and the corresponding stabilization approaches. Parameters in all cases are: viscosity coefficient µ = 0.3, source term a = 0.

In Figure 6 we report the numerical errors under mesh refinement in the surface height function on Ω⊥, surface velocities on Ω⊥, and the velocities on Ω. We observe that the error in the surface height case and the surface velocity case converges with order 1 for all the considered schemes, which is expected. The EE-Stabilized scheme is significantly more accurate compared to the EE-Unstabilized-W scheme.

In the surface velocity case, the error is comparable across all the schemes.

In the velocity case, the error decays with order 2 when the EE-Unstabilized-W scheme is used, but the error decay is of order 1 in all other cases. This implies that the stabilization term added to the Stokes problem is consistent up to the first order in each of the cases EE-Stabilized, SIE-FSSA, and EE-FSSA. The velocity error is significantly smaller in the EE-Unstabilized-W case compared to all the stabilized cases. Among those, the error is smallest in the EE-Stabilized case.

Among all the stabilized schemes, the EE-Stabilized scheme induced the smallest errors. The main point of interest in the free-surface flow simulations is the surface height, where EE-Stabilized was, in addition, more accurate than EE-Unstabilized-W.

5.1.3. Conservation of domain volume

Ideally, a numerical scheme for solving the free-surface coupled Stokes problem conserves the domain volume in time, i.e., the domain volume error is at the level of round-off errors (close to machine precision). In this experiment, we use N^x = 120 and N^y = 120 mesh elements in each direction and set the final simulation time to *t*ˆ = 4. In Figure (7), we display the domain volume error for different time step sizes. In the EE-Stabilized case, we observe that the domain volume is conserved up to round-off errors, which is in line with Proposition 7. The EE-FSSA scheme also conserves the domain volume. On the other hand, the SIE-FSSA scheme does not conserve the domain volume as observed from the figure.

FIG. 7. Domain Ω volume error (Newtonian fluid case) as a function of time for different choices of discretization in time and the corresponding stabilizations. We use a uniform mesh N^x = N^y = 120, viscosity µ = 0.3, source term a = 0.

5.2. Motivation for suppressing high-frequency oscillations due to unresolved gradients.

Even when the solution is stable in accordance with the stability estimate stated in 6, the free-surface function can in practice entail high-frequency oscillations, particularly when the mesh resolution is large. Theoretical reasoning is provided in the scope of Section 4. To suppress the oscillations, we use the edge stabilization term (4.1). In Figure 8 we visualize the free-surface solution computed from the stabilized explicit Euler scheme (4.26) at different times with a small time step (∆t = 0.05) and a fine mesh (N^x = N^y = 300), with and without edge stabilization. We observe that the edge stabilization term effectively filters out high-frequency oscillations.

5.3. Non-Newtonian fluid dynamics over a simplified Greenland geometry

In this section, we experimentally verify that our theory is correct for non-Newtonian fluids (nonlinear Stokes problem) on a domain with a noisy boundary shape. To this end, we consider an idealized ice sheet flow over a slice of the Greenland geometry (similar to the simplified Greenland test case in [28]). Domain Ω^h is represented using a simplified Greenland geometry as illustrated in Figure 1. The Greenland shape was extracted from the 3D Greenland data set BedMachine [24] and interpolated using cubic splines to form the functions h(x ^⊥) and b(x ^⊥). In our simulations, we use realistic ice parameters in SI units. Corresponding to the Greenland data we set Ω^⊥ = [−428675 m,489475 m]. In (2.4) we set density ρ = 910kgm−³ , gravity g = 9.82ms−² . In (2.2) we set p = 4/3 and viscosity coefficient

Newtonian fluid in a tank: surface height oscillations Unstabilized Edge stabilized 0.50 0.75 1.00 1.25 0.50 0.75 1.00 1.25 t = 0 t = 0.5 t = 0.9 t = 1.5 t = 4.0

FIG. 8. Surface height solutions (Newtonian fluid) at different t in the EE-Stabilized case without additional stabilization of the high-frequency oscillations (left plot), and with stabilization of the high-frequency oscillations using the edge stabilization term (4.1). We use a uniform mesh N^x = N^y = 300, viscosity µ = 0.3, source term a = 0.

1.0 0.5 0.0 0.5 1.0 0.25

(ice softness parameter) µ⁰ = 3.1688 · 10−²⁴ (Pa) −3 s −1 . In all cases, we use N^x = 300 and N^y = 20 mesh elements in horizontal and vertical directions, consider a homogeneous source term a = 0, use the Picard iteration to solve the nonlinear Stokes problem with a relative convergence tolerance 10−⁶ , and perform the simulation until final time *t*ˆ = 200 years.

In Figure 9 we display the domain shape and a reference solution.

1.0 0.5 0.0 0.5 1.0 0.25

FIG. 9. (a) Initial Greenland slice domain for the non-Newtonian fluid case. (b) free-surface height evolving in time from t = 0 years to t = 200 years.

5.3.1. Stability comparison of the time stabilization approaches.

To examine numerical stability we measure the normalized energy difference (5.1) in each time step. Results are given in Figure 10. First, the EE-Unstabilized scheme is unstable for all the considered time steps, confirming the result in Proposition 6. Second, EE-FSSA and SIE-FSSA schemes are unstable when using the largest considered time step ∆t = 50 years. The edge stabilization term (4.1) does not prevent instabilities in time. Third, the EE-Stabilized scheme is stable for all the considered time steps, which confirms our theoretical result in Proposition 8.

To visually inspect the stability behaviour of the numerical solutions, we plot the free-surface height solutions in Figure 11, for all the considered numerical schemes. In the EE-Unstabilized and SIE-FSSA cases, the chosen time step corresponds to the minimal value under which the solution becomes unstable. In the figure we observe that although the EE-Unstabilized solution is inaccurate at t = 200

Non-Newtonian fluid: Ice sheet dynamics over a slice of Greenland Normalized time-step energy difference

FIG. 10. Normalized consecutive time step energy difference (stability criterion) for non-Newtonian fluid case with source term a = 0 compared across three different choices of time discretization. We use a non-uniform mesh with N^x = 300 elements in the horizontal direction and N^y = 20 elements in the vertical direction.

years due to a large time step ∆t = 50 years when compared to the reference solution in Figure 9, it is the only example of a solution that does not entail instabilities.

Non-Newtonian fluid: Ice sheet dynamics over a slice of Greenland Free-surface height solutions with different stabilizations in time

FIG. 11. Numerical free-surface height solutions varying in time, when using schemes: EE-Unstabilized with ∆t = 0.5 years, EE-Stabilized with ∆t = 50 years, and SIE-FSSA with ∆t = 50 years. Edge stabilization term is added to all of the schemes. We use N^x = 300 mesh elements in the horizontal direction and N^y = 20 mesh elements in the vertical direction.

5.3.2. Conservation of domain volume

In Figure (12), we display the relative domain volume error for different time step sizes. In the EE-Stabilized case, we observe that the domain volume is conserved up to round-off errors, which confirms Proposition 7. The EE-FSSA scheme also conserves the domain volume. On the other hand, the SIE-FSSA scheme does not conserve the domain volume as observed from the figure.

Non-Newtonian fluid: Ice sheet dynamics over a slice of Greenland Relative domain volume error in time

FIG. 12. Domain Ω volume error (non-Newtonian fluid case) as a function of time for different choices of discretization in time and the corresponding stabilizations. We use a non-uniform mesh with N^x = 300 elements in the horizontal direction and N^y = 20 elements in the vertical direction.

6. Final remarks

We have investigated the stability and conservation properties of the coupled Stokes/free-surface system. The analysis covered Newtonian (linear Stokes) and non-Newtonian (nonlinear Stokes) fluids, both at the continuous and fully discrete levels.

On the continuous level, we derived energy estimates that couple the momentum balance and surface height evolution. Using these estimates, we examined implicit, semi-implicit, and explicit Euler time discretizations.

On the discrete level, we showed that the implicit method is unconditionally stable and conserves volume. The semi-implicit scheme, even when stabilized with FSSA, is only conditionally stable and does not conserve volume. The explicit Euler method conserves volume but is unstable. We theoretically found that the cause of the instabilities is related to the normal surface velocity norm. We proposed a stabilization term that accurately eliminates these instabilities and showed unconditional stability with respect to ∆t, while preserving the domain volume. The stabilization term is symmetric, simple to implement, and works without tuning parameters. Numerical tests confirmed all the theoretical results.

The stabilized explicit scheme (4.26) that we proposed offers a robust and efficient alternative to more expensive implicit solvers. This is the scheme we recommend using in practice.

Future work may entail: extending the newly proposed stabilization approach to higher-order time integrators, investigating positivity preservation of the free-surface height h, and developing adaptive time-stepping strategies based on the derived energy balances.

7. Acknowledgements

Josefin Ahlkrona and Igor Tominec were funded by the Swedish Research Council, grant number 2021- 04001. Lukas Lundgren was funded by the Swedish e-Science Research Centre (SeRC).

REFERENCES

• 1. Martin Alnaes, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie Rognes, and Garth Wells. The FEniCS project version 1.5. Archive of Numerical Software, 3(100), 2015.
• 2. Martin S. Alnaes, Anders Logg, Kristian B. Ølgaard, Marie E. Rognes, and Garth N. Wells. Unified form language: A domain-specific language for weak formulations of partial differential equations. ACM Trans. Math. Softw., 40(2):9:1–9:37, March 2014.
• 3. Martin S. Alnaes, Anders Logg, Kristian B. Ølgaard, Marie E. Rognes, and Garth N. Wells. Unified Form Language: A Domain-specific Language for Weak Formulations of Partial Differential Equations. ACM Trans. Math. Softw., 40(2):9:1–9:37, March 2014.
• 4. Miguel Andres-Mart ´ ´ınez, Jason P. Morgan, Marta Perez-Gussiny ´ e, and Lars R ´ upke. A new free-surface ¨ stabilization algorithm for geodynamical modelling: Theory and numerical tests. Physics of the Earth and Planetary Interiors, 246:41–51, 2015.
• 5. Emmanuel Audusse, Gabriel R. Barrenechea, Astrid Decoene, and Pierrick Quemar. Stability analysis of a finite element approximation for the Navier-Stokes equation with free surface. ESAIM Math. Model. Numer. Anal., 58(1):107–130, 2024.
• 6. L. Belenki, L. C. Berselli, L. Diening, and M. Ru˚ziˇ cka. On the finite element approximation of ˇ p-Stokes systems. SIAM J. Numer. Anal., 50(2):373–397, 2012.
• 7. Ed Bueler. Surface elevation errors in finite element stokes models for glacier evolution, 2024.
• 8. Erik Burman and Peter Hansbo. Edge stabilization for Galerkin approximations of convection-diffusionreaction problems. Comput. Methods Appl. Mech. Engrg., 193(15-16):1437–1453, 2004.
• 9. Gong Cheng, Per Lotstedt, and Lina von Sydow. Accurate and stable time stepping in ice sheet modeling. ¨ Journal of Computational Physics, 329:29 – 47, 2017.
• 10. Alexandre Ern and Jean-Luc Guermond. Weighting the edge stabilization. SIAM J. Numer. Anal., 51(3):1655– 1677, 2013.
• 11. Alexandre Ern and Jean-Luc Guermond. Finite elements I—Approximation and interpolation, volume 72 of Texts in Applied Mathematics. Springer, Cham, 2021.
• 12. O. Gagliardini, T. Zwinger, F. Gillet-Chaulet, and et al. Capabilities and performance of Elmer/Ice, a new generation ice-sheet model. Geosci. Model Dev., 6:1299–1318, 2013.
• 13. JW Glen. The creep of polycrystalline ice. Proc. Roy. Soc. Lond. A, 228(1175):519–538, 1955.
• 14. R. Greve and H. Blatter. Dynamics of Ice Sheets and Glaciers. Advances in Geophysical and Environmental Mechanics and Mathematics (AGEM² ), Springer, Berlin, 2009.
• 15. Edward M. Hinton, Andrew J. Hogg, and Herbert E. Huppert. Shallow free-surface Stokes flow around a corner. Philos. Trans. Roy. Soc. A, 378(2174):20190515, 10, 2020.
• 16. Adrian Hirn. Finite element approximation of singular power-law systems. Mathematics of Computation, 82(283):1247–1268, 2013.
• 17. Antony Jameson. Origins and further development of the Jameson–Schmidt–Turkel scheme. AIAA J., 55:1–23, 03 2017.
• 18. Boris J.P. Kaus, Hans Muhlhaus, and Dave A. May. A stabilization algorithm for geodynamic numerical ¨ simulations with a free surface. Physics of the Earth and Planetary Interiors, 181(1):12 – 20, 2010.
• 19. Dmitri Kuzmin and Joshua Vedral. Dissipation-based WENO stabilization of high-order finite element methods for scalar conservation laws. J. Comput. Phys., 487:Paper No. 112153, 20, 2023.
• 20. A. Lofgren, J. Ahlkrona, and C. Helanow. Increasing stable time-step sizes of the free-surface problem arising ¨ in ice-sheet simulations. Journal of Computational Physics: X, 16:100114, 2022.
• 21. Lukas Lundgren, Christian Helanow, Jonathan Wiskandt, Inga Monika Koszalka, and Josefin Ahlkrona. A potential energy conserving finite element method for turbulent variable density flow: application to glacierfjord circulation. Journal of Computational Physics, page 113981, 2025.

• 22. Andre L ´ ofgren, Thomas Zwinger, Peter Raback, Christian Helanow, and Josefin Ahlkrona. Increasing ¨ numerical stability of mountain valley glacier simulations: implementation and testing of free-surface

• stabilization in elmer/ice. EGUsphere, 2023:1–22, 2023.

• 23. Michael Manga and H. A. Stone. Buoyancy-driven interactions between two deformable viscous drops. J. Fluid Mech., 256:647–683, 1993.

• 24. M. Morlighem, C. N. Williams, E. Rignot, L. An, J. E. Arndt, J. L. Bamber, G. Catania, N. Chauche, J. A. ´ Dowdeswell, B. Dorschel, I. Fenty, K. Hogan, I. Howat, A. Hubbard, M. Jakobsson, T. M. Jordan, K. K. Kjeldsen, R. Millan, L. Mayer, J. Mouginot, B. P. Y. Noel, C. O'Cofaigh, S. Palmer, S. Rysgaard, H. Seroussi, ¨ M. J. Siegert, P. Slabon, F. Straneo, M. R. van den Broeke, W. Weinrebe, M. Wood, and K. B. Zinglersen. BedMachine v3: Complete bed topography and ocean bathymetry mapping of greenland from multibeam echo sounding combined with mass conservation. Geophysical Research Letters, 44(21):11,051–11,061, 2017.

• 25. Murtazo Nazarov. Convergence of a residual based artificial viscosity finite element method. Comput. Math. Appl., 65(4):616–626, 2013.

• 26. Ian Rose, Bruce Buffett, and Timo Heister. Stability and accuracy of free surface time integration in viscous flows. Physics of the Earth and Planetary Interiors, 262:90–100, January 2017.

• 27. Vidar Stiernstrom, Lukas Lundgren, Murtazo Nazarov, and Ken Mattsson. A residual-based artificial viscosity ¨ finite difference method for scalar conservation laws. J. Comput. Phys., 430:Paper No. 110100, 29, 2021.

• 28. Igor Tominec and Josefin Ahlkrona. Weak form shallow ice approximation models with an improved time step restriction. Journal of Glaciology, page 1–46, 2025.

• 29. Igor Tominec and Murtazo Nazarov. Residual viscosity stabilized RBF-FD methods for solving nonlinear conservation laws. J. Sci. Comput., 94(1):Paper No. 14, 31, 2023.

• 30. Junjie Wen and Murtazo Nazarov. An anisotropic nonlinear stabilization for finite element approximation of Vlasov–Poisson equations. J. Comput. Phys., 536:Paper No. 114079, 2025.

• 31. A. Wirbel and A. H. Jarosch. Inequality-constrained free-surface evolution in a full Stokes ice flow model (evolve glacier v1.1). Geoscientific Model Development, 13(12):6425–6445, 2020.