PaperlayerView exact revision on arXiv

Source-linked AI summary

A Discontinuous Galerkin discretization for the sea ice dynamics

Emma Lagracie, Thomas Richter

arXiv:2608.26810v1math.NA

TL;DR

Sea ice simulations must resolve sharp LKFs despite nonlinear rheology and discretization sensitivity. This paper develops a fully DG Hibler VP discretization and proves convergence of the continuous mEVP pseudo-time system; the DG method reproduces sharp LKFs, with second-order solutions more consistent across meshes.

  • Problem

    Sea ice numerics must accurately resolve sharp LKFs, whose number, orientation, and distribution are sensitive to spatial resolution and velocity discretization.

  • Method

    The paper uses a fully discontinuous LDG discretization for Hibler’s VP model and studies continuous pseudo-time mEVP dynamics using Lyapunov functionals.

  • Results

    The DG method reproduces sharp, numerous LKFs, while second-order solutions remain qualitatively consistent across grid types and resolutions; the continuous mEVP system converges to the VP limit.

  • Takeaways & Limitations

    Higher-order fully discontinuous discretization can improve qualitative robustness and produce sharp localized LKFs without requiring very fine meshes.

  • Takeaways & Limitations

    The LDG implementation has higher computational cost because discontinuous spaces introduce additional degrees of freedom and interface terms.

Abstract

from arXiv · show

Sea ice dynamics plays a crucial role in the Earth's climate system, making it an important component of weather and climate prediction models. Its numerical simulation remains challenging, however, as it exhibits complex mechanical behaviors due to a nonlinear ice rheology and interactions with external physical forcings. In particular, sea ice presents linear kinematic features (LKFs), i.e., narrow bands of intense deformation associated with processes such as lead opening or pressure-ridge formation. Accurately representing these features is necessary as they affect thermodynamics and ocean-atmosphere exchange. Yet their number, localization, and structure are highly sensitive to spatial resolution and to the chosen discretization of the sea ice velocity. In this work, we investigate the spatial discretization of Hibler's viscous plastic sea-ice model using a fully discontinuous Galerkin (DG) representation of all variables, including the velocity field. The ability of DG elements to represent discontinuities while having high order local polynomial approximation makes them well suited for resolving the complex ice deformation features, and insuring robustness towards mesh-induced numerical artefacts. We assess the fully DG method using an established sea ice dynamics benchmark and compare the obtained sea ice deformation with state-of-the-art discretizations. Additionally, we study the theoretical convergence of the modified Elastic-Viscous-Plastic (mEVP) formulation as a pseudo-time iterative solver for the viscous-plastic (VP) momentum equations. We prove the convergence of the underlying continuous pseudo-time dynamical system towards the viscous plastic formal limit, thereby providing a theoretical foundation for the use of mEVP as an iterative solver for Hibler's VP model.

1 Otto von Guericke Universit¨at, Magdeburg, Germany

The paper identifies sea ice dynamics, mEVP, the elastic-viscous-plastic model, and Local Discontinuous Galerkin as keywords.

  • Sea ice dynamics is a listed keyword.
  • mEVP is a listed keyword.
  • Local Discontinuous Galerkin is a listed keyword.

1 Introduction

The paper addresses numerical challenges in sea ice dynamics, especially the resolution of sharp LKFs, through a fully discontinuous Galerkin discretization and a convergence analysis of mEVP.

  • Sea ice models are important for climate and forecasting because sea ice influences energy budgets, ocean-atmosphere exchanges, and polar circulation.
  • LKFs are narrow, intense deformation bands whose faithful resolution is a central challenge in sea ice numerics.
  • The Hibler VP model is numerically difficult because nonlinear viscosities vary widely and produce nearly singular discrete operators.
  • Existing implicit solvers can be accurate but remain difficult to precondition, implement, tune, and scale as resolution increases.
  • The paper proves convergence of the underlying continuous pseudo-time mEVP system toward the VP limit, complementing earlier discrete numerical convergence results.
  • The paper presents a first fully DG discretization of Hibler’s VP model using the LDG framework, with discontinuous velocity and strain-rate representations.

2 The sea ice model

The model couples sea-ice momentum, nonlinear VP rheology, and transport equations for thickness and concentration, while mEVP introduces pseudo-time relaxation for iterative solution.

  • The VP model describes sea ice as a two-dimensional continuum with velocity, thickness, and concentration variables.
  • The momentum equation includes internal stress divergence, atmospheric wind stress, ocean drag, and Coriolis forcing.
  • Atmospheric and ocean stresses depend on relative velocities between the ice, wind, and ocean current.
  • The non-stress forces are grouped into F(u), combining wind stress, ocean drag, and Coriolis effects.
  • Thickness H and concentration A satisfy standard advective transport equations and are not the focus of this work.
  • The VP rheology defines strain from velocity and determines stress through nonlinear constitutive relations involving viscosities and yield-curve parameters.
  • mEVP introduces artificial elastic relaxation, promotes stress to an independent variable, and recovers the VP solution formally as subiterations converge.

