Efficiency constraints in applying numerical models to real problems force the use of a coarse discretization of the domain compared with the detailed scale required for the most adequate description of the system. Upscaling encompasses the methods that transfer small-scale information to the computational scale. The loss of small-scale information by upscaling aquifer properties to construct a numerical model largely modifies the true heterogeneous structure of the aquifer compromising the final predictions of solute transport. We present extensive Monte Carlo solute transport simulations in heterogeneous porous media to investigate (1) the impact of upscaling on the evolution of solute plumes and (2) the propagation of uncertainty through upscaling. In doing this we analyze the benefits of using enhanced block dispersion tensors in the advection-dispersion equation to compensate for the loss of information. It is shown that upscaling does not affect the median traveltimes of particles reaching a control location, but the tail of breakthrough curves is largely underestimated. It is also observed that normal upscaled transport models can largely underestimate the spreading of solute plumes even if block dispersivities are calculated as being representative of within-block heterogeneity. Discrepancies are mainly attributed to the mass transfer interaction between grid blocks of the numerical model and the use of a simplistic representation of the block hydraulic conductivity tensor. Using the small-perturbation stochastic approach, we found that mass transfer effects are still significant, even for block sizes of the numerical model greater than 30 correlation scales. Moreover, upscaling is shown to have a substantial impact on the uncertainty predicted by solute transport models, which is not compensated by the use of enhanced block dispersion tensors. For very heterogeneous media, simulation results showed the inability of the classical advection-dispersion equation to simultaneously reproduce the tail of the breakthrough curve and its uncertainty in the upscaled model. In those cases, results indicated that either a mass transfer model allowing solute mass exchange between high- and low-conductivity zones, encountered within a grid block of the numerical model, or the continuous time random walk model expressing the retardation of solute particles passing through low-conductivity regions inside a grid block by means of a transition time distribution function might suppose a better upscaled transport model.