In this paper, we firstly present a novel simple method based on a Picard integral type formulation for the nonlinear multi-dimensional variable coefficient fourth-order advection-dispersion equation with the time fractional derivative order $\alpha\in (1,2)$. A new unknown function $v(\mathbf{x},t)=\partial u(\mathbf{x},t)/\partial t$ is introduced and $u(\mathbf{x},t)$ is recovered using the trapezoidal formula. As a result of the variable $v(\mathbf{x},t)$ is introduced in each time step, the constraints of traditional plans considering the non-integer time situation of $u(\mathbf{x},t)$ are no longer considered. The stability and solvability are proved with  detailed  proofs and the precise description of error estimates is derived. Further, Chebyshev spectral collocation method supports accurate and efficient variable coefficient model with variable coefficients. Several numerical results are obtained and analyzed in multi-dimensional spatial domains and numerical convergence order is consistent with the theoretical value $3-\alpha$ order for different $\alpha$ under infinite norm.