Accurate numerical simulation of solute transport or multiphase flow through heterogeneous porous media is still a challenge. One reason is that in continuum models, sharp material interfaces are difficult to represent. The finite-element method allows representation of such discontinuities as element-to-element variations of materials; however, result variables are discretized on the finite-element nodes and when a complementary node-centered finite-volume mesh is used to model transport, this discretization smears out potential discontinuities in nodal concentration or saturation. To overcome this dilemma there are two options: either the nodes can be enriched at the interfaces by additional degrees of freedom or the model can be “exploded” along the interfaces. In the latter technique, interface nodes are multiplied so that they match in number the material domains that they join. We develop this latter approach, implementing a transport model that can evolve dependent variable discontinuities at material interfaces. For transport, we developed a new scheme that is higher order accurate in space, even across these discontinuities, and verify that it is consistent with the continuous method. This analysis shows that the new scheme removes a first-order dependence on mesh refinement observed for the continuous method so that coarser meshes can be used for material interfaces.