1. Introduction
The occurrence of drop-interface coalescence has been observed in a wide range of natural phenomena and industrial applications, such as rain/cloud formation (Raes et al. Reference Raes, Van Dingenen, Vignati, Wilson, Putaud, Seinfeld and Adams2000), atomisation (Villermaux Reference Villermaux2007) and also emulsification or de-emulsification processes (Ziegler & Wolf Reference Ziegler and Wolf2005). Over half a century ago, Charles & Mason (Reference Charles and Mason1960) observed coalescence in their ground-breaking experiments, and, ever since, researchers have been in constant pursuit of a better physical understanding of this phenomenon. However, it was not until the advent of high-speed imaging that it became possible for Thoroddsen & Takehara (Reference Thoroddsen and Takehara2000) to observe the self-similar coalescence cascade phenomena of a drop before its total coalescence. Since then, the significant interest in the field has led to a recent comprehensive review on the topic by Kavehpour (Reference Kavehpour2015), who concluded that further work is needed to understand the Marangoni effect during the drop-interface coalescence dynamics.
The dynamics commences with the drainage of the fluid between the drop and an interface separating this fluid from another bulk phase whose extent is typically much larger than the drop diameter. This drainage leads to the formation of a thin fluid layer between the drop and the interface. As the layer thickness decreases, van der Waals forces trigger its rupture with the generation of a hole which expands driven by capillarity. The hole expansion for drop coalescence has been widely studied by Eggers, Lister & Stone (Reference Eggers, Lister and Stone1999), Aarts et al. (Reference Aarts, Lekkerkerker, Guo, Wegdam and Bonn2005), Paulsen, Burton & Nagel (Reference Paulsen, Burton and Nagel2011), Paulsen et al. (Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012), Paulsen (Reference Paulsen2013), Paulsen et al. (Reference Paulsen, Carmigniani, Kannan, Burton and Nagel2014), Anthony et al. (Reference Anthony, Kamat, Thete, Munro, Lister, Harris and Basaran2017), Anthony, Harris & Basaran (Reference Anthony, Harris and Basaran2020) among others, concluding that the interfacial dynamics is solely governed by a balance between inertial, viscous and surface-tension forces, and subsequently the Ohnesorge number 
$Oh$ (i.e. ratio of viscous to capillary forces) is the most appropriate control parameter for this phenomenon. Different coalescence regimes have been identified depending on the order of magnitude of 
$Oh$: (i) the inertial regime (
$Oh \ll 1$), which is characterised by a nearly inviscid liquid, and the dynamics is surface-tension driven; (ii) the Stokes regime (
$Oh >1$) where the viscous forces play a major role in the interfacial dynamics; finally, (iii) an intermediate regime, the inertial-limited-viscous regime, which bridges the inertial and Stokes Regimes (i.e. no dominance by either viscosity or surface tension).
As pointed out earlier for intermediate values of 
$Oh$, Thoroddsen & Takehara (Reference Thoroddsen and Takehara2000) observed the so-called ‘coalescence cascade of a drop’ in which the drop coalescence leads to the generation of a smaller daughter droplet which results in a cascade of self-similar events until this successive coalescence process is completed. The process of formation of a daughter droplet is known as a ‘partial coalescence’ phenomenon, and its physical understanding came from the insightful experimental and numerical results of Blanchette & Bigioni (Reference Blanchette and Bigioni2006), who suggested that the occurrence of pinch-off depends solely on the competition between the vertical (inertia–viscous) and horizontal (capillary) pulls (the former aids the total coalescence and the latter the capillary breakup), rather than the mechanism of Rayleigh–Plateau instability. Additionally Blanchette & Bigioni (Reference Blanchette and Bigioni2006) have provided an extensive 
$Bo$-
$Oh$ phase diagram delineating the boundaries between partial and total coalescence. Here, 
$Bo$ is the Bond number which compares the importance of gravitational to surface-tension forces.
Importantly, Blanchette, Messio & Bush (Reference Blanchette, Messio and Bush2009), Thoroddsen et al. (Reference Thoroddsen, Qian, Etoh and Takehara2007) and Sun et al. (Reference Sun, Zhang, Che and Wang2018) have also considered situations in which there is a surface-tension mismatch between the drop and the interface triggering the generation of tangential Marangoni stresses in the plane of the common interface formed post-coalescence. In their experimental and numerical investigation of the coalescence of a water drop with an ethanol reservoir, they reported that Marangoni-induced flow leads to the ejection of an additional drop from its summit during its vertical stretching. Similarly, the generation of gradients of surface tension can also be triggered by the use of surfactants (Manikantan & Squires Reference Manikantan and Squires2020).
Current understanding of the coalescence dynamics in such surfactant-laden systems came from Dong, Weheliye & Angeli (Reference Dong, Weheliye and Angeli2019), who experimentally showed for the first time surfactant concentration profiles for systems characterised by high Bond numbers. Additionally, they suggested that surfactants have a strong effect on the interfacial dynamics inducing interfacial rupture (i.e. hole formation) in an off-axis location. Their interfacial concentration profiles agree qualitatively with the previous numerical work performed by Martin & Blanchette (Reference Martin and Blanchette2015). Finally, Shim & Stone (Reference Shim and Stone2017) suggested that the presence of surfactants decreases the air drainage time between the drop and interface (so-called, ‘damped-coalescence-cascade mechanism’).
The effect of surfactants on the thinning and pinch-off of liquids threads has been studied by Ambravaneswaran, Phillips & Basaran (Reference Ambravaneswaran, Phillips and Basaran2000); Craster, Matar & Papageorgiou (Reference Craster, Matar and Papageorgiou2002); Timmermans & Lister (Reference Timmermans and Lister2002) and Liao, Franses & Basaran (Reference Liao, Franses and Basaran2006) via linear stability analysis, one-dimensional models and full numerical simulations. This work concluded that surfactants are advected from the singularity point due to strong axial flow and do not modify the self-similar structure of the flow as breakup is approached; this structure, and associated scaling exponent, remain unaltered from the surfactant-free case (Eggers Reference Eggers1993; Brenner, Lister & Stone Reference Brenner, Lister and Stone1996). However, McGough & Basaran (Reference McGough and Basaran2006) and Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2018) have demonstrated that the addition of surfactants results in the formation of microthreads during thread thinning driven by Marangoni-induced flow near but not at the pinch point. Additionally, recent studies have examined the interplay between Marangoni stresses, capillarity and surface viscous effects (present at sufficiently high surfactant concentrations) via solution of the one-dimensional slender jet (Wee et al. Reference Wee, Wagoner, Kamat and Basaran2020) and full Stokes equations (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2021) reaching similar conclusions (i.e. surface rheological effects lead to a decrease in the rate of thinning). Recently, Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020) and Kamat et al. (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2020) showed the deleterious effect of Marangoni stresses on interfacial singularities in the context of the inhibition of the end-pinching mechanism during the capillary retraction of a liquid thread. Marangoni-induced flow results in the suppression of the stagnation point by flow reversal in the vicinity of the neck, and the higher generation of vorticity from the neck.
Although Martin & Blanchette (Reference Martin and Blanchette2015) have shown that the presence of surfactants is responsible for the inhibition of the partial coalescence event, there is still a lack of explanation of what causes this pinch-off inhibition. Moreover, the appreciation of the distribution of Marangoni stresses in crucial regions of the interface, and other significant insights into the flow fields close to the pinch-off, are also missing. The present study aims to clarify and answer these questions and to overcome all numerical difficulties presented in Martin & Blanchette (Reference Martin and Blanchette2015) by taking into account the nonlinear relation between the surfactant concentration and surface tension, which is of central importance in non-dilute systems. Additionally, we are able to explore parameter ranges corresponding to large density and viscosity contrasts, corresponding to air–water systems, without numerical difficulties, and to go beyond the surfactant elasticity range studied by Martin & Blanchette (Reference Martin and Blanchette2015).
The rest of this article is organised as follows: § 2 presents the governing equations, numerical set-up and the validation of the surfactant-free case against the experimental observations of Blanchette & Bigioni (Reference Blanchette and Bigioni2006). Section 3 provides a discussion of the results which are focused on the origin of the inhibition of the interfacial singularity, and a parametric study accounting for the strength of the Marangoni stress and sorption dynamics. Finally, concluding remarks are summarised in § 4.
2. Problem formulation and numerical method
With the purpose of studying the dynamics of interfacial coalescence in the presence of surfactants, we perform direct numerical simulations of the two-phase Navier–Stokes equations in a three-dimensional Cartesian domain 
$\boldsymbol {x} = (x, y, z )$ (see figure 1a). The treatment of the interface and its surface-tension forces is handled using a hybrid front-tracking/level-set technique, also known as the level contour reconstruction method (Shin & Juric Reference Shin and Juric2009; Shin, Chergui & Juric Reference Shin, Chergui and Juric2017), with surfactant transport being resolved both in the bulk and on the interface. More information on the numerical technique applied to surfactant transport can be found in the work of Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). Moreover, the dependence of the surface tension on the interfacial surfactant concentration is described by a nonlinear Langmuir equation of state (Muradoglu & Tryggvason Reference Muradoglu and Tryggvason2014; Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018).
Schematic representation of the flow configuration, and validation of the numerical procedure: (a) initial shape of the drop resting close to the interface, highlighting the computational domain of size 
$12R_o\times 12R_o\times 6R_o$ (not to scale) in a three-dimensional Cartesian domain, 
$\boldsymbol {x} = (x, y, z)$, with a resolution of 
$386^{3}$; (b) direct comparisons of our numerical predictions for a surfactant-free case (blue line) with experimental results reported by Blanchette & Bigioni (Reference Blanchette and Bigioni2006) for the post-coalescence dynamics of an ethanol drop in air prior to interfacial singularity formation with 
$Oh=0.011$ and 
$Bo=0.09$; also shown in red lines are the numerical solutions by Deka et al. (Reference Deka, Biswas, Sahu, Kulkarni and Dalal2019) for the same case.
2.1. Scaling
In what follows, all variables will be made dimensionless (represented by tildes) using
\begin{equation} \left.\begin{gathered} \tilde{\boldsymbol{x}}=\frac{\boldsymbol{x}}{R_o},\quad \tilde{t}=\frac{t}{t_c},\quad \tilde{\boldsymbol{u}}=\frac{\boldsymbol{u}} {U},\quad \tilde{p}=\frac{p}{\rho_l U^{2}},\\ \tilde{\sigma}=\frac{\sigma}{\sigma_s},\quad \tilde{\varGamma}=\frac{\varGamma}{\varGamma_\infty},\quad \tilde{C}=\frac{C}{C_\infty},\quad \widetilde{C_s}=\frac{C_s}{C_\infty}, \end{gathered}\right\} \end{equation}
where 
$t$, 
$\boldsymbol {u}$ and 
$p$ stand for time, velocity and pressure, respectively. The physical parameters correspond to the liquid density 
$\rho _l$, viscosity 
$\mu _l$, surface tension 
$\sigma$, surfactant-free surface tension 
$\sigma _s$ and gravitational acceleration, 
$g$; 
$t_c=\sqrt {\rho _l R_o^{3}/\sigma _s}$ is the capillary time scale and 
$R_o$ is the initial drop radius; hence, the velocity scale is 
$U=R_o/t_c= \sqrt {\sigma _s/(\rho _l R_o)}$. The interfacial surfactant concentration, 
$\varGamma$, is scaled on the saturation interfacial concentration, 
$\varGamma _{\infty }$, whereas the bulk and bulk sub-phase (the region immediately adjacent to the interface) surfactant concentrations given by 
$C$ and 
$C_s$, respectively, are scaled on the initial bulk surfactant concentration, 
$C_{\infty }$. As a result of the scaling in (2.1a–h), the dimensionless forms of the governing equations for the flow and the surfactant transport are respectively expressed as
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}=0, \end{gather}
\begin{gather} \tilde{\rho} \left(\frac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{t}}+\tilde{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{\boldsymbol{u}}\right) + \boldsymbol{\nabla} \tilde{p} ={-}Bo \boldsymbol{i}_z + Oh \boldsymbol{\nabla}\boldsymbol{\cdot} [\tilde{\mu} (\boldsymbol{\nabla} \tilde{\boldsymbol{u}} +\boldsymbol{\nabla} \tilde{\boldsymbol{u}}^\textrm{T})] \nonumber\\ \hspace{11.1pc} +\int_{\tilde{A}\widetilde{(t)}} \left(\tilde{\sigma} \tilde{\kappa} \boldsymbol{n} + \boldsymbol{\nabla}_s \tilde{\sigma} \right) \delta (\tilde{\boldsymbol{x}}-\tilde{\boldsymbol{x}}_{f})\mbox{d}\tilde{A}, \end{gather}
$$\begin{gather} \frac{\partial \tilde{C}} {\partial \tilde{t}}+\tilde{\boldsymbol{u}}\cdot \boldsymbol{\nabla} \tilde{C}= \frac{1}{Pe_b} \nabla^{2} \tilde{C}, \end{gather}$$
$$\begin{gather}\frac{\partial \tilde{\varGamma}}{\partial \tilde{t}}+\boldsymbol{\nabla}_s \boldsymbol{\cdot}(\tilde{\varGamma} \tilde{\boldsymbol{u}} _\textrm{t}) =\frac{1}{Pe_s} \nabla^{2}_s \tilde{\varGamma}+ Bi (k \widetilde{C_s} (1-\tilde{\varGamma})- \tilde{\varGamma}), \end{gather}$$
$$\begin{gather}\tilde{\sigma}=1 + \beta_s \ln{(1 -\tilde{\varGamma})}, \end{gather}$$
which correspond to the equations of mass and momentum conservation, the convective–diffusion equations for the surfactant bulk and interfacial concentrations and the nonlinear surfactant equation of state, respectively. Here, the density 
$\tilde {\rho }$ and viscosity 
$\tilde {\mu }$ are expressed by 
$\tilde {\rho }=\rho _g/\rho _l + (1 -\rho _g/\rho _l) \mathcal {H}(\tilde {\boldsymbol {x}},\tilde {t})$ and 
$\tilde {\mu }=\mu _g/\mu _l+ (1 -\mu _g/\mu _l) \mathcal {H}( \tilde {\boldsymbol {x}},\tilde {t})$ wherein 
$\mathcal {H}( \tilde {\boldsymbol {x}},\tilde {t})$ represents a smoothed Heaviside function, which is zero in the gas phase and unity in the liquid phase, where the subscript 
$g$ designates the gas phase; 
$\tilde {\boldsymbol {u}}_\textrm {{t}}= (\tilde {\boldsymbol {u}}_\textrm {{s}} \cdot \boldsymbol {t})\boldsymbol {t}$ represents the velocity vector tangential to the interface in which 
$\tilde {\boldsymbol {u}}_\textrm {{s}}$ corresponds to the interfacial velocity; 
$\kappa$ is twice the mean interface curvature calculated from the Lagrangian interface grid; 
$\boldsymbol {\nabla }_s=({\boldsymbol {I}}-\boldsymbol {n}\boldsymbol {n})\boldsymbol {\cdot }\boldsymbol {\nabla }$ stands for the surface gradient operator wherein 
$\boldsymbol {I}$ is the identity tensor and 
$\boldsymbol {n}$ is the outward-pointing unit normal to the interface; 
$\tilde {\boldsymbol {x}}_f$ is the parameterisation of the interface 
$\tilde {A} (\tilde {t})$; finally, 
$\delta$ represents a Dirac delta function that is non-zero when 
$\tilde {\boldsymbol {x}}=\tilde {\boldsymbol {x}}_f$ only. The numerical method used to solve the above equations is described in detail by Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018).
The dimensionless groups that appear in (2.2)–(2.6) are defined as
$$\begin{gather} Bo=\frac{\rho_l g R_o^{2}}{\sigma_s},\quad Oh=\frac{\mu_l}{\sqrt{\rho_l \sigma_s R}}, \end{gather}$$
$$\begin{gather}Bi=\frac{k_d R_o}{U},\quad k=\frac{k_a C_\infty }{k_d},\quad Pe_s=\frac{ U R_o}{D_s},\quad Pe_b=\frac{U R_o}{D_b},\quad \beta_s= \frac{\Re T \varGamma_\infty }{\sigma_s}, \end{gather}$$
where 
$Bo$ and 
$Oh$ are the Bond number (ratio of gravitational to capillary forces) and Ohnesorge number (ratio of viscous to surface tension forces), respectively. The surfactant elasticity parameter, 
$\beta _s$, measures the sensitivity of the surface tension to the surfactant concentration in which the parameter 
$\Re$ represents the thermodynamic ideal gas constant value 
$8.314$ J K
$^{-1}$ mol
$^{-1}$, and 
$T$ denotes temperature. The parameters 
$Pe_s$ and 
$Pe_b$ are the interfacial and bulk Péclet numbers that represent the ratio of convective to diffusive time scales in the plane of the interface and the bulk, respectively. The Biot number, 
$Bi$, stands for the ratio of characteristic desorptive to convective time scales. Finally, 
$k$ is the ratio of adsorption to desorption time scales where 
$k_a$ and 
$k_d$ refer to the surfactant adsorption and desorption coefficients, respectively.
At equilibrium, there is no surfactant exchange between the interface and the bulk, and the last term on the right-hand side of (2.5) reduces to the Langmuir adsorption isotherm
\begin{equation} \chi=\frac{\varGamma_{eq}}{\varGamma_{\infty}}=\frac{k}{(1+k)}, \end{equation} where 
$\chi$ stands for the fraction of the interface covered by adsorbed surfactant. Furthermore, the Marangoni stress, 
$\tilde {\tau }$, which appears in the third term on the right-hand side of (2.3), is expressed as a function of 
$\tilde {\varGamma }$ as follows:
\begin{equation} \tilde{\tau} \equiv \boldsymbol{\nabla}_s \tilde{\sigma} \cdot \boldsymbol{t} ={-}\frac{\beta_s}{1 -\tilde{\varGamma}}\boldsymbol{\nabla}_s\tilde{\varGamma} \cdot \boldsymbol{t}, \end{equation}
where 
$\boldsymbol {t}$ is the unit tangent to the interface. In all cases considered in the present study, the Marangoni time scale, 
$\mu R_o/ \Delta \sigma =O(10^{-4})$ s, as compared with the capillary and sorptive/desorptive time scales, which are of 
$O(10^{-3})$ and 
$O(10^{-3})-O(10^{-4})$ s, respectively; thus Marangoni stresses will play a crucial role in the coalescence phenomenon. Finally, the tildes are dropped henceforth with the understanding that, hereafter, all variables discussed are dimensionless unless stated otherwise.
2.2. Numerical set-up, validation and parameters
The numerical set-up closely follows the work done by Sun et al. (Reference Sun, Zhang, Che and Wang2018), and Martin & Blanchette (Reference Martin and Blanchette2015). Thus, the size of the dimensionless computational domain is chosen as 
$12R_o\times 12R_o\times 6R_o$, which is found to be sufficiently large to avoid the effect of artificial reflections from the boundaries. We define a radial component as 
$r=\sqrt {(x-x_o)^{2} + (y-y_o)^{2}}$ where 
$x_o$ and 
$y_o$ are the abscissa and ordinate drop position, respectively. Solutions are sought subject to Neumann boundary conditions on all variables at the lateral boundaries, 
$p=0$ at the top boundary 
$z=6R_o$ and no slip at the bottom 
$z=0$. At the interface, we impose 
$\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\nabla }\tilde {C} = -Bi Pe_b(k \tilde {C}_s(1 - \tilde {\varGamma }) - \tilde {\varGamma })$ as a condition on 
$\tilde {C}$ (we refer the reader to Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018) for more information). The initialisation of the interface corresponds to a spherical drop resting immediately above a horizontal flat interface before its interfacial rupture (e.g. all the velocities set to zero) where both drop and liquid pool are made up of the same liquid. Importantly, the drop is connected to the flat interface by a neck of radius 
$0.25R_o$ for the initialisation of the dynamics; a similar approach has been previously used by Blanchette & Bigioni (Reference Blanchette and Bigioni2006, Reference Blanchette and Bigioni2009) and Martin & Blanchette (Reference Martin and Blanchette2015). The assumption is based on the time scale associated with the retraction of the neck 
$t_{CR}=R_o/\sqrt {2\sigma _s/\rho _l \delta }$, which is too short to have an influence on the phenomenon.
Figure 1(b) highlights qualitative and quantitative validation of our numerical framework with results from the literature (Eggers Reference Eggers1993; Blanchette & Bigioni Reference Blanchette and Bigioni2006; Castrejón-Pita et al. Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, John and Basaran2015; Deka et al. Reference Deka, Biswas, Sahu, Kulkarni and Dalal2019). The numerical simulations have been benchmarked against the surfactant-free experimental results of Blanchette & Bigioni (Reference Blanchette and Bigioni2006) in terms of the temporal interfacial dynamics of the coalescence of an ethanol drop surrounded by air (displayed in figure 1b). Our numerical results are provided as snapshots of the interface location at times corresponding to those given by Blanchette & Bigioni (Reference Blanchette and Bigioni2006). We have also included the numerical predictions from Deka et al. (Reference Deka, Biswas, Sahu, Kulkarni and Dalal2019). Figure 1(b) demonstrates that our numerical framework is capable of predicting accurately the interfacial dynamics of the coalescence phenomenon for ‘clean’ interfaces. Since 
$Oh <<1$ (
$Oh=0.011$), during the early thinning stages, there is competition between the fluid inertia and the opposing capillary pressure which corresponds to the inertial regime. The dynamics is expected to transition from this regime (
$r_{min} \sim \tau ^{2/3}$) to the inertial–viscous regime (
$r_{min} \sim \tau$)) when the local Reynolds number drops to 
$Re_{local} \sim 1$, which occurs when 
$r_{min} \sim Oh^{2} \sim 10^{-4}$ (Notz, Chen & Basaran Reference Notz, Chen and Basaran2001). We note that the global three-dimensional nature of our numerical technique, taking into account the entire domain, makes it prohibitive computationally to reach the level of mesh refinement needed to capture all of the regime transitions accompanying the approach to pinch-off, which have been highlighted by the work of Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, John and Basaran2015).
The dimensionless quantities for the studied phenomenon are consistent with experimentally realisable systems. The Ohnesorge number was set to 
$Oh=2\times 10^{-2}$ because it allows for the observation of the interplay between the full range of dynamics as there is a competition between inertial, viscous and capillary forces. The chosen density and viscosity ratios, 
$\rho _g/ \rho _l =1.2 \times 10^{-3}$ and 
$\mu _g/\mu _l= 0.018$, respectively, are representative of an air–water system. The elasticity number 
$\beta _s$ depends on the interfacial concentration at saturation, 
$\varGamma _\infty$, which, in turn, is related to the critical micelle concentration (CMC) that is of 
$O(10^{-6})$ mol m
$^{-2}$. We have explored the range of 
$0.1<\beta _s<0.5$ which corresponds to 
$2.9 \times 10^{-6}<$ CMC 
$<1.4 \times 10 ^{-5}$ mol m
$^{-2}$. Typical values for the interfacial diffusion coefficient for surfactants such as sodium dodecyl sulphate (SDS), N-dodecyl-N,N-dimethylammonio-3-propane sulphonate and similar monomers in aqueous solution, are within the range of 
$10^{-12}< D_s< 10^{-8}$ m
$^{2}$ s
$^{-1}$ when 
$\varGamma$ is below the CMC (Joos, Bleys & Petre Reference Joos, Bleys and Petre1982; Siderius, Kehl & Leaist Reference Siderius, Kehl and Leaist2002); this range also covers phospholipid-based pulmonary surfactants, such as N-(7-nitrobenz-2-oxa-1,3-diazol-4-yl)-1,2-dihexadecanoyl-sn-glycero-3-phosphoch-oline (NBD-PC), which are considered effectively insoluble (Fallest et al. Reference Fallest, Lichtenberger, Fox and Daniels2010; Strickland, Shearer & Daniels Reference Strickland, Shearer and Daniels2015). Therefore, the interfacial Péclet number 
$Pe_s$ lies in the range 
$10^{3}< Pe_s<10^{6}$. Recently, Batchvarov et al. (Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020) suggested that the investigated interfacial dynamics reaches saturation above 
$Pe_s=100$; thus, the selected interfacial Péclet is set to 
$Pe_s=100$. In terms of the chosen bulk Péclet number, Agrawal & Neuman (Reference Agrawal and Neuman1988) suggested that the interfacial and bulk Péclet numbers are of the same order of magnitude; on this basis, hereafter, we set 
$Pe_b=Pe_s$. In summary, we have chosen the values of the surfactant-related parameters to ensure that all of the relevant physical processes associated with surfactant transport such as Marangoni stresses, surface/bulk diffusion and sorption kinetics are represented in the present study.
In terms of mesh resolution studies, we have ensured that our numerical simulations are mesh independent, and subsequently, for a resolution of 
$(386)^{3}$, the results do not change with decreasing cell size. We have also ensured that the liquid volume and surfactant mass conservation are satisfied with errors of under 
$10^{-3}\,\%$ and 
$10^{-2}\,\%$, respectively (see the Appendix for more information). Extensive mesh studies for surface-tension-driven phenomena using the same computational method have been published previously (Batchvarov et al. Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020). A discussion of the results is presented next.
3. Results
Following the good agreement between the surfactant-free coalescence simulation and the experimental results of Blanchette & Bigioni (Reference Blanchette and Bigioni2006), key surfactant effects will be investigated in this section. We first display our results related to the effect of insoluble surfactant, e.g. NBD-PC (Fallest et al. Reference Fallest, Lichtenberger, Fox and Daniels2010; Strickland et al. Reference Strickland, Shearer and Daniels2015), showing ultimately our insights regarding the surfactant-driven escape from a potential pinch-off singularity, such as the one depicted in figure 1(b). We then present the effect of the sorptive kinetics on the phenomenon through the use of soluble surfactants, e.g. SDS (Siderius et al. Reference Siderius, Kehl and Leaist2002). It is also worth mentioning that all surfactant simulations have been carried out until the neck has either pinched off or undergone reopening. Additionally, we provide with conclusive evidence that the neck reopening is driven by Marangoni stresses.
3.1. Insoluble surfactants
We start the discussion of the results by presenting the effect of the surface-active agents through the analysis of the elasticity parameters 
$\beta _s$ with 
$Oh=0.02$, 
$Bo=10^{-3}$, 
$Pe_s=100$ and 
$\varGamma _o=\varGamma _\infty /2$. Figure 2 shows the spatio-temporal interfacial dynamics for the surfactant-free and surfactant-laden cases as a function of the elasticity parameter. At the early stages of the dynamics, the neck expands as a result of the capillary retraction of the liquid bridge which separates the drop from the underlying liquid pool. The capillary retraction gives rise to the formation of capillary waves that travel upwards towards the drop summit. As pointed out by Blanchette & Bigioni (Reference Blanchette and Bigioni2006), the oscillations caused by the travelling capillary waves yield vertical stretching forming a nearly cylindrical drop, as shown in figure 2(c), before capillarity acts to drive the dynamics towards a more energy-favourable state by pulling on the sides of the drop. This capillary action leads to pinch-off of the liquid bridge via a singularity which culminates in the formation of a secondary droplet; this, in turn, follows a ‘cascade of coalescence events’ until the coalescence process is completed, as also shown by Thoroddsen & Takehara (Reference Thoroddsen and Takehara2000), Blanchette & Bigioni (Reference Blanchette and Bigioni2006), Blanchette & Bigioni (Reference Blanchette and Bigioni2009), Aryafar & Kavehpour (Reference Aryafar and Kavehpour2006) and Houssainy, Kabachek & Kavehpour (Reference Houssainy, Kabachek and Kavehpour2020). Similar phenomena in the surfactant-laden case are not the focus of the present work wherein we concentrate on elucidating the mechanisms by which the presence of surfactant leads to escape from singularity formation.
Effect of 
$\beta _s$ on the drop-interface coalescence dynamics for insoluble surfactants. Spatio-temporal evolution of the three-dimensional interface shape for surfactant-free, (a–d), and surfactant-laden coalescence for 
$\beta _s=0.1$, (e–h), 
$\beta _s=0.3$, (i–l) and 
$\beta _s=0.5$, (m–p). Here, the dimensionless parameters are 
$Oh=0.02$ and 
$Bo=10^{-3}$, and for the surfactant-laden cases, 
$Pe_s=100$ and 
$\varGamma _o=\varGamma _\infty /2$. The colour indicates the magnitude of 
$\varGamma$, and legend is shown in (e).
For all surfactant-laden cases, the generation of a secondary droplet is avoided even for the lower end of the elasticity parameter range. For 
$\beta _s=0.1$, the dynamics follows closely that of the surfactant-free case, where significant vertical stretching of the original droplet is observed (see figure 2e–h). At the point where surface tension is expected to dominate the narrowing of the neck, the presence of non-uniform surfactant concentration generates Marangoni stresses that change the outcome of the dynamics. By increasing 
$\beta _s$, surfactant redistribution along the interface is enhanced as displayed in figures 2(i–l) and 2(m–p), for 
$\beta _s=0.3$ and 
$\beta _s=0.5$, respectively. The surfactant concentration gradient, and associated transport, is seen to suppress the capillary waves and to limit the vertical stretch of the drop. For the region of 
$\beta _s>0.1$, Marangoni-induced flow inhibits the capillary singularity. However, a more exhaustive parametric study of the effect of the elasticity parameter is clearly worthwhile in order to identify the locus of 
$\beta _s$ for which Marangoni stresses do not result in neck reopening.
Figure 3 shows the immobilising effect brought about by the presence of surfactants through the analysis of the temporal dynamics of the maximum vertical stretch of the droplet, 
$z_{max}$, the neck radius, 
$r_{min}$, and the kinetic energy, 
$E_k=\int _V (\rho \boldsymbol {u}^{2}/2)\,{\textrm {d}}V$. Here, the 
$E_k$ values have been normalised by the surface energy 
$E_s=S \sigma _s$, where 
$S$ is the initial superficial area of the system. The evidence for damping of the upward drop oscillation can be seen in figure 3(a). Here, the increase in 
$\beta _s$ is seen to depress the maximum crest location of the drop. These observations confirm the expectations of Martin & Blanchette (Reference Martin and Blanchette2015) of suppression of the axial oscillation with an increase in the surfactant strength (though these authors were only able to run simulations for 
$\beta _s \leq 0.2$). Interestingly, the temporal evolution of 
$z_{max}$ exhibits a non-monotonic dependence on 
$\beta _s$, with the most suppressed crest being associated with the intermediate value of 
$\beta _s=0.3$. The physical explanation of this outcome will be provided in the discussion of figure 4 below. Furthermore, investigation of the temporal variation of the minimum neck radius, 
$r_{min}$ (see figure 3b), confirms neck reopening for all surfactant-laden cases, with a 50 % rise in the time associated with the onset of re-opening, 
$t_r$, corresponding to an increase in 
$\beta _s$ from 0.1 to 0.5. The non-monotonic dependence on 
$\beta _s$ is also exhibited by 
$t_r$: even though the longest delay in neck closure is observed for the highest 
$\beta _s$ studied, the 
$t_r$ value for 
$\beta _s=0.3$ is associated with the largest 
$r_{min}$. Further quantification of these physical phenomena is provided below. Additionally, figure 3(b) shows that the presence of surfactant rigidifies the interfacial dynamics by slowing down the neck growth at the early stages of the capillary-driven expansion.
Effect of 
$\beta _s$ on the temporal dynamics of the vertical extent of the drop (a), its minimum neck radius (b) and the system kinetic energy 
$E_k$ (c), for 
$Oh=0.02$, 
$Bo=10^{-3}$, 
$Pe_s=100$ and 
$\varGamma _o=\varGamma _\infty /2$.
(a–i) Effect of the elasticity parameter 
$\beta _s$ on the flow and surfactant concentration fields associated with the drop-interface coalescence phenomenon. Two-dimensional representation of the interface location, 
$\varGamma$, 
$\tau$ and the radial component of the interfacial velocity 
$u_{tr}$ are shown in (a–d) and (e–h) for 
$t=1.20$ and 
$t=1.68$, respectively. Note that the abscissa in (a,e) corresponds to the radial coordinate 
$r$, and in (b–d) and (f–h) to the arclength 
$s$. (i) Represents a magnified view of (h). The arrows in (g) indicate the directions of motion driven by the Marangoni stresses 
$\tau$; in (h,i), points P1, and S1 and S2 designate the peak in 
$u_{tr}$ and the stagnation points in the surfactant-free 
$u_{tr}$ profile, respectively. The diamond shapes in (i) show the location of the necks. The parameter values and the times for the radius are the same as in figure 3. (j,k) Effect of surfactants on the azimuthal vorticity 
$\omega _{\theta }$ for the surfactant-free (left panels), and the surfactant-laden cases (right panels), for 
$\beta _s=0.5$, at 
$t=1.20$, 
$t=1.68$, respectively. All other parameters remain unchanged from figure 3. The colour indicates the value of the azimuthal vorticity 
$\omega _{\theta }$. The arclength 
$s$ starts from the apex of the droplet.
Finally, inspection of the kinetic energy plots shows that the presence of surfactants induces a monotonically decreasing overall value of 
$E_k$ with 
$\beta _s$ over the range in time encompassing the creation of the cylindrically shaped drop (see figure 3c). This reduction of the kinetic energy is a result of the rigidification of the interface brought on by the tangential Marangoni stresses in agreement with Asaki, Thiessen & Marston (Reference Asaki, Thiessen and Marston1995). It is evident, however, that, for the surfactant-free case, 
$E_K$ decreases rapidly, as the drop breaks up via neck pinch-off, eventually dipping below those associated with 
$\beta _s=0.1$ and 
$\beta _s=0.3$.
The next part of the analysis focuses on the time evolution of a two-dimensional projection of the interfacial shape, 
$\tilde \varGamma$, 
$\tau$ and the radial component of the interfacial velocity, 
$u_{tr}$, presented in (a–i) of figure 4. We also show the interplay between the surface and the azimuthal component of the vorticity field (vorticity is defined as 
$\omega =\boldsymbol {\nabla } \times \boldsymbol {u}$), displayed in figures 4(j) and 4(k). In the surfactant-free case, it is seen from figure 4(d) that 
$u_{tr}<0$ and 
$u_{tr}>0$ upstream and downstream of the developing neck, which drives flow away from this region. The narrowing of the neck induces capillary-driven flow that leads to further neck thinning and the development of large peaks in 
$u_{tr}$, as shown in figure 4(e,h), which are typical of singularity formation. Close inspection of the 
$u_{tr}$ profile in figure 4(i) for the surfactant-free case reveals that it is characterised by the presence of a large velocity peak (P1) and two stagnation points (labelled S1 and S2) with the neck sandwiched between them. Over time, the inertio-capillary-induced flow ultimately culminates in interfacial breakup to form a daughter droplet. From the vorticity plots in figures 4(j) and 4(k), it is seen that, for the surfactant-free case, the vorticity generation is confined to the vicinity of the neck as the two stagnation points aid the fluid recirculation around the neck (so-called ‘vortex ring’, displayed in the left-panel of figure 4j). As time evolves, the interfacial curvature of the neck increases, and a large vorticity generation can be observed on the side of the bulk accompanying the eventual neck pinch-off, as depicted in the left-panel of figure 4(k). More information regarding the mechanisms which induce the generation of vorticity at the liquid–gas interface is provided below.
For the surfactant-laden cases, the accumulation of 
$\varGamma$ near the nascent neck can be seen in figure 4(b) thus giving rise to a local decrease of 
$\sigma$. The presence of 
$\varGamma$ gradients results in the generation of a large positive peak in the 
$\tau$ profile in the vicinity of the neck region, which is largest for the intermediate value of 
$\beta _s=0.3$, as shown in figure 4(c) (consistent with the non-monotonic response of the dynamics observed in figure 3). Upstream and downstream of the neck, 
$\tau > 0$ and 
$\tau <0$, respectively, which drives flow towards the drop summit and tail, reflected by 
$u_{tr}<0$ and 
$u_{tr}>0$, respectively. Although the overall shape of the 
$u_{tr}$ curve for the surfactant-free case is robust to the addition of insoluble surfactant, it is evident that the magnitude of 
$u_{tr}$ decreases with increasing 
$\beta _s$, particularly in the neck region; moreover, the oscillation in 
$u_{tr}$ in the surfactant-free case is damped out for 
$\beta _s>0$.
Additionally, by close inspection of the 
$u_{tr}$ plots for the surfactant-laden cases in figure 4(i), it becomes clear that only one stagnation point is present near the neck for 
$\beta _s=0.3$ and 
$0.5$; thus, 
$u_{tr}>0$ towards its tail. The Marangoni-induced flow has therefore led to the suppression of one of the stagnation points. Furthermore, by comparing the vorticity field pattern of the surfactant-free and the surfactant-laden cases, a change is observed as a result of the presence of surfactants, and the inhibition of a stagnation point. The generation of vorticity is also confined to the vicinity of the free surface; however, the ‘vortex ring’ no longer exists, as shown in the right panel of figure 4(j).
Further in time, when the escape of capillary singularity commences, we observe that vorticity is separated from the vicinity of the interface, and advected towards the bulk of liquid reservoir, which is consistent with the findings of Ananthakrishnan & Yeung (Reference Ananthakrishnan and Yeung1994), supporting the reopening of the neck, as displayed in the right panel of figure 4(k). The reason behind the vorticity separation is the inhibition of one of the stagnation points. This behaviour is similar to the phenomenon explained by Hoepffner & Paré (Reference Hoepffner and Paré2013) in terms of capillary retraction of surfactant-free viscous ligaments where they suggested that the advection of 
$\omega _{\theta }$ plays a crucial role in their escape from breakup. These observations also agree with the recent studies reported by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021) in terms of the escape of capillary singularity during the capillary retraction of a liquid thread, and the inhibition of jet-drop formation from bursting bubbles, respectively.
Next, we turn our attention to the role of the Marangoni stresses in the generation of vorticity at the gas–liquid interface; this is consistent with the work of Batchelor (Reference Batchelor1967), who concluded that vorticity in a homogeneous fluid is generated at the boundaries only. Several papers have discussed the generation of vorticity at a free surface (Longuet-Higgins Reference Longuet-Higgins1992; Cresswell & Morton Reference Cresswell and Morton1995; Peck & Sigurdson Reference Peck and Sigurdson1998; Lundgren & Koumoutsakos Reference Lundgren and Koumoutsakos1999; Brøns et al. Reference Brøns, Thompson, Leweke and Hourigan2014; Thoraval, Li & Thoroddsen Reference Thoraval, Li and Thoroddsen2016). Cresswell & Morton (Reference Cresswell and Morton1995) was the first to explain the formation of the vortex ring during the impact of a water drop on water surface. They have stated that the origin of the vorticity occurs as a boundary condition on the interface in order to simultaneously satisfy the stress-free boundary condition and irrotationality type of the flow. Once the vorticity is produced, it is diffused into a thin boundary layer and advected towards the pool. The formation of the vortex ring was presented experimentally by Dooley et al. (Reference Dooley, Warncke, Gharib and Tryggvason1997), who introduced surfactant into their system in order to lower the interfacial tension, but the vorticity in the vortex ring emerges much earlier than their sketches suggest.
Assuming that the interface behaves as a viscous free surface because of the small air-to-water viscosity and density ratios (similar assumptions have been made previously by Dooley et al. (Reference Dooley, Warncke, Gharib and Tryggvason1997) and Xia et al. (Reference Xia, He, Yu, Zhao and Zhang2017)). Kamat et al. (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2020) demonstrated that, for surface-tension-driven phenomena, vorticity generation depends solely on the interfacial boundary conditions when 
$Oh<<1$. As a result, the tangential stress at the interface is balanced by the surface-tension gradients, resulting in
\begin{equation} \boldsymbol{t} \cdot \boldsymbol{D} \cdot \boldsymbol{n} =\boldsymbol{t} \boldsymbol{\cdot}\boldsymbol{\nabla}_s \sigma, \end{equation}
in which, 
$\boldsymbol {D}$ represents the rate of deformation tensor (the symmetric part of the velocity gradient tensor). By further mathematical manipulation (Lundgren & Koumoutsakos Reference Lundgren and Koumoutsakos1999), the generation of vorticity at the free surface depends entirely on the velocity field, interfacial geometry and surface-tension gradients
\begin{equation} \omega= \omega_n+ \omega_t+ \omega_{\tau}={-}2 \frac{\partial \boldsymbol{u}\cdot \boldsymbol{n}}{\partial s}+ 2 \boldsymbol{u} \cdot \boldsymbol{t} \kappa + \boldsymbol{t} \boldsymbol{\cdot}\boldsymbol{\nabla}_s \sigma. \end{equation}
Similar results for 
$\omega _n$ and 
$\omega _t$ have been reported by Lundgren & Koumoutsakos (Reference Lundgren and Koumoutsakos1999) and Brøns et al. (Reference Brøns, Thompson, Leweke and Hourigan2014). The first two terms on the right-hand side of (3.2) correspond to the normal and tangential velocity-driven vorticity generation, respectively, whereas the last term is representative of the Marangoni stress vorticity contribution.
Figure 5 shows the vorticity distribution along the interface according to (3.2). By close inspection of the profiles, we observe that the vorticity generation at the interface is highly dominated by the interfacial curvature term, 
$\omega _{t}$. Moreover, in the vicinity of the neck we discover the existence of a positive peak in the Marangoni stress-driven vorticity production (i.e. 
$\omega _{\tau }$). The peak of 
$\omega _{\tau }$ has a different sign in comparison with 
$\omega _n$. The leading cause for this behaviour stems from the suppression of the stagnation points on both sides of the neck, as shown in figure 4(i). Ultimately, this analysis demonstrates the positive effect of surface-tension-driven vorticity generation on the neck reopening process. Additionally, this finding is in agreement with Kamat et al. (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2020), who concluded that the generation of vorticity arising from the presence of surfactants is generated over a time scale of similar magnitude to the capillary time scale.
Vorticity production, 
$\omega$, along the gas–liquid interface expressed in terms of the local normal and tangential components from the velocity field and the Marangoni stresses, represented by 
$\omega _n$, 
$\omega _t$ and 
$\omega _{\tau }$, respectively. The surfactant-laden case is characterised by 
$\beta _s=0.5$ at 
$t = 1.68$. All other parameters remain unchanged from figure 3. The diamond shape shows the location of the neck.
Finally, we aim to provide more conclusive evidence that the interfacial singularity inhibition is Marangoni driven rather than a result of the reduction of the surface tension (i.e. capillary pressure reduction). For this reason, we have performed an additional simulation in which Marangoni stresses have been suppressed (similar to what was done by Xu (Reference Xu2007) and Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2018)). Figure 6 reports the temporal evolution of the maximum axial position 
$z_{max}$, the neck radius 
$r_{min}$ and the kinetic energy 
$E_k$ for the surfactant-free and Marangoni-suppressed cases. Similar flow behaviours between the surfactant-free and Marangoni-suppressed cases are observed. The most remarkable finding is that, for the Marangoni-suppressed case, it is observed that the mean reduction of the surface tension does not prevent the horizontal collapse of the droplet (see figure 6b). The inspection of the kinetic energy plot shows that the Marangoni-suppressed and surfactant-free cases have almost identical behaviours (see figure 6c). Therefore, when Marangoni stresses are enabled fully, a change of the fate of the coalescence is observed via the reopening of the neck.
Demonstration that tangential Marangoni stresses are responsible for the inhibition of the interfacial singularity. Temporal evolution of the maximum vertical displacement of the interface, neck radius and kinetic energy; (a–c), respectively, for the surfactant-free, full-Marangoni 
$|\tau |>0$ and no-Marangoni cases 
$|\tau |=0$, for 
$Oh=0.02$, 
$Bo=10^{-3}$, 
$\beta _s=0.5$, 
$Pe_s=100$ and 
$\varGamma _o=0.5\varGamma _\infty$.
To conclude this section, we will turn our attention towards the role of Marangoni stress in the large-scale dynamics of the coalescence phenomenon (i.e. neck curvature). According to Alhareth & Thoroddsen (Reference Alhareth and Thoroddsen2020), when the axial curvature 
$\kappa _x$ overcomes the azimuthal curvature 
$\kappa _{\theta }$ this leads to a negative Laplace pressure preventing the capillary singularity (i.e. 
${\rm \Delta} p =\sigma ( \kappa _{\theta }-\kappa _x)$). This interplay between the curvatures can be seen as a result of the Marangoni-induced flow. A representation of the neck shapes with their respective curvatures for a surfactant-laden case can be seen in figure 7.
Spatio-temporal evolution of the drop neck for the surfactant-laden case characterised by 
$\beta _s=0.5$. We show a magnified view of the neck region, with 
$\kappa _{\theta }$ and 
$\kappa _x$ indicating the azimuthal and axial curvatures, respectively.
3.2. Soluble surfactants
In this subsection we present a discussion of the results associated with the effects of surfactant solubility and sorption kinetics, parameterised by 
$Bi$ and 
$k$, respectively. Unless stated otherwise, the parameters remain fixed to their ‘base’ values: 
$Oh=0.02$, 
$Bo=10^{-3}$, 
$\beta _s=0.5$ and 
$Pe_s=100$; the interfacial surfactant concentration is initialised using its equilibrium surfactant concentration, thus 
$\varGamma _o=\chi =k/(1+k)$. Once again, simulations are carried out until either neck pinch-off or reopening has been observed.
3.2.1. Effect of the Biot number, 
$Bi$
Figure 8 shows the effect of varying 
$Bi$ in the range 0.1–10 on the drop maximal vertical extent 
$z_{max}$, the neck radius 
$r_{min}$ and the kinetic energy 
$E_k$ with 
$k=1$; also shown are the curves associated with the insoluble surfactant and surfactant-free cases which respectively correspond to the 
$Bi \rightarrow 0$ and 
$Bi \rightarrow \infty$ (and/or 
$\beta _s \rightarrow 0$) limits. At the lower end of this range (e.g. 
$Bi =0.1$), the sorptive time scales are much larger than those associated with interfacial effects; therefore, the dynamics is dominated by capillarity and Marangoni stresses, and is therefore expected to be similar to that observed for the insoluble surfactant case. This is confirmed upon inspection of figure 8 as the curves associated with the 
$Bi=0.1$ case practically overlap with those generated for the insoluble surfactant case. For large 
$Bi$, the monomers desorb rapidly from the interface, which represents the case of a highly soluble surfactant characterised by a dynamics that is similar to that accompanying the surfactant-free case. This is also confirmed by comparing the curves associated with the surfactant-free and 
$Bi=10$ cases, the latter corresponding to the largest 
$Bi$ value studied.
Effect of 
$Bi$ on the temporal dynamics of the vertical strength of the drop (a), minimum neck radius (b) and kinetic energy (c), when 
$Oh=0.02$, 
$Bo=10^{-3}$, 
$Pe_s=100$, 
$\beta _s=0.5$, 
$k=1$ and 
$\varGamma _o=\chi$.
As depicted in figure 8(a), increasing the level of solubility leads to a decrease in the surfactant mass at the interface available to induce Marangoni stresses and, consequently, delays the retardation in the initialisation of the vertical stretch of the drop. Interestingly, the lowest 
$z_{max}$ is associated with the intermediate 
$Bi=1$ case, thus 
$z_{max}$ exhibits a non-monotonic dependence on the surfactant solubility. Turning our attention towards the effect of solubility on the neck size, 
$r_{min}$ displayed in figure 8(b), it is clearly seen that Marangoni-induced flow results in the inhibition of the capillary-driven singularity over the entire range of 
$Bi$ values studied. This effect becomes increasingly pronounced with decreasing 
$Bi$ and the 
$r_{min}$ vs 
$t$ profiles for the surfactant-laden cases are bounded between the insoluble and the surfactant-free cases. Finally, the 
$E_k$ profiles depicted in 8(c) show that decreasing 
$Bi$ leads to an overall reduction in the kinetic energy, which stems from the rigidifying effect of the Marangoni stresses; this is weakened by increasing the solubility and the enhanced surfactant desorption from the interface which leads to a reduction in 
$\tau$.
Evidence of surfactant-induced immobilisation with decreasing 
$Bi$ is further provided in figures 9(a,e) and 9(d,h), which depict the interface shape and 
$u_{tr}$, respectively; the rest of the parameters remain unaltered from figure 8. It is also clearly seen from figure 9(b,f) that higher surfactant desorption is observed as 
$Bi$ increases, driven by the mass transfer between the interface and bulk. Although the largest 
$\varGamma$ is associated with the smallest 
$Bi$ values, the largest gradients, and, therefore Marangoni stresses, 
$\tau$, are found for the intermediate Biot number, 
$Bi=1$, as shown in figure 9(c,g); this is consistent with the non-monotonic dependence of the vertical stretch of the drop on 
$Bi$, described in figure 8(a). Notably, a comparison of (b,f) and (c,g) of figure 9 reveals that, over time, the 
$\varGamma$ gradients become sharper, particularly for small and intermediate 
$Bi$, leading to an increase in 
$\tau$, as was also observed in the insoluble surfactant case. These Marangoni stresses counteract the direction of the inertio-capillary-induced flow, dampen the oscillations in 
$u_{tr}$ (see figure 9d,h) and act to prevent neck pinch-off where the efficacy is once more dependent on the magnitude of 
$Bi$. Finally, the same analysis regarding the role of surfactants in the inhibition of stagnation points, explained in § 3.1, can be extrapolated to the solubility parameter, 
$Bi$.
Effect of the solubility parameter 
$Bi$ on the flow and surfactant concentration fields associated with the drop-interface coalescence phenomenon. Two-dimensional representation of the interface location, 
$\varGamma$, 
$\tau$ and the radial component of the interfacial velocity 
$u_{tr}$ are shown in (a–d) and (e-h) for 
$t=1.20$ and 
$t=1.68$, respectively. Note that the abscissa in (a,e) corresponds to the radial coordinate 
$r$, and in (b–d) and (f–h) to the arclength 
$s$. Here, all other parameters remain unchanged from figure 8.
3.2.2. Effect of the adsorption parameter, $k$
Figure 10 shows the effect of the adsorption parameter, 
$k$, on the temporal dynamics of the drop-interface coalescence phenomenon through the analysis of the vertical stretch of the droplet 
$z_{max}$, the neck radius 
$r_{min}$ and the kinetic energy, 
$E_k$, for 
$Bi=1$, and 
$k=(0.01, 1, 5)$. We note that, as 
$k \rightarrow 0$, 
$\chi \rightarrow 0$, and this corresponds to vanishingly small equilibrium interfacial concentrations, which were used to initialise the simulations. In addition, from (2.5), 
$k \rightarrow 0$ implies that 
$\varGamma$ will remain small, and thus, in this limit, we expect the dynamics to be consistent with that associated with the surfactant-free case. For 
$k \gg 1$, on the other hand, the flow behaviour is similar to that observed in the insoluble surfactant case.
Effect of the adsorption parameter 
$k$ on the temporal dynamics of the vertical strength of the drop (a), minimum neck radius (b) and kinetic energy (c), when 
$Oh=0.02$, 
$Bo=10^{-3}$, 
$\beta _s=0.5$, 
$Pe_s=Pe_b=100$, 
$Bi=1$ and 
$\varGamma _o=\chi$.
We start the discussion of the effect of the 
$k$ parameter by analysing its effect on the vertical stretch of the droplet (shown in figure 10a). By inspection of the profiles, a monotonic response of the vertical stretch is observed with decreasing 
$k$ values (e.g. 
$k=0.1$) the dynamics is similar to that of the surfactant-free case characterised by neck formation and pinch-off, as shown in figure 11(a,e); this arises due to the increase in mass transfer from the interface to the bulk decreasing the interfacial concentration (see figure 11b,f) and reducing the magnitude of the Marangoni stress, which is maximal for 
$k=1$, as shown in figure 11(c,g). The trends highlighted in figure 10(a) are mirrored in figure 10(b), that displays the temporal dynamics of 
$r_{min}$, which also exhibits a monotonic dependence on 
$k$. Increasing 
$k$ alters the 
$u_{tr}$ profile in figure 11(d,h) in a similar manner to that observed upon increasing 
$\beta _s$ and/or decreasing 
$Bi$, as was shown previously in figures 4(d,h) and 9(d,h), respectively. As a result, it is seen clearly that the presence of surfactants alters the fate of the coalescence phenomenon as Marangoni-driven flow induces neck reopening. Finally, the 
$E_k$ plots shown in figure 10(c) support, once more, the high interfacial rigidification brought about by the presence of surfactants.
Effect of the adsorption parameter 
$k$ on the flow and surfactant concentration fields associated with the drop-interface coalescence phenomenon. Two-dimensional representation of the interface location, 
$\varGamma$, 
$\tau$ and the radial component of the interfacial velocity 
$u_{tr}$ are shown in (a–d) and (e–h) for 
$t=1.20$ and 
$t=1.68$, respectively. Note that the abscissa in (a,e) corresponds to the radial coordinate 
$r$, and in (b–d) and (f–h) to the arclength 
$s$. Here, all other parameters remain unchanged from figure 10.
4. Conclusions
A study of the effect of Marangoni-induced flow as a result of the presence of surfactants on the drop-interface coalescence was presented using a hybrid front-tracking/level-set method (Shin & Juric Reference Shin and Juric2009; Shin et al. Reference Shin, Chergui and Juric2017, Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). The surfactant transport equations were fully coupled to the Navier–Stokes equations in which the surface tension depends on the interfacial surfactant concentration through a nonlinear Langmuir equation of state. The numerical framework has been validated against the experimental work presented by Blanchette & Bigioni (Reference Blanchette and Bigioni2006) for the surfactant-free coalescence dynamics, and the inertio-viscous scaling laws regarding the temporal evolution of the neck towards its capillary singularity presented by Day, Hinch & Lister (Reference Day, Hinch and Lister1998) and Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, John and Basaran2015). We have selected a surfactant-free base case characterised by the dimensionless quantities 
$Oh=0.02$ and 
$Bo=10^{-3}$. The former parameter ensures a rich dynamics in the inertia–viscous–capillary flow regime, whereas the latter parameter ensures that gravity forces do not affect the dynamics of the system, which could mask effects related to the presence of surfactants.
In drop-interface coalescence in the presence of surfactants, we have demonstrated that Marangoni stresses drive motion from regions of high surfactant concentration (low surface tension) to low concentration (high tension) regions, resulting in retardation of the interfacial dynamics. This immobilising effect of the surfactants as a result of the Marangoni stresses is also observed via the strong reduction in the maximum stretching (by dampening the strength of the capillary waves which converge in the drop summit) of the drop and kinetic energy. We have also shown that the condition for the capillary singularity is the existence of two stagnation points close to the drop neck, which lead to the generation of vorticity in this area. In the presence of surfactants, Marangoni-induced flow suppresses one of the stagnation points, resulting in the advection of vorticity towards the liquid bulk and the reopening of the neck. This effect is strongest for insoluble surfactants and, for soluble surfactants, is maximal for an intermediate range of solubility and sorption kinetic parameter values.
Future directions are related to the performance of numerical simulations featuring three-dimensional behaviours occurring for large Bond numbers. Recently, the experimental work performed by Dong et al. (Reference Dong, Weheliye and Angeli2019) suggested that the presence of surfactants induces the rupture of the interface (i.e. hole formation) in an off-axis location at high Bond numbers. Thus, a fully three-dimensional retracting capillary wave will certainly affect the behaviour of the system rising a more complex coalescence dynamics, and constitute a fruitful area of future research. Additionally, the influence of an uneven initial surfactant distribution could result in a pull of a thin liquid film over the neck driven by Marangoni stresses (Thoroddsen et al. Reference Thoroddsen, Qian, Etoh and Takehara2007), and clearly deserves further study.
Funding
This work is supported by the Engineering and Physical Sciences Research Council, United Kingdom, through a studentship for CRCA in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (EP/L015579/1), and through the EPSRC MEMPHIS (EP/K003976/1) and PREMIERE (EP/T000414/1) Programme Grants. C.R.C.A. also acknowledges the funding and technical support from BP through the BP International Centre for Advanced Materials (BP-ICAM), which made this research possible. O.K.M. also acknowledges funding from PETRONAS and the Royal Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics. We also acknowledge HPC facilities provided by the Research Computing Service (RCS) of Imperial College London for the computing time. D.J. and J.C. acknowledge support through computing time at the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul Intensif) Grant 2020 A0082B06721. The numerical simulations were performed with code BLUE (Shin et al. Reference Shin, Chergui and Juric2017) and the visualisations have been generated using ParaView.
Declaration of interests
The authors report no conflict of interest.
Appendix. Mesh study
On account of showing mesh independent results, we have monitored the temporal variation of the liquid volume of the system for a resolution of 
$(386)^{3}$, which has been used throughout the entire study. Figure 12 shows the plot profiles for the surfactant-free and the surfactant-laden cases. It is evident that the numerical method is capable of capturing the rich interfacial dynamics with a conservation of volume of under 
$10 ^{-3}\%$. With respect to the accuracy of the surfactant equations, we refer to Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018), who carefully benchmarked the formulation and numerical implementation of the surface gradients of surfactant concentration and surface tension. For the studied phenomenon, we observed the conservation of surfactant mass under 
$10 ^{-2}\%$ for all the surfactant-laden cases. Additionally, extensive mesh studies for capillary phenomena, using the same numerical method, had been previously reported (Batchvarov et al. Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020, Reference Batchvarov, Kahouadji, Constante-Amores, Norões Gonçalves, Shin, Chergui, Juric, Craster and Matar2021; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020, Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021).
Relative variation of the liquid volume for the surfactant-free and full-Marangoni cases 
$|\tau |>0$, for 
$Oh=0.02$, 
$Bo=10^{-3}$, 
$\beta _s=0.5$, 
$Pe_s=100$ and 
$\varGamma _o=0.5\varGamma _\infty$.
