In this paper, the multiphase lattice Boltzmann flux solver (MLBFS), where the phase field model and the apparent liquid permeability model are built-in, is developed to simulate incompressible multiphase flows in fractal pore structure at the representative elementary volume scale. MLBFS takes advantage of the traditional Navier–Stokes solver (e.g., geometric flexibility and direct handling of complex boundary conditions) and lattice Boltzmann method (e.g., intrinsically kinetic nature, simplicity, and parallelism). It is easily applied to simulate multiphase flows transport in the porous medium with large density ratios and high Reynolds numbers. This study focuses on the fluid flow in fractal pore structures and provides an in-depth discussion of the effects of non-Newtonian index, fractal parameters, and density ratios on multiphase flow. The proposed model is validated with benchmark problems to test the applicability and reliability of the MLBFS in describing fluid flow in fractal pore structures with large density ratios and viscosity ratios. Simulation results show that the fractal parameters (i.e., fractal dimension, tortuous fractal dimension, porosity, and capillary radius ratio) can accurately characterize fractal pore structure and significantly affect the apparent liquid permeability. In addition, the flow rate increases with the fractal dimension and decreases with the tortuous fractal dimension, while both flow rate and apparent liquid permeability decrease as the capillary radius ratio. It is also noteworthy that the effect of nonlinear drag forces cannot be neglected for shear-thickened flows.