A parallel, semi-implicit numerical method is proposed for the study of naturally-fractured reservoir systems. It has proved to be computationally efficient in producing accurate numerical solutions for the dual-porosity model for immiscible, two-phase flow in such reservoirs. The method combines hybridized mixed finite elements, a new version of the modified method of characteristics, a sophisticated operator-splitting procedure for separating the pressure calculation in the fractures from that of the saturation, another operator splitting to handle the interaction of the matrix blocks and the fractures, and domain decomposition iterative procedures for both the pressure and the saturation in the fractures. It is easy to implement and permits moderately long time steps for the pressure and the saturation in the fractures by using short, inexpensive microsteps for the fracture saturation to treat the transport portion of the saturation equation. This talk is devoted to the formulation of the method and a discussion of numerical results for five-spot and vertical cross-section examples.