Drop encapsulation and bubble bursting in surfactant-laden flows in capillary channels

We present a parametric study of the unsteady phenomena associated with the flow of elongated gas bubbles travelling through liquid-filled square capillaries under high Weber number conditions. These conditions consistently induce the formation of a re-entrant jet at the back of the bubble that commonly gives way to a deep liquid cavity. Subsequent steps include pinch-off events in the cavity to generate one or multiple encapsulated drops which may coalesce, in conjunction with the bursting of the bubble-liquid interface by either the cavity or the drops. Some of these interfacial instabilities have previously been reported experimentally (Olbricht 1996) and numerically (Izbassarov&Muradoglu 2016) for liquid-liquid flow in microchannels. We carry out three-dimensional direct numerical simulations based on a hybrid interface-tracking/level-set method capable of accounting for the presence and dynamic exchange of surfactants between the liquid bulk phase and the liquid-gas interface. Our results indicate that the delicate interplay amongst inertia, capillarity, viscosity, surfactant adsorption/desorption kinetics, and Marangoni stresses has a dramatic influence over the non-axisymmetric morphological structures of the encapsulated drops-elongated bubble. This strong coupling also influences the pinch-off time, penetration depth of the cavity, and number, size, and velocity of the encapsulated drops across the bubble. The observed phenomena are summarised in three main morphological regimes based on surfactant-related parameters and dimensionless groups. A discussion of the flow regime maps is also provided.


Introduction
The problem of gas bubbles propagating in confined capillary channels has been the subject of extensive theoretical, experimental, and numerical investigation (see Etminan et al. (2021) and references therein).In addition to being a classic problem in fluid mechanics, as seen in the now-seminal works of Taylor (1960); Bretherton (1961); Schwartz et al. (1986), these flows are central to multiple applications in technology development and are found in several real-life systems.To name a few examples, gasliquid interfaces in the micro-and capillary scale are involved in two-phase coolers (Karayiannis & Mahmoud 2017), CO 2 sequestration and storage microdevices (Shim et al. 2014), obstructed pulmonary airways (Heil et al. 2008), and volcanic conduits (Carrigan & Eichelberger 1990), and have been proposed as a potential method for bacterial detachment and cleaning (Khodaparast et al. 2017).Since the early works mentioned above, extensive efforts have been undertaken to understand the phenomena dictating the bubble's behaviour along with the surrounding thin liquid film, largely through an idealisation of the operating conditions and the morphology of the bubble itself as a flat cylinder bounded by two semi-spherical menisci.It is well-known by experiments (Taylor 1960) and theory (Bretherton 1961) that the film thickness of the flat region scales as h ∼ Ca  Martinez & Udell 1990).As inertial effects become non-negligible, Han & Shikazono (2009) proposed a correction to the scaling exponents of Ca b and the introduction of the Weber number, W e b , in the scaling law expression.Particular attention has also been placed on exploring the influence of Ca b on the pressure drop ahead of the bubble (Wong et al. 1995), mapping the full bubble topological profiles (de Ryck 2002), analysing the undulatory structures that emerge in the vicinity of the bubble rear (Khodaparast et al. 2015;Magnini et al. 2017) as inertial effects begin to dominate, and predicting the final shape of the front meniscus (Giavedoni & Saita 1999).
Other researchers have delved into the problem by deviating from the traditional idealised assumptions of circular and straight channels, absence of surface-active materials, fully dominating capillary forces, and rigid spherical interfaces.These studies have elucidated critical differences between the steady-state and axisymmetrical features of idealised bubbles and those that adhere to more realistic scenarios.For square or rectangular channel cross sections, the experiments of Chen et al. (2015Chen et al. ( , 2016) ) and numerical simulations of Magnini & Matar (2020); Magnini et al. (2022) have uncovered significant angular and axial non-uniformities in the liquid film thickness, whereby the effects of the rigid walls tend to flatten the interface at the channel's centre-line (h min ∼ Ca b ), depleting it of the liquid phase, and inducing the formation of large lobes protruding towards the channel corners in the cross-section (h ∼ Ca 0.445 b ).Furthermore, non-circular channels promote an exponential thinning of the liquid film along the bubble length from the nose to the rear section, as reported in Lister et al. (2006); Magnini & Matar (2020).The effects of other non-ideal geometries have been examined by Sauzade & Cubaud (2018), who carried out an experimental campaign to study the distortion of a train of bubbles in a constricting-expanding microchannel under varying Ca b and bubble packing conditions.This study provided evidence for the occurrence of bubble width-to-length ratio hysteresis across the sinusoidal structures of the channel, and the attenuation of these effects with decreasing Ca b .
A wide array of investigations have focused on illustrating the effects of (in)soluble surfactants on liquid-gas confined flows, placing particular emphasis on their influence on liquid film thickness, bubble deformation, and pressure drop (see Olgac & Muradoglu (2013) and references).The literature largely agrees that, under typical surfactant characteristics and negligible inertia, the liquid film displays a thickening response to surfactants of the order of 1 to 4 2/3 times that predicted by Bretherton (1961)'s expression for uncontaminated interfaces.This thickening has been reported to vary with Ca b , surface elasticity, and across the bubble length due to Marangoni-related mechanisms, which act in opposition to liquid drag forces as they push the surfactants towards the bubble back.Nonetheless, Ghadiali & Gaver (2003)'s numerical results in the context of pulmonary airways suggest the possibility of liquid film thinning surfactant actions under conditions of large bulk Péclet number, P e c >> 1, and Ca b > 10 −1 due to a bifurcation and enhancement of Marangoni stresses, which drive fluid away from the thin film region and towards the bubble nose.In terms of pressure drop and bubble velocity, it has been shown that surfactants tend to increase the pressure drop across the dispersed phase (Ginley & Radke 1988;Luo et al. 2019) and significantly reduce the bubble's mobility along its length (Borhan & Mao 1992).Batchvarov et al. (2020) carried out a comprehensive computational exploration of multiple surfactant characteristics and dimensionless groups with the novelty of introducing inertial effects (Re = 443 − 728) to the analysis.A key finding of this work is the surfactant-induced dampening of the bubble's rear interfacial oscillatory structures, together with the evidence provided to affirm that this process is entirely Marangoni-driven rather than a consequence of localised lowered surface tension.Moreover, this investigation addressed the effect of low solubility surfactants on the bubble's dynamics, detecting a clear arrangement of the bubble's morphology into two separate regions; one fully covered in surfactants and of reduced mobility, extending from the bubble rear to a portion in the stream-wise direction that depends on P e c , and one region of diminished surfactant concentration that stretches to the bubble nose, characterised by mechanisms comparable to those of clean interfaces.
A noteworthy sequence of events transpires in the bubble as Ca b is increased to a critical value that strongly depends on its size and Re.As reported in a number of investigations (Abubakar & Matar 2022;Sauzade & Cubaud 2013;Taha & Cui 2006;Giavedoni & Saita 1999), the reduced dominance of capillary forces has a distorting impact on the curvature of the front and back menisci, dragging both interfaces towards the main flow direction and causing the bubble to adopt an elongated, 'bullet-like' shape.As inertial and viscous forces increase, the bubble's back undergoes a curvature inversion that allows the liquid phase to infiltrate into the bubble's domain, forming a small cavity.The fate of this cavity is highly influenced by a complex interplay of phenomena that might culminate in a continuous cavity growth in the main axial or radial direction, or in its capillary breakup to encapsulate and entrap liquid drops within the bubble.This type of encapsulation has been observed experimentally in the works of Goldsmith & Mason (1963); Olbricht & Kung (1992); Olbricht (1996) in liquid-liquid flows.Similar findings are reported in a limited number of numerical pursuits on multiphase confined flows, including the investigations of Tsai & Miksis (1994) and Izbassarov & Muradoglu (2016) for contracting/expanding capillaries in Newtonian and viscoelastic fluids, respectively, Pozrikidis (1992) for buoyancy-driven flows of viscous drops, Nath et al. (2017) for liquid drops in creeping flow, and Andredaki et al. (2021) and Atasi et al. (2018) for gas bubbles in the absence/presence of surfactants, respectively.The latter three studies have provided an examination of the specific conditions that induce the encapsulation phenomenon, a general elucidation of the mechanisms at play, and an overview of the encapsulation aftermath.Nath et al. (2017) have remarked that, in the inertialess limit and considering fully circular channels, there exists a critical Ca b beyond which drop breakup and entrapment events will occur.The authors found this critical value to decrease (increase) with the initial drop-to-channel radius (dispersed-to-continuous phase viscosity).Andredaki et al. (2021) and Atasi et al. (2018) have recently laid the groundwork for the formulation of general encapsulation and breakup regime maps for Re << 1, bounded by W e B and the relationship between inertia and pressure drop across the channel in the former, and Ca b and the relative importance of interface surfactant adsorption/desorption in the latter case.Despite these initial efforts, much remains unknown and unexplored about these encapsulation events.Notably, a quantitative and phenomenological description of the underlying dynamics that govern each step of the process is needed, in combination with a thorough inspection of the cavity characteristics.It is also crucial to address the multiple examples of interfacial singularities inherent to the system at hand, which have yet to be analysed from the perspective of the well-established problem of capillary breakup, from which multiple parallels naturally emerge between our system and others (e.g., contracting liquid ligaments and inkjet printing).Likewise, an analysis of the evolution of the drops beyond pinch-off is required as it is a topic much less understood than the steps that precede it (Martínez-Calvo et al. 2020) and could be crucial to potential applications (Izbassarov & Muradoglu 2016).Finally, and as will be seen in this paper, extending the set of surfactant parameters explored from those of Atasi et al. (2018) uncovers new encapsulation and breakup outcomes, which allows us to expand and enhance the original regime maps.The main objective of this paper is therefore to carry out an extensive characterisation of the cavity formation and (post) pinch-off dynamics, taking into account the interaction between surfactants and inertia, capillarity and viscosity in a non-circular channel geometry, which correspond to non-idealities of interest, as seen in the foregoing literature review.
The rest of the paper is organised as follows: §2 introduces our numerical approach, including the governing equations, boundary conditions, solution methods, computational setup, and a validation for our system.Our results are presented and discussed in §3, where we focus first on encapsulation in clean interfaces, on the effects of surfactants next, and finalise with a formulation of encapsulation regime maps that compile all behaviours observed.

