The viscous-plastic momentum equations for sea-ice dynamics are a model for the time evolution of the horizontal velocity field of a vertically homogeneous column of ice-pack continuum whose thickness and ice concentration vary in time and space. They constitute a highly nonlinear system of partial differential equations, and as such, both the mathematical analysis and numerical solution of its solution remain a challenge. To shed some light on this problem, we compared the performance of two known time discretization methods: the backward Euler (BE) and the Crank-Nicolson (CN) methods. Both methods are in theory unconditionally stable, but Euler's method is only first-order accurate while Crank-Nicolson is second-order. Centered finite differences were used for the spatial derivatives. This led to a nonlinear system of algebraic equations which was then solved using a Jacobian-free Newton-Krylov approach. First, the two methods were compared in terms of their ability to reproduce a synthetic solution to the viscous-plastic momentum equations when an artificial forcing is applied. The convergence of the two methods is assessed for a short time integration period of 12 hours and grid resolutions of $dx=5,10,20$ km and $dt=20,10,5$ minutes. While both methods showed an overall second-order convergence in space, only the BE method displayed the expected first-order convergence in timestep refinement. The CN's errors didn't decrease under time refinement, although they are consistently twice as small as the BE's errors. This behaviour persisted in longer 4-day runs, where the BE method displayed a more stable solution overall, especially at the coarse resolution with $dt=80$ minutes, where the CN failed to converge. The two methods are then compared in a 1-year climatic simulation with realistic wind and ocean forcings, using a grid resolution of 20 km where both methods performed fairly well. However, after 1 year the two solutions diverged somewhat significantly from one-another, providing hard evidence that climate models are sensitive to the underlying numerical methods.