Thickness Stabilization Scheme (TSS)

• TSS is a numerical stabilization method for the Shallow Shelf Approximation (SSA) that predicts next-step squared thickness to mitigate the explicit lag in the driving stress.
• The scheme enables significantly larger and stable timesteps in ice sheet simulations, with tests showing improvements in both low- and high-shear regimes.
• TSS retains a partitioned momentum–thickness workflow by modifying the driving stress term and incorporating semi-implicit thickness evolution with controlled artificial viscosity.

Searching arXiv for recent and directly relevant papers on “Thickness Stabilization Scheme” and related stabilization contexts. Thickness Stabilization Scheme (TSS) is a numerical stabilization scheme for the Shallow Shelf Approximation (SSA) in which an adapted, explicit Euler thickness evolution equation is inserted into the driving stress term, thereby treating that term implicitly (Dolling et al., 3 Oct 2025). It was proposed for partitioned SSA solvers, where the momentum equations are solved with the current thickness and the thickness is then advanced separately. In that setting, the explicit lag in the geometry-dependent driving force can trigger numerical instabilities when the time-step size is too large. TSS targets that specific coupling defect by predicting the next-step thickness contribution inside the SSA forcing, while retaining a partitioned momentum–thickness workflow. In the reported idealized low- and high-shear scenarios, the method enlarges the numerically stable timestep range substantially, including stable runs at $\Delta t = 50$–$100$ years and even $\Delta t = 10{,}000$ years, although the latter entails a reduction in accuracy (Dolling et al., 3 Oct 2025).

1. Definition and problem class

TSS is defined for SSA, one of the widely used vertically integrated models for large-scale Antarctic and Greenland ice sheet simulations. The motivating instability arises because the SSA driving stress depends on thickness, but standard partitioned time stepping evaluates that force with the old thickness $H_k$ and only afterwards updates the geometry. The paper identifies this explicit handling of the driving force as the source of the restrictive timestep limit (Dolling et al., 3 Oct 2025).

For floating ice in hydrostatic equilibrium, the driving term becomes proportional to $\nabla H^2$. This makes the momentum–geometry feedback particularly transparent: thickness changes alter the driving stress; the driving stress alters velocity; and the velocity then enters the thickness evolution equation. TSS addresses the lag in this loop without moving to a monolithic coupled nonlinear solve. The scheme therefore belongs to the class of semi-implicit stabilizations embedded in partitioned frameworks, rather than to fully implicit or purely explicit time integrators (Dolling et al., 3 Oct 2025).

The method is formulated for the right-hand-side driving stress in SSA, not for every thickness-dependent coefficient in the model. In particular, the paper focuses on correcting the explicit treatment of the forcing term rather than fully implicitly coupling all appearances of thickness. This design choice is central to the balance between increased numerical stability and modest per-step complexity (Dolling et al., 3 Oct 2025).

2. Governing equations and derivation

The starting point is the SSA momentum balance for floating ice shelves, together with the thickness evolution equation. In compact form, the vertically integrated momentum equation is written as

$\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$

where

$$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$

$H$ is ice thickness, and $\mathbf{T}$ is the vertically integrated stress tensor (Dolling et al., 3 Oct 2025).

The thickness equation is

$$\frac{\partial H}{\partial t} = - \nabla \cdot (H \mathbf{u}) + a_s - a_b,$$

with $100$0 the surface accumulation and $100$1 the melt rate (Dolling et al., 3 Oct 2025). Standard explicit Euler gives

$100$2

while the semi-implicit backward Euler form used for the actual geometry update is

$100$3

TSS is derived by multiplying the thickness equation by $100$4 and discretizing the resulting evolution equation for $100$5 with forward Euler. The paper gives the predicted next-step squared thickness as

$100$6

This predictor is then inserted into the momentum equation so that the forcing behaves as if it were treated implicitly (Dolling et al., 3 Oct 2025).

In weak form, the stabilized SSA equation becomes

$100$7

with user-controlled stabilization parameter $100$8, where $100$9 corresponds to no TSS and $\Delta t = 10{,}000$0 to full TSS (Dolling et al., 3 Oct 2025).

This construction does not alter the physical model itself. It modifies the numerical treatment of the driving force by embedding a first-order predictor for the geometry-dependent term inside the momentum balance. A plausible implication is that the method can be viewed as a one-step geometric predictor attached specifically to the forcing channel of the partitioned solver.

3. Algorithmic structure

