  Geometric flux-based Volume-of-Fluid (VOF) methods are widely considered
consistent in handling two-phase flows with high density ratios. However,
although the conservation of mass and momentum is consistent for two-phase
incompressible single-field Navier-Stokes equations without phase-change,
discretization may easily introduce small inconsistencies that result in very
large errors or catastrophic failure. We apply the consistency conditions
derived for the unstructured Level Set / Front Tracking method to flux-based
geometric VOF methods, and implement our discretization into the
plicRDF-isoAdvector geometrical VOF method. We find that computing the mass
flux by scaling the geometrically computed fluxed phase-specific volume
destroys the equivalence between the scaled volume fraction equation and the
mass conservation equation, depending on the choice for the temporal and
convective term discretization schemes. We propose two solutions. First, based
on the analysis of discretization errors, we suggest a consistent combination
of the temporal discretization scheme and the interpolation scheme for the
momentum convection term. Second, similar to our previous work on the
unstructured Level Set / Front Tracking method, we solve an auxiliary mass
conservation equation with a geometrical calculation of the face-centered
density. We prove the equivalence between these two approaches mathematically
and verify and validate their numerical stability for density ratios in the
range $[1,10^6]$ and viscosity ratios in the range $[10^2,10^5]$.
