Two-phase incompressible flow is often computed using a one-fluid representation, where density and viscosity vary sharply across the interface. This setting is challenging for the Navier–Stokes equations, and the difficulty becomes more pronounced for fully implicit solvers because pressure–velocity coupling yields ill-conditioned saddle-point systems with variable coefficients. The problem is further aggravated in Volume-of-Fluid (VOF) methods. Achieving a sharp interface while maintaining mass conservation is nontrivial, and many algebraic VOF advection updates contain a spurious source term proportional to the velocity divergence. When divergence of velocity is not zero, this term leads to systematic mass loss, which is often mitigated in practice by adding artificial compression. We present an algorithm designed to address the above challenges. The fluid dynamics step is solved implicitly using an SGS big-box Vanka smoother applied over the entire domain, providing an implicit pressure–velocity coupling and maintaining the discrete incompressibility constraint at each step. The phase fraction is then advanced explicitly using a CUIBS-type advection scheme. This combination avoids decoupling errors and suppresses the divergence-driven source term in algebraic VOF transport, enabling excellent mass conservation without any artificial compression terms, while keeping numerical diffusion low. A key component is the formulation of the SGS relaxation operator for the two-phase, variable-coefficient saddle-point system. We construct the big-box Vanka operator so that it preserves symmetry and positive definiteness while remaining compatible with the incompressibility constraint. The resulting method enables efficient simulations of challenging benchmarks, including three-dimensional rising bubbles and Rayleigh–Taylor instability.