The arbitrary Lagrangian-Eulerian (ALE) framework presented in [Sahin and Guventurk, An arbitrary Lagrangian-Eulerian framework with exact mass conservation for the numerical simulation of 2D rising bubble problem. Int J Numer Methods Eng. 2017;112:2110-2134] has been extended to simulate three-dimensional immiscible incompressible multiphase fluid flows. The div-stable side centered unstructured finite volume formulation is used for the discretization of the incompressible Navier-Stokes equations. A novel kinematic boundary condition that satisfies the local discrete geometric conservation law (DGCL) is implemented to ensure the mass conservation of each species at the machine precision. The pressure field is treated to be discontinuous across the interface with the discontinuous treatment of density and viscosity. Surface tension force is treated as a tangent force and discretized in a semi-implicit form. Two different approaches for the computation of unit normal vector have been implemented: the least squares biquadratic surface fitting (LSBSF) and the mean weighted by sine and edge length reciprocals (MWSELR). The resulting large system of algebraic equations is solved in a fully coupled manner in order to improve the time step restrictions. As a preconditioner, an approximate matrix factorization similar to that of the projection method is employed and the parallel algebraic multigrid solver BoomerAMG provided by the HYPRE library, which is accessed through the PETSc library, has been utilized for the scaled discrete Laplacian of pressure and the diagonal blocks of mesh deformation equations. The accuracy of the proposed method is initially validated on the static bubble problem, since the surface tension force is highly sensitive to the accurate evaluation of the unit normal vector and the inaccuracies significantly contribute to unphysical velocities, called parasitic currents. The calculations indicate that the parasitic currents can be reduced to machine precision for the MWSELR method. Then, the proposed approach is applied to the single bubble rising in a viscous quiescent liquid for both low and high density ratios. The calculations produce accurate predictions of the bubble shape, center of mass, rise velocity, and so forth. Furthermore, the mass of each species is conserved at machine precision and discontinuous pressure field is obtained in order to avoid errors due to the incompressibility restriction in the vicinity of liquid-liquid interfaces at large density and viscosity ratios.