  We present a Newton-Krylov solver for a viscous-plastic sea-ice model. This
constitutive relation is commonly used in climate models to describe the
material properties of sea ice. Due to the strong nonlinearity introduced by
the material law in the momentum equation, the development of fast, robust and
scalable solvers is still a substantial challenge. In this paper, we propose a
novel primal-dual Newton linearization for the implicitly-in-time discretized
momentum equation. Compared to existing methods, it converges faster and more
robustly with respect to mesh refinement, and thus enables numerically
converged sea-ice simulations at high resolutions. Combined with an algebraic
multigrid-preconditioned Krylov method for the linearized systems, which
contain strongly varying coefficients, the resulting solver scales well and can
be used in parallel. We present experiments for two challenging test problems
and study solver performance for problems with up to 8.4 million spatial
unknowns.
