Current sea ice models are based on a splitting in time between the momentum and continuity equations. As the ice strength is explicit when solving the momentum equation, this can lead to unrealistic ice stress gradients when using a large time step. As a consequence, noise develops in the numerical solution and these models can even become unstable at high resolution. To cure this problem, we have implemented an iterated IMplicit-EXplicit (IMEX) time integration method. This IMEX method was developed in the framework of an already implemented Jacobian-free Newton-Krylov solver. The basic idea of this iterated IMEX scheme is to move the explicit calculation of the sea ice thickness and concentration inside the Newton loop such that these fields evolve during the implicit integration. To obtain second-order accuracy in time, we have also modified the explicit time integration to a second-order Runge-Kutta approach and introduced a second-order backward difference method for the implicit integration of the momentum equation. These modifications to the code are minor and straightforward. A truncation error analysis and numerical experiments demonstrate that the approximate solution is second-order accurate in time. The new method permits to obtain the same accuracy as the splitting in time but by using a time step that is 10 times larger. Results show that the second-order scheme is more than five times more computationally efficient than the splitting in time approach for an equivalent level of error.