A multiphase extension of the δ-plus-SPH (smoothed particle hydrodynamics) model is developed for modeling non-Newtonian multiphase flow. A modified numerical diffusive term and special shifting treatment near the phase interface are introduced to the original δ-plus-SPH model to improve the accuracy and numerical stability of the weakly incompressible SPH model. The Herschel–Bulkley model is used to describe non-Newtonian fluids. A sub-particle term is added in the momentum equation based on a large eddy simulation. The graphic processing unit (GPU) acceleration technique is applied to increase the computational efficiency. Three test cases including, a static tank, Poiseuille flow, and submarine debris flow, are presented to assess the performance of the new multiphase SPH model. Comparisons with analytical solutions, experimental data, and previous numerical results indicate that the proposed SPH model can capture highly transient incompressible two-phase flows with consistent pressure across the interface.