Accurate modeling of sea ice dynamics is critical for predicting environmental variables, which in turn is important in applications such as navigating ice breaker ships, and has led to extensive research in both modeling and simulating sea ice dynamics. The most widely accepted model is the one based on the viscous-plastic formulation introduced by Hibler, which is intrinsically difficult to solve numerically due to highly nonlinear features. In particular, sea ice simulations often significantly differ from satellite observations. In this study we focus on improving the numerical accuracy of the viscous-plastic sea ice model. We explore the convergence properties for various numerical solutions of the sea ice model and in particular examine the poor convergence seen in existing numerical methods. To address these issues, we demonstrate that using higher order methods for solving conservation laws, such as the weighted essentially non-oscillatory (WENO) schemes, is critical for numerically solving viscous-plastic formulations whenever the solution is not smooth. Moreover, WENO yields higher order convergence for smooth solutions than standard central differencing does. Our numerical examples verify this, and in particular by using WENO, we are able to resolve the discontinuities in the sharp features of sea ice covers. We also propose an approach utilizing the idea of phase field method to develop a potential function method which naturally incorporates the physical restrictions of ice thickness and ice concentration in transport equations. Our approach results in modified transport equations with extra forcing terms coming from potential energy function, and has the advantage of not requiring any post-processing procedure that might introduce discontinuities and thus ruin the solution behavior.