Selection of a method for the numerical solution to the Stefan problem is considered, which permits one to calculate the ice boundary with accepted accuracy during a long period of time, and can be extended to multidimensional glaciation problems. A comparison of the different versions of continuous schemes with an exact analytical solution to the classical Stefan problem is presented. The criterion that permits estimation, under given assumptions, of the error behavior of calculating the phase boundary in numerical solution to the Stefan problem using a continuous calculation scheme is suggested. The advantage of this continuous calculation scheme for phase boundary calculation in problems of surface glaciation in seawater is demonstrated.