3 Convergence of the dynamical mEVP system

The continuous pseudo-time mEVP system is analyzed as a dynamical system whose equilibrium is the physical-time VP solution. Lyapunov arguments establish exponential convergence of velocity and stress under stated assumptions, while strain convergence requires additional coercivity information.

  • Dynamical-system formulation: The mEVP equations are interpreted as a pseudo-time explicit Euler discretization of a continuous dynamical system whose formal limit is the implicitly time-discretized VP momentum equation.Pseudo-time relaxation is considered over each fixed physical time step.
  • Variational structure: The VP rheology and forcing contributions derive from convex potentials, making the semi-discrete VP solution the unique minimizer of a strongly convex functional.The proof uses convexity of the VP potential, forcing potential, and combined functional.
  • Convergence results: The full stress variable and its temporal derivative converge with exponential rate, and the combined state converges in L2(Ω) × Hdiv(Ω) toward the VP solution.The divergence convergence follows from the velocity and stress estimates.
  • Limitations: Convergence of the strain tensor does not follow from velocity and stress convergence without additional assumptions because the VP potential lacks coercivity.A lower viscosity bound can yield boundedness and weak, then strong, L2 strain convergence.

4 A mixed form DG discretization for the sea ice model

The paper formulates sea-ice dynamics with an LDG mixed DG discretization for the mEVP system, using discontinuous velocity, stress, and strain approximations. It combines this spatial scheme with IMEX pseudo-time stepping, stabilization, tracer transport, and an element-coupled solution strategy.

  • 4.1 LDG scheme for the inner mEVP system: The mEVP momentum equation is discretized in mixed form with a Local Discontinuous Galerkin scheme.The formulation treats the system as first order and applies LDG numerical fluxes across element faces.
  • 4.1 LDG scheme for the inner mEVP system: Velocity and stress are represented in discontinuous polynomial spaces on triangular or quadrilateral meshes, with compatible trial spaces.The velocity space uses order n > 0, while the symmetric-matrix stress space uses order m ≥ 0.
  • 4.1 LDG scheme for the inner mEVP system: The elementwise weak formulation couples momentum, strain, and stress equations through numerical boundary fluxes.The formulation uses element scalar products, numerical fluxes on boundaries, and outward normals before summing over mesh elements.
  • 4.1 LDG scheme for the inner mEVP system: The LDG fluxes impose no-slip velocity and transparent no-jump stress-normal boundary conditions, with parameters (a, b) ∈ [0, 0.5) × (0, +∞).The boundary stress-normal flux is expressed using the stress, velocity, and flux parameters.
  • 4.1 LDG scheme for the inner mEVP system: For linear rheologies σVP(ε) = Lε with symmetric positive-definite L, an L2-energy estimate shows that the LDG discretization adds no instability.The result attributes this stability to compensation by the selected numerical fluxes.
  • 4.1 LDG scheme for the inner mEVP system: Optimal convergence of the primal variable requires velocity-jump stabilization b to scale as 1/h.In the sea-ice formulation, this stabilization also controls the auxiliary strain variable.
  • 4.2 Fully discontinuous sea-ice dynamics: The complete DG-IMEX mEVP scheme uses IMEX Euler discretization in pseudo-time together with implicit Euler discretization in physical time.The physical solution is advanced through pseudo-time states for velocity and stress.
  • 4.2 Fully discontinuous sea-ice dynamics: The mixed DG discretization is implicit, while forcings and VP rheology are explicit in pseudo-time; large viscosities require b ∝ 1/h scaled by maximal viscosity.The implementation sets a = 0.4 and b = 10^9/h, while a has negligible influence on the numerical solution.

5 Assessment on the sea ice benchmark

