This paper presents an explicit characteristics-based, conservative, finite-difference method for the simulation of incompressible multiphase flows. The method is based on the artificial compressibility concept, e,,tended to variable-density, and uses a time stretching procedure to relieve the acoustic constrain. We take advantage of the algorithmic simplicity and hyperbolicity provided by the artificial compressibility to develop a flow solver that is numerically robust, accurate and effective for massively parallel computations of incompressible multifluid flows. The resulting method, named Numerical Acoustic Relaxation or NAR, is a combination of the AC concept with the Level Set method for interface-capturing and the Ghost-Fluid method to compute flows with multiple, arbitrary density variation, free or stationary interfaces. In this paper we demonstrate convergence and accuracy of the solver by computing such standard test problems as the "Lid-Driven Cavity" and "Doubly Periodic Shear Layer". Competitiveness with approximate projection, vorticity stream function, pseudospectral, and Lattice Boltzmann Equation is also discussed. In addition, we demonstrate the interface-capturing features of NAR by means of the "simple" "Rayleigh-Taylor" and "Water Column Collapse" problems. (C) 2004 Elsevier Ltd. All rights reserved.