In the present study, we develop a three-dimensional two-way coupled Euler–Euler model to simulate the dilute suspensions of fine particles. In addition to the inter-phase drag term, commonly appearing in standard EE formulations, the model formulation includes inter-phase momentum exchange resulting from added mass, which is not negligible in solid–liquid systems. Moreover, through a two-phase pressure projection method, the present numerical model ensures that the incompressibility of the solid–liquid mixture is also taken into consideration. A series of numerical experiments on the particle-induced Rayleigh–Taylor (RT) instability is carried out to investigate bulk mixing attributable to the initial concentration of particles, covering a range of suspension from dilute to dense (O(0.001-0.05) in volume fraction). This study identifies deviations in the current two-phase simulations by comparing them with single-phase approximations. Our results indicate that the deviations are caused by non-equilibrium particle inertia and mixture incompressibility. In the dilute suspension, it is found that the non-equilibrium particle inertia enhances vertical motion of bubbles and spikes, resulting in a higher efficiency in vertical mixing, compared to the results from single-phase simulations. However, as initial concentration increases, the influence of mixture incompressibility becomes more pronounced and is able to induce a significant suppression of upward-moving motion of bubbles, which in turn decreases the efficiency of vertical mixing.