With the increased ability to model subsurface heterogeneity provided by stochastic simulation methods, there is a corresponding increase in the need for numerical flow models to accommodate the spatial variability of hydrodynamic parameters. These trends make it necessary to reexamine methods for calculating finite difference interblock conductivities originally developed for use in homogeneous media. This study describes a new method for calculating interblock conductivities better suited to the numerical modeling of unsaturated flow in heterogeneous media. An expression for interblock conductivity is derived by equating the finite difference flux between adjacent blocks with the exact analytical flux for one-dimensional steady unsaturated flow, assuming an exponential model for relative conductivity. This expression appears as a function of the saturated conductivities, alpha parameters, and matric potentials of the adjacent blocks and the potential at the interblock boundary. Since this last quantity is not directly featured in the discretized flow equation, it must be evaluated iteratively. The interblock conductivity calculation scheme is implemented within a block-centered seven-point finite difference model of steady unsaturated flow where nonlinearity is resolved using a Newton-Raphson iterative approach. In a series of numerical experiments this scheme is compared with the more familiar geometric averaging method. The proposed interblock conductivity scheme is shown to provide perfect agreement between numerical and analytical solutions for the one-dimensional case of steady infiltration into a perfectly layered medium. In problems of two-dimensional unsaturated flow in heterogeneous media, comparisons with the geometric averaging method show that the proposed scheme allows more rapid convergence of Newton-Raphson iterations and consideration of greater parameter spatial variability before the onset of divergence in these iterations. ¿ American Geophysical Union 1995