The solver workflow remains partitioned. At each time step, the momentum equation is solved first, using either the standard SSA weak form or the TSS-modified form, and the resulting velocity is then used in the thickness update (Dolling et al., 3 Oct 2025). The geometry evolution itself still uses the semi-implicit thickness equation

$\Delta t = 10{,}000$1

The paper is explicit that TSS is constructed with an explicit Euler time discretization and inserted into the SSA formulation, whereas the geometry evolution is solved using an implicit Euler time discretization.

The implementation described in the paper also applies artificial viscosity in the thickness solve:

$\Delta t = 10{,}000$2

with

$\Delta t = 10{,}000$3

and $\Delta t = 10{,}000$4 (Dolling et al., 3 Oct 2025). The paper also enforces a minimum thickness constraint,

$\Delta t = 10{,}000$5

This algorithmic organization preserves the conventional momentum-first, thickness-second operator split. TSS therefore does not replace the thickness equation or merge the stages into a single nonlinear system; it changes the momentum stage so that the explicit geometry lag in the driving stress is reduced. In practical terms, the method is intended for modelers who want larger usable time steps without abandoning a partitioned SSA–thickness architecture (Dolling et al., 3 Oct 2025).

4. Stability, accuracy, and timestep behavior

The analysis in the paper is empirical rather than a formal amplification-factor proof. Its conceptual claim is that instability originates in the explicit lagging of the driving force and that TSS enlarges the stable timestep range by approximately evaluating that force at the next time step (Dolling et al., 3 Oct 2025). The reported behavior is summarized below.

Regime	Unstabilized behavior reported	TSS behavior reported
Low shear	stable and physically realistic only up to $\Delta t = 10{,}000$6 years	$\Delta t = 10{,}000$7–$\Delta t = 10{,}000$8 years remain stable and accurate
High shear	only up to $\Delta t = 10{,}000$9 years	$H_k$0 and $H_k$1 years have much smaller errors
Very large timestep	—	$H_k$2 years is numerically stable, but inaccurate

In the low-shear case, the paper states that TSS allows $H_k$3 years to remain within 1% accuracy, whereas without TSS, $H_k$4 years already produces significant error (Dolling et al., 3 Oct 2025). The quantitative error table after 2000 years shows the trade-off clearly. At $H_k$5, thickness error is $H_k$6 without TSS and $H_k$7 with TSS; velocity error is $H_k$8 without TSS and $H_k$9 with TSS. At $\nabla H^2$0, thickness error is $\nabla H^2$1 without TSS and $\nabla H^2$2 with TSS; velocity error is $\nabla H^2$3 without TSS and $\nabla H^2$4 with TSS (Dolling et al., 3 Oct 2025).

The high-shear regime is more restrictive without stabilization and more favorable to TSS. At $\nabla H^2$5, thickness error is $\nabla H^2$6 without TSS and $\nabla H^2$7 with TSS; velocity error is $\nabla H^2$8 without TSS and $\nabla H^2$9 with TSS. At $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$0, thickness error is $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$1 without TSS and $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$2 with TSS; velocity error is $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$3 without TSS and $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$4 with TSS. The paper summarizes these differences by stating that at $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$5, thickness error is $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$6 times larger without TSS and velocity error $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$7 times larger, while at $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$8, thickness error is $\nabla \cdot \mathbf{T} = \rho' \nabla H^2,$9 times larger without TSS and velocity error $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$0 times larger (Dolling et al., 3 Oct 2025).

The paper also defines the update norm

$$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$1

using its decay toward zero as an indicator of approach to steady state. This is not presented as a formal proof of convergence, but as a diagnostic accompanying the time histories and domain-averaged errors (Dolling et al., 3 Oct 2025).

A recurring caveat is that stability is not accuracy. TSS can make extremely large time steps numerically stable, but this does not guarantee physically realistic results. The clearest example is $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$2 years, which remains stable with TSS yet exhibits large errors and does not approach the same steady-state thickness within the reported $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$3-year horizon (Dolling et al., 3 Oct 2025).

5. Numerical experiments, use cases, and limitations

The numerical experiments use an idealized floating ice shelf in a $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$4 km square domain with a central circular obstacle of radius $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$5 km at $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$6 km, representing a pinning point (Dolling et al., 3 Oct 2025). Boundary conditions include constant normal inflow velocity at $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$7, force balance at the outflow or calving front $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$8, and impenetrability on the lateral sides and cylinder boundary. Initial thickness is spatially uniform and fixed at the inflow boundary. This setup is specifically intended to create shear margins and strong geometry–flow coupling around the obstacle.

