Naturally fractured reservoirs can be seen as a set of low permeability matrix rock blocks and a high permeability network of fracture channels. This representation is conventionally described by a dual porosity model, which is the one used in the present work. More precisely, the porosities and absolute permeabilities at each point of a reservoir are considered to be those of two interpenetrating continua, the matrix and the fractures one. It is also assumed that the flow occurs in fractures only, i.e., the matrix permeability is equal to zero. Mass transfer between matrix and fractures is modeled by empirically determined transfer functions. Fracture permeabilities can differ in orders of magnitude, which results in very different flow velocities in different parts of the reservoir. This circumstance is advantageous for simulating the reservoir numerically with the streamline method [1,4,6,7,8,9,11,17,18,21]. Indeed, in the streamline method the transport part is solved along a set of one-dimensional streamlines, and the corresponding time step used for the update of the solution can be chosen independently for each streamline. Therefore, there is no global time step restriction, which is the case for finite difference simulation. This work is a further development of methods, proposed in [19,16] and . We are concerned with the efficient solution of saturation transport along streamlines. To this end, we propose to use an implicit method. This eliminates the CFL restriction on the time step which is inherent to explicit methods, at the price of possibly introducing more numerical diffusion. The numerical method can be further accelerated with the use of adaptive mesh refinement along streamlines . We show that the proposed technique can be useful when one needs to estimate the flow pattern in the reservoir quickly, e.g., for ranking of geological models.