while the $x$, $y$ and $z$ components of velocity are located at the faces. All spatial derivatives are approximated by standard second-order centered differences. The treatment of the free surface uses a hybrid Front-Tracking/Level-Set technique which defines the interface both by the Level-Set distance function field on the Eulerian grid as well as by triangles on the Lagrangian interface mesh.

The surface tension $\textbf{F}$ is implemented by the hybrid/compact formulation [@Shin07] $$\textbf{F} = \sigma \kappa_{H} \nabla I , \qquad
\kappa_{H}=\frac{\textbf{F}_{L}\cdot\textbf{N}}{\textbf{N}\cdot\textbf{N}}$$ where $\sigma$ is the surface tension coefficient and $\kappa_{H}$ is twice the mean interface curvature field calculated on the Eulerian grid, with $$\textbf{F}_{_{L}} = \int_{\Gamma(t)} \kappa_f  \textbf{n}_{f} \delta_f
\left( \textbf{x} - \textbf{x}_{f} \right) \mathrm{d}a ,
\qquad \textbf{N} = \int_{\Gamma(t)}  \textbf{n}_{f} \delta_f
\left( \textbf{x} - \textbf{x}_{f} \right) \mathrm{d}a$$ Here, $\textbf{x}_f$ is a parameterization of the time-dependent interface, $\Gamma(t)$, and $\delta_f \left( \textbf{x} -\textbf{x}_{f} \right) $ is a Dirac distribution that is non-zero only where $\textbf{x}=\textbf{x}_f$; $\textbf{n}_f$ stands for the unit normal vector to the interface and $ \mathrm{d}a$ is the area of an interface element; $\kappa_f$ is twice the mean interface curvature obtained on the Lagrangian interface. The geometric information, unit normal, $\textbf{n}_f$, and interface element length, $ \mathrm{d}a$ in $\textbf{N}$ 