## 2 Modelling
### Mathematical Formulation of CFD Domain  
For the CFD model in this study, the back-and-forth error compensation (BFECC) algorithm with the level set (LS) method was adopted in the physical modelling in Altair AcuSolve software. The integration of BFECC and LS improves the accuracy of a semi-Lagrangian integration for improved computational performance [39]. Various numerical studies have reported that the implementation of the BFECC algorithm results in second-order accuracy in the space and time domains and generates a relatively realistic animation of the two-phase flow [40; 41; 42; 43]. The CFD model numerically solves incompressible Navier-Stokes equations as follows:  
\[\rho\left[\partial_{t}\overset{\rightarrow}{V}+\overset{\rightarrow}{(V.\
abla)}\overset{\rightarrow}{V}\right]-\
abla.\sigma(\overset{\rightarrow} {V},p)=\overset{\rightarrow}{f} \tag{1}\]  
\[\
abla.\overset{\rightarrow}{V}=0 \tag{2}\]  
where \(\overset{\rightarrow}{V}\) and \(p\), respectively, represent the velocity and pressure profiles, \(\sigma\) and \(\rho\), respectively, denote the fluid stress tensor and density and \(\overset{\rightarrow}{f}\) is the summation of the bulk and surface forces per unit volume. In addition, \(\
abla.\overset{\rightarrow}{V}\) is defined as zero to physically represent the time rate of changing the volume of the flowing fluid element. In the case of fluids with a constant density despite the applied pressure, they adhere to the following equation:  
\[\sigma(\overset{\rightarrow}{V},p)=-p\mathbf{I}+2\mu\ D(\overset{\rightarrow}{V}) \tag{3}\]  
where \(\mu\) is the fluid absolute viscosity and \(D(\overset{\rightarrow}{V})\) represents the strain tensor rate of the fluid. In this work, both fluids satisfied the aforementioned equations. For the multiphase flow simulation, the numerical model was supported with a combination of predefined boundary and initial conditions at the fluid domain and the interface which are intended to solve Equations (1) and (2). In this numerical model, the no-slip boundary condition at the interface was considered to ensure that the defined pressure and velocity profiles were continuous along the interface. The LS method was represented by the Lipschitz scalar function \(\varnothing\), which is defined as a signed distance function that measures geometrical properties and tracks the transit topological changes in the interface of fluid regions \(\Omega_{i}\) and \(\Omega_{j}\) which share the same boundary \(\partial\Omega\) with an interface denoted by \(\Gamma_{ij}\). The signed distance function is defined as  
\[\varnothing(x)=\left\{\begin{array}{ll}d\big{(}x,\Gamma_{ij}\big{)}  
zero level implicitly representing the interface (\(\Gamma_{ij}\)). As the fluid regions \(\Omega_{i}\) and \(\Omega_{j}\) evolve with time and satisfy Navier-Stokes Equation (2), the LS function is governed by the following expression:  
\[\frac{\partial\varnothing}{\partial t}+\mathrm{V.}\
abla\varnothing=0 \tag{5}\]  
However, an auxiliary equation is required to solve Equation (5) at each time step because of the large variation and the difficulty of preserving the intrinsic property of \(|\
abla\varnothing|=1\) as a result of the velocity field in the second term (\(\mathrm{V.}\
abla\varnothing\)). Moreover, the curvature of the interface can develop irregularities on the LS function and thereby lead to dissipative errors. Therefore, the BFECC algorithm was introduced to compute the error compensation (\(e\)) from the LS function (\(\varnothing^{n}\)) at time \(t_{n}\) by evolving the reversible differential Equation (5) forward in time to measure \(\widetilde{\varnothing}^{n+1}\) at time \(t_{n+1}\) and backward in time to measure \(\varnothing^{n}1\) at time \(t_{n}\). As the same advection scheme was adopted for the forward and backward operations, the errors were expected to be similar in both steps.