Numerical modelling and non-dimensionalisation
We consider a horizontal liquid-filled square capillary of width and height D in a Cartesian three-dimensional domain, x = (x, y, z), with an elongated gas bubble propagating through its interior (see figure 1).Assuming incompressible flow, Newtonian fluids, and negligible gravitational effects, we perform direct numerical simulations (DNS) based on the two-phase Navier-Stokes equations and a hybrid front-tracking/level-set method to handle the interface and surface tension forces.This method, as formulated and implemented in Shin & Juric (2002), is coupled with the resolution of surfactant transport and exchange between the interface and the liquid phase bulk, as described in Shin et al. (2018).The relevant variables involved in the system have been rendered dimensionless (denoted by tildes) by using the scalings depicted in Eq. (2.1), where u, t, p, σ, Γ , C, and C s represent velocity, time, pressure, surface tension, interfacial surfactant concentration, bulk surfactant concentration, and bulk surfactant concentration in the vicinity of the interface, respectively.The physical parameters included in the scaling correspond to the width of the channel, D, the average inlet velocity of the liquid, U a , the density of the liquid phase, ρ l , the surface tension in a surfactant-free interface, σ s , the saturation interfacial concentration, Γ ∞ , and the initial bulk surfactant concentration, C ∞ , in line with the scaling proposed by Batchvarov et al.
(2020) for a similar system.In what follows, we refer to each variable by its name to refer to its dimensionless version, unless stated otherwise and add the subscript b to signify that the variable is reported at the bubble nose.The governing mass and momentum equations are written in dimensionless form and according to a 'one-fluid' formulation, as shown in Eq. (2.2)-(2.3), where I (x, t) is a smoothed Heaviside function that adopts the value of zero in the gas phase and unity in the liquid phase, ρ is the density and µ the viscosity of the fluids.
Here, we denote the liquid and gas phases with the subscripts l and g, respectively.
The normal and tangential components of the surface tension forces are represented by the last two terms on the right-hand-side (RHS) of Eq. (2.2), wherein κ corresponds to twice the mean interface curvature, n to a unit normal to the interface, δ x − xf to a Dirac delta function that is zero everywhere except for the interface (located at x = xf ), Ã( t) to the time-dependent interface area, and ∇ s to the surface gradient operator (Shin & Juric 2009;Shin et al. 2017).We adopt the convention that positive n vectors point outwards from the interface towards the liquid phase.Accordingly, positive values of κ describe a convex interface, as exemplified in figure 1.The mass conservation equations of surfactant species at the interface and bulk are given by Eq. (2.4)-(2.5)and the source term representing surfactant exchange between the interface and the bulk region immediately adjacent is given by Eq. (2.6): where ũt = (ũ s • t)t is the projection of the interface velocity vector, ũs , on the interface unit tangent, t.The dependence of surface tension on local interface surfactant concentration is represented by a non-linear equation of state derived from Langmuir adsorption isotherm, as expressed in Eq. (2.7) where β s = RT Γ ∞ /σ s is the surfactant elasticity number; R and T are the thermodynamic ideal gas constant and temperature, respectively.The non-linear equation of state described above is valid for very dilute systems in which Γ << Γ ∞ .As Γ increases, the equation of state yields unphysical negative values of surface tension.To circumvent this, a limiting value for σ, σ = 0.05, has been introduced, in accordance with Muradoglu & Tryggvason (2014) and Shin et al. (2018).The dimensionless groups that appear in the above equations and that characterise the system are defined as follows: where Re and Ca are the Reynolds and Capillary numbers, respectively (W e = CaRe, and Oh = Ca/Re).P e c and P e s are the bulk and interfacial Péclet numbers, which provide a measure of the importance of inertia relative to surfactant mass diffusion across the interface (modulated by a diffusivity D c ) or the bulk (modulated by a diffusivity D s ), respectively; Bi represents the Biot number, describing the competition between convective, t conv = D/U a , and interface surfactant desorption, t des = k −1 d , time scales.The Damköhler number, Da, provides a dimensionless measure of the initial interface saturation, and the number k reflects the interplay between the characteristic time scales of interface surfactant adsorption, t ad = (k a C ∞ ) −1 , and desorption.
For an exhaustive description of the discretisation schemes and numerical framework employed for solving the above mentioned governing equations, we refer the reader to Shin et al. (2018Shin et al. ( , 2017)), and to Shin & Juric (2009) for a comprehensive account of the hybrid front-tracking/level-set method, also called Level Contour Reconstruction Method (LCRM).Here, we present a short summary of the most relevant aspects.The LCRM considers the combination of a fixed structured Eulerian grid for the resolution of field equations and a moving and deforming Lagrangian grid for the interface, discretised via an unstructured triangular mesh.These interface elements are consequently advected through the integration of the Lagrangian equation dx f /dt = v, where v corresponds to the interface velocity, which is interpolated from the Eulerian grid.This integration is carried out with a second-order Runge-Kutta method.The spatial derivatives of the fields in the Eulerian grid are discretised through a standard cell-centred scheme for all terms, with the exception of the non-linear convective terms, for which a secondorder essentially non-oscillatory (ENO) procedure is used (Shu & Oshert 1989;Sussman et al. 1998).For the viscous term in the momentum equation, a second-order centred difference scheme is employed.With regard to the temporal terms, a second-order Gear method (Wang & Wen 2006) is adopted with an implicit time integration for the viscous terms.The time step is set to be adaptive, according to the following criterion: ∆t = min{∆t cap , ∆t vis , ∆t CFL , ∆t int }, where the time steps relate to viscosity, capillarity, Courant-Friedrichs-Lewy condition, and interface, respectively.The time-step related to the interface was found to be the limiting time-step in all of our simulations (∆t int ∼ O(10 −7 )s).

