1. Introduction
Over the last decade, pre-wetted surfaces and slippery liquid-infused porous substrate (SLIPS) surfaces have attracted significant academic and industrial attention due to their anti-icing, low friction and multifunctional properties (Wong et al. Reference Wong, Kang, Tang, Smythe, Hatton, Grinthal and Aizenberg2011; Anand et al. Reference Anand, Paxson, Dhiman, Smith and Varanasi2012; Rykaczewski et al. Reference Rykaczewski, Anand, Subramanyam and Varanasi2013; Smith et al. Reference Smith, Dhiman, Anand, Reza-Garduno, Cohen, McKinley and Varanasi2013; Subramanyam, Rykaczewski & Varanasi Reference Subramanyam, Rykaczewski and Varanasi2013; Solomon, Khalil & Varanasi Reference Solomon, Khalil and Varanasi2014; Hao et al. Reference Hao, Li, Liu, Zhou, Liu, Liu, Che, Zhou, Sun and Li2015; Park et al. Reference Park, Kim, Grinthal, He, Fox, Weaver and Aizenberg2016; Villegas et al. Reference Villegas, Zhang, Abu Jarad, Soleymani and Didar2019). These surfaces rely on a stable viscous fluid film to separate drops and other particles from the supporting solid, offering a relatively slippery, pressure-stable surface compared to their solid counterparts (Bohn & Federle Reference Bohn and Federle2004; Solomon et al. Reference Solomon, Subramanyam, Farnham, Khalil, Anand and Varanasi2016; Huang & Guo Reference Huang and Guo2019; Semprebon et al. Reference Semprebon, Sadullah, McHale and Kusumaatmaja2021). Over the last few years, wetting on soft materials has emerged as a topical field in fluid mechanics, where pronounced deformations of the soft substrate near the contact line greatly affect the spreading dynamics (Semprebon, McHale & Kusumaatmaja Reference Semprebon, McHale and Kusumaatmaja2017; Sadullah, Semprebon & Kusumaatmaja Reference Sadullah, Semprebon and Kusumaatmaja2018; McHale et al. Reference McHale, Orme, Wells and Ledesma-Aguilar2019, Reference McHale, Afify, Armstrong, Wells and Ledesma-Aguilar2022; Andreotti & Snoeijer Reference Andreotti and Snoeijer2020). A drop contacting a substrate coated with an immiscible fluid is one type of soft wetting problem that has received far less attention (Smith et al. Reference Smith, Dhiman, Anand, Reza-Garduno, Cohen, McKinley and Varanasi2013; Daniel et al. Reference Daniel, Timonen, Li, Velling and Aizenberg2017, Reference Daniel, Timonen, Li, Velling, Kreder, Tetreault and Aizenberg2018).
When a drop of an immiscible fluid comes in contact with a SLIPS surface, it first displays spreading-like behaviour (Carlson et al. Reference Carlson, Kim, Amberg and Stone2013) driven by the formation of a liquid bridge between the drop and the film, with flow that is driven by capillarity and limited by viscous friction. Here, the growth of the liquid bridge depends on the ratios between the material properties of the fluid phases, i.e. viscosity 
$\eta$, surface tension 
$\sigma$ and density 
$\rho$, as well as the ratio between the drop's radius 
$R$ and the film's thickness 
$H$. Cuttle et al. (Reference Cuttle, Thompson, Pihler-Puzović and Juel2021) decomposed the interfacial dynamics generated when an aqueous drop contacts an oil layer into different stages: first, the short-time spreading of the oil film on the aqueous drop, forming a radially expanding liquid bridge; then the complete engulfment of the drop by the oil film, followed by the drop being ‘pushed’ into the bath or film. For larger Bond numbers (ratio between gravity and surface tension forces), gravity can also affect the drop dynamics and its equilibrium shape. Depending on the spreading factor, finite contact angles can form between the drop and the film, creating a liquid lens (De Gennes, Brochard-Wyart & Quéré Reference De Gennes, Brochard-Wyart and Quéré2004; Kuiper & Hendriks Reference Kuiper and Hendriks2004; Cheng & Yeh Reference Cheng and Yeh2007; Hack et al. Reference Hack, Tewes, Xie, Datt, Harth, Harting and Snoeijer2020). Contrary to capillary spreading on a solid substrate (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009), the fluid film introduced by the SLIPS alleviates the stress singularity at the contact line while generating additional viscous dissipation that extends into the fluid film. A detailed description of the dynamics of the contact-line motion and the formation of a film on the drop for this soft wetting phenomenon is yet to be described, as it is challenging experimentally to measure the interfacial dynamics at these small time and length scales that are additionally complicated by the presence of three fluid phases. In this paper, we describe the details of the small-scale interfacial dynamics of the three-phase flow by using the lattice Boltzmann method (LBM) to solve the conservative phase field model, with simulation results compared directly to experiments. In addition, we address how the thickness of the coating film affects both the inertially limited and viscously limited drop engulfment dynamics, and describe where the dissipation is localized in these flows.
Post drop–film contact, the characteristic velocity of the radial growth of the liquid bridge is set by the capillary velocity 
$\sigma /\eta$. Thus it is natural to describe the drop's dynamics using the Ohnesorge number 
$Oh=\eta _o/\sqrt {\rho _o \sigma _{wo}R}$, where the subscript 
$o$ represents the oil film, and the subscript 
$w$ represents the aqueous drop. For thin films, 
$H/R \ll 1$, this radial growth has been shown to scale with an inertial time scale 
$t_\rho =\sqrt {\rho R^3/\sigma _{ao}}$ when 
$Oh\ll 1$, and with a viscous time scale 
$t_\eta =\eta _o R/\sigma _{ao}$ for 
$Oh\gg 1$, where the subscript 
$a$ represents air. The evolution of the liquid bridge's radial growth can then be found to follow the scaling law 
$r/R\sim (T/t_\rho )^{1/2}$ when 
$Oh\ll 1$, and 
$r/R\sim (T/t_\eta )^{7/20}$ when 
$Oh\gg 1$ (Carlson et al. Reference Carlson, Kim, Amberg and Stone2013). For thick films, 
$H/R\gg 1$, Cuttle et al. (Reference Cuttle, Thompson, Pihler-Puzović and Juel2021) described the long-time entrapment and transport into the bath. Despite these recent considerations, however, drop engulfment dynamics for moderate 
$H/R$ have been considered sparsely and are the focus of this work.
When an aqueous drop contacts an immiscible liquid film, the resulting interfacial dynamics can be interpreted as a combination of two widely studied phenomena, capillary-driven contact-line motion (De Gennes Reference De Gennes1985; De Gennes et al. Reference De Gennes, Brochard-Wyart and Quéré2004; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009; Snoeijer & Andreotti Reference Snoeijer and Andreotti2013) and miscible drop coalescence (Orme Reference Orme1997; Yao et al. Reference Yao, Maris, Pennington and Seidel2005; Frising, Noïk & Dalmazzone Reference Frising, Noïk and Dalmazzone2006). Both of these phenomena have been characterized in both the inertial and viscous regimes (Guido & Simeone Reference Guido and Simeone1998; Duchemin, Eggers & Josserand Reference Duchemin, Eggers and Josserand2003; Lee & Fischer Reference Lee and Fischer2006; Paulsen, Burton & Nagel Reference Paulsen, Burton and Nagel2011; Baroudi et al. Reference Baroudi, Nagel, Morris and Lee2015; Luo, Schiffbauer & Luo Reference Luo, Schiffbauer and Luo2016), as well as in more complex situations such as under electric fields (Mugele & Baret Reference Mugele and Baret2005) and with complex surface properties (Aarts et al. Reference Aarts, Lekkerkerker, Guo, Wegdam and Bonn2005; Burton & Taborek Reference Burton and Taborek2007; Kapur & Gaskell Reference Kapur and Gaskell2007; Mazutis & Griffiths Reference Mazutis and Griffiths2012).
During rigid solid surface spreading, the radius of the drop's contact line has been used as a proxy to characterize the spreading dynamics. For low-viscosity drops experiencing inertially limited contact-line motion, the spreading radius 
$r$ has been predicted to follow a power-law relation as 
$r/R\sim (T/t_{c})^{1/2}$, where 
$T$ is the time, and 
$t_{c}$ is the characteristic time scale (Biance, Clanet & Quéré Reference Biance, Clanet and Quéré2004; Courbin et al. Reference Courbin, Bird, Reyssat and Stone2009; Carlson, Do-Quang & Amberg Reference Carlson, Do-Quang and Amberg2011; Carlson et al. Reference Carlson, Kim, Amberg and Stone2013; Eddi, Winkels & Snoeijer Reference Eddi, Winkels and Snoeijer2013). For viscously limited spreading on perfectly wetting substrates, analogies to viscous coalescence have been made, motivated by the curvature between the drop and its precursor film (Bird, Mandre & Stone Reference Bird, Mandre and Stone2008; Eddi et al. Reference Eddi, Winkels and Snoeijer2013). In non-perfectly wetting systems, the spreading radius 
$r$ has been found to follow a power-law-like behaviour 
$r \sim t^\beta$, where 
$\beta$ has been shown to vary with the surface wettability at short time (Bird et al. Reference Bird, Mandre and Stone2008; Winkels et al. Reference Winkels, Weijs, Eddi and Snoeijer2012; Mitra & Mitra Reference Mitra and Mitra2016), before transitioning to a Tanner’s-law-like spreading behaviour at long time (Tanner Reference Tanner1979). These differences in the short-time power-law behaviour have been shown to be a function of contact-line friction, which itself is a function of surface wettability, and recent scalings have been shown to collapse these differences (Carlson, Bellani & Amberg Reference Carlson, Bellani and Amberg2012).
During drop–drop coalescence, much focus has been placed on the shape and evolution of the neck radius 
$r$ that forms as the two miscible drops touch, as it gives a way to characterize the coalescence dynamics and identify the different flow regimes (Eggers, Lister & Stone Reference Eggers, Lister and Stone1999; Aarts et al. Reference Aarts, Lekkerkerker, Guo, Wegdam and Bonn2005). In inertially limited systems, the neck radius 
$r$ expands rapidly post coalescence, and the curvature 
$\kappa \sim R/r^2$ generates a capillary pressure 
$p_\sigma \sim \sigma \kappa \sim \sigma R/r^2$. Tracing the streamlines from the drop interior to the interface, it follows then from the Bernoulli equation that 
$p_\sigma =({\rho }/{2})u_r^2$, where 
$u_r=\textrm {d}r/\textrm {d}t$ is the neck velocity. Consequently, the evolution of spreading radius in time can be written as 
$r/R\sim (T/t_\rho )^{1/2}$ (Eggers et al. Reference Eggers, Lister and Stone1999), where 
$t_\rho =\sqrt {\rho R^3/\sigma _{ao}}$ is the inertial time scale. In viscously limited systems, the momentum equation simplifies to the Stokes equation 
$\eta \,\nabla ^2u=\boldsymbol {\nabla } p_\sigma$, and a viscous time scale defines the dynamics, 
$t_{\eta }=\eta R / \sigma _{ao}$. Here, the evolution of spreading radius in time can be written as 
$r/R\sim (T/t_\eta )^{\beta }$, where 
$\beta$ is determined by the viscous dissipation in the liquid. Simulations and theoretical analysis of two equal-sized drops coalescing indicates 
$\beta = 0.5-1$, with potential for logarithmic correction at early times (Hopper Reference Hopper1984; Eggers et al. Reference Eggers, Lister and Stone1999; Paulsen et al. Reference Paulsen, Carmigniani, Kannan, Burton and Nagel2014).
In recent years, simulations of contact-line spreading and drop coalescence have helped to reveal additional dynamics in both the viscous and inertial regimes that are otherwise challenging or nearly impossible to observe experimentally (Lee & Fischer Reference Lee and Fischer2006; Paulsen et al. Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012; Baroudi et al. Reference Baroudi, Nagel, Morris and Lee2015; Xia, He & Zhang Reference Xia, He and Zhang2019; Baroudi & Lee Reference Baroudi and Lee2020; Anthony, Harris & Basaran Reference Anthony, Harris and Basaran2020; Anthony et al. Reference Anthony, Kamat, Thete, Munro, Lister, Harris and Basaran2017). Paulsen et al. (Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012), for example, used simulations to argue that in the early time of drop coalescence, a third regime beyond the Stokes and inertial regimes exists, where inertia, viscosity and surface tension forces are all important. Baroudi et al. (Reference Baroudi, Nagel, Morris and Lee2015) employed the LBM to model the coalescence of two drops in both a saturated vapour phase and a non-condensable gas phase, showing good agreement with experiments (Yao et al. Reference Yao, Maris, Pennington and Seidel2005), and Xia et al. (Reference Xia, He and Zhang2019) was able to resolve the transition region between the inertially limited and viscously limited drop–drop coalescence regimes. More recently, Baroudi & Lee (Reference Baroudi and Lee2020) used the Cahn–Hilliard equation to show the effect of interfacial mass transport during inertial drop spreading on solids with varying wetting properties. Most of these previous works, however, are restricted to systems with 
$Oh\ll 1$ or 
$Oh\sim 1$, as simulations for 
$Oh\gg 1$ have proved too costly. In this work, we address this gap by considering not only systems with 
$Oh\ll 1$ and 
$Oh\sim 1$, but also systems with 
$Oh \gg 1$.
2. Methodology
2.1. Experimental set-up
Figure 1(a) describes the experimental set-up. Here, a glycerol/water pendant drop of radius 
$R=1$ mm comes in contact with a quiescent silicone oil film of thickness 
$H$. The radial expansion of the neck of the liquid bridge formed as the oil spreads on the aqueous drop during engulfment was imaged using a Photron Fastcam SA5 with a Nikon 200 mm f/4 AF-D Macro lens and two teleconverters (Nikon 
$2\times$ and 
$1.7\times$) at 3000 fps. This neck radius 
$r$ was tracked in Matlab using partial area subpixel edge detection techniques affording resolutions below 
$6\ \mathrm {\mu }{\rm m}\ {\rm pixel}^{-1}$ (Trujillo-Pino et al. Reference Trujillo-Pino, Krissian, Alemán-Flores and Santana-Cedrés2013; Trujillo-Pino Reference Trujillo-Pino2019). A complete list of the materials used is included in A.2.
Figure 1. (a) Experimental set-up. Water/glycerol pendant drops of radius 
$R\approx 1$ mm and viscosity 
$\eta _w=[0.035, 0.068, 0.154, 0.705]\ \textrm {Pa}\ \textrm {s}$ are brought into contact with a silicone oil film of thickness 
$H$ and viscosity 
$\eta _o = [0.33, 0.73, 1.5, 7.6]\ \textrm {Pa}\ \textrm {s}$. Viscosity ratios are held at 
$\eta _w/\eta _o = 1/10$. The radius of the post-contact liquid bridge 
$r$ is quantified using high-speed photography. (b) Initialization of the three-phase flow simulations. The red phase represents the aqueous drop with density 
$\rho _w$, viscosity 
$\eta _w$ and radius 
$R$. The blue phase represents the oil film with density 
$\rho _o$, viscosity 
$\eta _o$ and film thickness 
$H$. The colourless background represents air with density 
$\rho _a$ and viscosity 
$\eta _a$.
The drop's initial radius 
$R$ is held constant, and 
$Oh$ is varied by changing the oil film viscosity 
$\eta _o=[0.33, 0.73, 1.5, 7.6]\ \textrm {Pa}\ {\rm s}$. This range of oil viscosities was created by mixing 20 cSt and 1000 St silicone oils. Drops consisted of deionized (DI) water and glycerol mixtures with viscosities 
$\eta _w=[0.035, 0.068, 0.154, 0.705]\ {\rm Pa}\ {\rm s}$ chosen to produce drop–film viscosity ratios 
$\eta _w/\eta _o\approx 0.1$. All viscosities were characterized on an Anton-Paar 702 Rheometer with parallel plate geometry from shear rates 
$0.1$–
$100$ s
$^{-1}$, and all mixtures were confirmed to behave as Newtonian fluids.
Film thickness ratios were varied as 
$H/R = [0.2\unicode{x2013}4.4]$ by varying the film's thickness 
$H$. For 
$H < 1$ mm, a glass slide cleaned with DI water and isopropyl alcohol was placed on a laboratory balance, and a prescribed mass of silicone oil was deposited on the slide. The slide was then placed on a level surface, and the oil was allowed to spread fully until it became pinned by the slide's corners, a phenomenon that occurs even for fluids that wet a surface perfectly (Oliver, Huh & Mason Reference Oliver, Huh and Mason1977), forming a level oil–air interface of known thickness (Brochard-Wyart et al. Reference Brochard-Wyart, Hervet, Redon and Rondelez1991; Noblin, Buguin & Brochard-Wyart Reference Noblin, Buguin and Brochard-Wyart2004). For 
$H \geqslant 1$ mm, glass slides were glued together with a hot-glue gun to form a rectangular well. The well depth was controlled using machinery shims. A prescribed mass of oil was again added using a laboratory balance. To avoid complications in side view imaging from the concave meniscus formed by the oil on the glass wall, care was taken to ensure that the desired film thickness was slightly larger than the well depth, so imaging would occur on a level convex meniscus not obstructed by the oil's contact line.
The surface tension between the oil phase and the glycerol/water phase was measured to be 
$\sigma _{ow} \approx 0.02\ \textrm {N}\ \textrm {m}^{-1}$ by comparing the change in contact angle 
$\alpha$ of a sessile glycerol/water drop on a clean glass slide that is surrounded by either air or silicone oil. The Young–Dupré equation was then applied, assuming that the surface tension between the air and the glass is 
$0$ and that 
$\alpha = 0^\circ$ for oil on glass.
2.2. Simulation set-up and parameters
Numerical simulations of the drop engulfment process were performed in a rectangular domain; see figure 1. A water/glycerol drop with radius 
$R$, density 
$\rho _w$ and viscosity 
$\eta _w$ is placed in contact with an oil film. The centre of the drop is placed at 
$H+R+\delta$, where 
$H$ is the film thickness, 
$\delta$ is the diffuse interface thickness, 
$\rho _o$ is the density, and 
$\eta _o$ is the oil's viscosity. Air, with density 
$\rho _ a$ and viscosity 
$\eta _ a$, makes up the background phase. Table 1 provides the relative magnitudes of the parameters and the range of parameters explored in the simulations and experiments. Densities and viscosities are fixed in the simulations: 
$\rho _w/\rho _a=\rho _o/\rho _a=100$ and 
$\eta _o/\eta _a=1000$, 
$\eta _w/\eta _a=100$. The principal interest is to characterize the viscous and inertial engulfment dynamics for different ratios 
$H/R$, which is achieved by exploring the parameter space 
$Oh=[0.01\unicode{x2013}40]$ and 
$H/R=[0.03\unicode{x2013}4]$; see table 1.
Table 1. Physical properties of the liquids used in the ternary flow system.
We consider a ternary flow system that is composed of a water/glycerol drop, an oil film and air, where the spreading coefficients determine the equilibrium contact angles (Pannacci et al. Reference Pannacci, Bruus, Bartolo, Etchart, Lockhart, Hennequin, Willaime and Tabeling2008; Guzowski et al. Reference Guzowski, Korczyk, Jakiela and Garstecki2012). The spreading factor can be calculated as 
$S_o=\sigma _{wa}-\sigma _{wo}-\sigma _{ao}$. Since the surface tension of oil is typically lower than that of water/glycerol, the silicone oil will perfectly wet the water/glycerol surface. Thus we consider only spreading coefficients 
$S_o\geqslant 0$, when the oil fully engulfs the water/glycerol drop at equilibrium.
Inherent to the phase field model are diffuse interfaces with finite thicknesses. Here, we ensure that the interface thickness is small by using a Cahn number 
$Cn=\delta /R\leqslant 0.05$, which achieves the sharp interface limit by performed convergence tests shown in § A.3. The convergence tests also indicate that the oil film thickness does not affect the results, even for 
$H/R=0.1$. Since we are interested in drops with radii smaller than the capillary length, we are justified in neglecting the gravitational force in our simulations. From the simulations, we extract the fluid velocity fields, the interfacial shapes, the position of the apparent contact line, and the drop's centre of mass calculated as 
$C_m=\sum \phi _w y/({\rm \pi} R^2 )$, where 
$\phi _w$ denotes the order parameter of the water droplet, and 
$y$ is the vertical displacement scaled by the drop diameter 
$2R$.
A detailed comparison between two-dimensional (2-D) and three-dimensional (3-D) simulations is conducted in § A.1, which shows that there are only minor differences between the two. For the 3-D simulations, the initial movement of the contact line is also restricted by another curvature, which makes the evolution slower at short times than in the 2-D case. Nevertheless, from the comparison, it does not affect the power law for the contact liquid bridge growth. Thus we focus the article on 2-D simulations using a 
$D2Q9$ lattice, as these give a good qualitative description of the flow.
Although numerous validation tests of the conservative phase field method have been conducted, we admit that the method still has drawbacks. For example, simulations of large interface deformations in the inertial regime may produce non-physical ghost droplets of oil along the water droplet's interface near the water/oil/air contact line. The reason for these ghost droplets might be pressure oscillations due to the inertial regime's low viscosities and large capillary forces. In addition, the initial conditions are such that there is a small overlap between the droplet and film interfaces, which causes inaccuracy in the surface tension calculation. As the interfaces evolve, this overlap will quickly disappear. Other possibilities include the enforcement of a constant mobility, which is a non-physical parameter in the conservative phase field method (Sun & Beckermann Reference Sun and Beckermann2007). In order to mitigate non-physical droplets that are not present in the data presented here, we consider an oil spreading factor 
$S_o=0$ in the inertial regime in our simulations, and use the experimental value 
$S_o=0.5\sigma _{wa}$ in the viscous regime; see table 1.
2.3. Conservative phase field lattice Boltzmann method
The velocity–pressure LBM (Baroudi & Lee Reference Baroudi and Lee2021) is utilized to solve the mass and momentum transport equations, and the conservative phase field LBM is employed to capture the interfaces in order to simulate the three-phase engulfment process in two and three dimensions. In the following equations, the order parameter 
$0\leqslant \phi _i\leqslant 1$ is used to separate different phases, where 
$\phi _i=1$ represents the bulk region of component 
$i$, and 
$0<\phi _i<1$ represents the interface region. In our simulation, we have three order parameters: 
$\phi _o$ for the oil film, 
$\phi _w$ for the water droplet, and 
$\phi _a$ for the air. Since the summation of three order parameters is 
$\sum _i\phi _i=1$, we here solve for only 
$\phi _w$ and 
$\phi _o$, and obtain 
$\phi _a$ from 
$\phi _a=1-\phi _o-\phi _w$. For order parameter 
$\phi _ i$, the discrete Boltzmann equation can be expressed as
\begin{equation} \left(\frac{\partial}{\partial t}+\boldsymbol{e}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}\right)h_\alpha^i={-}\frac{1}{\lambda_\phi}\,(h_\alpha^i-h_\alpha^{i,eq})+{S_{\alpha}^i}, \end{equation}
where 
$h_\alpha ^i$ represents the particle distribution function for order parameter 
$\phi _i$ along the 
$\boldsymbol {e}_\alpha$ direction. Here, 
$\boldsymbol {e}_\alpha$ denotes the discrete velocity vector in a 
$D2Q9$ or 
$D3Q27$ lattice (Lee & Lin Reference Lee and Lin2005), and 
$\lambda _\phi$ is the relaxation time. The equilibrium distribution function 
$h_\alpha ^{i,eq}$ can be calculated as
\begin{equation} h_\alpha^{i,eq}=\phi_i\varGamma_\alpha, \end{equation}
where 
$\varGamma _\alpha$ is the Maxwell–Boltzmann equilibrium distribution function presented by the method of Hermite polynomials and discretized by Gauss quadrature (He & Luo Reference He and Luo1997), and 
$\varGamma _\alpha$ takes the form
\begin{equation} \varGamma_\alpha=t_\alpha\left[1+\left(\frac{\boldsymbol{e}_\alpha \boldsymbol{\cdot}\boldsymbol{u}}{c_s^2}+\frac{\left(\boldsymbol{e}_\alpha \boldsymbol{\cdot}\boldsymbol{u}\right)^2}{2c_s^4}-\frac{\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{u}}{2c_s^2}\right)\right], \end{equation}
where 
$t_\alpha$ is the weight, 
$c_s$ is the speed of sound, and 
$\boldsymbol {u}$ is the macroscopic velocity vector of the flow field (He & Luo Reference He and Luo1997). The last term, 
$S_\alpha ^i$, is the intermolecular forcing term
\begin{equation} {S}_\alpha^i=\varGamma_\alpha(\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left(\frac{4 }{\delta}\,\frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\,\phi_i(1-\phi_i)- \frac{\phi_i^2}{\sum_{j=1}^3\phi_j^2}\sum_{j=1}^3\frac{4 }{\delta}\,\frac{\boldsymbol{\nabla}\phi_j}{|\boldsymbol{\nabla}\phi_j|}\,\phi_j(1-\phi_j)\right), \end{equation}
where 
$\delta$ denotes the finite interface thickness.
The discrete Boltzmann equation for mass and momentum transport can be expressed as
\begin{equation} \left(\frac{\partial }{\partial t}+\boldsymbol{e}_\alpha \boldsymbol{\cdot}\boldsymbol{\nabla}\right)g_\alpha ={-}\frac{1}{\lambda}\,(g_\alpha-g_\alpha^{eq})+F_\alpha, \end{equation}
where 
$g_\alpha$ represents the particle distribution function for the pressure, with the function
\begin{equation} g_\alpha^{eq}=t_\alpha\bar{p}+\varGamma_\alpha c_s^2 -t_\alpha c_s^2, \end{equation}
where 
$\bar {p}=P/\rho$. Here, 
$P$ represents the dynamic pressure, and 
$\rho$ is the mixture density. Also, 
$F_\alpha$ represents the forcing term
\begin{align} F_\alpha&={-}\varGamma_\alpha (\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left(\frac{1}{\rho}\,\boldsymbol{\nabla} P\right)+\varGamma_\alpha(0)\,(\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left(\boldsymbol{\nabla} \bar{p}\right)\nonumber\\ &\quad +\varGamma_\alpha (\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left[\frac{\nu}{\rho}\,( \boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{\rm T})\, \boldsymbol{\nabla}\rho+\frac{1}{\rho}\,\boldsymbol{F}_s\right], \end{align}
where 
$\nu$ is the kinematic viscosity, and 
$\boldsymbol {F}_s$ is the surface tension force. In our simulations, the continuum surface tension force model is employed (Kim Reference Kim2005):
\begin{equation} \boldsymbol{F}_{s}={-}\frac{3\delta}{2}\sum_i\sigma_i\,\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\right) |\boldsymbol{\nabla}\phi_i|\,\boldsymbol{\nabla}\phi_i, \end{equation}
where 
$\sigma _i=(\sigma _{ij}+\sigma _{ik}-\sigma _{jk})/2$. The gradient terms of the macroscopic values that appear above can be calculated by the second-order isotropic finite difference method (Lee & Lin Reference Lee and Lin2005).
The governing equations recovered through the Chapman–Enskog expansion are
$$\begin{gather} \frac{\partial \bar{p}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \bar{p}+ c_s^2\, \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$
$$\begin{gather}\frac{\partial
\boldsymbol{u}}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{uu})={-}\frac{1}{\rho}\,\boldsymbol{\nabla}
P+\frac{1}{\rho}\,\boldsymbol{\nabla}\boldsymbol{\cdot}
\eta\left(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^ {\rm T}\right)\nonumber\\
\quad -\frac{3\delta}{2\rho}\sum_{i=1}^{3}\sigma_i\,
\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\right)
|\boldsymbol{\nabla}\phi_i|\,\boldsymbol{\nabla}\phi_i,
\end{gather}$$
$$\begin{gather}\frac{\partial\phi_i}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\phi_i\boldsymbol{u})=
\boldsymbol{\nabla}\boldsymbol{\cdot}M\left(\vphantom{\frac{\phi_i^2}{\sum_{j=1}^3\phi_j^2}\sum_{j=1}^3}\boldsymbol{\nabla}\phi_i-\frac{4}{\delta}\,
\frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\,\phi_i(1-\phi_i)\right.\nonumber\\
\quad \left.+\frac{\phi_i^2}{\sum_{j=1}^3\phi_j^2}\sum_{j=1}^3\frac{4}{\delta}\,\frac{
\boldsymbol{\nabla}\phi_j}{|\boldsymbol{\nabla}\phi_j|}\,\phi_j(1-\phi_j)\right),
\end{gather}$$and the density and viscosity defined for the three phases are
$$\begin{gather} \rho=\sum_{i=1}^3\rho_i\phi_i, \end{gather}$$
$$\begin{gather}\eta=\sum_{i=1}^3\eta_i\phi_i. \end{gather}$$At the boundaries of the domain, we apply symmetric boundary conditions to the left and right boundaries, while the no-slip boundary condition is applied to the top and bottom boundaries. The details of the LBM derivation and Chapman–Enskog expansion can be found in Zhao & Lee (Reference Zhao and Lee2023).
3. Results and discussion
3.1. Inertial-capillary drop engulfment 
$Oh<1$
Figure 2 shows the 2-D simulation results of an aqueous drop contacting an oil film in the inertial-capillary regime, 
$Oh \ll 1$, for varying film thicknesses 
$H/R$. Immediately post-contact, the interface near the film–drop–air contact line deforms, creating a liquid bridge with a large curvature tangential to the wall. This curvature drives the radial expansion of the liquid bridge's neck. As the neck expands, the drop becomes engulfed by the oil film and either is supported from below by the wall with an interface reminiscent of a spherical cap, or sinks into the film, adopting a circular interface.
Figure 2. Three snapshots in time of simulations of an aqueous drop being engulfed, in each case, by a thin 
$H/R=0.1$ oil film on the left, and a thick 
$H/R=4$ oil film on the right, in the inertial regime, 
$Oh=0.07$. The contours indicate the three different phases: air (colourless), droplet phase (red) and oil film (blue). The inset shows the zoomed-in contours of the diffuse interface for three order parameters. Black dashed lines indicate the initial oil film thickness. Post drop–film contact, the drop deforms and oscillates, and the drop is quickly engulfed by the oil. As the film's thickness increases, capillary waves are formed on the film's interface, affecting the drop's dynamics.
Evident from figure 2 with 
$Oh=0.07$ are the striking morphological variations post drop–film contact in the inertial-capillary regime. Here, the contours of the three different liquid phases are shown for 
$H/R=0.1$ and 
$4$. We notice that in the short time 
$T/t_{\rho }=0.08$, the drop's interfacial shapes for each film thickness 
$H$ appear identical, and a large radial velocity is generated near the three-phase line as the oil spreads onto the drop. This radial velocity is due to the region of large curvature generated by the oil film as it spreads on the drop, and the subsequent capillary force pulls the drop in the direction tangential to the wall, reminiscent of a drop spreading on a solid surface. This rapid morphological deformation of the drop causes its interface to deform significantly and oscillate, a result of the lack of viscous resistance to motion in the inertial-capillary regime. In a relatively short time, the oil completely engulfs the drop for the thick film, inducing a capillary force that pushes the drop towards the substrate. A squeeze flow in the gap between the drop interface and the wall is generated, with a build up of pressure that is naturally a function of the gap height (figure 2b). The difference in the flow-induced pressure underneath the drop significantly affects the drop's interfacial dynamics for 
$T/t_{\rho }=0.8\unicode{x2013}1.6$. For thin films, 
$H/R=0.1$, the drop extends along the wall, and for thick films, 
$H/R=4$, the drop is transported into the film with capillary waves formed on the film's interface. For thin films at long times, not shown here, the drop adopts a spherical cap shape, completely covered by the oil, with a wedge formed at the drop's most radial point (Style & Dufresne Reference Style and Dufresne2012; Park et al. Reference Park, Weon, Lee, Lee, Kim and Je2014). In contrast, for thick films, the drop is transported into the film and relaxes into a circular drop below a flat oil interface. As there are no buoyancy effects in the simulations, only viscous resistance will halt the drop's motion, which in the inertial regime is small.
3.1.1. Droplet centre of mass motion and liquid bridge radial growth
To characterize the drop's engulfment dynamics, we extract the position of the centre of mass 
$C_m$ in the wall-normal direction relative to the position of the initial film–air interface and the initial drop radius 
$R$ from the simulations. Here, the centre of mass 
$C_m$ is used as a measure for the drop's motion into the oil film. In figure 3, we compare the evolution of the drop's centre of mass 
$C_m$ for oil films of different thicknesses. Initially, there appears to exist a short-time regime that is fairly independent of the film height for 
$T/t_\rho <2$. This follows a similar curve and indicates where it is the outward radial growth of the liquid bridge formed between the drop and the film that drives the displacement of the drop's centre of mass. At long times, the position of the centre of mass 
$C_m$ is highly affected by the ratio 
$H/R$, with larger film thickness 
$H$ allowing for larger wall-normal displacements.
Figure 3. The scaled centre of mass (
$C_m/2R$) position (wall-normal) of the drop as a function of time post-contact with an oil film of varying thickness (
$H/R=[0.1\unicode{x2013}4]$), for 
$Oh=0.07$. The large inertia causes an oscillatory motion of the drop's centre of mass for thin films with oscillations damped when 
$H/R\approx 4$.
A wave-like motion of 
$C_m$ appears for 
$H/R\leqslant 1$ in figure 3, which is formed by an inertial-capillary wave caused by the imbalance between the drop's initial shape and its equilibrium state (Ding et al. Reference Ding, Li, Zhang, Sui, Spelt and Thoroddsen2012). The inertial-capillary wave motion can also be seen in figure 2 when 
$H/R=0.1$ as the drop oscillates up and down between times 
$T/t_{\rho }=0.8$ and 
$T/t_{\rho }=1.6$. This wave-like motion becomes less pronounced as the film thickness increases, and when the film is sufficiently thick, 
$H/R>4$, no inertial-capillary waves are formed. The pressure built up in the film by the squeeze flow between the drop and the wall ultimately arrests the drop at different positions in the film as a function of 
$H/R$.
To describe the dynamics, we also extracted the radius of the neck of the growing liquid bridge 
$r$ from the simulations. In the inertial-capillary regime, the log-log representation of the neck radius 
$r$ for 
$Oh=[0.01\unicode{x2013}0.07]$ and 
$H/R = 1$ is shown in figure 4(a), and for 
$Oh=0.07$ and 
$H/R= [0.1\unicode{x2013}4]$ in figure 4(b). Here, we find that the outward radial growth of the liquid bridge is insensitive to the film height. Further, the simulations of figure 4 indicate that 
$r$ has a power-law dependence with time as 
$r/R\sim ({T/t_\rho })^{1/2}$, which is consistent with miscible inertial drop–film coalescence and inertially limited engulfment of an aqueous drop on a very thin oil film (Eggers et al. Reference Eggers, Lister and Stone1999; Carlson et al. Reference Carlson, Kim, Amberg and Stone2013). The differences between the simulations for various film thicknesses 
$H/R=[0.1\unicode{x2013}4]$ appear to be present only at long times and are associated with the drop being transported by capillarity into the film.
Figure 4. The radius of the growing liquid bridge 
$r/R$ plotted against time 
$T/t_{\rho }=[0\unicode{x2013}0.2]$ for (a) 
$Oh=[0.01\unicode{x2013}0.07]$ and 
$H/R=1$, and (b) 
$Oh=0.07$ and 
$H/R=[0.1\unicode{x2013}4]$ in logarithmic axis. The inset shows the plot for 
$Oh=0.07$ with different oil film thicknesses (linear axis). The oil film thickness effect is very small. Reference triangles indicate slope 
$\beta =0.5$. Regardless of 
$Oh$ and 
$H/R$, it can be seen that 
$r/R\sim (T/t_\rho )^{1/2}$, analogous to the coalescence of two drops in the inertial regime.
3.2. Visco-capillary drop engulfment, $Oh>1$
Next, we move on to study the less understood visco-capillary regime 
$Oh>1$. Figure 5 shows the 2-D simulation results of an aqueous drop contacting an oil film in the visco-capillary regime for 
$Oh = 5.27$ and two different film thicknesses, 
$H/R=0.1$ and 
$H/R=0.4$. Post drop–film contact, the velocity vectors mostly appear near the contact region, indicating that at short times, the region of high curvature formed upon contact drives the radial growth of the liquid bridge. Compared to the 
$Oh=0.07$ case, the drop's apparent spreading is slower, and the engulfment process is free of any inertial-capillary waves. Despite the fact that the capillary driving force is the same, any waves are effectively damped by the viscosity (Cuttle et al. Reference Cuttle, Thompson, Pihler-Puzović and Juel2021). Comparing the simulation for different film thicknesses at long times shows that the drop can be fully engulfed for thick films, whereas when the oil film is very thin, the velocity inside the film becomes small and the simulations are unable to reach the fully engulfed state expected at equilibrium.
Figure 5. Simulations of an aqueous drop being engulfed, in each of three cases, by a thin 
$H/R=0.1$ oil film on the left, and a thick 
$H/R=4$ oil film on the right, in the viscous regime, 
$Oh=5.27$. The inset shows the diffuse interface contours for three order parameters. Black dashed lines indicate the initial oil film thickness 
$H$. Immediately post-contact, the drop and film form a liquid bridge whose neck grows radially. It can be seen that regardless of film thickness, the radial growth of the neck occurs at the same rate, with similar drop morphologies at each time step.
3.2.1. Droplet centre of mass motion and liquid bridge radial growth
In figure 6, we compare the evolution of the drop's centre of mass 
$C_m$ in the visco-capillary regime (
$Oh>1$) for oil films of different thicknesses. As in the inertial-capillary regime, the visco-capillary regime suggests a region of universal short-time dynamics that is fairly independent of the film height for 
$T/t_\eta <0.7$, suggesting that the growth of the liquid bridge formed between the drop and the film is driving the displacement of the drop's centre of mass. At long times, the position of the centre of mass 
$C_m$ is highly affected by the ratio 
$H/R$, with larger film thicknesses 
$H/R$ allowing for larger wall-normal displacements. Variations in the thickness of the film show a diminishing effect on the engulfment process other than for the final position of the centre of mass 
$C_m$. Here, no oscillations of the drop's centre of mass are present. After complete engulfment, the drop is arrested by the viscous force and reaches equilibrium.
Figure 6. The scaled centre of mass (
$C_m/2R$) position (wall-normal) of the drop as a function of time post-contact with an oil film of varying thickness (
$H/R=[0.1\unicode{x2013}4]$), for 
$Oh=5.27$. Liquid bridge growth ends at 
$T/t_\eta =0.6$. As 
$H/R$ increases, the drop is able to sink further into the film.
Figure 7 shows the log-log representation of the growth of the liquid bridge's radius 
$r$ for 
$Oh=[1\unicode{x2013}40]$. In figure 7(a), we see that the radius 
$r$ has a power-law dependence in time with an exponent that decreases with increasing 
$Oh$. This variation in the power-law exponent is a consequence of being in a transitional regime 
$Oh\approx 1$, where there can still be some inertial effects (Paulsen et al. Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012). However, as 
$Oh$ continues to increase, the exponent of the power law appears to approach a constant value as 
$r/R\sim (T/t_\eta )^{2/5}$, which is comparable to the previous observations for aqueous drop engulfment by a thin viscous oil film (Carlson et al. Reference Carlson, Kim, Amberg and Stone2013). As in the inertial regime, we can see that the short-time dynamics is insensitive to the pre-wetted film height, as shown in figure 7(b), with a power-law dependence 
$r/R\sim (T/t_\eta )^{3/5}$ for 
$Oh = 10.54$.
Figure 7. The radius of the growing liquid bridge 
$r/R$ against time 
$T/t_\eta =[0\unicode{x2013}1]$ for (a) 
$Oh=[ 2.35 , 5.27, 10.54 , 20, 30, 40]$ with film thickness 
$H/R=0.2$, and (b) 
$Oh=10.54$ with various film thicknesses 
$H/R=[0.1\unicode{x2013}4]$. Reference triangles highlight the power-law 
$\beta$ dependence with time. In (a), 
$\beta$ converges to 
$2/5$ from 
$1$ as 
$Oh$ increases from 
$Oh=2.35$ to 
$Oh=40$. In (b), 
$\beta$ is invariant to the film thickness, 
$H/R = [0.1\unicode{x2013}4]$, for constant 
$Oh$.
To better understand the invariance of the power-law dependence of the neck radius in time on the film thickness, we compare the dissipation over time when 
$Oh=5.27$ for different 
$H/R$. Here, the viscous dissipation 
$\varGamma _\eta$ of the whole simulation domain can be formulated as
\begin{equation} \varGamma_\eta=\int_{\varOmega} \eta(\phi)\,\boldsymbol{\nabla} u:(\boldsymbol{\nabla} u +\boldsymbol{\nabla}^{\rm T} u)\,{\rm d}\varOmega. \end{equation}
It can then be scaled by the surface energy 
$E_s=\sigma _{ow} U R$, where the capillary speed 
$U=\sigma _{ow}/\eta _w$ is considered as the reference velocity (Carlson et al. Reference Carlson, Do-Quang and Amberg2011). The dimensionless form of the viscous dissipation is denoted 
$\varGamma _\eta ^*=\varGamma _\eta /E_s$.
Figure 8 shows the viscous dissipation per volume for thin 
$H/R = 0.1$ and thick 
$H/R = 4$ films at times in the range 
$T/t_{\eta } = [0.04 \unicode{x2013}0.6]$ for 
$Oh = 5.27$. It is clear that during the initial rapid radial expansion of the liquid bridge, the majority of the dissipation is centred around a small area in the region of high curvature near the triple contact point, regardless of the film thickness. This invariance to film thickness in the profile of the viscous dissipation supports the prior observation of the invariance in the power-law exponent for varying film thicknesses; see figure 7(b).
Figure 8. Contour plot of the viscous dissipation 
$\varGamma _\eta ^*$ of a thin film 
$H/R=0.1$ (a), and a thick film 
$H/R=4$ (b), for 
$T/t_\eta =0.04$, 
$0.2$, 
$0.6$ and 
$Oh=5.27$. The three black lines show 
$\phi _i=0.5$ for the three order parameters. The density and viscosity are interpolated by the linear functions given in (2.12) and (2.13) near the contact line. The viscous dissipation is localized at the contact line for thin and thick oil films, and gradually decreases.
Figure 9(a) shows the rate of the total viscous dissipation for varying film thicknesses 
$H/R = [0.1\unicode{x2013}4]$ for 
$Oh = 5.27$. Here, we see that variations in the film thickness for 
$T/t_{\eta }<1$ have little effect on the system's total viscous dissipation, with deviations mainly when 
$T/t_{\eta }>1$ as the viscosity of the oil is larger than that of the air, supporting again the invariance of the power-law exponent of figure 7(b) to film thickness. Additionally, the first peak ends at 
$T/t_\eta \approx 1$, which is consistent with the maximum radial growth of the liquid bridge. It can be seen that a thinner film will overall have less viscous dissipation, indicating that the flow structure at the contact line is smaller than the pre-wetted film thickness. Figure 9(b) shows the rate of the scaled viscous dissipation for various 
$Oh$ with 
$H/R = 0.2$. When 
$Oh$ is large, i.e. 
$Oh>10$, we see that the rate of viscous dissipation over time converges to a single trend and becomes insensitive to 
$Oh$, parallelling the convergence of the power-law exponent to 
$2/5$ in figure 7(a). When the Ohnesorge number is in the transition regime, i.e. 
$Oh=[1\unicode{x2013}10]$, the rate of dissipation during liquid bridge growth varies with 
$Oh$, resulting in the varying power-law exponents of figure 7(a).
Figure 9. (a) The rate of viscous dissipation as a function of time for 
$Oh=5.27$ with 
$H/R=[0.1\unicode{x2013}4]$. (b) The rate of viscous dissipation as a function of time for 
$Oh=[2.35\unicode{x2013}40]$ and 
$H/R = 0.2$. (a) The viscous dissipation first rapidly increases as the post-contact liquid bridge formed between the aqueous drop and the oil film grows due to the large initial curvature near the three-phase line and gradually dissipates over time as viscosity diffuses the energy. The film thickness shows little effect on the rate of dissipation for a constant 
$Oh$. (b) The viscous dissipation is depressed for 
$Oh\approx 1$ and converges to a similar trend as 
$Oh$ approaches 
$40$. This variation in dissipation over time for increasing 
$Oh$ accounts for the changing power-law evolution of the growing liquid bridge from 
$\beta =4/5$ to 
$\beta =2/5$.
Figure 10 shows the variation of the drop's centre of mass against the radius of the growing liquid bridge in the visco-capillary regime for 
$Oh = 5.27$. It can be seen that the initial evolution of the liquid bridge's radius, defined as rapid first-stage engulfment in Cuttle et al. (Reference Cuttle, Thompson, Pihler-Puzović and Juel2021), is quite similar for different film thicknesses. At later times, however, when 
$T/t_\eta >0.6$ and 
$r/R>1$, the differences in the engulfment dynamics for varying film thicknesses become obvious, with the centre of mass being pulled into the film further and faster as the film thickness increases. The large curvature of the film's interface coupled with the large viscosity and localized viscous dissipation explain the invariance to film thickness during this initial regime, 
$T/t_\eta <0.6$. This invariance to film thickness during the initial regime further explains the invariance in the mass centre evolution to film thickness, as it is the rapid radial expansion of the growing liquid bridge that initially pulls the drop's centre of mass towards the film. As time progresses, however, the pressure difference becomes the only effect driving engulfment, allowing the centre of mass to be pulled further into the film as the film thickness increases.
Figure 10. The centre of mass (wall-normal direction) 
$C_m/2R$ plotted against the radius of the growing liquid bridge 
$r/R$ for 
$Oh=5.27$ and 
$H/R=[0.1\unicode{x2013}4]$. Immediately post-contact during the initial rapid radial growth of the drop–film liquid bridge 
$T/t_\eta <0.6$, the position of the drop's centre of mass with respect to 
$r/R$ is quite similar for varying film thicknesses 
$H/R$. After 
$T/t_\eta =0.6$, the liquid bridge has reached its maximum extent, and the drop begins sinking into the film. Here, the effect of the film thickness becomes apparent as it affects the pressure difference between the air and the oil film.
3.3. Experimental considerations
In order to validate the simulation results for visco-capillary drop engulfment from § 3.2, we consider experiments where 
$Oh>1$. Figure 11 shows the results of a glycerol/water drop of radius 
$R \approx 1$ mm contacting a viscous oil film of thickness 
$H/R = 0.2$. It can be seen qualitatively from figure 11 for 
$Oh > 10$ that the 2-D simulations capture well the shape of the liquid–gas interface during the initial rapid radial growth of the post-contact liquid bridge. Discrepancies between the simulations and the experiments may arise, though, from either the experimental drop's attachment to the needle tip or the simulation's lower dimensionality. Figures 11(a,b) show an obvious connection between the drop and the needle, which does influence the interfacial shape of the upper region of the drop. This influence, however, affects primarily the upper portion of the drop, leaving the region of high curvature near the three-phase line unaffected. Additionally, in § A.1, we perform the associated 3-D simulations to test the dimensional dependency of the simulation, and show that the initial radial growth of the liquid bridge can be reduced to a 2-D problem. With these considerations in mind, we can compare meaningfully the results of the experiments with the 2-D simulations shown previously.
Figure 11. Comparison between the simulation and the experiment for 
$H/R=0.2$. Simulation: 
$Oh =40$. Experiment: 
$Oh = 57$, 
$\eta _o = 7.6\ {\rm Pa}\ {\rm s}$, 
$\eta _w = 0.705\ {\rm Pa}\ {\rm s}$. Snapshots are taken at similar radii 
$r/R$. The simulation captures well the morphology of the drop's interface during the initial post-contact radial expansion of the liquid bridge.
Figures 12(a,b) show the scaled liquid bridge radius 
$r/R$ against the scaled time 
$T/t_\eta$ extracted from the experiments for varying 
$Oh=[2.51\unicode{x2013}57]$ and varying film thicknesses 
$H/R=[0.2\unicode{x2013}4.43]$, respectively. When the film thickness is fixed at 
$H/R=0.2$, it can be seen from figure 12(a) that increases in the Ohnesorge number 
$Oh$ lead to decreases in the power-law exponent from 
$4/5$ to 
$2/5$, mirroring the simulation results shown previously in figure 7(a). Exponents are fitted from the data for 
$r/R$ in the range 
$0.4$–
$1$. Additionally, figure 12(b) shows the effect of increasing film thicknesses 
$H/R=[0.2\unicode{x2013}4]$ for 
$Oh \approx 11$, and we see that the evolution of the liquid bridge's radius can be described by the same power-law exponent, 
$r \sim t^{3/5}$, as shown by the simulations of figure 7(b). We find that varying 
$H/R$ has no effect on the evolution of the growing liquid bridge, supporting the previous conclusion that the radial growth is affected only by the small region of high curvature near the three-phase line.
Figure 12. Experimentally measured liquid bridge radii 
$r/R$ against time 
$T/t_\eta$ for immiscible glycerol/water drops (
$\blacktriangle,\blacktriangleleft,\blacktriangledown,\bullet,\blacksquare$) and miscible silicone oil drops (
$\blacklozenge$) coalescing with oil films for (a) 
$Oh = [2.5\unicode{x2013}57]$ with fixed film thickness 
$H/R = 0.2$, and (b) 
$Oh \approx 11$ with 
$H/R = [0.2\unicode{x2013}4.4]$, shown in logarithmic axes. Drop–film viscosity ratios are held at 
$\eta _w/\eta _o=1/10$. (a) Increases in 
$Oh$ for constant 
$H/R$ are shown to decrease the power-law slope from 
$4/5$ to 
$2/5$, following a similar convergence trend as shown in the simulations. The decrease in velocity of 
$r/R$ in time for 
$Oh=57$ near 
$T/t_\eta = 0.07$ is thought be a result of the needle tip influencing the curvature of the pre-detached drop. (b) Increases in 
$H/R$ for constant 
$Oh$ are shown to have no effect on the power-law slope, as predicted from the simulations. Comparison with miscible oil drop and oil film coalescence for similar viscosities as in the immiscible case recovers the same power-law slope 
$\beta \approx 3/5$. Here, the local curvature of the oil film near the contact point is thought to drive the liquid bridge growth, limited only by viscous dissipation.
Figure 12(b) also shows the results of a miscible oil drop coalescing with an oil film. Interestingly, we see that when the viscosity ratio between the drop and the film is the same as in the immiscible case, the temporal evolution of the liquid bridge's radius formed during this miscible coalescence evolves with the same power-law exponent. As was shown previously by Aarts et al. (Reference Aarts, Lekkerkerker, Guo, Wegdam and Bonn2005) for viscously limited drop–film coalescence, the coalescence velocity is set by the capillary velocity 
$\sigma /\eta$, driven by the interfacial curvature and limited by the viscosity. Here, the temporal collapse of the interfacial dynamics for the miscible drop–film coalescence and immiscible drop–film engulfment can be explained if the relevant capillary velocity is 
$\sigma _{oa}/\eta _o$, and the limiting viscosities of the film, the drop and the air are comparable between the miscible and immiscible cases. As the surface tension of the oil is at least half that of the drop in the immiscible case, this collapse shows that it is the curvature of the oil film in the neck region and subsequently the film's surface tension that is relevant to the radial growth of the liquid bridge.
4. Conclusion
Measuring the flow field of a liquid spreading on another liquid is experimentally challenging, whereas via simulations, the internal velocity field can be extracted easily, rendering trivial calculations of the viscous dissipation as well as the extraction of the drop's centre of mass, the interfacial boundaries, and the location of the three-phase line. In this paper, we utilized both simulations and experiments to investigate the effect of the film thickness and the Ohnesorge number 
$Oh$ on the growth of the post-contact drop–film liquid bridge and the engulfment of the aqueous drop by the oil film in both the viscously limited 
$Oh>1$ and inertially limited 
$Oh \ll 1$ regimes. Here, we simulated the three-phase flow by using the conservative phase field equations (Chiu & Lin Reference Chiu and Lin2011; Geier, Fakhari & Lee Reference Geier, Fakhari and Lee2015) solved by the LBM (He & Luo Reference He and Luo1997; Mohamad Reference Mohamad2011), which has recently been proposed to have improved mass conservative properties (Zheng et al. Reference Zheng, Lee, Guo and Rumschitzki2014; Geier et al. Reference Geier, Fakhari and Lee2015; Aihara, Takaki & Takada Reference Aihara, Takaki and Takada2019; Baroudi & Lee Reference Baroudi and Lee2021).
In the inertial-capillary regime, 
$Oh<1$, we found that post drop–film contact for film thicknesses 
$H/R<1$, the drop's centre of mass would oscillate during engulfment, and for film thicknesses 
$H/R>1$, inertial-capillary waves on the surface of the oil film would be produced. In the visco-capillary regime, 
$Oh>1$, we found that decreases in the film thickness led to increases in the time required for the drop to become fully engulfed.
During the post-contact rapid radial growth of the drop–film liquid bridge, we found that regardless of 
$Oh$ or film thickness, the growth of the radius of the post-contact liquid bridge 
$r$ evolved in time with a power-law scaling 
$r/R\sim (T/t_\rho )^{1/2}$ in the inertial-capillary regime, consistent with prior work for drops contacting thin oil films (Carlson et al. Reference Carlson, Kim, Amberg and Stone2013). In the visco-capillary regime, we similarly saw that the power-law growth of the liquid bridge was invariant to the film thickness; however, variations in 
$Oh$ for 
$Oh = [1\unicode{x2013}10]$ varied the power-law evolution of the liquid bridge's neck 
$r/R\sim (T/t_\eta )^{\beta }$, with 
$\beta$ ranging from 
$4/5$ for 
$Oh \approx 1$ to 
$2/5$ for 
$Oh>10$. We rationalize that the invariance to film thickness on this power-law growth is a result of the viscous dissipation being initially concentrated locally at the contact line near the growing liquid bridge.
We compared our simulation results to experiments in the visco-capillary regime, validating both the qualitative shape of the drop post-contact as well as the power-law evolution of the growing liquid bridge for a similar range of 
$Oh$. Interestingly, we found that the miscible coalescence between an oil drop and an oil film with the same viscosity ratio as the immiscible aqueous drop and oil film case recovered the same temporal evolution of the liquid bridge's growth. We rationalize that this collapse of the miscible and immiscible cases is a consequence of the similar curvatures of the oil film near the growing liquid bridge coupled with the similar resistances to growth of the liquid phase's viscosities.
Wetting on soft materials, such as the SLIPS surfaces considered here, is becoming increasingly relevant to the fabrication of self-healing, low-hysteresis, water-repellent surfaces. However, drop engulfment as well as many other soft wetting phenomena is far from trivial when considered experimentally. In this work, we have shown that simulations utilizing the LBM are able to capture accurately the drop engulfment dynamics during soft wetting on SLIPS surfaces, opening the door for future research probing increasingly complicated systems with various rheological properties and geometries, free from the challenges of physical observation.
Funding
We acknowledge the financial support of the Research Council of Norway through the program NANO2021 (project no. 301138) and the PIRE project ‘Multi-scale, Multi-phase Phenomena in Complex Fluids for the Energy Industries’, founded by the Research Council of Norway and the National Science Foundation of USA under award no. 1743794. This research was supported, in part, under National Science Foundation Grants CNS-0958379, CNS-0855217, ACI-1126113 and OEC-2215760 (2022), and the City University of New York High Performance Computing Center at the College of Staten Island. The computations were also performed on resources provided by Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway.
Declaration of interests
The authors report no conflict of interest.
Appendix A
A.1. Comparison between 2-D and 3-D simulations at drop engulfment dynamics
The short-time contact-line dynamics on solid or coalescence of drops has been shown to be described by 2-D theory (Eggers et al. Reference Eggers, Lister and Stone1999; Eddi et al. Reference Eddi, Winkels and Snoeijer2013). To verify that our 2-D simulations capture well the drop dynamics for this soft wetting problem, we have conducted both 2-D and 3-D simulations of drops on the thin oil film with 
$H/R=0.03$, and compare their results.
Three-dimensional simulations of the drop evolution on the film with 
$H/R=0.03$ for 
$Oh=0.2$ and 
$Oh=4.8$ are shown in figures 13 and 14. The log-log representations of 2-D and 3-D drop spreading radius are shown in figures 15(a) and 15(b), which show only a very minor effect of the 3-D flow, and both predict the radius to follow 
$r/R\sim (T/t_c)^{1/2}$.
Figure 13. Three-dimensional simulation of the drop dynamics (
$H=0.03R$) at (a) 
$T/t_\rho =0.01$,(b) 
$T/t_\rho =0.05$, (c) 
$T/t_\rho =0.1$ and (d) 
$T/t_\rho =0.15$ for 
$Oh=0.2$.
Figure 14. Three-dimensional simulation of the drop dynamics (
$H/R=0.03$) at (a) 
$T/t_\eta =0.01$,(b) 
$T/t_\eta =0.05$, (c) 
$T/t_\eta =0.1$ and (d) 
$T/t_\eta =0.15$ for 
$Oh=4.8$.
Figure 15. Log-log representation of the water drop's apparent spreading radius with oil film thickness 
$H/R=0.03$ for (a) 
$Oh=0.2$ and (b) 
$Oh=4.8$. Both simulation results indicate the power law 
$r/R\sim (T/t_c)^{1/2}$.
A.2. Materials
Glass slides (Avantor, Article no. 631-1550P, 
$76\ {\rm mm}\times 26\ {\rm mm}$), 1000 St silicone oil VIS-RT100K (European Article no. 99515.271), 20 cSt silicone oil (CAS no. 63148-62-9) and Glycerine (CAS no. 56-81-5) were purchased from VWR International (Oslo, Norway), and 30 gauge syringe tips (Weller, Mfr. Part no. KDS3012P) were purchased from RS components (Oslo, Norway).
A.3. Convergence test
We conduct simulations for the inertial regime and the viscous regime. We use 
$Cn=0.02$ in simulations of the inertial regime where 
$Oh<1$ for which the interface thickness is 
$\delta /\Delta x=4$, and 
$R/\Delta x=200$. In the viscous regime for 
$Oh>1$, 
$Cn=0.04$, 
$\delta /\Delta x=4$ and 
$R/\Delta x=100$ are used. The reasons why we choose a larger 
$Cn$ for the viscous regime are: (i) the shape of the droplet deforms a lot in the inertial regime, for which we need to increase the resolution to collect a correct result; (ii) the simulations of the viscous regime have a significant time consumption, for which we anticipate dozens of times higher time consumption as we increase resolution. We present here the convergence tests for each regime's simulation.
In figure 16, we present the comparison between 
$Cn=0.01$ and 
$Cn=0.02$ for 
$Oh=0.07$. In the test, we fix the interface thickness as 
$\delta /\Delta x =4$, and increase the radius of the droplet from 
$R/\Delta x=200$ to 
$R/\Delta x=400$. Since we increase the radius of the droplet but keep the same interface thickness, we change the initial set-up to keep the simulation consistent. This modification will induce a time shift between two simulations. From the comparison, the evolutions of the contact lines are highly overlapped.
Figure 16. Temporal evolution of the water/oil/air interfacial profiles for 
$Cn=0. 02$ (black) and 
$Cn=0.01$ (red) with 
$Oh=0.07$. The interfaces are defined as the contours of the three different components 
$i$ when 
$\phi _i=0.5$. The times shown below each panel are for 
$Cn=0.01$, and the time shift between the two simulations is 
$0.02t_\rho$, which means that for 
$Cn=0.02$, the times are 
$T/t_\rho =[0.04,0.12,0.22,0.32]$, respectively.
We then show the convergence test for 
$Oh=10.54$ between 
$Cn=0.02$ and 
$Cn=0.04$ in figure 17. As with the convergence test for the inertial regime, we change the distance between the droplet and the film surface. The simulation results are remarkably similar.
Figure 17. Temporal evolution of the water/oil/air interfacial profiles for 
$Cn=0.04$ (black) and 
$Cn=0.02$ (red) with 
$Oh=10.54$. The interfaces are defined as the contours of three different components 
$i$ when 
$\phi _i=0.5$. The times shown below each panel are for 
$Cn=0.02$, and the time shift between the two simulations is 
$0.0114t_\eta$, which means that for 
$Cn=0.04$, the times are 
$T/t_\eta =[0.05,0.21,0.41,0.61]$, respectively.
We further consider the effect of the initial condition for the simulation. In figure 18, we compare three different initial conditions: the centre of mass is placed at 
$C_m=H+R+\delta$ (blue), 
$C_m=H+R+0.5\delta$ (black) and 
$C_m=H+R$ (red). We notice that the closer we initialize the droplet to the film surface, the quicker the spreading process. Nevertheless, this modification of the initial condition does not affect the scaling of the spreading process.
Figure 18. Comparison of three different initial conditions, 
$C_m=H+R+\delta$ (blue), 
$C_m=H+R+0.5\delta$ (black) and 
$C_m=H+R$ (red), for 
$Cn=0.04$ and 
$Oh=10.54$. The times shown are for 
$C_m=H+R$. The time shift is 
$0.02t_\eta$, which means that for 
$C_m=H+R+0.5\delta$, the exact times are 
$T/t_\eta =[0.06,0.22,0.42]$, and for 
$C_m=H+R+\delta$, the exact times are 
$T/t_\eta =[0.08,0.24,0.44]$.
Based on the above convergence tests, the difference for the inertial regime between 
$Cn=0.01$ and 
$Cn=0.02$, and for the viscous regime between 
$Cn=0.02$ and 
$Cn=0.04$, is negligible when we consider a small time shift. This time shift is induced only by the initial condition and will not affect the analysis in the paper.
Open access
