We present a numerical methodology for three-dimensional large-scale simulations of two-fluid flow for generalised Newtonian fluids exhibiting non-Newtonian behaviour such as a non-zero yield stress and power-law dependency on strain-rate. The incompressible continuity and Cauchy momentum equations, along with appropriate rheological models, are solved using a computational framework initially developed at Lawrence Berkeley National Laboratory. The solver uses second-order Godunov method- ology for the advective terms and semi-implicit diffusion in the context of an approximate projection method to evolve the system in time. We have extended the algorithm to en- able the simulation of Herschel-Bulkley fluids by means of a mathematical regularisation of the constitutive equation which describes the fluid rheology. Additionally, interfaces between fluids with different properties are treated using a passively advected indicator function. The performance of the software is validated for two-dimensional displacement flow and tested on a three-dimensional viscoplastic dambreak.