Abstract: An operator-splitting algorithm for the two-dimensional convection-dispersion equation is developed. The governing equations are split into two successive initial value problems, which include a pure convection problem and a pure dispersion problem. For any arbitrary triangle calculation grids, suitable algorithms will be applied to resolve the problems at each grid point according operator characteristics. In the case of the convection problem, the method of characteristics can be applied. While for the dispersion problem, a semi-implicit finite element method will be employed. The issue concerning interpolation in convection problem will be discussed in detail, and a cubic interpolation method will be proposed to reduce the numerical damping effect. To overcome the numerical oscillation problem in high-order interpolation, it is required to keep the continuity of function and its partial derivatives. The testing results show that the numerical simulations of this study agree well with analytical solutions. The proposed operator-splitting algorithm can significantly reduce the numerical damping effect and overcome the numerical oscillation problem in high-order interpolation when solving pure convective transport problems.