A numerical stabilization scheme for the shallow shelf approximation

TL;DR

The paper addresses the stability bottleneck of depth-averaged SSA ice-sheet simulations by introducing the Thickness Stabilization Scheme (TSS), which mimics an implicit driving-stress treatment without solving a full 3D velocity field. TSS modifies the SSA weak form by replacing the current thickness term with a predicted next-step thickness, controlled by a stabilization parameter $\theta$, and is complemented by an artificial-viscosity term to suppress spurious oscillations. Numerical experiments in a 2D domain with a central obstacle demonstrate that TSS dramatically increases the largest stable time step from the order of decades to up to $10^4$ years, while maintaining acceptable accuracy in many regimes; low-shear cases can achieve $\Delta t$ of 100 years with errors under 1% in $H$ and $|\mathbf{u}|$, whereas high-shear cases show order-of-magnitude improvements in error when using TSS. The method reduces computational cost and can be extended to other reduced-order or vertically-integrated frameworks, offering a practical pathway to efficient, large-scale ice-sheet simulations and potentially to other shallow-fluid problems coupled to geometry evolution.

Abstract

We present the Thickness Stabilization Scheme (TSS), a numerical stabilization scheme suitable for the Shallow Shelf Approximation (SSA), one of the most widely-used models for large-scale Antarctic and Greenland ice sheet simulations. The TSS is constructed by inserting an adapted, explicit Euler thickness evolution equation into the driving stress term, thereby treating the term implicitly. We investigate the applicability of TSS across low- and high-shear idealized scenarios, by altering the inflow velocity and initial ice thickness. TSS demonstrates an increase in the numerical stability of SSA, allowing large time-step sizes of dt = 50-100 years to remain numerically stable and accurate, while time-step sizes of dt > 5 in high-shear simulations show significant error without TSS. Remarkably, a time step size as large as dt = 10 000 years is numerically stable with TSS, albeit with a reduction in accuracy. TSS offers greater flexibility for ice-sheet modeling by allowing the re-allocation of computational resources. This method is applicable not only to ice-sheet modeling, including in coupled frameworks, but also to other vertically-integrated computational fluid dynamics problems that couple momentum and geometry evolution equations.

Figures (14)

• Figure 1: The flow of ice in a marine-terminating domain. The upper ice surface, $z_s(x, y, t)$, represents the ice-atmosphere interface, the lower ice surface, $z_b(x, y, t)$, represents the ice-bed or ice-ocean interface, and $b(x, y)$ represents the bed elevation. The ice sheet flows with contact to the bedrock below ($z_b = b$) before beginning to float ($z_b > b$) at the grounding line to form an ice shelf. The otherwise floating ice can come into contact with anomalies in the bedrock to form pinning points referred to as ice rises or ice rumples, where $z_b = b$. The ice thickness is the vertical distance between the upper and lower ice surfaces, i.e. $H(x, y, t) = z_s(x, y, t) - z_b(x, y, t)$.
• Figure 2: A schematic of the 2D model domain, $\Omega$, used in the simulations. The square domain represents a floating ice shelf, with boundaries indicated by $\Gamma_\mathrm{in}$ (inflow), $\Gamma_\mathrm{out}$ (calving front), $\Gamma_\mathrm{left}$ and $\Gamma_\mathrm{right}$ (lateral sides), and $\Gamma_\mathrm{cyl}$ (cylindrical obstacle). Cross sections for simulation analysis at $x=15$ km (blue), $x=35$ km (yellow), and $x=50$ km (green) are indicated.
• Figure 3: The $H_0 = 300$ m, $u_{y,0} = 300$ m/yr reference simulation at $2\,000$ years with a time-step size of $\Delta t = 0.5$ yrs and no TSS ($\theta = 0$). The panels show (a) the velocity streamlines around the obstacle with the color showing the velocity magnitude, and (b) the ice thickness field, $H$.
• Figure 4: The final ice thickness, $H$, after $2\,000$ years of simulation with a time-step size $\Delta t = 40$ years for the $H_0=300$ m, $u_{y,0}=300$ m/yr initial conditions. The panels show (a) the ice thickness with TSS ($\theta = 1$), (b) the ice thickness without TSS ($\theta = 0$), and (c) the difference between the simulation with TSS and the simulation without TSS. The gray circle centered at $(x, y) = (50, 50)$ km indicates the location of the obstacle.
• Figure 5: The absolute velocity, $|\mathbf{u}|$(m/yr), after $2\,000$ years of simulation with a time-step size of $\Delta t = 40$ years for the $H_0=300$ m, $u_{y,0}=300$ m/yr initial conditions. The panels show (a) the velocity field with TSS ($\theta = 1$), (b) the velocity field without TSS ($\theta = 0$), (c) the difference between the simulation with TSS and the simulation without TSS. The gray circle centered at $(x, y) = (50, 50)$ km shows the location of the obstacle.