The fractional advection-reaction-diffusion equation plays a key role in describing the processes of multiple species transported by a fluid. Different numerical methods have been proposed for the case of fixed-order derivatives, while there are no such methods for the generalization of variable-order cases. In this paper, a numerical treatment is given to solve a variable-order model with time fractional derivative defined in the Atangana-Baleanu-Caputo sense. By using shifted Gegenbauer cardinal function, this approach is based on the application of spectral collocation method and operator matrices. Then the desired problem is transformed into solving a nonlinear system, which can greatly simplifies the solution process. Numerical experiments are presented to illustrate the effectiveness and accuracy of the proposed method.