In this paper, we consider the discrete Hahn polynomials fH n g and investigate their application for numerical solutions of the time fractional variable order mobile-immobile advection-dispersion model which is advantageous for modeling dynamical systems. This paper presented the operational matrix of derivative of discrete Hahn polynomials. The main advantage of approximating a continuous function by Hahn polynomials is that they have a spectral accuracy in the interval [0, N]. Furthermore, for computing the coefficients of the expansion uðxÞ ¼ P 1 n¼0 c n H n ðxÞ, we have to only compute a summation and the calculation of coefficients is exact. Also an upper bound for the error of the presented method, with equidistant nodes, is investigated. Illustrative examples are provided to show the accuracy and efficiency of the presented method. Using a small number of Hahn polynomials, significant results are achieved which are compared to other methods.