A numerical simulation model is developed for characterizing three dimensional multiphase fluid flow phenomena in naturally fractured reservoirs. This new model treats the reservoir as a double porosity medium consisting of heterogeneous, isotropic primary rock matrix, and an anisotropic, heterogeneous fracture system. The statistical distribution in space and orientation of the fracture system is manifest in terms of a fracture permeability tensor. Such an approach provides a model which permits flow through the rock matrix and the fractures to be described simultaneously. A fluid interaction term is utilized to couple the flow between the rock matrix and the fractures. Gravity effects, capillary forces, and the fluid velocity in the fractures are also taken into account in the simulator model.