## 2 Colour Gradient Lattice Boltzmann Model for Two-Phase Flow
### Single-Phase Collision Operator  
To achieve high accuracy and good stability, we implemented the collision step in accordance with the two-relaxation-time (TRT) model [28; 29; 30]--a special multiple relaxation time (MRT) model with only two relaxation rates [31]. With the TRT model, the single-phase collision operator is expressed as  
\[\left(\Omega_{i}^{k}\right)^{(1)}=-(\mathbf{M}^{-1}\mathbf{S}\mathbf{M})_{ij} (f_{j}^{k}-f_{j}^{k,eq}), \tag{4}\]  
where \(\mathbf{M}\) is the transformation matrix and is given by  
\[\mathbf{M}=\begin{pmatrix}1&1&1&1&1&1&1&1\\ -4&-1&-1&-1&-1&2&2&2\\ 4&-2&-2&-2&-2&1&1&1&1\\ 0&1&0&-1&0&1&-1&-1&1\\ 0&-2&0&2&0&1&-1&-1&1\\ 0&0&1&0&-1&1&1&-1&-1\\ 0&0&-2&0&2&1&1&-1&-1\\ 0&1&-1&1&-1&0&0&0&0\\ 0&0&0&0&0&1&-1&1&-1\end{pmatrix}. \tag{5}\]  
The diagonal relaxation matrix \(\mathbf{S}\) is defined as  
\[\mathbf{S}=(1,s_{e},s_{e},1,s_{q},1,s_{q},0,s_{q},s_{v},s_{v}). \tag{6}\]  
For the conserved moments of density and momentum, the relaxation rates were set to \(1\). \(s_{e}\) and \(s_{v}\) are related to the bulk and shear viscosities, and \(s_{e}\) and \(s_{q}\) are free parameters. Following the recommendations of [PERSON] et al. [32], these relaxation rates were taken as\[s_{e}=s_{e}=s_{v}=\frac{1}{\tau},\quad s_{q}=8\frac{(2-s_{v})}{(8-s_{v})}, \tag{7}\]  
where \(\tau\) is the dimensionless relaxation time. The mass conservation for each fluid distribution function and the total momentum conservation require  
\[\rho_{k}=\underset{i}{\sum}f_{i}^{k}=\underset{i}{\sum}f_{i}^{k,eq},\quad\rho \boldsymbol{u}=\underset{i}{\sum}\underset{k}{\sum}f_{i}^{k}\boldsymbol{c}_{i} =\underset{i}{\sum}\underset{k}{\sum}f_{i}^{k,eq}\boldsymbol{c}_{i}, \tag{8}\]  
where \(\rho\) is the total density, and \(\rho_{k}\) is the density of fluid \(k\), and \(\boldsymbol{u}\) is the local fluid velocity.  
The equilibrium distribution function (\(f_{i}^{k,eq}(\boldsymbol{x},t)\)) is a truncated Taylor series expansion up to the second order in Mach number of the Maxwell-Boltzmann equilibrium distribution function, which was chosen to satisfy Equation (8):  
\[f_{i}^{k,eq}(\rho_{k},\boldsymbol{u})=\rho_{k}W_{i}\left[1+\frac{\boldsymbol{ c}_{i}\cdot\boldsymbol{u}}{c_{s}^{2}}+\frac{\left(\boldsymbol{c}_{i}\cdot \boldsymbol{u}\right)^{2}}{2c_{s}^{4}}-\frac{\boldsymbol{u}^{2}}{2c_{s}^{2} }\right], \tag{9}\]  
in which, \(c_{s}=1/\sqrt{3}\) is the lattice sound speed, and the weight coefficients are \(W_{0}=4/9\), \(W_{1-4}=1/9\) and \(W_{5-8}=1/36\) for the D2Q9 lattice model. As we focused on the oil-water two-phase systems, the density ratio of both fluids was set to be unity.  
To ensure a constant viscosity across the interface when both fluids were of different viscosities, the following harmonic mean [33] was employed to determine the viscosity of the fluid mixture:  
\[\frac{1}{\eta}=\frac{1+\rho^{N}}{2\eta_{R}}+\frac{1-\rho^{N}}{2\eta_{B}}, \tag{10}\]  
where \(\eta_{k}\) (\(k=R\) or \(B\)) is the dynamic viscosity of fluid \(k\), and \(\eta\) is the dynamic viscosity of the fluid mixture which is related to the dimensionless relaxation time \(\tau\) by \(\eta=c_{s}^{2}\rho(\tau-0.5)\delta t\); \(\rho^{N}\) is the indicator function used to distinguish different fluids, and it is defined as  
\[\rho^{N}(\boldsymbol{x},t)=\frac{\rho_{R}(\boldsymbol{x},t)-\rho_{B}( \boldsymbol{x},t)}{\rho_{R}(\boldsymbol{x},t)+\rho_{B}(\boldsymbol{x},t)}. \tag{11}\]