Simulation setup and validation
As depicted in figure 1, our simulation domain consists of a three-dimensional rectangular box of dimensions 27D × D × D, corresponding to the x (channel length), y (channel height), and z (channel width)-directions, respectively.The value of the channel length was selected to ensure the attainment of a steady-state bubble propagation for W e < 10, following the results of Batchvarov et al. (2020) and Magnini & Matar (2020).Under conditions of W e > 10, highly unsteady encapsulation/bursting phenomena are detected in the system.The channel length chosen allowed for the observation of all relevant phenomena, including the formation of the liquid cavity, drop encapsulations and coalescence events, and bubble bursting, for all conditions tested.The elongated bubble Figure 1: (not to scale) Schematic of simulation setup, re-entrant liquid cavity at the bubble back, and encapsulated drop.The direction of outwards pointing normal vectors, liquid film thickness, h(x), and curvature sign according to the adopted convention, are also highlighted.The location of the Cartesian axis represents the outset of the geometrical domain (x = 0, y = 0, z = 0) was initialised in quiescent conditions close to the channel inlet and at the centre-line of the y and z-directions as a horizontal cylinder of length 3.12D with two hemispherical caps of diameter 0.94D on each end.A fully-developed liquid velocity profile was imposed at the channel inlet with a Neumann condition for pressure (∂p/∂n = 0).All walls were treated as no-slip boundaries, and Neumann conditions were imposed for all variables at the outlet.The bubble-liquid interface was initially covered in surfactant with uniform concentration Γ 0 , determined from an initial equilibrium state between surfactant adsorption and desorption (last term on the RHS of Eq. (2.4) equal to zero), which, assuming C s = C ∞ initially, reads as: (2.9) This set-up closely resembles the approaches followed in previous numerical investigations for similar systems (Atasi et al. 2018;Nath et al. 2017;Luo et al. 2019).
We explore the effect of the dimensionless numbers and flow parameters that characterise the system (see Eq. (2.8)) considering the following ranges: Re = 100 − 443, Ca = 0.0089 − 0.0693 (W e = 3.94 − 30.70,Oh = 4.48 × 10 −3 − 1.25 × 10 −2 ), P e c,s = 100, Bi = 0.10 − 1, Da = 0.01 − 1, k = 0.10 − 10.00, and β s = 0.25 − 1.For all cases simulated, the viscosity and density ratio between the phases is kept constant and in line with representative values for water and air: ρ l /ρ g ∼ O(10 3 ) and µ l /µ g ∼ O(10 2 ).We define a 'base' case from which our parametric sweep was carried out with the following conditions: Re = 443, Ca = 0.0693, P e c,s = 100, Bi = 0.10, Da = 1, k = 0.10, and β s = 0.50.Unless stated otherwise, all results correspond to simulations under these conditions.The selection of the testing ranges was based on common values encountered in bubbly flow systems in the presence of surfactants, as reported in Atasi et al. (2018), as well as values of W e that would allow a significant disruption of the bubble back to generate the aforementioned re-entrant cavity and consequent drop encapsulation (Giavedoni & Saita 1999;Andredaki et al. 2021).For typical ionic and non-ionic surfactants in water, such as sodium dodecyl sulphate (SDS), N-dodecyl-N,N-dimethylammonio-3-propane sulfonate, and Triton X-100 (TX100), the value of Γ ∞ ranges between O(10 −6 − 10 −5 ) mol/m 2 and D s between O(10 −12 − 10 −8 ) m 2 /s, which results in P e s = O(10 3 − 10 6 ) (Kalli et al. 2022;Constante-Amores et al. 2021c).As noted previously in Constante-Amores et al. (2021a) and Batchvarov et al. (2020), saturation of the system's dynamics is reached above P e s ∼ 100.P e c was set to be equal to P e s , following the study of Agrawal & Neuman (1988).
Based on the work of Kamat et al. (2018) and Constante-Amores et al. (2021a) and our range of characterising Oh, we define an inertial-capillary or Rayleigh time scale as t in−cap = ρ l D 3 /σ s for Oh << 1 and a Marangoni time scale as t m = (µ l D)/∆σ, where ∆σ = σ s − 0.05σ s is a measurement of the surface tension gradients.The ranges for the characteristic time scales of the system are estimated to be: 10 −3 − 10 −2 )s, and t m ∼  O(10 −4 )s (for surfactant-laden cases).These time scales ensure that sorptive/desorptive, inertial, capillary, and Marangoni phenomena will play an important role in the dynamics of the system.
Following numerous studies making use of the numerical methods described above for inertial and capillary phenomena (see Constante-Amores et al.On account of ensuring this resolution produces mesh-independent results, we have performed a mesh independence study in surfactant-free and surfactant-laden conditions, detailed in Appendix A along with other resolution considerations.The mesh size selected also allowed to capture accurately the thin liquid film surrounding the bubble, which, according to the works of Gupta et al. (2009) and Magnini & Matar (2020), requires a minimum of 5-10 computational cells covering its domain for a correct development of the liquid velocity profile.In our simulations, we have ensured the placement of at least 5-6 cells in the liquid film section for the cases that exhibit encapsulation (W e ∼ 30), for which h(x)/D ∼ 5 −2 (see figure 11 of Magnini & Matar (2020)).This condition is further achieved in the surfactant-laden simulations, which develop a thicker liquid film.We make the note, however, that the size of the entire domain, as well as the comprehensive sweep over dimensionless numbers presented, make it prohibitively expensive in terms of computation to achieve the refinement levels reported in works such as Castrejón-Pita et al. (2015); Li & Sprittles (2016) to fully capture all regime transitions leading to the interfacial singularities of pinch-off.
The numerical framework and simulation setup presented herein have previously been carefully validated in the context of elongated bubbles in circular capillary channels.For this validation, we refer to the work of Batchvarov et al. (2020), in which the code was benchmarked against the well-known correlation of Han & Shikazono (2009) for steady-state liquid film thickness in conditions of non-negligible inertia.For Ca ∼ O(10 −2 ) and Re ∼ O(10 2 ), deviations of up to 10% were reported, which are within the uncertainty of the correlation.A supplementary validation for square capillary channels is conducted using the numerical data and conditions described in Magnini & Matar (2020) for surfactant-free scenarios under non-negligible inertia.Satisfactory quantitative agreement between our results and those of Magnini & Matar (2020) is achieved in terms of bubble-to-liquid velocity ratio, U b /U l , and gas area fraction, (see figure 2(a)-(b)), with maximum deviations of 2.7% and 3.1%, respectively.We underline that the cases that do not have a direct counterpart from the above mentioned reference (i.e., Ca = 0.0089, Ca = 0.0377, and Ca = 0.068 for Re = 443) adequately follow the qualitative trend of increasing U b /U l and decreasing with Ca and Re.
Figures 2(c)-(d) exemplify the well-known bubble morphological regimes found in noncircular channels, as well as our model's correct capturing of the non-axisymmetric bubble shapes that arise below a threshold Ca. (Kolb & Cerro 1991;Hazel & Heil 2002;Chen et al. 2015).In accordance with Ferrari et al. (2018), the corner flow effects exhibited in Ca 0.01 promote the development of a saddle-shaped interface nearby the channel centre lines (see outermost series in figures 2(c)-(d)), which brings about the formation of four lobes projecting towards the channel corners.As portrayed in figure 2(c), these lobes are located furthest away from the channel centre-point (highest r/R h ) and will continue to grow as the saddle deepens (r d > r c ) with decreasing Ca.The presence of the saddle shifts the minimum liquid film thickness, h min , from the centre-line of the channel to a position, l y , that varies as a function of Ca.We compute l y /R h = 0.2 for the case highlighted in the figure, which matches remarkably well with the value reported by Magnini & Matar (2020) for W e << 1 (l y /R h ≈ 0.195, see figure 7(b) of the reference).The high pressure zones at the saddles showcased in figure 2(d) are evidence of significant capillary flow drainage away from the thin-liquid film region, as reported in Magnini & Matar (2020).
The other morphological regimes represented in figures 2(c)-(d) correspond to an additional non-axisymmetric cross-section regime for 0.01 Ca 0.05, characterised by r d > r c and the absence of the saddle (l y = 0, see Ca = 0.02 in the figure), and a fully symmetric bubble topology (r d = r c , see Ca = 0.1 in the figure) for the final regime.We draw attention to the almost perfect alignment of these regimes and the ranges of Ca in which they occur between our numerical predictions and those of the literature, notwithstanding that the regimes documented in the literature were observed under conditions of vanishing inertia.

Results and discussion
Following the description of the numerical methods, simulation setup, and validation, the results of our parametric study are presented in three main parts.The first ( §3.1) describes and analyses the underlying mechanisms responsible for the changes in curvature of the interface at the bubble tail and nose observed for increasing Ca in surfactant-free cases.These changes in curvature result in a re-entrant cavity and subsequent drop encapsulation, which are illustrated in detail.§3.2 expands on these concepts by incorporating into our analysis of encapsulation the effect of several different surfactant parameters related to adsorption/desorption kinetics.Finally, §3.3 explores the multiple non-axisymmetric structures that develop at the highly-deformed bubble tail after encapsulation for cases where adsorption phenomena are dominant.

Bubble dynamics in surfactant-free cases: tail undulations, curvature disruptions, and encapsulations
We start our discussion by showcasing the bubble topology at a two-dimensional projection in the (x,y) plane at the channel centre-line (z = 0.5), as well as the normalised mean curvature, k, and the normal, τn , and tangential, τt , components of the viscous stress along the interface (see figures 3(a)-(b)).Here, these quantities are defined as k = k/k c , τn = τ n /(D/U a ), and τt = τ t /(D/U a ), respectively, where rate of deformation tensor), and k c corresponds to the curvature of a sphere with a volume equal to the initial bubble volume.Under the sign conventions adopted, τn > 0 denotes stresses exerted on the interface towards the liquid phase.Conforming to previous observations, such as those by Taha & Cui (2006); Giavedoni & Saita (1999); Sauzade & Cubaud (2013); Abubakar & Matar (2022), which uncover the prominent impact of Ca on bubble shape, our plots reveal a loss of sphericity (typical of fully dominant capillary forces) at the bubble front and back with increasing Ca, in conjunction with an overall axial bubble elongation.The interface at the bubble front sharpens and expands, whilst the bubble back becomes progressively flattered with Ca, adopting a 'bullet-like' shape (κ f ront > κback ).This process continues until a curvature inversion emerging from ỹ = 0.5 is detected for Ca = 0.0693, which eventually leads to a drop encapsulation.
Detailing the behaviour of κ for Ca = 0.0693 along the interface path going from the cavity vertex to the bubble nose, a zone of negative curvature is first observed at 3.8 x 4.1, followed by a positive curvature meniscus at 3.3 x 3.8.This particular curvature sign inversion is a manifestation of the early stages that precede pinch-off and drop encapsulation, in which the fragmenting liquid filament undergoes the formation of a discernible pinching neck with a vanishing radius.This process, and its connection to the well-established scaling laws, are examined further along in this paper.Exiting the interior of the liquid cavity towards the region adjacent to the liquid film (3.3 x 7.0), a non-uniform low curvature zone that encompasses the majority of the bubble is identified, later culminating in the high-curvature bubble nose seen for x > 7.0.Contrasting this evolution to that of the two lowest Ca, it can be noticed that the spatial location that marks the beginning of the bubble nose (sudden curvature rise from the thin-liquid film region, highlighted in gray in figure 3(a)) shifts to lower x for decreasing Ca.Notable curvature inversions at the bubble tail are also identified for the low Ca cases, expressed in the undulatory structures emphasized in figure 3(b).These interfacial waves, previously reported in Magnini et al. (2017); Edvinsson & Irandoust (1996); Giavedoni & Saita (1999), moderately diminish in frequency and amplitude with Ca at the same Re, as illustrated in the figure .The dynamics of the bubble-liquid interface and its curvature gradients are readily mapped onto the τn,t spatial profiles of figure 3. Pertaining to the trailing undulatory structures (shown in yellow and green for Ca = 0.0089 and Ca = 0.0377, respectively), it is found that these are first characterised by a sudden surge in τn that terminates in a local maximum, which pushes the interface towards the channel walls and thus creates the crests of the interfacial waves.These peaks are promptly followed by a rapid τn decrease approaching a local minimum that contributes to the formation of the wave troughs by pulling the interface back towards the bubble.Concurrently, the tangential stresses feature inverse patterns of local maxima/minima that promote the sharpening of the crests/troughs.A substantial drop in τn,t , mirrored by κ, is observed by departing from the rear undulatory region to the high-pressure central section of the bubble.τn,t remain quasi-constant and close to zero throughout this section until the reappearance of local maxima/minima marks the onset of the bubble nose.Note the generally higher |τ t | exerted on the bubble nose at ỹ = 0.5 for decreasing Ca, partly contributing to a lower curvature (more spherical) bubble front.Similar remarks about the close relationship between curvature, capillarity, and interfacial stresses have been made by Atasi et al. (2018).
Figure 4(a) shows a three-dimensional representation of the interface at t = 0.89, coloured by the x-component velocity relative to the bubble nose, ũref,x = (u x − U b )/U b .Figure 4(b) depicts the azimuthal component of the normalised vorticity field, ωz = ω z /(D/U a ) (top half), and the flow topology parameter, Q (bottom half) as contour plots over an (x, y) plane for the three Ca studied.Here, Q, which measures the influence of the bubble-liquid interface on the flow field and relates this information to the dominant flow type, is defined as follows (Soligo et al. 2020)  The effects of this are seen in the progressive shrinkage of purely rotational flow areas and their substitution by shear and elongational flows for high Ca.These larger areas of pure elongation/shear as opposed to rotation provide partial evidence to explain the axially longer bubbles that result from higher Ca.
To further investigate this point, we refer to the streamlines depicted in figure 4(c) for the two limiting Ca, where a few key differentiating features between the two cases are unveiled.From left to right we first recognise two stagnation points near the channel walls and almost reaching the bubble rear for Ca = 0.0089.These stagnation points are also observed for Ca = 0.0693, albeit detached from the interface and closer to the channel centre.The trailing interface undulations identified for low Ca are emphasised by the vortical streamlines displayed, markedly absent in high Ca.Moving along the x-axis, we arrive at the core section of the bubble, which features a pair of central vortices covering much of the regions nearby the interface.Tying in the above analysed rotation patterns and dominating flow phenomena at the bubble nose to the arrangement of the streamlines, we point to the two vortices that develop at the bubble nose in Ca = 0.0089, and that have vanished in Ca = 0.0693.Lastly, we note that the streamlines in the liquid phase ahead of the bubble represent the occurrence of a 'bypassing effect' at high Ca, in which the streamlines crossing the thin liquid layer travel in the main flow direction without recirculating back to the bubble, as seen at low Ca.The above described features of the streamlines surrounding the bubble, and how these are influenced by Ca, are consistent with the findings of Taha & Cui (2006); Abadie et al. (2013); Atasi et al. (2018).
We now proceed to detail the encapsulation process by first illustrating the numerous flow singularities produced by high Ca conditions in uncontaminated interfaces.In figure 5(a) we plot the temporal evolution of the maximum normalised velocity in x-direction, M ax ũx .This evolution reveals five major local peaks in Ca = 0.0693 in the interval 0 < t < 3.5, which materialise in (i) the previously analysed curvature inversion at the bubble back, (ii)-(iii) two pinch-offs that give rise to a 'satellite' and a primary drop, (iv) an abrupt retraction and subsequent expansion of the former, and (v) a coalescence event between the two drops (see figure 5(b)).A first description of the mechanisms by which these interfacial singularities come about is shown in figure 5(c), where a cavity deepening in the positive x-direction takes place at 0.44 < t < 1.11.The fate of this liquid cavity can be understood through the hybrid lens of the well-known problems of liquid filament breakup surrounded by a gas phase (Eggers & Villermaux 2008;Castrejón-Pita et al. 2012) and liquid dripping faucets (Ambravaneswaran et al. 2004).Our liquid cavity at t < 1.11 can be thought of as a quasi-cylindrical inviscid (Oh = 0.0125) liquid filament with a low ratio between its length, 2L cavity , and radius, R cavity .As in the systems studied by Wang et al. (2019a), the radius of our cavity is not small enough for viscous drag effects from the gas phase to be significant (R cavity ∼ O( 10 5(c), we observe a halt in cavity axial growth as capillary pressure builds at the bubble rear, inducing a vertical expansion towards the channel centre-line (ỹ = 0.5), creating an essentially flat bubble back at t = 1.99 and a nascent pinching neck characterised by a high |ω z | ring.The interface areas covered by the ring are further pulled in the counter-flow direction, prompting an additional curvature inversion at the bubble rear and the first pinch-off thereupon (see t = 2.22).This capillary singularity promotes a swift and opposite-direction recoil of the bubble back and the newly-entrapped liquid structure, thus forming the globoid-shaped bubble rear observed in t = 2.39 − 2.63 (notice the change in direction of rotation from negative to positive ωz at the pinching neck between t = 2.22 and t = 2.30).The rapidly forming rear pinching neck witnesses the incipience of an additional capillary neck at the interior of the bubble (see highlight and high ωz ring in t > 2.22).This concludes in a second pinch-off and a precipitous recoil into a satellite drop (see t = 2.49).Henceforth, we will refer to the first of these pinch-off events as 'back pinch-off', and to the second as 'interior pinch-off'; the times for these events are denoted as tp−o,bk and tp−o,int , respectively.These phenomena of successive back and interior pinch-offs in surfactant-free cases are in agreement with the numerical observations of Andredaki et al. (2021).
The post second pinch-off dynamics are driven by rich rotation mechanisms around both the satellite and primary drops.The two pinch-off events create a small contracting liquid filament, characterised by local Oh ∼ 0.05 and L ligament /R ligament ∼ 2.62 and similar to those studied by Notz & Basaran (2004); Driessen et al. (2013), and others.This small aspect ratio, in conjunction with the strong contracting kinetic energy derived from the second pinch-off, overcome capillary pressure forces to avoid an additional pinchoff instance in the small ligament.Observing the vorticity contours and the highlights at t = 2.39 − 2.49, an example of a vortex ring detachment can be noticed at the forming neck of the filament, where the vorticity boundary layer detaches and advects towards the centre of the bulbous region, pushing flow away from the neck and promoting the escape from pinch-off (Hoepffner & Paré 2013;Wang et al. 2019a;Constante-Amores et al. 2020).As seen in t = 2.57, succeeding the small liquid filament retraction in the counterflow direction ( t = 2.49) we spot a counter expansion (see positive ωz at t = 2.57) that generates the bullet-like satellite drop depicted at t = 2.63 and that eventually coalesces with the primary drop.The process of incorporation of the satellite into the primary drop exhibits the progressive growth of the finite liquid bridge (Eggers et al. 1999), driven by the two opposite-signed high |ω x | rings in its vicinity (see t = 3.10) as well as complex dynamics of interfacial waves across the newly-formed drop (see t = 3.19 − 3.37).

Influence of surfactant parameters on cavity formation and encapsulation dynamics
Having examined in detail the different encapsulation phenomena observed for high Ca in surfactant-free interfaces, we continue our analysis to explore the influence of key surfactant dimensionless groups.We begin by highlighting the influence of Marangoni stresses on the cavity shape, followed by describing the relative effect of β s and Bi.Thereafter, we analyse the competition between the rates of interface surfactant adsorption and desorption, represented by k and Da simultaneously, as well as Re.

Marangoni stresses, surface elasticity, and Biot number
To begin our analysis, we elucidate the behaviour of a surfactant-laden system under the conditions ascribed to the 'base' case (refer to §2.2) and compare it to that of surfactant-free interfaces.In order to isolate the effects of lower surface tension from those arising from Marangoni stresses on the encapsulation dynamics, we have set up an additional case, denoted by |τ m | = 0, in which we have suppressed the last term on the RHS of Eq. (2.2) to inhibit Marangoni stresses while simultaneously allowing the reduction of surface tension.Figures 6(a)-(b) show the shape of the liquid cavity and a few other characteristics along its interface at t = 2.22.For the sake of clarity, we have duplicated the plots of interfacial shape on both column panels.Examination of the plots uncovers an almost surfactant-free interface at the liquid cavity for |τ m | = 0, which exhibits a significant drop in Γ going from the edge to the interior of the cavity in the stream-wise direction and a slight surge after surpassing the incipient pinching neck.This rapid loss of surfactant for |τ m | = 0, as opposed to the Marangoni-supported case, |τ m | > 0, is explained by noticing that the operating conditions of the base case promote surfactant desorption (large Da and small k) as Γ → 1 (see Eq. (2.6)) as Figure 6(c) portrays a key difference in the encapsulation behaviour between the three cases, where the succesion of back pinch-off first and interior pinch-off second observed for clean interfaces is reversed in the presence of surfactants and the overall pinchoff dynamics are delayed (see the figure caption for the pinch-off times).This delay  6(a), where these tangential stresses counteract capillary forces at both the back and interior pinch-off regions.These observations are in complete agreement with the conclusions of Kamat et al. (2020Kamat et al. ( , 2018)); Constante-Amores et al. (2020) for liquid threads, although we note that in our system the Marangoni forces exerted on the interface are not sufficient to escape pinch-off completely.The inversion of back and interior pinch-off times in contaminated interfaces can be attributed to the noticeably non-uniform surfactant distribution along the liquid cavity, which displays its highest values in the regions neighbouring the capillary neck at the back of the cavity, lowering surface tension and providing additional disruption to neck closure in comparison to the interior cavity neck.
For the first part of our parametric investigation, we study the effect of β s on encapsulation and the post pinch-off dynamics.As figure 7 printing, increasing surfactant strength gradually increases the speed of cavity formation and decreases the rate of neck thinning at both the back and interior of the cavity (see locations of cavity nose in figure 7(a) and pinch-off times in the figure caption).This delay is materialised in the larger |κ| peaks seen for β s = 0.25 at the incipient pinching necks, as well as the higher |τ m | that arise at their vicinity, seeking to counteract capillary draining (Ambravaneswaran & Basaran 1999).Figure 7(b) illustrates the progression of the primary encapsulated drop as it travels across the bubble domain in the form of its volume-averaged x-velocity, u x,d , normalised by U b (top), interfacial area, A d , normalised by its area after pinch-off, A d,0 (middle), and ratio between drop length in the x and ydirections, a/b (bottom).The encapsulated drops exhibit three common characteristics irrespective of the presence of surfactants; namely, a higher stream-wise velocity than the bubble nose, highlighting the potential for a scenario in which the bubble nose interface is ruptured by the drop, an overall decrease in interfacial area, and dampening a/b oscillations in time, which culminate in virtually spherical drops.This bubble nose rupture is further analysed in §3.3.Under the conditions considered, the addition of a strong surfactant (β s = 1) decreases the rate of surface area reduction and significantly increases drop velocity relative to the bubble nose.The non-monotonic behaviour of u x,d /U b with regard to a clean interface and weak surfactants (β s = 0.25 − 0.50) can be explained by observing in figure 7(c) the spatial distribution of Γ , σ and τm in the encapsulated drop at t = 7.82.Despite the largest | τ m | opposing drop propagation in the stream-wise direction for the strongest surfactant, it is seen that these effects are fully countered by the lower surface tension across the drop's domain, yielding a higher drop velocity than the surfactant-free case.In contrast, the weak surfactant cases feature an almost clean interface with surface tension values very close to the clean case, but with non-zero Marangoni effects that induce a delay in drop propagation.
The morphological deformations that characterise drop behaviour after pinch-off, as depicted in figure 7(b)(bottom), demonstrate the presence of a semi-periodic oscillating pattern monotonically influenced by surface elasticity.In the inset of this figure, a higher deviation from a fully spherical drop can be appreciated for the surfactant-free case in both the major crests and troughs of the oscillatory structures.Considering these results, we draw parallels with the study of Wang et al. (2019b) in an ink-jet system, whose main results suggest that the extent of deformation decreases with Ca (after a critical Ca has been reached), in line with our numerical observations for increasing β s .Estimating the frequency of a/b oscillations through the classic expression of Rayleigh (1879) (f ∼ σ/ρ c r 3 d , where r d is the drop's spherical equivalent radius) at the latest time recorded ( t = 13.29),we obtain f ∼ 1.49, 1.55, 1.55, 1.57 for surfactant-free and β s = 0.25, 0.50, 1, respectively, implying a dominating effect of lower drop radius over lower surface tension on the frequency of deformation cycles for increasing β s .
The effects arising from altering Bi in the range 0.01 − 1 are depicted in the plots of figure 8.An interface covered by a highly soluble surfactant (i.e., high Bi) undergoes rapid rates of mass exchange between the interface and bulk as a result of its comparatively low (high) characteristic desorptive (inertial) time.This is demonstrated by noticing the significantly lower Γ distributed along the liquid cavity for the highest Bi considered, which, barring the surfactant-free case (Bi → ∞), exhibits the slowest rate of cavity formation and the fastest pinch-off events (see location of cavity nose in figure 8(a)(top) and figure caption for the pinch-off times), in alignment with preceding investigations (see Milliken & Leal (1994)).Similar to our previous remarks, the signature |τ m |, |κ|, and capillary pressure peaks (not shown) that signal the onset of capillary neck closure Figure 7: Effect of β s on cavity formation and encapsulation.(a) Cavity shape, surfactant interfacial concentration (left axis, continuous line) and surface tension (right axis, dotted line), Marangoni stresses, and interface curvature at t = 2.22.(b) Temporal evolution of primary encapsulated drop velocity normalised by bubble velocity, primary drop interfacial area normalised by area after pinch-off, and ratio between drop length in the x and y directions.Data points computed for all cases every 10 −3 s.(c) Primary encapsulated drop shape, surfactant interfacial concentration (left axis, continuous line), and surface tension (right axis, dotted line) and Marangoni stresses at t = 7.82.Plots on two-dimensional projection in x-y plane (z = 0.5) for (a,c).tp−o,bk = 2.28, 2.52, 2.73, 2.90 and tp−o,int = 2.45, 2.64, 3.02, 3.94 for surfactant-free, and β s = 0.25, 0.50, 1, respectively.All other parameters remain unchanged from those specified in §2.2 for the base case and the advent of pinch-off are already developing in a noticeable manner for Bi = 1, in contrast to Bi = 0.01 − 0.10.
The post pinch-off temporal evolution of the primary drop, as depicted in figure 8(b), reveals similar trends to those referenced in the above analysis of β s in terms of drop shrinkage, its attenuation in surfactant-laden cases, and semi-periodic drop size oscillations.Nonetheless, we observe the incidence of non-monotonic dynamics with regard to Bi and the surfactant-free case.In particular, the highest values of u x,d /U b appear to occur for the most soluble surfactant, followed by the surfactant-free case, leaving the less soluble cases with the lowest u x,d /U b .The disruption of the direct relationship between higher pinch-off times and higher u x,d /U b seen for β s suggests a complex interplay between Marangoni stresses, surface tension, pinch-off times, and drop velocity with varying Bi.The spatial locations of the encapsulated drops at t = 7.82 (figure 8(c) (top)) uncover an additional layer of non-monotonic behaviour, where the overall fastest (slowest) drop corresponds to Bi = 1 (surfactant-free), despite these two cases being the closest in terms of pinch-off times and local σ.The large velocities exhibited by Bi = 1 in comparison to lower Bi and the surfactant-free case is explained by noticing its virtually clean interface and thus the complete absence of Marangoni stresses to counteract its axial motion at later times, while also having a low enough σ at the earlier times to increase the rate of cavity formation and delay pinch-off.

Effect of surfactant adsorption/desorption kinetics (Da and k) and Re
We now proceed to examine the influence of directly contrasting the characteristic times of adsorption and desorption, along with the initial interface saturation through varying k and Da simultaneously.This concurrent variation allows us to maintain all other parameters constant.Figure 9(a) portrays a monotonic response of the cavity depth and Γ with k and Da, where the rapid rates of surfactant adsorption inherent to the high k (low Da) cases are materialised in the relatively large and comparatively constant Γ distributions along the cavity.In line with our previous results, general reductions in local surface tension at the bubble rear, brought about by higher k or β s and lower Bi conditions, promote a faster and deeper infiltration of liquid into the bubble domain.In figure 9(b) we record the time until the first pinch-off instance at the bubble interior as a function of k for the entire parameter space explored.A retardation of the capillary instabilities that lead to end-pinching is seen as k is increased in the range 0 (surfactantfree) < k < 1, consistent with our previous remarks about the combined effect of σ and Marangoni stresses on opposing drop encapsulation.This behaviour is mirrored by the τm and κ peaks exposed in the vicinity of the nascent pinching necks for the lowest k cases displayed in figure 9(a) (k = 0.1 − 1.0), while the oscillations seen for k = 10 correspond to capillary waves in the liquid cavity, as reported by Constante-Amores et al. (2020) for long liquid ligaments.
Similar to our previous analyses of β s , Bi, and Marangoni-free case for k = 0.10, it is found that increasing k in the range k < 0.10 results in a sequence inversion of the two main pinch-off events ( tp−o,bk > tp−o,int ) identified in clean interfaces ( tp−o,bk < tp−o,int , see caption of figure 9).Interestingly, it is also found that for all cases above k = 0.25, the large accumulation of surfactants nearby the bubble rear provides a strong resistance to capillary pressure buildup in those regions, resulting in the partial elimination of the 'back pinch-off' event until bubble bursting events occur.This will be further discussed in the following section.Another notable feature to highlight from figure 9(b) is the sharp decrease in tp−o,int above the critical value of k = 1.This discontinuity in the function tp−o,int vs. k, in conjunction with our observations about the elimination of the back pinch-off suggest a partition of the encapsulation dynamics into two major regions according to the rates of adsorption/desorption: i) the above described region for k < 1, in which the presence of surfactants exerts a stabilising effect, leading to a delay in neck thinning in relation to a fully clean interface; and, ii) a region of rapid cavity formation in the axial direction alongside a significant acceleration of the first interior pinch-off event for k > 1.The first region is further segregated by the occurrence (k < 0.25, region i.a) or elimination (k > 0.25, region i.b) of the back pinch-off.An explanation for the marked transition from the first to the second region in terms of thinning rates is as follows.A general increase of surfactant coverage due to high rates of adsorption introduces high All other parameters remain unchanged from those specified in §2.2 for the base case rates of viscous and inertial deformation to the interface, as noted by Bazhlekov et al. (2006), hence, increasing the axial velocity of the cavity and its aspect ratio.Recalling the observations of Wang et al. (2019a) and Ambravaneswaran & Basaran (1999) in contacting/expanding liquid filaments, we show in figure 9(c) how in cases above the critical k, L cavity /R cavity surpasses a given threshold value to enter the beak-up regime and accelerate end-pinching, even if the local values of surface tension are lower than those of smaller k conditions.
Moving on to the final exploration of the dimensionless space of the system, we test the influence of varying Re in the range 100 − 443 (W e = 6.93 − 30.70) and maintaining all other characterising numbers of Eq. (2.8) constant: Ca = 0.0693, P e c,s = 100, Bi = 0.1,  (bottom).The values of R cavity were calculated at the midpoint between the cavity nose and the start of the cavity.Plots on two-dimensional projection in x-y plane (z = 0.5).tp−o,int = 2.45, 2.64, 2.93, 2.96, 2.05, 2.00, 2.03 for surfactant-free and k = 0.10, 0.25, 1, 2.50, 5.00, 10.00 (Da = 1, 0.40, 0.10, 0.04, 0.02, 0.01), respectively.A few cases of the full range of k have been omitted in (a) for an easier visualisation.All other parameters remain unchanged from those specified in §2.2 for the base case Da = 0.1, k = 1, and β s = 0.5.Figure 10 demonstrates the pronounced effects of Re on cavity formation and the overall surfactant dynamics, where the presence of comparatively lower inertial forces of reduced Re and constant Ca (lower W e) appears to oppose the complete development of the liquid cavity and the process of drop entertainment, whilst still allowing the rear curvature inversion described for large Ca (see §3.1).The spatial profiles of Γ show a cavity largely covered in a uniform manner by surfactants for Re = 100, manifested in the very low τm exerted and the virtually planar κ of a significant portion of its cavity.This contrasting behaviour against the highly non-uniform and smaller Γ describing high Re provides evidence for arguing that, in the locality of the cavity, inertio-capillary forces are greatly overcome by viscous-related forces (higher Oh), strongly opposing liquid infiltration, necking, and ultimately pinching.We therefore highlight the aggregate effect of both Ca and Re on the successful entrapment of liquid drops within the bubble domain, which necessitates an initial bubble rear inversion, a sufficiently fast cavity development, and the incidence of capillary pressure buildup that leads to pinch-off.

Bubble bursting, back deformations, and regime maps
We finalise our discussion by providing a more detailed account of the post pinchoff dynamics, focusing specifically on the process of bubble bursting via the liquid cavity/entrapped drops and the temporal fate of the bubble back and accompanying trailing structures.In figure 11(a)-(b) we summarise the three main bursting behaviours observed in the β s -k and Bi-k (Da) spaces, leaving all other conditions specified in the baseline case unchanged.As captured in the figure, a first regime (I), distinctly bounded by a surfactant-free interface and k below the critical value defined in the previous section, is identified for all β s and Bi tested.This regime adheres to a behaviour mode of slow infiltration of liquid, monotonic delay of the first interior pinch-off event with increasing (decreasing) β s (Bi), and most importantly, eventual restoration of the interfacial morphology of the bubble back, depicted in figure 11(c, top).This restoration ensures a pseudo-stable entrapment of one or multiple drops of varying size and Γ within the bubble domain, only disrupted by a potential rupture of the bubble nose by the first drop.We note that this rupture was not observed for any of the cases encompassed within (I) given the channel length employed in our simulations and the significant elongation of the bubble across the channel.The stability of the encapsulated drops-bubble compound in regime I poses a potential avenue for future research.
A second regime (II) emerges as k approaches the critical value, where a much faster and deeper liquid cavity develops across the bubble and, ensuing from diminished capillary forces, further delays to pinch-off are recorded.These delays lead to the entrapment of the first drop in regions approaching the bubble nose in the axial direction, promptly leading to its rupture, as illustrated at t = 5.32 in figure 11(c, middle).Following a series of end-pinching events, the liquid cavity undergoes several instances of escapes of pinch-off, in conjunction with a surge of surface capillary waves across its domain, and a radial growth that gives rise to a large bulbous end (see t = 8.86).These mechanisms are closely related to the intricate phenomena stemming from retracting liquid filaments, as addressed in Anthony et al. (2019); Constante-Amores et al. (2020).This substantial growth causes the liquid cavity to swiftly burst the interface near the channel walls, thereby fully separating the original bubble into two individual entities.Notably, the fate of the newly-formed small entity located further down the streamwise direction mirrors that of the original bubble as new instances of encapsulation are observed.These bursting mechanisms are wholly in line with the observations of Atasi et al. (2018) for close-to-spherical bubbles under similar conditions of k and Ca (see figure 16 of this reference) and the reports of Izbassarov & Muradoglu (2016) for contraction/expansion channels with a viscoelastic liquid phase.Following this busting, capillary forces will attempt to minimise surface area by closing the openings of the newly-formed bubbles.
The third regime (III) is confined to regions of above critical k with moderate β s and 0.01 < Bi < 1 and closely resembles the mechanisms of II, with the pinch-off times being the key differentiating characteristic.In III, the fast expansion of the liquid cavity expedites the first pinch-off with respect to the surfactant-free case (see previous subsection).Although the location of the first entrapment is much farther away from the bubble nose than in II, the encapsulated drops exhibit a comparatively faster velocity (e.g., u x,d /U b = 1.27 and u x,d /U b = 1.47 at t = 4.87 for the cases in II and III emphasised in figure 11(a), respectively), also leading to bubble nose ruptures by the drops and the consequent radial rupture by the bulbous end of the cavity mentioned above (see t = 10.63).We now briefly draw attention to the multitude of non-axisymmetric structures that unfold in II and III at the trailing ends of the original bubble after the rupturing sequence.Instead of the restoration of the bubble back seen in I, the liquid drag and prevailing adsorption-induced accumulation of surfactant promotes the thinning and elongation of the bubble back in the counter-flow direction, as highlighted in figure 11(c, bottom), whereupon rich interfacial dynamics transpire, including numerous bubble bursting and coalescence instances.These intricate dynamics at the trailing structures of the bubbles are consistent with the numerical results of Nath et al. (2017) (e.g, figure 15 of this reference).

Concluding remarks
This work presents an in-depth characterisation of the highly unsteady mechanisms describing elongated bubbles flowing across liquid-filled channels via numerical simulations.These mechanisms arise in opposition to the well-known steady-state phenomena of Taylor bubbles under the traditional assumptions of negligible inertial/viscous effects, symmetric cross-section channels, and the absence of surface-active agents.As reported in limited instances (see Atasi et al. (2018) and Sauzade & Cubaud (2013)) and confirmed here, interesting dynamics emerge as W e → O(10).Notably, the concomitant reduced capillary effects of high Ca (and W e) induce a loss of sphericity at the bubble's front and back, triggering an elongation in the flow direction of the former and a flattening and subsequent curvature inversion of the latter.The underlying mechanisms by which this curvature inversion evolves into a liquid cavity infiltrating the bubble, how this cavity collapses into small liquid drops entrapped in the bubble domain, and the overall role of surfactants in these processes have been analysed in detail.We emphasise in this study the multiple parallels that can be readily drawn between the various interfacial features, singularities, and intricate interplay amongst the multitude of physical mechanisms that ensue in ours and other well-known systems involving capillary liquid breakup, such as inkjet printing, two-phase microfluidics, and contracting liquid filaments.
By performing a systematic sweep of a number of surfactant parameters and characterising dimensionless groups, we have elucidated the coupled effect of lower surface tension and Marangoni stresses, brought about by the surfactants, and their interactions with inertia and viscosity.Two distinct modes of cavity breakup/drop encapsulation were found for comparatively low rates of surfactant adsorption (low k): one closure of the liquid cavity at the bubble's rear and one end-pinching breakup mode at the interior of the cavity.Our simulations have shown that the combined and individual effects of both surface tension and Marangoni stresses induce a delayed response in terms of pinch-off times, as well as an inversion of the sequence of breakup modes with respect to clean interfaces.Provided that k is low enough, these effects are maintained for increasingly strong (large β s ) and less soluble (small Bi) surfactants.For these cases, complex post pinch-off non-monotonic behaviours were observed in terms of drop size, semi-periodic deformations across the bubble, and velocity relative to the bubble nose, suggesting noticeable effects from Marangoni stresses exerted on the encapsulated drops.Under the baseline conditions considered, a critical value of k = 1 represents a phenomenological limit from which an increase of k represents an acceleration of the second mode of cavity breakup when compared to an uncontaminated interface due to the fast and deep liquid infiltration in the bubble, promoting capillary breakup.A summarising map with three well-defined regimes of behaviour was constructed in the β s -k and Bi-k (Da) spaces, where regime I, limited by the critical k, allows for a semi-stable drop entrapment due to the eventual closure and restoration of the bubble back in the aftermath of pinch-off.Conversely, in regimes II (delayed pinch-off) and III (expedited pinch-off), the deep infiltration of the cavity that precedes end-pinching, together with the high relative velocity of the entrapped drops are responsible for the subsequent bursting of the bubble nose.These two regimes, in addition, display a multitude of rich interfacial events following the cavity end-pinching, including escapes from pinch-off, radial growth, and an eventual bubble bursting in the vicinity of the channel walls.
To the best of our knowledge, this is the first thorough and detailed characterisation of the phenomenon of drop encapsulation in moving bubbles in the presence of surfactants.There are, however, a few avenues worthy of further pursuit in this topic.For instance, an examination of the influence of surfactant parameters, β s for example, on the critical k value.Based on our numerical results, it is likely that this critical k will decrease with increasing β s , but additional evidence is required.A deeper exploration of the stability of the encapsulated drops in regime I may also be of interest to determine if there exists a set of conditions that allow for a perpetual entrapment of these drops, and the potential applications of such systems in the context of emulsification, as remarked previously by Izbassarov & Muradoglu (2016).
Declaration of Interests.The authors report no conflict of interest.

Figure 2 :
Figure 2: Validation for square capillary channels and W e > 1.(a) Bubble nose velocity relative to average liquid phase velocity.(b) Gas phase area fraction.Each marker corresponds to a different Ca, where the colour red represents our numerical results and the remaining colours Magnini & Matar (2020)'s results.(c) Distance from channel centre to any point of interface normalised by the channel's hydraulic radius, R h .(d) Normalised relative pressure, pref = (p − p b )/ρ l U 2 a .Results in (b)-(d) are measured at a cross-sectional plane normal to the stream-wise direction located at a distance of 5.5D behind the bubble nose, as reported in Magnini & Matar (2020).From outermost to innermost series: Ca = 0.0089, Ca = 0.02, Ca = 0.0377, and Ca = 0.1 for Re = 443−500 in (c) and (d) = D : D, O 2 = O : O, and O = (∇u − ∇u T )/2 corresponds to the rate of rotation tensor.Insights into the velocity distribution along the interface, and how this relates to the overall bubble motion across the channel can be gained from the three-dimensional structures displayed in figure4(a).These images show that, although the bubble nose velocity increases in magnitude with Ca (and W e) due to increasingly prominent inertia, a large portion of the interface for Ca = 0.0693 remains notably slower than the bubble nose.In particular, a quasi-uniform spatial distribution corresponding to u ref,x < 0 is observed in the regions adjacent to the thin liquid film.In contrast, Ca = 0.0089 exhibits a sequential pattern of large positive u ref,x at the diagonal crosssection lobes and negative u ref,x at the saddle-like, high-pressure regions (refer to figures 2(c)-(d)).These differences in velocity distribution are attributed to corner flow effects, which drive the cross-sectional interfacial non-uniformity at low Ca, leading to a very thin film (h(x)/D ∼ 3.1 × 10 −3 ) surrounding the saddles with u ref,x < 0 and a thicker

Figure 3 :
Figure 3: Effect of Ca on bubble characteristics for surfactant-free cases at t = 0.89.(a) Full bubble length in the axial direction.(b) Zoom over the trailing undulatory structures at low Ca (3.25 < x < 4.5).From top to bottom the plots correspond to bubble-liquid interface shape, normalised mean curvature, k, normalised normal component of viscous stress on the interface, τn , and normalised tangential component of viscous stress on the interface, τt , respectively.Dotted line for Ca = 0.0693 highlights the re-entrant cavity.All other parameters remain unchanged from those specified in §2.2 for the base case −4 )m > µ c µ d /ρ c σ s ) nor is our local W e characterising the gas phase (W e d = W e c ρ d /ρ c ∼ O(10 −2 ), W e c = ρ c U 2 cavity R cavity /σ s ) large enough for gas inertial effects to be influential (W e d > 0.2) (van Hoeve et al. 2010; Lin & Reitz 1998).Considering this parameter space, we are within the 'dripping' or 'jetting' regimes, as defined in van Hoeve et al. (2010).

