This paper examines the effect of pore scale heterogeneity on the shape of breakthrough curves for a reactive solute through the application of a pore network model. The particular objective was to determine whether pore scale heterogeneity, in the absence of kinetic effects, could contribute to the increased dispersion that is characteristic of measured BTCs for reactive solutes. The pore network was generated stochastically, while flow and solute transport were simulated deterministically. A lognormal pore size distribution was assigned, and reactions were represented by a distribution coefficient that was defined on a unit surface area basis and was assumed to be constant everywhere in the model. Flow and transport within the network were simulated by finite element and particle tracking methods, respectively. The nonreactive solute exhibited almost ideal behavior, while the reactive solute showed a considerable degree of nonideality, with skewed and dispersed shapes. The results of a sensitivity analysis showed greater nonideality in reactive BTCs for networks with broader pore size distribution, shorter travel distance, and greater distribution coefficient. Unlike the predictions of kinetic models, BTCs were found to be independent of average pore water velocity. Although it is a highly idealized representation of a porous medium, the network model gave symmetrical BTCs for nonreactive solutes, lending credibility to its ability to represent solute transport in a qualitatively correct manner. Therefore, the fact that the model produced nonideal BTCs for reactive solutes, in the absence of kinetic process, provides considerable support for pore scale heterogeneity in retardation as a cause of nonideality. ¿ American Geophysical Union 1995