*3.1. Lattice Boltzmann Equations*

The hydrodynamics of droplets passing through porous metal foam was simulated using LBM. The lattice Boltzmann simulation includes two steps. First, the collision between fluid particles is calculated, and the second, the particles are streamed. In the case of two immiscible fluids, the fluid distribution function for each fluid is described as the following:

$$f\_i^k(\mathbf{x} + \mathbf{e}\_i \delta\_{l\prime}, t + \mathbf{e}\_i \delta\_{l\prime}) - f\_i^k(\mathbf{x}, t) = -\frac{1}{\tau\_k} [f\_i^k(\mathbf{x}, t) - f\_i^{k, \text{eq}}(\mathbf{x}, t)] \tag{1}$$

where *k* stands for different particles, **x** is the discrete lattice location, **e***<sup>i</sup>* is the discrete velocity direction for the D3Q19 lattice structure, δ*<sup>t</sup>* is time step and is chosen as 1, τ*<sup>k</sup>* is relaxation parameter and *f <sup>k</sup> <sup>i</sup>* (**x**, *t*) and *f k*,*eq <sup>i</sup>* (**x**, *<sup>t</sup>*) are the number density distribution function and local equilibrium function, respectively. *f k*,*eq <sup>i</sup>* (**x**, *<sup>t</sup>*) is calculated for the LBM D3Q19 model as

$$f\_i^{k,eq}(\mathbf{x},t) = \omega\_i n\_k \left[ 1 + \frac{\mathbf{e}\_i \cdot \mathbf{u}\_k^{eq}}{c\_s^2} + \frac{\left(\mathbf{e}\_i \cdot \mathbf{u}\_k^{eq}\right)^2}{4c\_s^4} + \frac{\mathbf{u}\_k^2}{4c\_s^4} \right] \tag{2}$$

$$\mathbf{e}\_{l} = \mathbf{c} \begin{pmatrix} 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & 1 & -1 & -1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & -1 & -1 & 1 & 0 & 0 & 0 & 1 & -1 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & 1 \end{pmatrix} \tag{3}$$

$$\omega\_{i} = \begin{cases} \frac{1}{3}, \mathbf{e}\_{i} = 0 \\\\ \frac{1}{18}, \mathbf{e}\_{i} = \mathfrak{c}^{2} \\\\ \frac{1}{36}, \mathbf{e}\_{i} = 2\mathfrak{c}^{2} \end{cases} \tag{4}$$

where *nk* is the number density, and *c<sup>k</sup> <sup>s</sup>* is the sound velocity of the *k*th component that satisfies the formula *c<sup>k</sup> <sup>s</sup>* = 1/3.

The mass density and velocity of the *k*th component can be obtained:

$$
\rho\_k = m\_k n\_k = m\_k \sum\_i f\_i^k \tag{5}
$$

$$
\rho\_k u\_k = m\_k \sum\_i \mathbf{e}\_i f\_i^k \tag{6}
$$

where *mk* is the molecular mass of each component. In the Shan-Chen scheme, the interaction force is incorporated into the model by the equal velocity **u***eq <sup>k</sup>* in the equilibrium equation. The **<sup>u</sup>***eq <sup>k</sup>* is determined by the relation

$$\mathbf{u}\_k^{eq} = \mathbf{u}' + \frac{\tau\_k F\_k}{\rho\_k} \tag{7}$$

where **u** is the common velocity, which satisfies the equation

$$\mathbf{u}' = \left(\sum\_{k} \frac{\mu\_k \rho\_k}{\tau\_k}\right) \Big/ \left(\sum\_{k} \frac{\rho\_k}{\tau\_k}\right) \tag{8}$$

The total force acting on each component *Fk* includes the fluid–fluid interaction *Ff*, the fluid-solid interaction *Fads* and the external force *Fe*

$$F\_k = F\_f + F\_{ads} + F\_e \tag{9}$$

The SC model assumes that the interaction between particles is nonlocal. This means that the interaction only considers the nearest-neighbors **x** = **x** + **e***i*. The fluid–fluid interaction acting on the *k*th component at site *x* is

$$F\_{f,k} = -\Psi\_k(\mathbf{x}) \sum\_{\mathbf{x}'} \sum\_k G\_{k\overline{k}}(\mathbf{x}, \mathbf{x}') \Psi\_{\overline{k}}(\mathbf{x}') (\mathbf{x}' - \mathbf{x}) \tag{10}$$

The effective mass of *k*th component Ψ*k*(**x**) is related to the number density

$$\rho\_k: \Psi\_k(\mathbf{x}) = \rho\_{k0} \left[ 1 - \exp\left(-\frac{\rho\_k}{\rho\_{k0}}\right) \right] \tag{11}$$

*Gkk* **x**, **x** is Green's function,

$$G\_{k\overline{k}}(\mathbf{x}, \mathbf{x}') = \begin{cases} \begin{array}{c} \mathcal{g}\_{k\overline{k}'} \left| \mathbf{x} - \mathbf{x}' \right| = 1 \\ \mathcal{g}\_{k\overline{k}'} / 2, \left| \mathbf{x} - \mathbf{x}' \right| = \sqrt{2} \\ 0, otherwise \end{array} \tag{12}$$

where *gkk* controls the interaction strength between different particles, while its sign determines whether the interaction is attraction (negative) or repulsion (positive), and it is sufficient to only involve the nearest neighbor interactions. The interaction *Fads* for the *k*th component at the fluid/solid interface is defined as

$$F\_{\rm ads,h} = -\Psi\_k(\mathbf{x}) \sum\_k G\_{\rm w,k} \mathbf{s}\left(\mathbf{x}'\right) \left(\mathbf{x}' - \mathbf{x}\right) \tag{13}$$

where *s* **x** is the switch parameter which is 1 when **x** is in the pore space of the solid lattice, otherwise it is 0. *Gw*,*<sup>k</sup>* determines the interaction's type and magnitude between fluid and solid wall. Different wettability can be obtained by adjusting *Gw*,*k*. Non-wetting fluid and wetting fluid are obtained when *Gw*,*<sup>k</sup>* > 0 and *Gw*,*<sup>k</sup>* < 0, respectively. *Gw*,*<sup>k</sup>* is obtained by comparing the measured contact angles and the simulation results.
