<p>This paper aims to present the spectral collocation method-based approximation for the numerical simulation of the variable-order space-time fractional advection–diffusion equation in the two-dimensional domain. We employ the shifted Chebyshev polynomial as an orthogonal polynomial, and fractional order is defined according to Caputo’s definition. Error and convergence studies of the present approach are also provided. We validate and examine the effectiveness of the proposed method by applying it to a number of numerical situations, some of which involve a non-smooth solution. We discuss an application of atmospheric pollution distribution with fractional derivatives to illustrate the significance of variable order over constant order in such models.</p>