The benchmark evaluates LDG sea-ice simulations across mesh types, resolutions, approximation orders, and comparisons with established discretizations. Results show strong first-order grid dependence, while second-order LDG substantially reduces it and produces qualitatively improved LKF patterns.

  • Benchmark setup: The benchmark solves sea-ice dynamics on a 512 km × 512 km square for 2 days with prescribed wind and ocean forcings.Homogeneous no-slip velocity conditions and specified initial velocity, height, and concentration are used.
  • Discretizations and comparisons: LDG simulations use first- or second-order velocity DG spaces on triangular and quadrilateral meshes, with compatible stress and strain spaces.The implementation uses Δt = 6 min, α = β = 1, Δτ = 1000, and 400 mEVP subiterations.
  • Comparison with established methods: LDG produces a similar number of LKFs to ICON and Gascoigne CR-Q0 at both first and second order.The compared nonconforming elements also permit velocity discontinuities and generate numerous sharp LKFs.
  • Grid dependence: At first order, triangular and quadrilateral meshes produce substantially different LKF orientations and locations across mesh resolutions.Solutions cluster into distinct triangular- and quadrilateral-mesh patterns, with same-grid nonconforming comparisons showing similar behavior.
  • Grid dependence: Total LKF number and cumulative length may be unreliable comparison metrics for first-order approximations because numerical artifacts can dominate, especially on triangular meshes.Artifacts are classified as LKFs that disappear when approximation order and mesh resolution increase.
  • Higher-order behavior: Second-order LDG LKF patterns remain broadly consistent across grid types, reducing mesh dependence and improving qualitative accuracy even on coarse meshes.Quadrilateral solutions are smoother than same-order triangular solutions but require larger strain spaces, increasing computational cost.

6 Discussion

The paper establishes continuous-system convergence for mEVP and assesses a fully discontinuous Galerkin discretization for Hibler’s viscous-plastic model. The numerical results show sharper, more robust LKFs with second-order LDG, while highlighting mesh and computational-cost limitations.

  • mEVP convergence: The continuous pseudo-time mEVP system converges toward the viscous-plastic limit at fixed physical time, independently of spatial discretization.This complements earlier stability and discrete-convergence results and theoretically justifies mEVP as a pseudo-implicit VP solver.
  • Fully discontinuous Galerkin: The fully LDG formulation represents every variable, including velocity, in discontinuous finite-element spaces and is evaluated on a standard sea-ice dynamics benchmark.The method targets Hibler’s viscous-plastic sea-ice model.
  • Benchmark assessment: The benchmark reproduces sharp and numerous LKFs, although the fully discontinuous method incurs higher computational cost.This establishes the method’s qualitative ability to capture localized deformation while identifying an efficiency trade-off.
  • Mesh and approximation effects: For first-order velocity approximation, triangular versus quadrilateral grids strongly changes LKF number, orientation, and spatial distribution.Similar patterns across methods on the same grid indicate that some structures may be mesh-oriented rather than physically determined.
  • Mesh and approximation effects: Second-order LDG solutions remain qualitatively consistent across grid types and resolutions, producing sharp localized LKFs on relatively coarse meshes.Increasing approximation order can remove features interpreted as numerical artifacts and improve qualitative robustness without very fine meshes.
  • Limitations and future work: The LDG implementation requires optimization because discontinuous spaces add degrees of freedom and interface terms relative to continuous formulations.Future work includes integrating LDG into neXtSIM-DG alongside high-order tracer-advection methods.

A Appendix

The appendix develops energy estimates supporting stability analysis of the semi-discretized mEVP LDG scheme. The proof combines pointwise inequalities, forcing-term bounds, Bregman-divergence estimates, and a positive-definite rheology assumption.

  • Auxiliary inequalities: The appendix obtains additional lower and upper bounds through pointwise algebra, integration, and bounds for a cubic polynomial on [0, 1].These intermediate inequalities support the final energy estimates.
  • Convexity estimate: A Bregman-divergence representation is inserted into the energy estimate to derive a differential inequality leading to the stated energy bound.The divergence is defined from ψ by Dψ(v, v*) = ψ(v) − ψ(v*) − ∇ψ(v*)·(v − v*).
  • Stability proposition: The stability proposition assumes a symmetric positive-definite rheology matrix L, so (., L^-1.) defines a scalar product.Under this assumption, the semi-discretized mEVP LDG system satisfies an energy estimate.
  • Energy estimate: The discrete energy is Eh = 1/2∥uh∥2 + 1/2(σh, L^-1σh), combining velocity and stress contributions.The scalar product is taken in the discrete velocity or stress spaces.
  • Energy estimate: The proof tests the semi-discrete equations with stress, velocity, and transformed stress functions before summing them to obtain an estimate.The selected test functions are τh = σh, vh = uh, and th = L^-1σh.
  • Forcing and cancellation: The Coriolis term cancels in the energy estimate because (k⃗ × u)·u = 0, while forcing terms are bounded using generalized Young’s inequality.The argument separately treats contributions involving ocean velocity and air forcing.
Loading 2608.26810v1…