We present a fully second order IMplicit/EXplicit (IMEX) time integration technique for solving incompressible multi-phase flow problems. A typical incompressible multi-phase flow model consists of the Navier-Stokes equations plus an interface dynamics equaiton (e.g, the level set equation). Our IMEX strategy is applied to such a model in the following manner. The hyperbolic terms of the Navier-Stokes equations together with the interface dynamics equation are solved explicitly (Explicit Block) making use of the well-understood explicit numerical schemes. On the other hand, the non-hyperbolic (stiff) parts of the flow equations are solved implicitly (Implicit Block) within the framework of the Jacobian-Free Newton Krylov (JFNK) method. In our algorithm implementation, the explicit block is embedded in the implicit block in a way that it is always part of the non-linear functions evaluation. In this way, there exists a continuous interaction between the implicit and explicit algorithm blocks meaning that the improved solutions (in terms of time accuracy) at each nonlinear iteration are immediately felt by the explicit block and the improved explicit solutions are readily available to form the next set of nonlinear residuals. This continuous interaction between the two algorithm blocks results in an implicitly balanced algorithm in that all the non-linearities due to coupling of different time terms are converged with the desired numerical time accuracy. In other words, we obtain a self-consistent IMEX method that eliminates the possible order reductions in time convergence that is quite common in certain types of nonlinearly coupled systems. We remark that an incompressible multi-phase flow model can be a highly non-linearly coupled system with the involvement of very stiff surface tension source terms. These kinds of flow problems are difficult to tackle numerically. In other words, highly non-linear terms may remain unconverged leading to time inaccuracies or time order reductions to the first order even though the overall numerical scheme is designed as high order (second order or higher). These and few more issues are addressed in this paper. We have numerically tested our newly proposed scheme by

solving several multi-phase flow settings such as an air bubble rising in water, a Rayleigh-Taylor instability problem that is initiated by placing a heavy fluid on top of a lighter one, and a droplet problem in which a water droplet hits the pool of water. Our numerical results indicate that we have achieved the second order time accuracy without any order reductions. Moreover, the interfaces

between the fluids are captured reasonably well.