## 1 Introduction  
The simultaneous flow of two immiscible fluids in porous media is often found in nature and industrial processes such as non-aqueous phase liquid (NAPL) contamination of groundwater, carbon dioxide injection and storage, and enhanced oil recovery (EOR) operations. For the spontaneous imbibition of oil, carbon dioxide or water, which is immiscible to the oil phase, is usually imbibed into permeable rock formations which are characterized by low porosity and low in situ permeability. The displacement process is affected by matrix properties such as the heterogeneity, wettability, and fluid properties of each phase [1]. Experimental efforts have been made to understand interactions between fluids and rock during fluid displacement [2, 3]. Theoretical works usually use simplified pore geometries which allow theoretical solutions to the microscopic flow patterns [4]. As a complementary method to the experimental and theoretical approaches, accurate and reliable numerical tools are needed to understand the underlying mechanisms of multiphase displacement. Conventional macro-scale simulations have achieved great success in solving the continuity, momentum, and species balance equations which rely on the empirical estimation of relative permeabilities. This approach, however, may not be able to capture effects related to micro-scale structure in a multiphase system. On the contrary, pore-scale simulation can provide local detailed information about the flow field. The previous works on pore-scale simulation of immiscible two-phase flows in porous media [5; 6; 7; 8; 9; 10] have been shown to be capable of capturing the effects of various factors, e.g., capillary number, viscosity ratio, surface wettability, and media heterogeneity. However, due to the inherent complexity related to interface evolution, heterogeneous geometry, dominant capillarity, and moving contact line, the computational simulation of multiphase displacement in realistic porous media still remains a research challenge. As a mesoscopic method, the lattice Boltzmann method (LBM) has been widely accepted as a useful alternative for simulating multiphase flows, in particular, with the advantage of dealing with complex geometries such as porous media [11]. A number of multiphase LBM models have been proposed, e.g., the colour gradient model [12; 13], the pseudo-potential model [14], the free energy model [15] and the mean-field model [16]. The colour gradient model is known to be the first multiphase model. Later, a new multiphase LBM model was proposed by [PERSON] and [PERSON] [14], where the interparticle forces (whether attractive or repulsive) mimicking the molecular potentials are introduced by a modified equilibrium velocity in the equilibrium distribution functions. The free energy model is thermodynamically consistent, and it considers a generalised equilibrium distribution function that includes a nonideal pressure tensor term. The original free energy model suffers from a lack of Galilean invariance, but later works have mitigated this problem [17; 18]. The mean-field model describes interfacial dynamics based on the mean-field theory, and uses a pressure distribution function instead of its density counterpart to reduce discretization errors in the calculation of density gradient and to improve the numerical stability for variable density. Extensive reviews of these LBM multiphase models can be found in previous research [8; 19]. The viscosity contrast is large while the density difference is small for the water or carbon dioxide involved oil recovery process. However, most LBM models can only access a limited range of viscosity ratios for dynamic fluid displacement. [PERSON] et al. [20] performed a systematic comparison of three popular multiphase models, i.e., the pseudo-potential model, the free energy model, and the colour gradient model, for flow simulations in porous media.