In this work, we present an efficient implementation of Eulerian Total Variation Diminishing methods for simulating solute transport in heterogeneous porous media using GPU power. The study focuses on solving the advection-diffusion-dispersion equation in 2D and 3D domains with high-order temporal schemes optimized for GPUs in C++/CUDA. A parallel GPU-based strategy is proposed for generating lognormal permeability fields to enable large-scale Monte Carlo simulations. The methods are applied to compute longitudinal and transverse macrodispersion coefficients under various transport scenarios, including pure advection, advection-diffusion, and advection-dispersion. By averaging over 100 permeability realisations , we explore the dimensional effects on dispersion behavior. Previous studies have shown that transverse macrodispersion vanishes in 2D but increases steeply with heterogeneity in 3D, linked to flow-line braiding and velocity correlations. These results are particularly relevant for media with non-smooth or locally anisotropic conductivity, challenging assumptions of isotropic hydraulic conductivity in smooth media, which may not capture key transport mechanisms. Our approach demonstrates efficient computation in domains of up to 134.5 million cells on a single GPU, with a detailed comparison of explicit and implicit methods. These results emphasise the capability of fully Eulerian GPU-based algorithms to achieve high-resolution macrodispersion analyses, offering valuable insights into solute transport mechanisms in highly heterogeneous media and paving the way for further applications in complex geological systems.