Figure 4 :
Figure 4: Effect of Ca on velocity characteristics for surfactant-free cases at t = 0.89.(a) Three-dimensional relative velocity in x-direction, ũref,x = (u x − U b )/U b .U b /U a = 1.0679, 1.2284, 1.3272 for Ca = 0.0089, 0.0377, 0.0693, respectively.(b) Normalised azimuthal vorticity (top) and flow topology parameter (bottom).(c) Relative bubble velocity in x-direction and streamlines drawn in the bubble-tip reference frame.Contour plots on two-dimensional projection in x-y plane (z = 0.5) in (b) and (c).All other parameters remain unchanged from those specified in §2.2 for the base case

Figure 5 :
Figure 5: Temporal evolution of interface characteristics in surfactant-free cases (a) Maximum dimensionless velocity in x-direction (left) and zoom on interfacial singularities found for Ca = 0.0693 (right).(b) Detailed evolution of encapsulation process for Ca = 0.0693.Contour plots of normalised azimuthal vorticity (top) and three-dimensional bubble shape (bottom).Refer to the Supplementary material for an animation of the encapsulation process
(a) depicts, and in line with the studies of Antonopoulou et al. (2021); Zhong et al. (2018) in the context of ink-jet

Figure 8 :
Figure 8: Effect of Bi on cavity formation and encapsulation.(a) Cavity shape, surfactant interfacial concentration (left axis, continuous line) and surface tension (right axis, dashed line), Marangoni stresses, and interface curvature at t = 2.22.(b) Temporal evolution of primary encapsulated drop velocity normalised by bubble velocity, primary drop interfacial area normalised by area after pinch-off, and ratio between drop length in the x and y directions.Data points computed for all cases every 10 −3 s.(c) Primary encapsulated drop shape, surfactant interfacial concentration (left axis, continuous line), surface tension (right axis, dashed line), and Marangoni stresses at t = 7.82.Plots on twodimensional projection in x − y plane (z = 0.5) for (a,c).tp−o,bk = 2.28, 2.48, 2.73, 2.80 and tp−o,int = 2.45, 2.75, 3.02, 3.21 for surfactant-free and Bi = 1, 0.10, 0.01, respectively.All other parameters remain unchanged from those specified in §2.2 for the base case