Two principal regimes are studied. The low-shear case uses $$\rho' = \frac{1}{2}\rho_i g \bigg(1-\frac{\rho_i}{\rho_o} \bigg),$$9 m and $H$0 m/yr, with tested time steps $H$1 years. The high-shear case uses $H$2 m and $H$3 m/yr, with tested time steps $H$4 years. Additional TSS-only runs in the low-shear setup use $H$5 years up to $H$6 years of simulated time (Dolling et al., 3 Oct 2025).

The practical significance of TSS is tied to timestep flexibility rather than to a cheaper single step. The paper notes that runtime with and without TSS is similar for a fixed $H$7, so the benefit comes from safely using larger timesteps. One example states that $H$8 takes about 40 seconds, roughly $H$9 faster than $\mathbf{T}$0 for the same $\mathbf{T}$1-year simulation (Dolling et al., 3 Oct 2025). This makes TSS relevant when computational resources must be reallocated toward higher spatial resolution, longer integrations, larger ensembles, or coupled workflows.

The method is nevertheless bounded by the assumptions of the paper. The presented formulation is limited to floating ice shelves and does not include basal friction. SSA itself neglects vertical shear, so TSS improves the numerical properties of SSA but does not remove SSA model-form limitations. The experiments are idealized rather than fully realistic, and the implementation uses first-order finite elements together with artificial viscosity, so the reported timestep thresholds should not be read as mesh-independent universal limits (Dolling et al., 3 Oct 2025).

The paper explicitly mentions potential relevance beyond the specific setup, including DIVA, Blatter-Pattyn-type approximations, ISCAL, and coupled lower/higher-order model hierarchies. It also suggests applicability to other vertically integrated computational fluid dynamics problems that couple momentum and geometry evolution equations (Dolling et al., 3 Oct 2025). This suggests a broader interpretation of TSS as a stabilization pattern for lagged geometry-dependent forcing in partitioned reduced-order flow models.

6. Terminology, related stabilization ideas, and common ambiguities

The named method “Thickness Stabilization Scheme” appears explicitly in the SSA context just described (Dolling et al., 3 Oct 2025). In adjacent literatures, however, closely related ideas exist without that name. In Cahn–Hilliard numerics, for example, stabilization is often organized around the interface thickness parameter $\mathbf{T}$2, but the relevant schemes are not called TSS. “Convergence Analysis of an Unconditionally Energy Stable Linear Crank-Nicolson Scheme for the Cahn-Hilliard Equation” proposes a stabilized linear Crank–Nicolson scheme with two linear stabilization terms and emphasizes that $\mathbf{T}$3 is the interface thickness parameter; the paper is directly concerned with thin-interface stiffness and with error prefactors that depend polynomially rather than exponentially on $\mathbf{T}$4, but it does not define a Thickness Stabilization Scheme (Wang et al., 2017). “An Energy Stable Linear Diffusive Crank-Nicolson Scheme for the Cahn-Hilliard Gradient Flow” develops a related stabilized semi-implicit diffusive Crank–Nicolson method, again with rigorous energy stability and polynomial-in-$\mathbf{T}$5 error behavior, but again not under the name TSS (Wang et al., 2020).

A second nearby usage appears in codimensional contact simulation. “Unlocking Thickness Modeling for Codimensional Contact Simulation” is not phrased in terms of TSS, but its Barrier-Filtering method functions as a thickness stabilization mechanism by reducing the effective barrier activation distance for local stencil-neighbor pairs with $\mathbf{T}$6, thereby removing contact-locking artifacts while retaining IPC-style non-intersection robustness (Gomez-Nogales et al., 11 Oct 2025). This suggests that “thickness stabilization” can refer either to stabilization of a thickness-dependent forcing term, as in SSA, or to stabilization of thickness-aware contact processing, as in codimensional mechanics.

The acronym TSS is also not unique across arXiv. In “Tele-Operated Oropharyngeal Swab (TOOS) Robot Enabled by TSS Soft Hand for Safe and Effective COVID-19 OP Sampling,” TSS stands for “TStone Soft,” referring to a soft robotic hand rather than to any stabilization scheme (Chen et al., 2021). This distinction is important because the SSA method is a numerical stabilization of momentum–thickness coupling, whereas the robotics usage is a named hardware architecture.

In that broader terminological landscape, Thickness Stabilization Scheme most precisely denotes the SSA method that predicts the next-step squared thickness $\mathbf{T}$7 and inserts that prediction into the driving stress, thereby mimicking implicit treatment of the geometry-dependent forcing while preserving a partitioned solve structure (Dolling et al., 3 Oct 2025).