Numerical solution of the Hibler model of sea-ice dynamics: Comparing backward Euler and Crank-Nicolson methods

Salim Bensassi; Boualem Khouider; Clint Seinen; M'hamed Kesri; Salim Bensassi; Boualem Khouider; Clint Seinen; M'hamed Kesri

doi:10.3934/math.2026329

AIMS Mathematics

2026, Volume 11, Issue 3: 7980-8013. doi: 10.3934/math.2026329

Previous Article Next Article

Research article

Numerical solution of the Hibler model of sea-ice dynamics: Comparing backward Euler and Crank-Nicolson methods

1.
Dynamical Systems Laboratory, Department of Analysis, Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 Bab Ezzouar, 16111 - Algiers, Algeria
2.
Mathematics and Statistics, University of Victoria, P.O. Box 1700 STN CSC, Victoria, Canada

Received: 30 September 2025 Revised: 03 March 2026 Accepted: 05 March 2026 Published: 25 March 2026
MSC : 35Q74, 35Q86, 65M06

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 $ dx = 5, 10, 20 $ km and $ dt = 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 $ dt = 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.
- viscous-plastic sea-ice equations,
- Jacobian-free Newton-Krylov,
- numerical solutions,
- Crank-Nicolson,
- backward Euler,
- numerical stability
Citation: Salim Bensassi, Boualem Khouider, Clint Seinen, M'hamed Kesri. Numerical solution of the Hibler model of sea-ice dynamics: Comparing backward Euler and Crank-Nicolson methods[J]. AIMS Mathematics, 2026, 11(3): 7980-8013. doi: 10.3934/math.2026329

Related Papers:

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 $ dx = 5, 10, 20 $ km and $ dt = 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 $ dt = 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.

References

