Advection–diffusion equations play a crucial role in modeling the transport of heat, mass, pollutants, and contaminants in porous and heterogeneous media. These models arise in many practical applications, including environmental engineering, hydrology, geophysics, and material sciences. However, classical integer-order advection–diffusion models are often unable to accurately describe anomalous transport processes observed in complex systems. Such systems frequently exhibit long-range interactions, nonlocal dynamics, and memory effects. As a result, the transport behavior may deviate from classical Fickian diffusion and display sub-diffusive or super-diffusive characteristics that cannot be captured using traditional formulations. To address these limitations, this study investigates a variable-order fractional advection–diffusion integro-differential equation that incorporates fractional derivatives in both time and space under appropriate initial and boundary conditions. Unlike many existing studies that primarily focus on constant-order fractional models, the present formulation allows the fractional order to vary with time and space. This feature provides greater flexibility in modeling evolving memory effects and spatial heterogeneity of the medium. A finite difference numerical scheme is developed to approximate the variable-order fractional derivatives together with the convolution-type integral term. The proposed method achieves an accuracy of order (2−α ) on a nonuniform temporal mesh and reduces to first-order accuracy on a uniform mesh. Furthermore, the stability of the numerical scheme is established through theoretical analysis, ensuring the reliability of the computational results. To illustrate the effectiveness of the proposed approach, numerical simulations are carried out for three representative cases: constant fractional order (α=constant), time-dependent order α(t), and fully variable order α(x,t). The numerical results demonstrate good agreement with the analytical solutions and confirm that the proposed method effectively captures complex anomalous transport dynamics with high flexibility and accuracy.