The simulation of Naturally Fractured Reservoirs (NFRs) has always been a challenging topic from the time of their discovery. One of the main problems in their simulation is calculating the matrix-fracture transfer which governs the dynamic behavior of the reservoir. This work is a thorough investigation of the transfer term. It will be shown that the conventionally practiced method to calculate the potential difference for displacement of one fluid by another is unreliable; however, when the saturation distribution is homogeneous (e.g. for single-phase expansion drive or solution gas dive), the conventional transfer equation can be reliably used. In 2004, Heinemann[44] suggested that recovery curves could be used to calculate the matrix-fracture transfer term. This research focuses on this idea contributing to its theoretical foundation and practical implementation. A workflow for using recovery curves for simulation of NFRs is (1) field investigations, (2) generating the recovery curves for the matrix blocks, (3) lumping (weighted-averaging of) matrix block recovery curves to generate the simulation cell recovery curves, (4) using the recovery curves in the mathematical model and finally, (5) utilizing recovery curves in the simulation model. New concepts are devised and different tools and methods that were necessary are designed and developed which include: “Single Matrix Block Analysis” method and software tools are developed that allow detailed study of matrix blocks with different fluid and rock properties under various initial and boundary conditions. Each investigation produces a recovery curve and average saturation versus time plots. “Matrix Block Classes” are defined that allow scaling the recovery curve of matrix blocks with different porosities, permeabilities, and even different shapes by use of a theoretically-derived “dimensionless time” scaling factor. “Recovery Curve Regions” are defined on the full field. For each such region and distinct drive mechanism a “lumped recovery curves” represents the overall recovery from all the different matrix blocks in the simulation cell. “Using Recovery Curves instead of the Matrix-Fracture Transfer Function”. Several challenges existed for each of these steps. They were successfully overcome and the developed methodology was implemented in an industrial reservoir simulator. The new concept can be easily implemented in any mature dual-porosity reservoir simulator by making only minor extensions to it. Thus they can also benefit from the additional accuracy of the recovery curve method without reducing the industrial practicability of the simulators.