## 1 Introduction
There are many LB models proposed for solving multiphase flow, including:  
* Color gradient model based on [PERSON] and [PERSON] lattice gas model [32; 33];
* Shan-Chen (pseudopotential) model [34];
* free energy model [35; 36]; and
* phase field based LB [37; 38]. In the color gradient model, each phase is assigned a color, and for each color, a distribution function is introduced to calculate the interaction between different phases. The interface is determined using the color gradient. In the free energy model, two particle distribution functions are utilized. One distribution function calculates an order parameter, which defines each phase, and the other distribution function is responsible for the predicted velocity without considering a pressure gradient. To find the exact pressure and velocity, another Poisson-type equation must be solved in each iteration. In the phase field method, in addition to the fluid flow equation, another equation is solved to determine the interface location. The fluid flow distribution function for the local density and momentum is transformed to the mean density and momentum to improve stability, and then a cohesion force is introduced in the mean field flow. The He-Chen-Zhang (HCZ) model is a variation of the phase field-based multiphase method, in which two sets of distribution functions are used to recover the Navier-Stokes (NS) and Cahn-Hilliard (CH) equations [29; 36; 37; 39; 40; 41]. The Shan-Chen model and all of its variations are among the most used models because the interface is not introduced as a boundary condition and no extra equation is required to track the interface between different phases [29]. However, the original Shan-Chen model has several issues, including, but not limited to: Large spurious current, limitations regarding simulating fluids with a high density and viscosity ratio, and avoidance of coalescences of all bubbles [28]. [PERSON] et al. [28] classified the ways to improve the original Shan-Chen model. The improvements consist of the following: Using a realistic equation of state [42], using an interaction force with a higher order of isotropy [43], and introducing a mid-range repulsive force [44]. The phase field LB model has been successfully used to simulate many problems of multiphase flow [45; 46; 47; 48]. Most of the researchers have solved the CH equation by another set of LB distribution functions. However, solving the CH equation by the LB technique has drawbacks that may compromise the results. First, selecting the relaxation parameter for the CH equation is not apparent. Many articles use the same fluid flow relaxation time for the CH equation. This selection is not physically acceptable since different values of relaxation time may lead to different results. Next, the LB model cannot entirely recover the CH equation, but the recovered CH equation always has higher order terms. It is also hard in the LB model to discretize the convection term in the CH equation by upwind schemes since the streaming process implies a central discretization. Finally, it is hard to use higher order explicit time stepping for the LB method since LB models usually use second-order explicit time stepping, which makes the time step size significantly small for the CH equation [49; 30; 41]. In this article, the CH equation is solved by a weighted essentially non-oscillatory (WENO) scheme for the convection term and a third order Runge-Kutta time stepping scheme for the time marching. There are few articles in which the cellular automaton (CA) and LB methods were coupled to simulate bubble dynamics during dendrite growth; two of them are Refs. [7; 50]. In these cases, the model predicted gas-liquid-solid interactions correctly and the evolution of bubbles during solidification. However, in both articles, the original Shan-Chen model was used, hence suffering from all the problems mentioned earlier. Most importantly, they did not discuss the high spurious current generated by the Shan-Chen model.