A two-level factored implicit scheme is considered for solving a two-dimensional unsteady advection-dispersion equation with spatio-temporal coefficients and source terms subjected to suitable initial and boundary conditions. The approach reduces multi-dimensional problems into pieces of one-dimensional subproblems and then solves tridiagonal systems of linear equations. The computational cost of the algorithm becomes cheaper and makes the method more attractive. Furthermore, the two-level approach is unconditionally stable, temporal second-order accurate and spatial fourth-order convergent. The developed numerical scheme is faster and more efficient than a broad range of methods widely studied in the literature for the considered initial-boundary value problem. The stability of the proposed procedure is analyzed in the $L^{\mathrm{\infty }}(t_0,T_f;L^2)$-norm whereas the convergence rate of the algorithm is numerically analyzed using the $L^2(t_0,T_f;L^2)$-norm. Numerical examples are provided to verify the theoretical result.