When modeling transport of chemicals or solute in realistic large-scale subsurface systems, numerical issues are a serious concern, even with the continual progress made over the past few decades in both simulation algorithms and computer hardware. The problem becomes even more difficult when dealing with chemical transport in a vadose zone or multiphase flow system using coarse, multidimensional regular or irregular grids because of the known effects of numerical dispersion associated with moving plume fronts. We have investigated several total variation diminishing (TVD) or flux-limiter schemes by implementing and testing them in the T2R3D code, one of the TOUGH2 family of codes. The objectives of this paper are (i) to investigate the possibility of applying these TVD schemes, using multidimensional irregular unstructured grids, and (ii) to help select more accurate spatial averaging methods for simulating chemical transport, given a numerical grid or spatial discretization. We present an application example to show that such TVD schemes can effectively reduce numerical dispersion.