Oil production, enhanced oil recovery, waste disposal, and CO2 storage applications in naturally fractured oil, gas, coalbed methane, geothermal reservoirs, and aquifers are characteristically controlled by an interaction between matrix and fracture. The correct estimation of the relative permeabilities for matrix-fracture interaction is essential in the performance analysis of such reservoirs. Conventional relative permeability measurement techniques are not suitable for this type of processes as the driving force is capillary rather than viscous. An alternative to these techniques is a pore scale modeling of the process. In this study, matrix-fracture interaction by capillary (spontaneous) imbibition was numerically simulated using the Lattice Boltzmann Method (LBM). The classical LBM algorithm was modified to add the effects of capillary characteristics such as wettability and interfacial tension. The model was validated using experiments on two - dimensional sand pack models, where the strongly water-wet model saturated with oil was exposed to water to displace oil by capillary interactions. Further, LBM simulations were used to investigate the critical parameters that have impacts on relative permeabilities such as different wettabilities, matrix boundary conditions that cause co-current interaction, and gravity (vertical and horizontal interaction). Finally, the LBM results were used to generate relative permeability curves by incorporating the algorithms based on single – phase normalization techniques. The LBM images were used to quantify the saturation values. The effects of different parameters listed above on the end points and the shape of the relative permeability curves as well as the residual oil saturation were identified. The results and observations provide qualitative and quantitative data which can be used in modeling studies for naturally fractured oil, gas, and coalbed methane reservoirs.