A methodology for incorporating subgrid effects in coarse-scale numerical models of flow in heterogeneous porous media is presented. The method proceeds by upscaling the deterministic fine-grid permeability description and then solving the pressure equation over the coarse grid to obtain coarse-scale velocities. The coarse-grid saturation equation is formed through a volume average of the fine-scale equations and includes terms involving both the average and fluctuating components of the velocity field. The terms involving the fluctuating components are subgrid effects that appear as length- and time-dependent dispersivities. A simplified model for the coarse-scale dispersivity, in terms of these subgrid velocity fluctuations, is proposed, and a numerical scheme based on this model is implemented. Results using the new method are presented for a variety of two-dimensional heterogeneous systems characterized by moderate permeability correlation length in the dominant (horizontal) flow direction and small correlation length in the vertical direction. The new method is shown to provide much more accurate results than comparable coarse-grid models that do not contain the subgrid treatment. Extensions of the overall methodology to handle more general systems are discussed. ¿ 2000 American Geophysical Union