Figure 9 :
Figure 9: Effect of k and Da on cavity formation and encapsulation.(a) Cavity shape, surfactant interfacial concentration, Marangoni stresses, and interface curvature at t = 2.22.(b) Interior pinch-off time vs. k and Da, and schematic of the three regions that divide thinning behaviour.(c) Contour plots of normalised vorticity in z-direction and comparison of cavity shape for k = 0.1, Da = 1.0 (top) and k = 10.0,Da = 0.01 (bottom).The values of R cavity were calculated at the midpoint between the cavity nose and the start of the cavity.Plots on two-dimensional projection in x-y plane (z = 0.5).tp−o,int = 2.45, 2.64, 2.93, 2.96, 2.05, 2.00, 2.03 for surfactant-free and k = 0.10, 0.25, 1, 2.50, 5.00, 10.00 (Da = 1, 0.40, 0.10, 0.04, 0.02, 0.01), respectively.A few cases of the full range of k have been omitted in (a) for an easier visualisation.All other parameters remain unchanged from those specified in §2.2 for the base case

Figure 11 :
Figure 11: Encapsulation and bursting regime maps in the (a) β s -k and (b) Bi-k spaces.Bi = 0.10 for (a) and β s = 0.50 for (b).The pinch-off times reported correspond to those of the first interior pinch-off event observed.(c) Contour plots of normalised vorticity in z-direction and comparison for the three regimes identified.All other parameters remain unchanged from those specified in §2.2 for the base case

Figure 12 :
Figure 12: Results of mesh independence study for surfactant-free and surfactant-laden cases

Table 1 :
Additional features of the mesh resolutions tested