Improving the Jacobian free Newton-Krylov method for the viscous-plastic sea ice momentum equation

Sea ice plays a central role in regulating Earth's radiative budget because of its high albedo effect, and its melting level during the summer season is considered to be an important index of global warming. Contemporary earth system models (ESM) utilize complex dynamical models for the pack ice to account for the variations of sea ice cover and its feedback on the earth system. There is a wide consensus in the climate modeling community that the ice pack is most accurately modeled as a viscous-plastic flow with a highly nonlinear rheology consisting of an elliptic-yield-curve constitutive law. However, sea ice dynamics remains one of the most uncertain factors in the ESM's ability to address the climate change problem. The difficulty in accurately and efficiently solving numerically the associated highly nonlinear partial differential equations is believed to be a big contributor to this uncertainty. This work builds on recent efforts to construct fast and accurate numerical schemes for the viscous-plastic sea ice equations, based on a Jacobian free Newton-Krylov (JFNK) algorithm. Here, we propose to improve on the JFNK approach by using a fully second order Crank-Nicolson-type method to discretize the sea ice momentum equations (SIME) instead of the previously used first order backward Euler. More importantly, we improve on the Jacobian free approximation by expressing the derivatives of the least cumbersome and linear terms in the discretized SIME functional in closed form and use a second order Gateaux derivative approximation for the remaining terms, instead of using a first order approximation of the Gateaux derivative for the whole Jacobian matrix. Numerical tests performed on a synthetic exact solution for an augmented set of equations demonstrated that the new scheme is indeed second order accurate and the second order approximation of the Jacobian matrix was revealed to be crucial for the convergence of the nonlinear solver. One of the main difficulties in the JFNK approach resides in deciding on a stopping criterion for the Newton iterations. Our tests show that iterating beyond a certain level can in fact deteriorate the solution and prevent convergence. To overcome this issue, we suggest to use a conditional termination strategy by stopping the iterations as soon as the residual starts to increase. The resulting gain in efficiency overshadows any gains in accuracy when requiring formal convergence. (C) 2017 Elsevier B.V. All rights reserved.