A computational recovery of multi-phase intrusion was discussed with the modified multi-relaxation time Lattice Boltzmann method (MRT-LBM). Originally proposed dual-matrix computation is developed to address the different phase separation and interface tracking for the multi-phase problem. A comprehensive validation is performed with the previously theorized observation of the mercury-water system. Results show that the dual-matrix computation is feasible to provide converged output under narrowed density difference down to 18%. The wetting and non-wetting behaviour resulted from form solid–fluid interaction is realized with arbitrary boundaries, in which the contact variance is up to 4.14%. The linear relations described by Laplace's law and Washburn's equation were three-dimensionally recovered with determination coefficients of 96.34% and 94.19%, respectively. A third fluid intrusion status of partial-intrusion is captured in addition to complete-intrusion and non-occupation in porous boundary, demonstrating the advanced function of the phase-separation and interface tracking in problems with further increased heterogeneity.