Imbibition is a commonly encountered multiphase problem in various fields, and exact prediction of imbibition processes is a key issue for better understanding capillary flow in heterogeneous porous media. In this work, a numerical framework for describing imbibition processes in porous media with material heterogeneity is proposed to track the moving wetting front with the help of a partially saturated region at the front vicinity. A new interface treatment, named the interface integral method, is developed here, combined with which the proposed numerical model provides a complete framework for imbibition problems. After validation of the current model with existing experimental results of one-dimensional imbibition, simulations on a series of two-dimensional cases are analysed with the presences of multiple porous phases. The simulations presented here not only demonstrate the suitability of the numerical framework on complex domains but also present its feasibility and potential for further engineering applications involving imbibition in heterogeneous media.