## 2 Models and Methods
### Level Set Method  
The level set method is a numerical technique used to represent and track moving interfaces between two fluids. In this method, the material interface is represented as a zero level isosurface of the function \(\varphi\) (\(x\), \(t\)), which is a level set function. The level set function has a positive value outside the liquid bridge, a zero-value interface between gas and liquid and a negative value inside the liquid bridge. By solving a convection equation using the level set function, the motion of the interface can be predicted and tracked over time. This method was first introduced by [PERSON] and [PERSON] [45] and has since become widely used in various computational applications related to fluid dynamics and material sciences. With respect to the level set function \(\varphi\) (\(x\), \(t\)), the movement of the interface is forecasted by solving the subsequent convection equation:  
\[\varphi_{t}+\mathbf{V}\cdot\
abla\varphi=0 \tag{12}\]  
To ensure the level set function \(\varphi\) retains its status as a distance function, it is necessary to perform a process known as level set function reinitialization. This process involves updating the level set function to conform to the expected properties of a distance function [46]:  
\[\left\{ \begin{array}{l}\varphi_{t}=sign(\varphi_{0})(1-|\
abla\varphi|)\\ \varphi(X,0)=\varphi_{0}(X)\end{array} \right. \tag{13}\]  
where \(|\
abla\varphi|=\sqrt{\varphi_{x}^{2}+\varphi_{x}^{2}}\). To simplify the solution, the \(sign(\varphi_{0})\) needs to be smoothed as [14]:  
\[sign(\varphi_{0})=\varphi/\sqrt{\varphi^{2}+\epsilon^{2}} \tag{14}\]  
Despite the completion of the re-initialization process of the level set function, the area conservation cannot be reached yet due to the numerical dissipation. Therefore, it is recommended to solve the equation for the compensation of the area until it reaches a stable state [47]:  
\[\frac{\partial\varphi}{\partial\overline{t}}+[1-A(t)/A_{0}]f(c)|\
abla\varphi |=0 \tag{15}\]  
where \(A(t)\) represents the liquid bridge area at time \(t\) as determined by the level set function, \(A_{0}\) represents the initial area of the liquid bridge. The function \(f(c)\) is a representation of the area constraint and can be regarded as a function of the local curvature. It is subject to variation with respect to the parameters \(h\) and \(n\), which is expressed as:  
\[f(c)<1.0+h^{\prime}\gamma^{n} \tag{16}\]  
The values \(h^{\prime}=0\) and \(n=0\) are denoted, as they result in an accelerated convergence of the numerical procedure. The criteria for convergence can be formulated as follows:  
\[\left|A_{0}-A\left(\overline{t}\right)\right|<1.0\times 10^{-4} \tag{17}\]  
By employing this method, the total mass is conserved with a high level of accuracy, not deviating beyond 0.01%. To mitigate against numerical instability on the interface, particularly in the event of a significant density ratio existing within the two-phase system, it becomes necessary to smooth the physical properties in close proximity to the interface. The interface is assumed to comprise of a thin layer that possesses a thickness of 2\(\alpha\) and maintains continuity without abrupt alterations concerning physical properties. To aid in this process, both the Heaviside and Dirac functions are introduced, and their definitions are listed as follows [48]:  
\[H_{\alpha}(\varphi)=\left\{ \begin{array}{l}0,\;if\;\varphi<-\alpha\\ [1+\varphi/\alpha+\sin(\pi\varphi/\alpha)]/2,\;if\;|\varphi|\leq\alpha\\ 1,\;if\;\varphi>\alpha\end{array} \right. \tag{18}\]  
\[\delta_{\alpha}(\varphi)=\left\{ \begin{array}{l}(1+\cos(\pi\varphi/\alpha) )/2\alpha,\;if\;|\varphi|<\alpha\\ 0,\;if\;|\varphi|\geq\alpha\end{array} \right. \tag{19}\]  
The density \(\rho\) and dynamic viscosity \(\mu\) can be exhibited as follows:  
\[\left\{ \begin{array}{l}\rho(\varphi)=\rho_{l}+\left(\rho_{g}-\rho_{l}\right) H_{\alpha}(\varphi)\\ \mu(\varphi)=\mu_{l}+\left(\mu_{g}-\mu_{l}\right)H_{\alpha}(\varphi)\end{array} \right.