This paper deals with a two-level factored Crank-Nicolson method in an approximate solution of two-dimensional evolutionary advection-diffusion equation with time dependent dispersion coefficients and sink/source terms subjects to appropriate initial and boundary conditions. The procedure consists of reducing problems in many space variables into a sequence of one-dimensional subproblems and then find the solution of tridiagonal linear systems of equations. This considerably reduces the computational cost of the algorithm. Furthermore, the proposed approach is fast and efficient: unconditionally stable, temporal second order accurate and fourth order convergent in space and it improves a large class of numerical schemes widely studied in the literature for the considered problem. The stability of the new method is deeply analyzed using the $L^{\infty}(t_{0},T_{f};L^{2})$-norm whereas the convergence rate of the scheme is numerically obtained in the $L^{2}$-norm. A broad range of numerical experiments are presented and critically discussed.