A computationally efficient numerical method is developed for solving sea ice momentum equations that employ a nonlinear viscous-plastic rheology. The method is based on an alternating direction implicit (ADI) technique that involves a direct solution of the momentum equations. This method is therefore more computationally efficient than those employing an iterative procedure in solving the equations. The ADI method for modeling sea ice dynamics is dynamically consistent since it rapidly approaches a viscous-plastic solution described by the sea ice rheology. With different model configurations of varying spatial resolutions and decreasing time step intervals the ADI method converges to the same viscous-plastic solution as another numerical method that uses a line successive relaxation procedure to solve the ice momentum equations. This indicates that the ADI method is also numerically consistent. The approximateness of numerical solutions of sea ice, resulting from coarse model resolutions in time, is addressed. It is found that a significant bias, up to 10% or more, in the solution is likely to occur for a typical but coarse time step interval. This indicates that an assessment of the numerically created bias from a crude time integration may be necessary when model data comparisons are performed. In addition, suggestions are given for selecting appropriate time step intervals to enhance numerical accuracy in model applications. ¿ 2000 American Geophysical Union