Abstract The accurate representation of the interactions between matrix and fracture in the dual porosity model is crucial to the proper prediction of recovery from fractured petroleum reservoirs. In the standard dual porosity model, matrix-fracture transfer is represented by a single term in the mass balance equations, and the matrix-matrix flow is assumed to be negligible. This was found to be inadequate when gravity effects are important. Various methods were reported in the recent literature which account for the gravity drainage process. These include the gravity segregated model, the sub-domain method, and the dual permeability model. In this paper, we will examine the gas-oil gravity drainage process and evaluate the different approaches for solving this problem. Introduction The dual porosity approach has become the most widely accepted technique for the simulation of naturally fractured reservoirs. In this approach, the reservoir is represented by two collocated continua. They are referred to as fracture and matrix. The fracture continuum represents the interconnected network of fractures and/or vugs which constitute the primary conduits for fluid flow. The matrix continuum represents the pore volume in the rock which comprises the majority of the storage in the reservoir. Because the permeability in the matrix is substantially lower than that in the fracture, earlier models usually assume that the matrix does not contribute to the global mass transfer and act primarily as source and sink terms to the fracture. An early account of a full-field study of a model of this type can be found in Kazemi et al.(1). The fractured continuum is discretized into computational grid cells. The matrix continuum is assumed to be comprised of matrix blocks which are separated spatially by fractures. Strictly speaking, these matrix blocks can be highly irregular, and their sizes are not constant throughout the reservoir. To make the analysis manageable, it is typically assumed that the matrix blocks can be represented by regular geometries such as rectangular prisms(1,2) or cylinders(3). The dimensions of these matrix blocks can be variable throughout the reservoir and are primarily a function of the fracture spacing, orientation, and width. One computational grid cell may contain several matrix blocks. It is usually not feasible to insist on having one matrix block per computational cell. Because the matrix block is the primary unit for matrix-fracture exchange calculation, it is usually assumed that all the matrix blocks within one computational cell are geometrically and dynamically similar. Hence, the mass transfer between the matrix and the fracture in one computational grid is just the mass transfer between one matrix block and its surrounding fractures multiplied by the number of matrix blocks in the fracture cell. The adequacy of this assumption will be examined in later sections. The matrix-fracture transfer is usually represented by a single term with a geometric factor in the mass balance equations. The mass balance finite-difference equations have been described in several papers on the subject and will not be repeated here (see, for example, Reference 6).