Multi-phase flow in carbonate reservoirs, which hold about half of the world's remaining oil reserves, is strongly influenced by fractures present in the geological formations. Fractures are often the main flow conduits, leaving most of the oil behind in the low permeability rock matrix, and cause early water breakthrough. An accurate characterization of fracture flow and fluid exchange between fracture and matrix is needed to forecast oil recovery and optimise production in fracture-dominated and fracture-assisted reservoirs. Dual-porosity simulations are traditionally used to model naturally fractured reservoirs. However, the classical dual-porosity models miss some key physics of fracture - matrix fluid exchange: (1) They tend to simplify mass transfer due to spontaneous imbibition since they typically cannot include arbitrary petro-physical properties due to e.g. wettabilities, viscosity ratios, initial water content, etc. (2) They cannot account for the fact that matrix blocks within a single simulation grid cell have various sizes and permeabilities, giving rise to different transfer rates, which should be captured by a distribution of transfer functions in each grid cell. In this paper we present a novel multi-rate dual-porosity model. It is based on an unstructured finite element - finite volume technique, which solves the governing equations for two-phase flow fully implicitly. This allows us to represent complex large-scale geological structures (e.g., non-orthogonal faults and fracture corridors) accurately while small-scale diffuse fractures are modeled with an improved dual-porosity approach. Until recently, calculating the mass transfer between fractures and matrix blocks due to spontaneous imbibition for arbitrary petro-physical and fluid properties exactly was not possible because a general and exact transfer rate for arbitrary petro-physical and fluid properties was not known. In our new dual-porosity model we compute fracture-matrix transfer using the only known analytical and general solution of the Darcy equation including capillarity. This provides us with a generalised transfer function for arbitrary wettability, viscosity ratios, rock types, initial water content and boundary conditions. Each reservoir simulation grid block can contain multiple of these generalised transfer functions to model different matrix permeabilities and/or matrix block sizes present at the sub-grid scale. Using a series of proof-of-concept simulations, we show that the difference in oil recovery using a standard single-rate dual-porosity model and our new multi-rate dual-porosity model cannot be neglected. This demonstrates that our proposed model with the generalised transfer function predicts oil recovery more accurately compared to a classical dual-porosity model.