The shallow ice approximation (SIA) model in strong form is commonly used for inferring the flow dynamics of grounded ice sheets. The solution to the SIA model is a closed-form expression for the velocity field. When that velocity field is used to advance the ice surface in time, the time steps have to take small values due to quadratic scaling in the horizontal mesh size. In this paper, we write the SIA model in weak form and add in the free-surface stabilization algorithm (FSSA) terms. We find numerically that the time-step restriction scaling is improved from quadratic to linear, but only for large horizontal mesh sizes. We then modify the weak form by adding the initially neglected normal stress terms. This allows for a linear time-step restriction across the whole range of horizontal mesh sizes, leading to improve efficiency. Theoretical analysis demonstrates that the inclusion of FSSA stabilization terms transitions the explicit time-stepping treatment of second derivative surface terms to an implicit approach. Moreover, a computational cost analysis, combined with numerical results on stability and accuracy, advocates for preferring the SIA models written in weak form over the standard SIA model.