Existing continuous schemes for computing the specific discharge in heterogeneous porous media result in inaccuracy in simulating the local specific discharge field at and in the vicinity of an interface between subdomains of different material properties. The domain decomposition method used in multilayered systems and a consistent scheme in a randomly heterogeneous domain are applicable only to some particular cases, such as when the material interfaces are parallel to an axis. Along a material interface, the law of refraction of streamlines, related to both the tangential and the normal specific discharge components, should be satisfied. By incorporating these interface conditions into the weak formulation of Darcy's law, two consistent numerical schemes have been developed. These schemes can be applied to accurately and efficiently simulate the specific discharge field in a heterogeneous domain with different values of the ratio of hydraulic conductivities between adjacent elements. In particular, they provide an easy way for treating aquifers with semipervious layers as a single system. Numerical simulations have been used to demonstrate the capability and accuracy of the proposed schemes. Comparison between the specific discharge fields computed by the proposed schemes and by other techniques shows that the new schemes are superior to the traditional ones. ¿ 2001 American Geophysical Union