@article{Bensassi2026, 
author = {Salim Bensassi and Boualem Khouider and Clint Seinen and M'hamed Kesri},
title = {Numerical solution of the Hibler model of sea-ice dynamics: Comparing backward Euler and Crank-Nicolson methods},
year = {2026},
journal = {AIMS Mathematics},
volume = {11},
number = {3},
pages = {7980-8013},
keywords = {viscous-plastic sea-ice equations, Jacobian-free Newton-Krylov, numerical solutions, Crank-Nicolson, backward Euler, numerical stability},
url = {https://www.sciopen.com/article/10.3934/math.2026329},
doi = {10.3934/math.2026329},
abstract = {The viscous-plastic momentum equations for sea-ice dynamics are a model for the time evolution of the horizontal velocity field of a vertically homogeneous column of ice-pack continuum whose thickness and ice concentration vary in time and space. They constitute a highly nonlinear system of partial differential equations, and as such, both the mathematical analysis and numerical solution of its solution remain a challenge. To shed some light on this problem, we compared the performance of two known time discretization methods: the backward Euler (BE) and the Crank-Nicolson (CN) methods. Both methods are in theory unconditionally stable, but Euler's method is only first-order accurate while Crank-Nicolson is second-order. Centered finite differences were used for the spatial derivatives. This led to a nonlinear system of algebraic equations which was then solved using a Jacobian-free Newton-Krylov approach. First, the two methods were compared in terms of their ability to reproduce a synthetic solution to the viscous-plastic momentum equations when an artificial forcing is applied. The convergence of the two methods is assessed for a short time integration period of 12 hours and grid resolutions of    d  x  =  5  ,  10  ,  20 km and    d  t  =  20  ,  10  ,  5 minutes. While both methods showed an overall second-order convergence in space, only the BE method displayed the expected first-order convergence in timestep refinement. The CN's errors didn't decrease under time refinement, although they are consistently twice as small as the BE's errors. This behaviour persisted in longer 4-day runs, where the BE method displayed a more stable solution overall, especially at the coarse resolution with    d  t  =  80 minutes, where the CN failed to converge. The two methods are then compared in a 1-year climatic simulation with realistic wind and ocean forcings, using a grid resolution of 20 km where both methods performed fairly well. However, after 1 year the two solutions diverged somewhat significantly from one-another, providing hard evidence that climate models are sensitive to the underlying numerical methods.}
}