An improved approach for handling boundaries, interfaces, and continuous depth dependence with the elastic parabolic equation is derived and benchmarked. The approach is applied to model the propagation of Rayleigh and Stoneley waves. Depending on the choice of dependent variables, the operator in the elastic wave equation may not factor or the treatment of interfaces may be difficult. These problems are resolved by using a formulation in terms of the vertical displacement and the range derivative of the horizontal displacement. These quantities are continuous across horizontal interfaces, which permits the use of Galerkin's method to discretize in depth. This implementation extends the capability of the elastic parabolic equation to handle arbitrary depth dependence and should lead to improvements for range-dependent problems.