A computationally efficient numerical method for the solution of nonlinear sea ice dynamics models employing viscous-plastic rheologies is presented. The method is based on a semi-implicit decoupling of the x and y ice momentum equations into a form having better convergence properties than the coupled equations. While this decoupled form also speeds up solutions employing point relaxation methods, a line successive overrelaxation technique combined with a tridiagonal matrix solver procedure was found to converge particularly rapidly. The procedure is also applicable to the ice dynamics equations in orthogonal curvilinear coordinates which are given in explicit form for the special case of spherical coordinates. ¿ 1997 American Geophysical Union