The development of a computationally efficient high-resolution viscous-plastic sea ice model

This thesis presents the development of a high-resolution viscous-plastic (VP) sea ice model. Because of the fine mesh and the size of the domain, an efficient and parallelizable numerical scheme is desirable. In a first step, we have implemented the nonlinear solver used in existing VP models (referred to as the standard solver). It is based on a linear solver and an outer loop (OL) iteration. For the linear solver, we introduced the preconditioned Generalized Minimum RESidual (pGMRES) method. The preconditioner is a line successive overrelaxation solver (SOR). When compared to the SOR and the line SOR (LSOR) methods, two solvers commonly used in the sea ice modeling community, pGMRES increases the computational efficiency by a factor of 16 and 3 respectively. For pGMRES, the symmetry of the system matrix is not a prerequisite. The Coriolis term and the off-diagonal part of the water drag can then be treated implicitly. Theoretical and simulation results show that this implicit treatment eliminates a numerical instability present with an explicit treatment. During this research, we have also observed that the approximate nonlinear solution converges slowly with the number of OL iterations. Furthermore, simulation results reveal: the existence of multiple solutions and occasional convergence failures of the nonlinear solver. For a time step comparable to the forcing time scale, a few OL iterations lead to errors in the velocity field that are of the same order of magnitude as the mean drift. The slow convergence is an issue at all spatial resolutions but is more severe as the grid is refined. It is attributed in part to the standard VP formulation that leads to a momentum equation that is not continuously differentiable. To obtain a smooth formulation, we replaced the standard viscous coefficient expression with capping by a hyperbolic tangent function. This provides a unique solution and reduces the computational time and failure rate. To further improve the computational efficiency, the Jacobian-free Newton-Krylov (JFNK) method was implemented. Results show that JFNK is 3.5-6.6 times faster than the standard solver (depending on the resolution and termination criterion). As the standard solver, the JFNK method exhibits a relatively small number of convergence failures. A globalization "trick", which involves a switch in the linearization approach, eliminates the failures for the standard solver. Globalization approaches such as the line search did not improve the failure rate of the JFNK method. Future work will focus on the implementation of an effective globalization approach for JFNK.