[1]	M. S. Arif, K. Abodayeh, Y. Nawaz, Numerical modeling of mixed convective nanofluid flow with fractal stochastic heat and mass transfer using finite differences, Front. Energy Res., 12 (2024), 1373079. https://doi.org/10.3389/fenrg.2024.1373079 doi: 10.3389/fenrg.2024.1373079
[2]	M. S. Arif, K. Abodayeh, Y. Nawaz, A two-stage computational approach for stochastic Darcy-forchheimer non-Newtonian flows, Front. Phys., 13 (2025), 1533252. https://doi.org/10.3389/fphy.2025.1533252 doi: 10.3389/fphy.2025.1533252
[3]	W. J. Campbell, The wind-driven circulation of ice and water in a polar ocean, J. Geophys. Res., 70 (1965), 3279–3301. https://doi.org/10.1029/JZ070i014p03279 doi: 10.1029/JZ070i014p03279
[4]	T. F. Chan, K. R. Jackson, Nonlinearly preconditioned Krylov subspace methods for discrete Newton algorithms, SIAM J. Sci. Stat. Comput., 5 (1984), 533–542. https://doi.org/10.1137/0905039 doi: 10.1137/0905039
[5]	S. Chatta, B. Khouider, Well posedness of the regularized-hibler model of sea-ice dynamics, ESS Open Archive, 2025. https://doi.org/10.22541/essoar.175309004.43816040/v1
[6]	S. Chatta, B. Khouider, M. Kesri, Linear well posedness of regularized equations of sea-ice dynamics, J. Math. Phys., 64 (2023), 051504. https://doi.org/10.1063/5.0152991 doi: 10.1063/5.0152991
[7]	M. D. Coon, G. A. Maykut, R. S. Pritchard, D. A. Rothrock, A. S. Thorndike, Modeling the pack ice as an elastic–plastic material, AIDJEX Bulletin, 24 (1974), 1–105.
[8]	D. R. Durran, Numerical methods for fluid dynamics: With applications to geophysics, Springer, 32 (2010). https://doi.org/10.1007/978-1-4419-6412-0
[9]	D. L Feltham, Sea ice rheology, Annu. Rev. Fluid Mech., 40 (2008), 91–112. https://doi.org/10.1146/annurev.fluid.40.111406.102151
[10]	T. C. Fisher, M. H. Carpenter, J. Nordström, N. K. Yamaleev, C. Swanson, Discretely conservative finite-difference formulations for nonlinear conservation laws in split form: Theory and boundary conditions, J. Comput. Phys., 234 (2013), 353–375. https://doi.org/10.1016/j.jcp.2012.09.026 doi: 10.1016/j.jcp.2012.09.026
[11]	O. Guba, J. Lorenz, D. Sulsky, On well-posedness of the viscous–plastic sea ice model, J. Phys. Oceanogr., 43 (2013), 2185–2199. https://doi.org/10.1175/JPO-D-13-014.1 doi: 10.1175/JPO-D-13-014.1
[12]	W. D. Hibler, A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9 (1979), 815–846.
[13]	W. D. Hibler III, A viscous sea ice law as a stochastic average of plasticity, J. Geophys. Res., 82 (1977), 3932–3938. https://doi.org/10.1029/JC082i027p03932 doi: 10.1029/JC082i027p03932
[14]	E. C. Hunke, W. H. Lipscomb, A. K. Turner, N. Jeffery, S. Elliott, CICE: The Los Alamos sea ice model documentation and software user's manual, Version 5.1, Los Alamos National Laboratory, 2015.
[15]	E. C. Hunke, J. K. Dukowicz, An elastic-viscous-plastic model for sea ice dynamics, J. Phys. Oceanogr., 27 (1997), 1849–1867.
[16]	E. C. Hunke, Y. Zhang, A comparison of sea ice dynamics models at high resolution, Mon. Weather Rev., 127 (1999), 396–408.
[17]	E. C. Hunke, Viscous–plastic sea ice dynamics with the evp model: Linearization issues, J. Comput. Phys., 170 (2001), 18–38. https://doi.org/10.1006/jcph.2001.6710 doi: 10.1006/jcph.2001.6710
[18]	O. B. Johnson, L. I. Oluwaseun, Crank-nicolson and modified crank-nicolson scheme for one dimensional parabolic equation, Int. J. Appl. Math. Theor. Phys., 6 (2020), 35–40. https://doi.org/10.11648/j.ijamtp.20200603.11 doi: 10.11648/j.ijamtp.20200603.11
[19]	S. Juricke, T. Jung, Influence of stochastic sea ice parametrization on climate and the role of atmosphere–sea ice–ocean interaction, Philos. Trans. A Math. Phys. Eng. Sci., 372 (2014), 20130283. https://doi.org/10.1098/rsta.2013.0283 doi: 10.1098/rsta.2013.0283
[20]	J. F. Lemieux, B. Tremblay, Numerical convergence of viscous-plastic sea ice models, J. Geophys. Res. Oceans, 114 (2009), C05009. https://doi.org/10.1029/2008JC005017 doi: 10.1029/2008JC005017
[21]	J. F. Lemieux, B. Tremblay, J. Sedláček, P. Tupper, S. Thomas, D. Huard, et al., Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton–Krylov method, J. Comput. Phys., 229 (2010), 2840–2852. https://doi.org/10.1016/j.jcp.2009.12.011 doi: 10.1016/j.jcp.2009.12.011
[22]	J. F. Lemieux, D. A. Knoll, B. Tremblay, D. M. Holland, M. Losch, A comparison of the Jacobian-free Newton–Krylov method and the EVP model for solving the sea ice momentum equation with a viscous-plastic formulation: A serial algorithm study, J. Comput. Phys., 231 (2012), 5926–5944. https://doi.org/10.1016/j.jcp.2012.05.024 doi: 10.1016/j.jcp.2012.05.024
[23]	R. J. LeVeque, Finite volume methods for hyperbolic problems, Cambridge University Press, 2012. https://doi.org/10.1017/CBO9780511791253
[24]	M. Losch, S. Danilov, On solving the momentum equations of dynamic sea ice models with implicit solvers and the elastic–viscous–plastic technique, Ocean Model., 41 (2012), 42–52. https://doi.org/10.1016/j.ocemod.2011.10.002 doi: 10.1016/j.ocemod.2011.10.002
[25]	C. Mehlmann, S. Danilov, M. Losch, J. F. Lemieux, N. Hutter, T. Richter, et al., Simulating linear kinematic features in viscous-plastic sea ice models on quadrilateral and triangular grids with different variable staggering, J. Adv. Model. Earth Syst., 13 (2021), e2021MS002523. https://doi.org/10.1029/2021MS002523 doi: 10.1029/2021MS002523
[26]	C. Mehlmann, T. Richter, A modified global Newton solver for viscous-plastic sea ice models, Ocean Model., 116 (2017), 96–107. https://doi.org/10.1016/j.ocemod.2017.06.001 doi: 10.1016/j.ocemod.2017.06.001
[27]	Á. Nagy, I. Omle, H. Kareem, E. Kovács, I. F. Barna, G. Bognar, Stable, explicit, leapfrog-hopscotch algorithms for the diffusion equation, Computation, 9 (2021), 92. https://doi.org/10.3390/computation9080092 doi: 10.3390/computation9080092
[28]	H. Nessyahu, E. Tadmor, Non-oscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87 (1990), 408–463. https://doi.org/10.1016/0021-9991(90)90260-8 doi: 10.1016/0021-9991(90)90260-8
[29]	E. Ólason, G. Boutin, A. Korosov, P. Rampal, T. Williams, M. Kimmritz, et al., A new brittle rheology and numerical framework for large-scale sea-ice models, J. Adv. Model. Earth Syst., 14 (2022), e2021MS002685. https://doi.org/10.1029/2021MS002685. doi: 10.1029/2021MS002685
[30]	E. Ólason, G. Boutin, T. Williams, A. Korosov, H. Regan, J. Rheinlænder, et al., The next generation sea-ice model nextsim, version 2, EGUsphere, 2025. https://doi.org/10.5194/egusphere-2024-3521
[31]	C. Seinen, A fast and efficient solver for viscous-plastic sea ice dynamics, Master's thesis, Canada: The University of Victoria, 2017.
[32]	C. Seinen, B. Khouider, Improving the Jacobian free Newton–Krylov method for the viscous–plastic sea ice momentum equation, Phys. D, 376-377 (2018), 78–93. https://doi.org/10.1016/j.physd.2017.09.005 doi: 10.1016/j.physd.2017.09.005
[33]	M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys., 114 (1994), 146–159. https://doi.org/10.1006/jcph.1994.1155 doi: 10.1006/jcph.1994.1155
[34]	B. Tremblay, L. A. Mysak, Modeling sea ice as a granular material, including the dilatancy effect, J. Phys. Oceanogr., 27 (1997), 2342–2360.
[35]	M. Yaremchuk, G. Panteleev, On the jacobian approximation in sea ice models with viscous-plastic rheology, Ocean Model., 177 (2022), 102078.
[36]	J. Zhang, W. D. Hibler III, On an efficient numerical method for modeling sea ice dynamics, J. Geophys. Res. Oceans, 102 (1997), 8691–8702. https://doi.org/10.1029/96JC03744 doi: 10.1029/96JC03744

Reader Comments

Your name:*

Email:*
© 2026 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)