Abstract This paper provides a review of common approaches for simulation of naturally fractured reservoirs and a new model formulation that is more amenable to the utilization of detailed geologic information from deterministic models, multipoint statistical simulations (MPS), and discrete fracture network (DFN) models. Unlike the common dual-porosity and dual-permeability models, the new model considers flow in several sets of fractures (that is, micro, macro, and mega fracture levels) and the wells can intercept any class of fractures. Matrix flow is also included in the formulation. A dual-mesh-computing algorithm is used to capture the major orientation of fractures in the flow network. The algorithm consists of the conventional five-point discretization for the coarse grid and a special nine-point discretization scheme for the fine grid. The size of the coefficient matrix for the discretization scheme can be reduced because of the dominance of vertical flow drainage in fractured reservoirs.