Transport in the subsurface is strongly characterized by the heterogeneity of the porous medium. As the subsurface is strongly influenced by human activity, obtaining an accurate description of mass transport in such media remains a crucial task. Tracer transport in porous media has been extensively treated in the literature. However, considering the many pollutants in industry and agriculture that enter the subsurface further research is ongoing. In addition, for many pollutants (nuclear waste, microplastics, pharmaceuticals, etc.), contamination remains critical at very low concentrations and cannot be ignored. Here, heavy tailing is a major issue particularly in water remediation. Further understanding of the coupling between transport, diffusion, and exchange kinetics (due to stagnant zones or adsorption), including tailing effects, remains a crucial task in the actual environmental context. Dynamics of spreading and the temporal spatial extension of the pollutants but also the time to reach a critical location (e.g. aquifer, well) strongly depend on the underlying heterogeneity of the permeability field. Moreover, in the case of solute retention, transport is also strongly impacted by local exchange kinetics that depend on the local aquifer properties. Consequently, exchange (retention) times are expected to be spatially heterogeneous. In this work, we simulate transport in a two-dimensional heterogeneous medium under spatially varying permeability and mobile-immobile mass transfer parameters. Equations are solved using a Lattice-Boltzmann TRT algorithm. We assume the following relation between the local permeability K and the local retention time 1/τ: 1/τ ≈K^γ. Taking into account this relation, we investigate the impact of the Damköhler number (Da, ratio of the retention and convection time scales), the disorder of the permeability field and the value of the exponent of the coupling function (γ) on the spatial evolution of the concentration field and the breakthrough curve. We show that, depending on the parameters (Da, γ, etc.), we can observe normal or anomalous dispersion, which are characterized by power-law tails of solute breakthrough curves and non-linear evolution of the spatial variance of the solute distribution. These behaviors are upscaled using a continuous time random walk approach based on a spatial Markov model for particle velocities that couples advective-dispersive transport and heterogeneous mass transfer through a compound Poisson process. The model can be parameterized by the distributions of permeabilities and exchange parameters. Coupling between the heterogeneous permeability field and the local mass transfer properties can strongly influence transport and explain the experimentally observed non-Gaussian behavior.