Palas Kumar Farsoiya

Stéphane Popinet

Howard A Stone

Luc Deike
email: ldeike@princeton.edu

Palas Kumar Farsoiya

A coupled Volume of Fluid -Phase Field method for direct numerical simulation of insoluble surfactant-laden interfacial flows and application to rising bubbles

I. INTRODUCTION A. General context

Surfactants -a contraction of the term surf ace-active agents -are found in a broad spectrum of applications owing to their unique ability to lower surface tension and facilitate interactions between different phases. Found in everyday products, surfactants are essential components in detergents and cleaning agents, enabling emulsification and removal of oils and dirt from surfaces. For example, surfactants are widely used in personal care items and cosmetics as they contribute to foaming and emulsifying properties, while in pharmaceuticals they aid in drug delivery and solubilization. Also, surfactants play key roles in agriculture as adjuvants in pesticide formulations, and in the oil and gas industry, they are crucial for enhanced oil recovery. They find application in paints, coatings, textiles, and even firefighting foams. Additionally, surfactants are integral to biomedical applications, environmental remediation, construction, metal processing, and various industrial processes, illustrating their versatile and indispensable use across diverse sectors [START_REF] Rosen | Surfactants and Interfacial Phenomena[END_REF][START_REF] Myers | Surfactant Science and Technology[END_REF].

The ubiquity of the liquid-air interface in both nature and industry is evident [START_REF] Lohse | Bubble puzzles: from fundamentals to applications[END_REF][START_REF] Lohse | Fundamental fluid dynamics challenges in inkjet printing[END_REF][START_REF] Deike | Mass transfer at the ocean-atmosphere interface: The role of wave breaking, droplets, and bubbles[END_REF]. Its manifestations range from rising air bubbles in water, rain droplets, sloshing in tanks, and breaking waves in the ocean to fuel injection and sprays. Among multiple examples of the influence of surfactants in fluid mechanics problems, recent investigations have shown delayed coalescence of droplets [START_REF] Bruning | Delayed coalescence of surfactant containing sessile droplets[END_REF] of different surface tension, and how the presence of surfactants delayed the break-up of droplets in inkjet printing and prevented the formation of satellite droplets [START_REF] Antonopoulou | Effect of surfactants on jet break-up in drop-ondemand inkjet printing[END_REF]. Separately, bubbles forming on a liquid surface coalesce more at low surfactant concentrations and less at high concentrations [START_REF] Néel | Collective bursting of free-surface bubbles, and the role of surface contamination[END_REF], while also affecting mass transfer at the ocean-atmosphere interface [START_REF] Bell | Dimethylsulfide gas transfer coefficients from algal blooms in the southern ocean[END_REF]. Similarly, spray generation from bursting bubbles shows non-monotonic behavior as the surfactant concentration is increased [START_REF] Néel | Role of contamination in optimal droplet production by collective bubble bursting[END_REF], which affects the generation of jet drops [START_REF] Pierre | Influence of surfactant concentration on drop production by bubble bursting[END_REF] as well as film drops [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF]. While surfactant effects are known to be important at the scale of the capillary length (typically mm to cm), the dynamics induced by Marangoni stresses due to surfactant gradients can affect large-scale phenomena such as waves, as demonstrated by laboratory experiments [START_REF] Liu | The effects of surfactants on spilling breaking waves[END_REF][START_REF] Liu | An experimental study of surfactant effects on spilling breakers[END_REF][START_REF] Liu | Weakly breaking waves in the presence of surfactant micelles[END_REF] on spilling breakers and plunging breakers (Erinin et al. [START_REF] Erinin | The effects of surfactants on plunging breakers[END_REF]). A historical exposition by Lohse [START_REF] Lohse | Surfactants on troubled waters[END_REF] on surfactant effects on waves provides a timeline from the initial critical observation by Franklin et al. [START_REF] Franklin | Of the stilling of waves by means of oil. Extracted from Sundry Letters between Benjamin Franklin[END_REF] of smooth wakes of ships that were releasing oil to Pockels [START_REF] Pockels | Surface tension[END_REF] experiments on measuring surface tension, and highlights the need for improved numerical methods to describe quantitatively the role of surfactants on interfacial fluid dynamics.

B. Surfactant effect on bubbles

Manikantan and Squires [START_REF] Manikantan | Surfactant dynamics: hidden variables controlling fluid flows[END_REF] discuss how surfactants can introduce hidden variables that control fluid flows and have effects that are still not well understood. In their work, they review some classic problems such as bubbles rising in liquids. Interfaces with high surface tension like air/water, liquid metals/air, and aqueous liquids/non-polar liquids are especially susceptible to significant interfacial effects even for trace amounts of surface-active contaminants (contaminants hereafter). In the context of bubbles and drops, the theoretical analyses started by Frumkin [START_REF] Frumkin | On surfactants and interfacial motion[END_REF] and Levich [START_REF] Levich | Physicochemical Hydrodynamics[END_REF] to investigate contaminant effects assume creeping flow and soluble surfactants. Experiments by Savic [START_REF] Savic | Circulation and distortion of liquid drops falling through a viscous medium[END_REF] on surfactant-laden rising bubbles analyzed the case of the insoluble surfactant and observed increased drag. The increased drag is explained through a stagnant cap model, which results from the collection of surfactants at the back of the bubble, where they cause a surface tension gradient, and hence Marangoni stress, tangentially directed towards the top of the bubble. This force opposes the motion of the fluid and creates a zone where the tangential velocity vanishes and an immobile interface is formed, i.e. the stagnant cap. The creeping flow problem for a spherical bubble is formulated in terms of an infinite set of algebraic equations for the coefficients of a series solution. Savic [START_REF] Savic | Circulation and distortion of liquid drops falling through a viscous medium[END_REF] truncated this series after six terms and Davis and Acrivos [START_REF] Davis | The influence of surfactants on the creeping motion of bubbles[END_REF] improved the approximate solution by retaining 150 terms. For a foundational understanding of the motion of surface-active bubbles and drops, Harper [START_REF] Harper | The motion of bubbles and drops through liquids[END_REF] offers an early comprehensive literature survey and theoretical analysis.

Complementary to experimental work, there has been a wide range of numerical studies, with different methods summarized in Table 1, leading to the discussion of various physical effects of surfactant on bubble dynamics. Bel Fdhila and Duineveld [START_REF] Fdhila | The effect of surfactant on the rise of a spherical bubble at high reynolds and peclet numbers[END_REF] conducted experiments and performed a numerical analysis of the surfactant effect on rising bubbles using the stagnant cap model. Atasi et al. [START_REF] Atasi | Influence of soluble surfactants and deformation on the dynamics of centered bubbles in cylindrical microchannels[END_REF] discuss the reduction of rise velocities of spherical bubble and deformation as surfactant effects are increased using Direct numerical simulations (DNS, solving the full Navier-Stokes equations), followed by Atasi et al. [START_REF] Atasi | Lifetime of surface bubbles in surfactant solutions[END_REF] investigating surfactant effect on bubble lifetime at the liquid-air surface. Examples of three-dimensional rising bubble in soluble surfactant solutions are discussed by Muradoglu and Tryggvason [29] while Constante-Amores et al. [START_REF] Constante-Amores | Dynamics of a surfactant-laden bubble bursting through an interface[END_REF] investigated surfactant-laden bursting bubbles among many others. However, a systematic numerical investigation quantifying the effect of insoluble surfactants on a rising unconfined bubble, including three-dimensional configurations with large bubble deformations and unstable paths has not been systematically performed.

C. The need for a versatile numerical framework

The experimental approaches provide a picture of the blended effect on complex interfacial phenomena coupled with the Marangoni effect. Interfacial fluid mechanics problems with surfactants are also dependent on initial conditions that are sometimes difficult to control in experiments. The surfactant concentration on interfaces is difficult to measure and generally experimental data is provided in terms of the bulk concentration of surfactants. The tracking of how much of the surfactants are adsorbed at the interface dynamically is challenging. These challenges are reviewed in Manikantan and Squires [START_REF] Manikantan | Surfactant dynamics: hidden variables controlling fluid flows[END_REF] and highlighted in the recent focus paper by Lohse [START_REF] Lohse | Surfactants on troubled waters[END_REF] and call for the development and release of versatile, open-source numerical frameworks to simulate interfacial flows with surfactants.

Numerical simulations can provide a comprehensive picture of surfactant transport and distribution, and its effect on hydrodynamics, as a result of the ability to control initial conditions and determine all flow field data. Also, numerical simulations provide the opportunity to switch on and off some effects such as adsorption-desorption to switch from a soluble to an insoluble surfactant or include a thermal gradient.

There are many numerical studies and methods considering the Marangoni effects due to thermal gradients, e.g. [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF][START_REF] Seric | Direct numerical simulation of variable surface tension flows using a volume-of-fluid method[END_REF][START_REF] Tripathi | Motion of an air bubble under the action of thermocapillary and buoyancy forces[END_REF][START_REF] Mialhe | An extended model for the direct numerical simulation of droplet evaporation. influence of the marangoni convection on leidenfrost droplet[END_REF]. However, in the presence of surfactant, the Marangoni effects need to be coupled with surfactant transport equations on the interface. Transporting surfactants on a moving interface has been challenging in Eulerian numerical frameworks, and is the focus of the present work.

Recent efforts have investigated how to simulate surfactant-laden flows using different approaches, which we summarize in Table I. The different representations of interfaces provide access to different kinds of physical problems that can be investigated. Stone and Leal [START_REF] Stone | The effects of surfactants on drop deformation and breakup[END_REF] developed a numerical procedure based on the boundary element method to simulate the effect of insoluble surfactant on drop breakup. Further, Li and Pozrikidis [START_REF] Li | The effect of surfactants on drop deformation and on the rheology of dilute emulsions in stokes flow[END_REF] and Yon and Pozrikidis [START_REF] Yon | A finite-volume/boundary-element method for flow past interfaces in the presence of surfactants, with application to shear flow past a viscous drop[END_REF] studied surfactant effects on 3D Stokes flows using the boundary element method. Jan [38] used the Front-Tracking approach to study the effect of surfactants on rising bubbles. The numerical methods for soluble surfactants using front-tracking-based methods were further developed for two-dimensional axisymmetric flows by Muradoglu and Tryggvason [39] and extended to three-dimensional flows [START_REF] Muradoglu | Simulations of soluble surfactants in 3d multiphase flow[END_REF]. Another full three-dimensional front-tracking-based method was developed by De Jesus et al. [START_REF] De Jesus | A 3d front-tracking approach for simulation of a two-phase fluid with insoluble surfactant[END_REF] using the Leibniz Transport Formula for the time derivative of integrals over moving surfaces previously used by Dziuk and Elliott [START_REF] Dziuk | Finite elements on evolving surfaces[END_REF] and Lenz et al. [START_REF] Lenz | A convergent finite volume scheme for diffusion on evolving surfaces[END_REF]. A hybrid approach was used by Ceniceros [START_REF] Ceniceros | The effects of surfactants on the formation and evolution of capillary waves[END_REF] coupling the Front Tracking and Level Set method to study the effect of surfactants on capillary waves. Similarly, Shin et al. [START_REF] Shin | A hybrid interface tracking-level set technique for multiphase flow with soluble surfactant[END_REF] developed a procedure by coupling level-set front-tracking method (adapted from [START_REF] Muradoglu | A front-tracking method for computation of interfacial flows with soluble surfactants[END_REF]) to simulate two-phase flows (in the code BLUE), using a stationary Eulerian structured mesh for velocity and pressure and tracking the interface with an unstructured mesh. Level-set approaches were developed by Xu and Zhao [START_REF] Xu | An Eulerian formulation for solving partial differential equations along a moving interface[END_REF] to solve surfactant transport on a moving interface and were later coupled to the Navier-Stokes equations [START_REF] Xu | A level-set method for interfacial flows with surfactant[END_REF]. Moreover, Piedfert et al. [START_REF] Piedfert | Numerical simulations of a rising drop with shape oscillations in the presence of surfactants[END_REF] and Atasi et al. [START_REF] Atasi | Influence of soluble surfactants and deformation on the dynamics of centered bubbles in cylindrical microchannels[END_REF] used similar approaches for the code DIVA (Dynamics of Interfaces for Vaporization and Atomization) and JADIM respectively to investigate axisymmetric two-dimensional surfactant-laden flows. The diffuse-interface methods also provide an alternative approach to tracking interfaces. Many numerical methods have been developed using the free energy approach by Van der Sman and Van der Graaf [START_REF] Van Der Sman | Diffuse interface model of surfactant adsorption onto flat and droplet interfaces[END_REF], Yun et al. [START_REF] Yun | A new phase-field model for a water-oil-surfactant system[END_REF] and the Cahn-Hillard-based diffused interface framework by Teigen et al. [START_REF] Teigen | A diffuse-interface approach for modeling transport, diffusion and adsorption/desorption of material quantities on a deformable interface[END_REF]. More recently Jain [START_REF] Jain | A model for transport of interface-confined scalars and insoluble surfactants in two-phase flows[END_REF] developed transport of surfactant or interface confined scalars for second-order phase-field methods. The Volume of Fluid (VoF) method for interface tracking provides a mass-conserving approach to solving two-phase Navier-Stokes equations. Renardy et al. [START_REF] Renardy | A new volume-of-fluid formulation for surfactants and simulations of drop deformation under shear at a low viscosity ratio[END_REF] extended this numerical framework for a linear equation of state, i.e. relationship between surface tension and surfactant concentration, using the continuum surface stress method to include surfactant effects. The surfactant transport was achieved by advecting a scalar between two layers of VoF-tracked interfaces. James and Lowengrub [START_REF] James | A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant[END_REF] provided a VoF-based surfactant-induced Marangoni flow for axisymmetric flows which supported an arbitrary equation of state for surface tension. However, the extension of the geometric area advection three-dimensional, adaptive mesh, and parallelization is complex. Within OpenFOAM, the Arbitrary Lagrangian-Eulerian (ALE) interface tracking approach has been used to simulate soluble surfactant-laden flows by Dieter-Kissling et al. [START_REF] Dieter-Kissling | Direct numerical simulation of droplet formation processes under the influence of soluble surfactant mixtures[END_REF] and Lakshmanan and Ehrhard [START_REF] Lakshmanan | Marangoni effects caused by contaminants adsorbed on bubble surfaces[END_REF]. Within the algebraic VoF interface tracking, the authors used an immersed boundary approach for a modified level-set interface representation, similarly to Antritter et al. [START_REF] Antritter | Spreading of micrometer-sized droplets under the influence of insoluble and soluble surfactants: A numerical study[END_REF]. In this context, there is a need for a geometric VoF-based method for three-dimensional two-phase flows, accounting for surfactant transport and Marangoni stresses. Such a method will be well suited for free-surface flows where distinct phases, separated by a sharp interface [START_REF] Popinet | Numerical models of surface tension[END_REF], and would allow investigating problems involving interface reconnection, such as break-up and breaking waves [START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF][START_REF] Mostert | High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production[END_REF]. Moreover, the use of adaptive mesh refinement is essential for resolving multiscale problems in multiphase flows [START_REF] Mostert | High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. Proposing such methods in a fully open-source format is also desirable for the community.

In the present study, we address this gap, by adapting a method developed by Jain [START_REF] Jain | A model for transport of interface-confined scalars and insoluble surfactants in two-phase flows[END_REF] for an interface confined scalar represented by a phase field and incorporate it within the geometric VoF framework Basilisk [START_REF] Popinet | Numerical models of surface tension[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. The present numerical approach uses volume fractions to represent the interface and solve momentum equations. The phase field is only used to transport the surfactant and is computed from a signed distance from the interface [START_REF] Limare | A hybrid level-set/embedded boundary method applied to solidification-melt problems[END_REF]. As a signed distance calculation based on volume fractions is an expensive step, we transport the phase field as developed by [START_REF] Jain | Accurate conservative phase-field method for simulation of two-phase flows[END_REF] and compute the signed distance at a suitable frequency. The Marangoni force is implemented based on continuum surface force modeling proposed by Seric et al. [START_REF] Seric | Direct numerical simulation of variable surface tension flows using a volume-of-fluid method[END_REF] and improved by Tripathi and Sahu [START_REF] Tripathi | Motion of an air bubble under the action of thermocapillary and buoyancy forces[END_REF]. The framework is implemented in Basilisk (Popinet and collaborators [65]), which is a fully open-source code, with well-documented test cases, supports adaptive mesh refinement (AMR), is highly scalable, and can be applied to a wide range of interfacial fluid mechanics problems. Here, we illustrate our methods on the classic problem of the rising bubble in an otherwise quiescent fluid and investigate the role of surfactant and Marangoni stresses. We demonstrate the effect of insoluble surfactant on an axisymmetric and 3D bubble rising in quiescent liquid providing quantitative and qualitative features of this complex problem.

D. Outline

Section II presents the numerical framework, the set of equations to be considered, the implementation of surfactant transport, Marangoni stresses, and several validation cases. Section III illustrates the use of the numerical framework on the problem of a rising bubble in a quiescent fluid in the presence of surfactant. We present an extensive parameter sweep for the axisymmetric configuration, demonstrating the slow down of the bubble due to Marangoni effects, and illustrate the 3D capabilities of the numerical approach by showing how Marangoni effects can induce a transition in the bubble dynamics from ellipsoidal to zigzag motion. Section IV presents conclusions and prospects.

II. NUMERICAL FRAMEWORK A. The Navier-Stokes and surfactant transport equations

We solve the three-dimensional, incompressible, two-phase Navier-Stokes equations using the open-source solver Basilisk [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]65,[START_REF] Popinet | A quadtree-adaptive multigrid solver for the serre-green-naghdi equations[END_REF],

∂u ∂t + ∇ • (uu) = 1 ρ -∇p + ∇ • (µ(∇u + ∇u T )) + γ ρ κnδ s + δ s ∇ s γ (1)
∇ • u = 0, (2) ∂T ∂t + u • ∇T = 0. ( 3
)
where u, p, µ, ρ, γ, κ, and T are, respectively, the velocity, pressure, viscosity, density, surface tension coefficient, curvature, unit normal to the interface, n = ∇T /|∇T |, and volume fraction fields. Also, δ s is the surface delta function (with dimensions 1/length).

The Basilisk solver has been validated extensively for complex interfacial flows including drop and bubble break-up, waves, and bursting bubbles, though without Marangoni forces [START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF][START_REF] Mostert | High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production[END_REF][START_REF] Farsoiya | Axisymmetric viscous interfacial oscillations-theory and simulations[END_REF][START_REF] Farsoiya | Free and Forced, Linearised Interfacial Waves, Theses[END_REF][START_REF] Berny | Role of all jet drops in mass transfer from bursting bubbles[END_REF][START_REF] Farsoiya | Bubble-mediated transfer of dilute gas in turbulence[END_REF][START_REF] Farsoiya | Direct numerical simulations of bubble-mediated gas transfer and dissolution in quiescent and turbulent flows[END_REF][START_REF] Farsoiya | Role of viscosity in turbulent drop break-up[END_REF]. It uses the projection method to compute the velocity and pressure and the geometric VoF method for the evolution of the interface between two immiscible fluids [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]. The Piecewise Linear Interface Calculation (PLIC) geometric interface and flux reconstruction ensures a sharp representation of the interface [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]. The surface tension force is modeled by a continuum surface force model (Brackbill et al. [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]) and uses generalized height functions for curvature by Popinet [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. The surface gradient of interfacial tension in the last term in Eq. ( 1) is modeled as proposed by Seric et al. [START_REF] Seric | Direct numerical simulation of variable surface tension flows using a volume-of-fluid method[END_REF] and the surface delta function has a different implementation than the static surface tension discussed in Popinet [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF].

The time-dependent advective-diffusion equation for surfactants with surface concentration Γ(x, t) along a deforming interface is given by (Scriven [START_REF] Scriven | Dynamics of a fluid interface equation of motion for newtonian surface fluids[END_REF], Stone [START_REF] Stone | A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming interface[END_REF])

∂Γ ∂t + ∇ s • (Γu s ) + Γ(∇ s • n)(u • n) = D s ∇ 2 s Γ (4)
where ∇ s = (I -nn) • ∇ is the surface gradient operator and D s is the surfactant diffusivity along the interface.

Multiphase flows are usually characterized by a set of non-dimensional numbers, together with viscosity and density ratios between the two phases. Flows with with significant viscous and capillary effects can be characterized by the Ohnesorge number, Oh ≡ µ/ √ ρσL, with L a characteristics length (for example the bubble/drop size or the film thickness).

The introduction of surfactant transport and induced stresses leads to additional non-dimensional numbers. The interfacial surfactant transport on the interface can be characterized by the Péclet number Pe s ≡ U L/D s , with U a characteristic velocity. The surfactants alter the surface tension and are generally described by an equation of state or isotherm [START_REF] Manikantan | Surfactant dynamics: hidden variables controlling fluid flows[END_REF]. A linear isotherm relates surface tension as a linear function of surfactant concentration,

γ(x, t) = γ 0 1 -β Γ(x, t) Γ 0 , (5)
with β is a non-dimensional number characterizing the strength of the surface tension variation with a change in surfactant concentration. Other isotherms can be used within the present numerical framework, and have been proposed to describe experimental results on surfactants tested in the laboratory and summarised by Manikantan and Squires [START_REF] Manikantan | Surfactant dynamics: hidden variables controlling fluid flows[END_REF]. The effect of variable surface tension introduces another parameter, often called the Marangoni number.

The Marangoni number will compare the surface tension gradient (leading to Marangoni stresses) and surfactant transport, which can be either surface diffusivity or an advection velocity. Different definitions of the Marangoni numbers have therefore been proposed, depending on the regimes of interest, as discussed in the review from Manikantan and Squires [START_REF] Manikantan | Surfactant dynamics: hidden variables controlling fluid flows[END_REF]. In the present study, we consider moderate to high inertial effects (or Reynolds number), equivalent to low Oh number regimes where surfactant transport is dominated by inertia and P e ≫ 1 so that the Marangoni number is defined through a characteristic velocity, Ma = βγ 0 L/(U µ).

B. Numerical implementation of surfactant transport and Marangoni stresses

We use an approach recently developed by Jain [START_REF] Jain | A model for transport of interface-confined scalars and insoluble surfactants in two-phase flows[END_REF] for interface-confined scalars for diffuse interface/phase-field methods which was not coupled to a flow solver. For the transport of insoluble surfactants having volumetric concentration c, an advection-diffusion equation is solved,

∂c ∂t + ∇ • (uc) = ∇ • D ∇c - 2(0.5 -ϕ)nc ϵ , (6)
where c = Γξ and ξ = ϕ(1 -ϕ)/ϵ. The second term in the parenthesis at the right-hand side of the above equation is an anti-diffusion term that keeps the scalar c concentrated at the interface. The phase field ϕ is computed from the signed distance χ of VoF field T ( [START_REF] Limare | A hybrid level-set/embedded boundary method applied to solidification-melt problems[END_REF]), as indicated in Fig. 1(a,b,c),

ϕ = 1 2 1 -tanh χ 2ϵ , (7)
where ϵ defines the thickness of phase-field ϕ and is kept at 0.75∆x min in the present study. Maintaining the positivity of the scalar concentration field is a requirement that should be consistently met throughout the simulation. Jain and Mani [START_REF] Jain | A computational model for transport of immiscible scalars in two-phase flows[END_REF] derived a criterion sufficient to maintain the positivity of scalar concentration c at all times. The discretized Eq. ( 6) is constrained by the positivity criterion ∆x ≤ 2D s /(|u| max + D s /ϵ) where ∆x, |u| max are the grid size and maximum interfacial velocity respectively. This criterion is however prohibitive even at moderate Péclet numbers. In the present study, we verify that the magnitude of negative values remains at most four orders less than the maximum surfactant concentration for the grid-converged solutions. Moreover, the explicit time discretization constrains the time step size ∆t ≤ (∆x) 2 /2N d D s , where N d is the number of dimensions. 6)). The signed-distance field χ is used to reinitialize the phase field ϕ (see Eq. ( 7)).

As an illustration, the computed fields f (x) in one dimension x, where the interface is at x = 0 between fluid 1 and fluid 2, are shown in Fig. 1. The volume fraction field in an interfacial cell has a value 0 < T < 1 (shown as a black dot), with 0 and 1 assigned to the bulk phases. The field χ is the signed distance from the interface y = x line, shown in blue in Fig. 1, while the phase field ϕ varies as the hyperbolic tangent 0.5(1 -tanh(χ/(2ϵ))) and the concentration field c varies as the hyperbolic secant squared sech 2 (χ/(2ϵ))) [START_REF] Jain | A model for transport of interface-confined scalars and insoluble surfactants in two-phase flows[END_REF]. Computing the signed distance from VoF is an expensive step ( [START_REF] Limare | A hybrid level-set/embedded boundary method applied to solidification-melt problems[END_REF]), especially at high resolutions and for three-dimensional flows. To reduce the computational cost we evaluate Eq. 7 at a suitable frequency and for intermediate steps we use an accurate conservative phase-field/diffuse interface model (ACDI) proposed by [START_REF] Jain | Accurate conservative phase-field method for simulation of two-phase flows[END_REF],

∂ϕ ∂t + ∇ • (uϕ) = ∇ • ζ ϵ∇ϕ - 1 4 1 -tanh 2 ψ 2ϵ ∇ψ |∇ψ| , ( 8
)
where ζ is the velocity scale parameter and should be greater than |u| max . ψ is an auxiliary signed-distance-like variable ψ = ϵ log[(ϕ + ε)/(1 -ϕ + ε)] where ε = 10 -6 . The use of the phase field simplifies surfactant transport, which involves the calculation of the surface gradient ∇ s Γ (see Eq. 4) [START_REF] Jain | A model for transport of interface-confined scalars and insoluble surfactants in two-phase flows[END_REF].

C. Validation of surfactant transport and Marangoni forces

Here we report the validation of the numerical framework for the transport of the surfactant at the interface coupled with Marangoni forces (Eqs. 1 to 4). First, the surfactant transport Eq. ( 4) is tested below for interface expansion, advection, and diffusion. Surfactant transport is restricted to the interface and quantitatively agrees with the analytical solutions. Second, the Marangoni force calculation is tested against a widely used analytical solution by [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF]. A test for topological changes is also provided. To ensure full reproducibility of the results, each validation case is accompanied by the corresponding documented source code on the Basilisk website [START_REF] Farsoiya | Expanding Circle Test[END_REF].

Expansion test

Both circle(2D) and sphere(3D) of diameter d b are initialized at the center of the computational domain of size L. The third term on the left-hand side in Eq. ( 4), Γ(∇ s • n)(u • n) accounts for the transport of surfactants due to the expansion of the interface [START_REF] Stone | A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming interface[END_REF]. For a given velocity field u(r, θ, ϕ) = Krê r , Γ = Γ/Γ 0 , ũ = u/(Kd b ), t = Kt, P e = Kd 2 b /D s Eq. ( 4) reduces to , where m(t) = V c(t) dv, (i) Second-order accuracy as the resolution is increased (number of grid points per diameter). These results can be reproduced using the code in [START_REF] Farsoiya | Expanding Circle Test[END_REF].

∂ Γ ∂ t + Γ( ∇s • n)(ũ • n) = 1 P e ∇2 Γ. (9)
For a uniform surfactant concentration ∇s Γ = 0 and P e > 0, Eq. ( 9) reduces to

∂ Γ ∂ t + Γ( ∇s • n)(ũ • n) = 0, (10)
where ( ∇s • n) and (ũ • n) are, respectively, the total curvature and the normal velocity at the interface. Integrating Eq. [START_REF] Néel | Role of contamination in optimal droplet production by collective bubble bursting[END_REF] gives the analytical expression Γ( t) = e -t and Γ( t) = e -2 t for circle(2D) and sphere(3D) respectively. where m(t) = V c(t) dv. (f) First-order accuracy of Γ as the resolution is increased. These results can be reproduced using the code in [START_REF] Farsoiya | Rotating Circle Test[END_REF].

We performed the numerical test for an initial diameter d b /L = 0.4, initial uniform concentration Γ(0) = 1, P e = 16 and for three uniform resolutions, i.e., ∆x = d b /12.8, d b /25.6, d b /51.2. Fig. 2(a,d),(b,e),(c,f) show the initial ( t = 0), intermediate ( t = 0.35), and final ( t = 0.7) surfactant concentration on the expanding circle and sphere for the highest resolution. Fig. 2(g) shows the time evolution of surfactant concentration for three different resolutions. With more than twenty grid points per diameter good agreement with the analytical solution (red curve) is obtained. Fig. 2(h) shows the evolution with time of the relative error on surfactant mass: Ẽ = (m(t) -(m(0)))/m(0) < 10 -6 where m(t) = V c(t) dv. Second-order accuracy can be seen in Fig. 2(i) for the error computed relative to the analytical solution with respect to the number of grid cells per unit diameter d b /∆x.

Advection-Diffusion test

The advection term in Eq. ( 4) is constrained by the cell Péclet number [START_REF] Jain | A model for transport of interface-confined scalars and insoluble surfactants in two-phase flows[END_REF]. Hence, advection and diffusion both are validated in the following single test. We The squares show the mean velocities of axisymmetric and 3D droplets respectively, converging towards the theoretical prediction [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF]. These results can be reproduced using the code in [START_REF] Farsoiya | Migrating Bubble Test[END_REF].

Droplet migration in a thermal gradient

We test the Marangoni force calculation by considering a droplet (density and viscosity ratio being unity) migrating due to a thermal gradient. The test has been used extensively in earlier studies by [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF][START_REF] Seric | Direct numerical simulation of variable surface tension flows using a volume-of-fluid method[END_REF][START_REF] Tripathi | Motion of an air bubble under the action of thermocapillary and buoyancy forces[END_REF][START_REF] Herrmann | Thermocapillary motion of deformable drops and bubbles[END_REF] to validate the Marangoni force calculation but does not consider surfactant transport at the interface. The temperature gradient causes a surface tension variation at the interface of the drop and induces thermo-capillary motion. The surface tension coefficient is considered a function of temperature,

γ = γ 0 + γ T (T -T 0 ). ( 12
)
An expression for the migration velocity of a spherical droplet in a constant thermal gradient for small Re ≡ ρ bulk U d b /µ bulk and Ca = µ bulk U/γ 0 is given by [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF]

U drop = - γ T d b ∇T µ bulk (6 + 9µ drop /µ bulk ) (13)
The numerical test is performed for a droplet of diameter d b , domain L = 7.5d b , Re = 0.132, Ca = 0.066, viscosity and density ratio unity. Fig. 4(a) and (b) show the axisymmetric two-dimensional and three-dimensional droplets migrating in a thermal gradient respectively. The droplet migrates towards warmer regions due to the non-equilibrium force created by variable surface tension. The vectors show the direction of the local velocity. Fig. 4(c) shows the time evolution of migration velocity center of mass of the droplet. At the lowest grid resolution of 16 cells per diameter of the droplet, the migration velocity oscillates around the value predicted by Eq. 13. Increasing the resolution from 32 to 64 cells per diameter of the droplet provides good agreement with the analytical expression. Figure 5 shows the surfactant concentration and the interface at different times. multiple topological changes (break-up and coalescence, the surfactant remains on the interface and its total mass is conserved. These results can be reproduced using the code in [START_REF] Farsoiya | Topological Test[END_REF].

the surfactant transport involving topological changes. The cylindrical fluid thread predictably undergoes Plateau-Rayleigh instability as kd b /2 < 1. Figures 5(a-f) show the interface and surfactant concentration at different times non-dimensionalized by ω c = γk 3 /ρ b . The mass of the surfactant is calculated on the interfaces of the initial thread and the individual droplets m = c(t) dv/ c(0) dv and shown in Figure 5(g). In Figure 5(b) tω c = 0.94, a thinning of the thread occurs before the breakup where the surfactant is moved away from the central part of the thread. The breakup is a topological change where the thread interface is broken into three parts and shown in Figure 5(c) at tω c = 1.1. The surfactant concentration remains restricted to the interface without any singular values, and mass is conserved. Figure 5(g) shows that the mass gets from the thread (red squares) to be distributed during the breakup between the top (red circle) and bottom (black circle) drops as well as in a satellite drop (blue diamonds). In Figure 5(d) tω c = 1.57, the subsequent retraction of the central droplet further breaks into three droplets (Black triangles in Figure 5(g)) before coalescing in Figure 5(e) at tω c = 1.6, with redistribution of surfactant well captured in both cases.

Advantages and Limitations of numerical framework for the surfactant transport

The numerical framework is thus validated for the transport of surfactants, the Marangoni force calculation, and the topological changes. The method is at least first-order accurate and a simple implementation avoids the calculation of surface derivatives on a moving fluid-fluid interface for surfactant transport. Therefore, it is adaptable to various other numerical frameworks and efficient in investigating two and three-dimensional two-phase flow problems. The surfactant transport is achieved by an anti-diffusion term and resolution requirements become prohibitive for high Péclet numbers. In the next section, we investigate a bubble laden with insoluble surfactant.

III. EFFECT OF SURFACTANT ON RISING BUBBLES A. Problem statement

Experimental data summarized by Clift et al. [START_REF] Clift | Bubbles, drops, and particles[END_REF] from Grace [START_REF] Grace | Shapes and velocities of single drops and bubbles moving freely through immisicible liquids[END_REF] provide an empirical prediction of a bubble's terminal velocity and terminal shape given a set of two independent parameters, Ga = ρ l d b √ gd b /µ l and Bo = ρ l gd 2 b /γ 0 , where d b is the diameter of the undeformed bubble, compares viscous and gravity and gravity and surface tension effects respectively. The bubble's terminal velocity can then be expressed as a non-dimensional number, usually either expressed as a Reynolds number or a drag coefficient, which is a function of the other two non-dimensional numbers Bo and Ga.

For small Bo and Ga, slow-moving spherical bubbles (creeping flow), Hadamard [START_REF] Hadamard | Slow permanent movement of a liquid and viscous sphere in a viscous liquid[END_REF] and Rybczynski [START_REF] Rybczynski | On the progressive motion of a liquid sphere in a viscous medium[END_REF] independently found the solution for the terminal velocity. However, the Hadamard-Rybczynski theory of internal circulation and terminal velocity is not commonly observed. The explanation provided by Frumkin [START_REF] Frumkin | On surfactants and interfacial motion[END_REF] and Levich [START_REF] Levich | Physicochemical Hydrodynamics[END_REF] for the absence of internal circulation is due to surfactants, which tend to accumulate at the interface reducing the surface tension and result in a gradient that causes tangential stress and so retards the surface motion. Direct Numerical Simulations (DNS) of a weakly soluble surfactant-laden axisymmetric bubble were carried out by Cuenot et al. [START_REF] Cuenot | The effects of slightly soluble surfactants on the flow around a spherical bubble[END_REF] confirming the stagnant cap model and drag reduction, and consistent with the dependence on stagnant cap angle similar to the theory by Sadhal and Johnson [START_REF] Sadhal | Stokes flow past bubbles and drops partially coated with thin films. part 1. stagnant cap of surfactant film-exact solution[END_REF]. Later, the effect of adsorption/desorption of soluble surfactant was investigated by Lakshmanan and Ehrhard [START_REF] Lakshmanan | Marangoni effects caused by contaminants adsorbed on bubble surfaces[END_REF] for axisymmetric bubbles.

Apart from terminal velocity and terminal shape, a phenomenon that has been an active area of research is the transition from a straight to non-straight (either helicoidal or zigzag) path of a rising bubble (traced to the original observations of Leonardo da Vinci by Krishna and Baten [START_REF] Krishna | Simulating the motion of gas bubbles in a liquid[END_REF]; see also the discussion of "Leonardo's paradox" by Prosperetti [START_REF] Prosperetti | Bubbles[END_REF]). Studies performing stability analysis on fixed spheroids Mougin and Magnaudet [START_REF] Mougin | Path instability of a rising bubble[END_REF] or oblate spheroids Tchoufag et al. [START_REF] Tchoufag | Linear instability of the path of a freely rising spheroidal bubble[END_REF] demonstrated the existence of a wake instability behind the particle causing a helical or zigzag path. Separately, stability analysis of deformable spheroids by Zhou and Dušek [START_REF] Zhou | Marginal stability curve of a deformable bubble[END_REF], Bonnefis [START_REF] Bonnefis | Study of Wake, Shape and Trajectory Instabilities of Bubbles Using a Global Linear Stability Approach[END_REF], Herrada and Eggers [START_REF] Herrada | Path instability of an air bubble rising in water[END_REF] and Bonnefis et al. [START_REF] Bonnefis | Path instability of deformable bubbles rising in newtonian liquids: A linear study[END_REF] revealed that the instability due to the interplay of the bubble deformation and the bubble flow can also be responsible for the bubble path transitioning from a rectilinear to oscillatory trajectory. Such instability when the fixed point becomes unstable with a periodic orbit is known as the Hopf bifurcation. Direct numerical simulations by Cano-Lozano et al. [START_REF] Cano-Lozano | Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability[END_REF] have found regions of parameters, especially for small Bo and high Ga, where the helical or zigzag path is neither explained by the wake instability nor by the bubble deformation. These observations and analyses point to the bubbles showing a behavior somewhere between solid and fluid.

In the present study, we explore the effect of surfactants on a rising bubble in both axisymmetric and threedimensional configurations. Figure 6(a) shows the qualitative features of a rising bubble laden with insoluble surfactant. The surface-restricted transport concentrates the surfactants at the rear of the bubble and creates a gradient of surface tension resulting in Marangoni forces in the direction of this gradient ∇ s γ, slowing down the bubble.

The effect of surfactants transported along the interface therefore introduces a non-dimensional parameter, the Marangoni number, comparing surface tension gradient effects and surfactant advection. Here we consider the characteristic velocity √ gd b and define Ma = βγ 0 / √ gd b µ l , where β is the Gibbs elasticity. The Gibbs elasticity is related to the surfactant isotherm provided in Eq. ( 5). Next, using our numerical framework we explore the role of Ma on the bubble rise velocity and demonstrate that several classical results can be retrieved and quantified.

B. Approach: Direct Numerical Simulations

We perform 3D and axisymmetric DNS of a single buoyantly rising bubble (subscript b) of undeformed diameter d b in a denser quiescent liquid (subscript l). The domain uses periodic boundary conditions and is of size L = 20d b . The computational domain is resolved using adaptive mesh refinement (as illustrated in figure 6(b)) by having 25-100 grid cells per diameter of the bubble, with the high-resolution grid being focused on the bubble interface and the boundary layer around it. The chosen adaptivity criteria allow us to properly resolve the influence of the bubble's wake on its dynamics while making sure that given the periodic boundary conditions, the bubble is not influenced by the previous wake. shows a snapshot from one simulation, where the bubble interface is shown by the solid black curve, and illustrates the adaptive grid and the formation of the stagnation cap, where the surfactant (in color, with red regions corresponding to high surfactant concentration and blue regions of low concentration) is concentrated at the bottom of the bubble. The adaptive mesh refinement (AMR) around the bubble shown for L11 ∆x = L/2 11 keeps computational cost low as our simulations at L11 involve ≈ 1.4 × 10 4 number of grid points in 2D(axisymmetric) and ≈ 4 × 10 6 in 3D. Comparing such resolution on the bubble to a uniformly-refined domain, the number of cells would be (2 11 ) 2 ≈ 4.2 × 10 6 and (2 11 ) 3 ≈ 8.6 × 10 9 for 2D (axisymmetric) and 3D respectively, demonstrating the interest of AMR for such flows. Grid convergence is discussed in Appendix A.

We consider the case where surfactant transport is convection dominated so that the Péclet number is large; here P e ≡ d b √ gd b /D s = 100 (using √ gd b as the characteristic advection velocity). We consider large density and viscosity ratios ρ l /ρ b = 1000 and µ l /µ b = 100.

We limit the axisymmetric simulations in Table II to the region where the rising bubble remains spheroidal and the path is rectilinear, as observed by Cano-Lozano et al. [START_REF] Cano-Lozano | Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability[END_REF]. Three-dimensional simulations consider cases showing helicoidal and zigzag trajectories as discussed in the numerical work from Cano-Lozano et al. [START_REF] Cano-Lozano | Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability[END_REF] as well as in the experiments from Tagawa et al. [START_REF] Tagawa | Surfactant effect on path instability of a rising bubble[END_REF].

As discussed above, we consider a linear isotherm, Eq. ( 5). Experimental observations have shown that the minimum surface tension achieved is surfactant-specific [START_REF] Giribabu | Adsorption of nonionic surfactants at fluid-fluid interfaces: Importance in the coalescence of bubbles and drops[END_REF]. The maximum packing concentration of surfactants limits a further decrease in the surface tension coefficient γ ∞ /γ 0 = 0.2 -0.6 [START_REF] Erinin | The effects of surfactants on plunging breakers[END_REF][START_REF] Kosior | Initial degree of detaching bubble adsorption coverage and the kinetics of dynamic adsorption layer formation[END_REF]. In the present study, the transport of surfactants is restricted to the interface and no desorption and adsorption has been considered. An approach of considering a minimum threshold value of surface tension has been used in independent numerical studies by [START_REF] Atasi | Lifetime of surface bubbles in surfactant solutions[END_REF][START_REF] Atasi | Presence of surfactants controls the stability of bubble chains in carbonated drinks[END_REF] to avoid negative values. In the present study, we found that the sharp threshold value produces a surface tension gradient which is always grid dependent. We choose a smooth isotherm γ(x, t) = γ 0 (1 -tanh[βΓ(x, t)/Γ 0 ]) that behaves as Eq. ( 5) for βΓ(x, t)/Γ 0 → 0 and approaches zero as βΓ(x, t)/Γ 0 → ∞. We have validated that this isotherm behaves as Eq. ( 5), remains grid independent, and does not affect the bubble rise velocity and trajectory. We now present a systematic study of bubbles rising in clean and contaminated conditions, i.e., surfactants are present at the interface, for a wide range of non-dimensional parameters, as summarized in Table 2, varying Ga from 5 to 100 and keeping Bo= 0.5 while Ma varies from 0 (clean) to 20. As discussed in the literature, these conditions correspond to a straight path and the axisymmetric geometry is appropriate.

Figure 7 shows the terminal state from the axisymmetric simulations of the rising bubble at various Ma for Ga = 100 and Bo = 0.5 and illustrates the effect of surfactant and Marangoni stresses on rising bubbles; the effect of increasing Ma number is evident. The streamlines and vorticity field are plotted from the reference frame of the rising bubble and illustrate the flow structure, while the surfactant concentration is also shown. In particular, the clean bubble corresponding to Ma = 0 is shown in Figure 7(a). We observe that the high vorticity region is inside the bubble due to the internal circulation of lower-viscosity gas driven by the liquid outside and streamlines show that the flow is, qualitatively, similar to a potential flow, outside the bubble. The other snapshots (fig. 7b,c,d) show, for increasing Marangoni numbers, the flow features when the terminal velocity is reached once surfactants have been introduced. Figure 7(b) shows Ma = 0.2, with surfactant concentrated at the rear of the bubble and a stagnant cap formation is visible at the rear of the bubble due to the surface immobilization, itself a consequence of the Marangoni force acting in the direction opposite to the surfactant flow. The effect of surfactants is therefore similar to a no-slip condition, which arises at the stagnant cap and a flow separation can be seen from the ends of the stagnant cap. The streamlines show the first indications of wake formation behind the bubble. As the Marangoni number increases from 0.2 to 2 (figure 7(c,d)), the stagnant cap and wake expand, and there is a change in the bubble shape, which gets more spherical as the bubble. In the figure panels, the streamlines are colored with the magnitude of velocity, which confirms the reduction in internal circulation as the surface is immobilized due to the Marangoni force.

Next, we report in Figure 8 the interfacial velocity in the reference frame of the bubble, for the same cases as in figure 7, as a function of the polar angle θ measured from the nose of the bubble (as defined in figure 6(a)). The velocity is decomposed into a normal component and a tangential component (the normal component being 0 in the frame of reference of the bubble). In the case without surfactant (Fig. 8a,Ma=0), the interfacial velocity tangential component increases and reaches a maximum around a polar angle θ = π/2 and then decreases to zero at the bottom point of the bubble which is the stagnation point for a clean bubble. These features are qualitatively consistent with the potential flow solution. The stagnant cap angle for this case is thus θ 0 = π.

As Ma increases, e.g., Ma= 0.2 in Figure 8(b), the tangential component of the bubble velocity goes to zero before the bottom point of the bubble, so that the stagnant cap angle is 3π/4 < θ 0 < π, and reaches π/2 < θ 0 < 3π/4 for Ma = 2 (Fig. 7(c) and Fig. 8(c)) and θ 0 ≈ 0 for Ma = 20 (Fig. 7(d) and Fig. 8(d)) where the stagnant cap covers a larger surface area of the bubble and, correspondingly, the surface is immobilized. These effects produce stronger wakes behind the bubble and increase the drag and hence decrease the rise velocity. Also, the shape of the bubble is changed from oblate to spherical as Ma is increased. The stagnant cap angle θ 0 corresponds to the location where the tangential velocity u s reaches zero. For θ > θ 0 the surfactant accumulates, as indicated by the red dashed line in Figure 8(b-d), which shows the surfactant concentration as a function of θ. The presence of surfactant for θ < θ 0 is due to the competition between the diffusive transport and advective transport, i.e., D s |∇Γ| > u s Γ or P e s = u s d b /D s < d b |∇Γ| /Γ. As P e s → ∞, the surfactant concentration at the stagnant cap angle Γ(θ 0 ) → 0.

The clean bubble case Ma = 0 reaches a terminal velocity u c b (∞), which is reduced when surfactants are introduced and as the strength of the Marangoni stresses is increased (increasing Ma). Figure 9 shows the time evolution of the rise velocity normalized by the terminal velocity of the clean bubble, u b (t)/(u c b (∞)), for (a) Ga = 5 and (b) Ga = 100 with Bo = 0.5 and increasing Marangoni numbers. Similar trends are observed for all Ga and Bo tested, and results are grid converged as demonstrated in Appendix A. Increasing Ma leads to a decrease in the rise velocity at all times, and the transient regime is longer when Ma is increased. The transient regime is affected by the transport of the surfactant on the interface of the bubble, and in contaminated cases, the terminal velocity is only achieved once the surfactant concentration along the bubble has reached a steady state (e.g., see Figure 7). The terminal rise velocity for both values of Ga decreases as Ma is increased, as seen in Figures 9a and9b. A larger Ma is required for Ga = 100 to reach a similar reduction in velocity, as a higher Marangoni force is required to counteract the force of buoyancy. At sufficiently high Ma, the reduction saturates, with a similar reduction in velocity for Ma=10 and 20 for all Ga. The terminal velocity, normalized by the terminal velocity of a clean bubble u b (∞)/u c b (∞), with increasing Ma for a range of Ga = 5 -100 (and Bo = 0.5) is reported in Figure 9(c). The sensitivity to the Marangoni number on the terminal velocity decreases as the Galileo number is increased, in the sense that a higher Ma is necessary to start observing a decrease in rise velocity. For example, we observe a reduction of 20% in rise velocity at Ma=0.1 for Ga=5, while Ma=1 is necessary at Ga=50. However, in all cases, a decrease in rise velocity is observed for Ma > 1, and the reduction in rise velocity plateaus for Ma > 5. The rise velocity reduction for Ma ≫ 1 increases slightly from about 50% to 60% when increasing Ga. The transition from a clean to fully contaminated terminal rise velocity occurs over a very sharp gradient in Ma number at high Ga while a smooth transition is observed at low Ga. The behavior described in Figure 9c when increasing the Marangoni number is in good agreement with the observations reported in experiments by Bel Fdhila and Duineveld [START_REF] Fdhila | The effect of surfactant on the rise of a spherical bubble at high reynolds and peclet numbers[END_REF], where a plot similar to ours is presented with the concentration of surfactant used for the x-axis instead of the Ma number. A quantitative comparison of the effect of surfactants on the rise velocity with experimental data is challenging because we would need to be able to link the bulk concentration reported experimentally to the bubble surface concentration (which is challenging to measure), while also requiring experimental data on the specific isotherm from the experiment (the isotherm being surfactant-specific; see the discussion in Manikantan and Squires [START_REF] Manikantan | Surfactant dynamics: hidden variables controlling fluid flows[END_REF]). Indeed, the details of the curve shown in figure 9c are expected to be sensitive to the details of the surface isotherm, which vary with the type of surfactant. However, the general behavior of retardation of rise velocity is consistent.

Next, the normalized drag coefficient is defined as CD = (C D -C rigid )/(C rigid -C clean ), where C rigid and C clean are the drag coefficients for the bubble having stagnant cap angle θ 0 = 0 and θ 0 = π/2 respectively (clean and fully contaminated cases) and the standard drag coefficient is calculated as

C D = 4gd b /(3(u b (∞)) 2
). The corresponding results for Ga = 5 and Ga = 100 are shown in Figure 9(d). For high Ma, the stagnant cap angle is close to zero, the bubble has a completely immobilized interface, and the drag coefficient corresponds to a rigid sphere. Similarly, for low Ma, the interface is mobile and the drag coefficient corresponds to that for clean bubbles. The inset in Figure 9(d shows the variation of the stagnant cap angle θ 0 as the Marangoni number increases for Ga = 5, 100, which effectively providing a way to go from Marangoni number to the stagnant cap (for different Ga numbers) and finally to the drag coefficient.

Even though the equivalent Reynolds number is not necessarily low, the dependence of the drag coefficient as a function of the stagnant cap angle follows a trend predicted for creeping flows by Sadhal and Johnson [START_REF] Sadhal | Stokes flow past bubbles and drops partially coated with thin films. part 1. stagnant cap of surfactant film-exact solution[END_REF]. Other independent numerical studies by Cuenot et al. [START_REF] Cuenot | The effects of slightly soluble surfactants on the flow around a spherical bubble[END_REF] and Kentheswaran et al. [START_REF] Kentheswaran | Direct numerical simulation of gas-liquid mass transfer around a spherical contaminated bubble in the stagnant-cap regime[END_REF] have found similar trends for drag coefficient as a function of stagnant cap angle θ 0 . This suggests that a simple parameterization for the drag coefficient (or reduction in rise speed) could be proposed as a function of Ma.

D. Three-dimensional motion of bubbles: spiral and zigzag trajectories

Three-dimensional simulations of surfactant-laden bubbles are necessary to explore the transitions from rectilinear (straight path) to complex paths (helicoidal or zigzag motion). The complex physics of the non-straights paths of rising bubbles is an active area of research as discussed earlier [START_REF] Bonnefis | Study of Wake, Shape and Trajectory Instabilities of Bubbles Using a Global Linear Stability Approach[END_REF][START_REF] Herrada | Path instability of an air bubble rising in water[END_REF][START_REF] Bonnefis | Path instability of deformable bubbles rising in newtonian liquids: A linear study[END_REF]. Investigations using numerical simulations provide control over the effect of surfactants complementing experiments [START_REF] Tagawa | Surfactant effect on path instability of a rising bubble[END_REF] and could help better understand the various regimes between straight, zigzag, and helical paths, as well as the mechanisms behind each dynamical regime. Here we consider three-dimensional simulations for a set of parameters from Cano-Lozano et al. [START_REF] Cano-Lozano | Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability[END_REF] where they have observed path instability of a rising bubble in clean conditions. We validated the trajectories and the rise velocities of the bubble with their simulations, in terms of terminal rise velocity and oscillation frequency.

We consider a case from Cano-Lozano et al.

[101] Ga = 100.25, Bo = 10 and displays a spiral trajectory as shown in Fig. 10(a) in clean conditions Ma = 0, as shown by the black curve. When surfactant effects are introduced, here increasing the Marangoni number to Ma = 1, we observe a transition in trajectory from spiral to zigzag. Such transition from spiral to zigzag is also observed by Tagawa et al. [START_REF] Tagawa | Surfactant effect on path instability of a rising bubble[END_REF] (their experiment being conducted for Ga = 280 and Bo = 0.5) as the surfactant concentration is increased in the experiments. These comparisons suggest that surfactant effects at relatively high Ga and various Bo introduce a change in the non-linear path, from spiral to zigzag trajectories. That behavior might be relatively generic since it was observed with a linear isotherm in our numerical configuration while the surfactant used in Tagawa et al. [START_REF] Tagawa | Surfactant effect on path instability of a rising bubble[END_REF] are 1-Pentenol and Triton X-100.

The qualitative comparison is encouraging as it concurs with the experimental results of Tagawa et al. [START_REF] Tagawa | Surfactant effect on path instability of a rising bubble[END_REF] who explored the effect of surfactants on bubble trajectories. The results call for a systematic three-dimensional study of the role of the Marangoni number on bubble trajectories in the non-linear regime, as well as a discussion on the role of the specific surfactant isotherm.

The streamlines calculated from the reference frame of the bubble for Ga = 100.25, Bo = 10 are shown in Fig. 10(b,c). The curved streamlines around the clean bubble in Fig. 10(b) are due to the periodic flapping deformation of the bubble which lets the liquid move towards the other side of the bubble. The periodic deformation also validates the stability analysis on rising deformable bubbles by Bonnefis [START_REF] Bonnefis | Study of Wake, Shape and Trajectory Instabilities of Bubbles Using a Global Linear Stability Approach[END_REF], Herrada and Eggers [START_REF] Herrada | Path instability of an air bubble rising in water[END_REF], Bonnefis et al. [START_REF] Bonnefis | Path instability of deformable bubbles rising in newtonian liquids: A linear study[END_REF]. The surface immobilization due to the presence of surfactants, causes a wake formation behind the bubble possibly resulting in a zigzag trajectory.

IV. CONCLUSION

The present study provides a numerical framework to simulate insoluble surfactant-laden 2D and 3D interfacial flows. We provide a hybrid method in which the momentum and surfactant transport equations, respectively, use the VoF and phase-field representation of the interface. The surfactant transport and Marangoni force calculation are tested against standard cases. The method conserves the mass of the surfactant and can be used with arbitrary interfacial isotherms or equations of state for the variation of surface tension. We also demonstrate that the method can handle surfactant transport with topological changes produced by the Plateau-Rayleigh instability. The implementation supports parallel computing and adaptive mesh refinement, is fully open-source, and scripts are provided to reproduce the test cases and results. The surfactant transport equation in the present framework requires high resolution as the surface Péclet number increases and cannot be applied to strictly advective transport. Despite these limitations, the method can be applied to study various interfacial flows, such as drop and bubble dynamics, waves, or thin film motions.

In the present work, we apply the numerical framework for the case of surfactant-laden rising bubbles, considering an idealized linear isotherm. The strength of the surfactant effects is then determined by the Marangoni number, and the behavior when increasing the Marangoni number is akin to increasing the concentration of surfactants in experiments. We retrieve the qualitative and quantitative features of the effect of surfactants previously described in the literature, theoretically and experimentally. In the axisymmetric simulations, as the Marangoni number is increased, we observe the formation of a stagnant cap and immobilization of the interface, as well as the change of the bubble geometry to a more spherical shape at high Marangoni number. As a result, we observe and characterize the decrease of the terminal velocity as a function of Marangoni number. The slowing down of the rise velocity can be summarized in terms of the variation of the drag coefficient as a function of the stagnant cap angle and good agreement is observed when compared to the theoretical expression derived by Sadhal and Johnson [START_REF] Sadhal | Stokes flow past bubbles and drops partially coated with thin films. part 1. stagnant cap of surfactant film-exact solution[END_REF] for creeping flow, even though our simulations are performed at relatively high Reynolds or Galileo numbers. This suggests that by parameterizing the dependence of the stagnant cap angle on the Marangoni number, a prediction for the drag coefficient of the bubble in the presence of surfactant for a wide range of conditions should be possible.

We also perform three-dimensional simulations of rising bubbles in the presence of surfactant. We consider conditions studied by Cano-Lozano et al. [START_REF] Cano-Lozano | Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability[END_REF] of a bubble, with a clean interface, which follows a helicoidal trajectory. In the presence of surfactant, i.e., for a finite Marangoni number, we obtain a change in the bubble trajectory from helicoidal to zigzag. Such observations are in agreement with the experimental study by Tagawa et al. [START_REF] Tagawa | Surfactant effect on path instability of a rising bubble[END_REF] for a similar range of parameters. Further analysis of surfactants affecting the path instability is of great interest as discussed by Herrada and Eggers [START_REF] Herrada | Path instability of an air bubble rising in water[END_REF] and we believe the present open-source numerical framework will be beneficial to the community. Finally, experimental and numerical studies should also be performed, with measurements of the specific surfactant isotherm, to quantitatively link the Marangoni number and the experimental conditions, as discussed and suggested in Erinin et al. [START_REF] Erinin | The effects of surfactants on plunging breakers[END_REF] and Lohse [START_REF] Lohse | Surfactants on troubled waters[END_REF].

FIG. 1 .

 1 FIG. 1. Schematic representing the VoF field T , the signed-distance field χ, the phase-field ϕ, and the volumetric surfactant concentration c. The momentum equations use the volume fraction T to update the fluid properties. The surfactant concentration c transport is achieved by solving for the phase-field ϕ (see Eq. (6)). The signed-distance field χ is used to reinitialize the phase field ϕ (see Eq. (7)).

FIG. 2 .

 2 FIG.2. Normalized surfactant concentration Γ on an expanding circle(2D) and sphere(3D) (sliced to show the cross-section) having initial diameter d b /L = 0.4 (a,d) at times t = 0, (b,e) t = 0.35, and (c,f) t = 0.7. (g) Concentration time evolution in 2D and 3D, with a comparison of the numerics (symbols) and theory (solid and dashed). (h) Relative error Ẽ = (m(t) -(m(0)))/m(0) of surfactants mass conservation remains below 10 -6 , where m(t) = V c(t) dv, (i) Second-order accuracy as the resolution is increased (number of grid points per diameter). These results can be reproduced using the code in[START_REF] Farsoiya | Expanding Circle Test[END_REF].

FIG. 3 .

 3 FIG. 3. Surfactant advection test for a rotating circle. Normalized surfactants concentration Γ having initial diameter d b /L = 0.4 at (a) tω/(2π) = 0 and (b) tω/(2π) = 1.6 (c) tω/(2π) = 3.2. (d) Concentration time evolution at the zenith of the circle (θ = π/2) compared to the analytical solution. (e) The relative error Ẽ = (m(0)-(m(t)))/m(0) on the total mass of surfactants remains below 10 -6 where m(t) = V c(t) dv. (f) First-order accuracy of Γ as the resolution is increased. These results can be reproduced using the code in[START_REF] Farsoiya | Rotating Circle Test[END_REF].

  consider a circle of size d b /L = 0.4, initial surfactant concentration as Γ(θ, 0) = 2 + sin θ with the flow field representing solid body rotation u(r, θ) = d b ω/2 e θ . The analytical solution is Γ(θ, t) = 2 + sin ( t + θ)e -4 t/P e (11) where t = ωt and P e = d 2 b ω/D s The test is performed for a Peclet number of P e = 160, and three uniform resolutions, i.e., ∆x = d b /12.8, d b /25.6, d b /51.2.

Figure 3 (

 3 a), (b), and (c) show the initial, intermediate, and final (tω/2π = 0.4, 1.6 and 3.2) surfactant concentrations on the rotating circle for the highest resolution.

Figure 3 (

 3 d) shows the time evolution of surfactant concentration for three different resolutions and having more than twenty grid points per diameter of the circle provides good agreement with the analytical solution (red curve).

  Figure 3(e) shows the time evolution of the relative error on surfactant mass, which remains negligible, as Ẽ = (m(0) -(m(t)))/m(0) < 10 -6 . First-order accuracy can be seen in Fig.3(f) for the error on surfactant concentration computed relative to the analytical solution with respect to the number of grid cells per unit diameter d b /∆x.

FIG. 4 .

 4 FIG. 4. Migration of a droplet in a thermal gradient (no buoyancy) for Re = 0.132, Ca = 0.066, ∆x = d b /68 (a) Axisymmetric droplet interface shown by the white curve, velocity vectors show the direction of migration towards the temperature gradient shown in color. (b) 3D droplet (sliced in half) showing the velocity vectors, with magnitude color-coded. (c) Migration velocity as a function of time. The squares show the mean velocities of axisymmetric and 3D droplets respectively, converging towards the theoretical prediction[START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF]. These results can be reproduced using the code in[START_REF] Farsoiya | Migrating Bubble Test[END_REF].

4 .

 4 Surfactant transport with topological changes in Plateau-Rayleigh instabiltiy One critical phenomenon for interfacial flows is topological change including like breakup and coalescence. We test the present framework for mass conservation during topological change, by considering the classic Plateau-Rayleigh instability, occurring in circular jets. We initialize an axisymmetric cylindrical thread having diameter d b , density ρ b , and viscosity µ b perturbed axially with wavenumber k in a fluid having density ρ l and viscosity µ l . The fluid properties are set to ρ b /ρ l = 1, µ b /µ l = 100, Oh = µ b / √ ρ b γd b = 0.01 and Pe s = γ/(kρ b )/D s = 10 where γ and D s is the surface tension coefficient between the fluids and surfactant diffusivity respectively. Such configuration is prone to the Plateau-Rayleigh instability.

Figure 5 (

 5 a) shows the thread is initially perturbed as, 2r(z)/d b = 1 -2a p /d b cos(kz) with wavenumber kd b = 0.2π, and amplitude a p /d b = 0.2.A non-uniform surfactant distribution is initialized as Γ = Γ/Γ 0 = 1 + Γ p /Γ 0 cos(kz) similarly to the interface with wavenumber k and amplitude Γ p /Γ 0 = 0.5. The surface tension is kept constant throughout as we test only

FIG. 5 .

 5 FIG. 5. Time stills showing surfactant concentration(Blue-Red) and interface (Black curve) of axisymmetric thread undergoing Plateau-Rayleigh instability. (a) tωc = 0 The initial condition of the perturbed thread (b) tωc = 0.94 Thinning of the thread before breakup (c) tωc = 1.1 surfactant remains at the interface after the breakup (d) tω = 1.57 Subsequent breakup of satellite droplet (e) tω = 1.6 Coalescence of satellite droplets (f) tωc = 1.73 Oscillations before equilibrium state. (g) During themultiple topological changes (break-up and coalescence, the surfactant remains on the interface and its total mass is conserved. These results can be reproduced using the code in[START_REF] Farsoiya | Topological Test[END_REF].

Figure 5 (

 5 f) tω c = 1.73 shows the oscillations of droplets before reaching equilibrium. The net sum of distributed mass among all the droplets is shown as a black dashed line in Figure5(g) showing good conservation of mass.

FIG. 6 .

 6 FIG. 6. (a) Schematic of a rising bubble at the terminal state in the presence of surfactant. The stagnant cap angle θ0 is measured from the nose of the bubble. Figure adapted from [86]. (b) Still image from the present study, illustrating the bubble interface (solid black line), the surfactant concentration (blue colors, with dark blue indicating higher concentration) and the formation of the stagnant cap at the rear of the bubble. The adaptive mesh refinement (AMR) around the bubble is shown for L11 ∆x = L/2 11 , with the number of cells of ≈ 1.4 × 10 4 yielding a significantly smaller number of grid points than a simulation with a uniform grid for which the number of cells would be (2 11 ) 2 ≈ 4.2 × 10 6 .

  Three different resolutions are considered, with smallest grid size ∆x being L9 ∆x ≡ L/2 9 ≡ d b /25, L10 ∆x ≡ L/2 10 ≡ d b /50, L11 ∆x ≡ L/2 11 ≡ d b /100.

Figure 6 (

 6 b)

FIG. 7 .

 7 FIG. 7. Bubble shape, dynamics and surfactant concentration at steady (terminal velocity) state for Ga = 100, Bo = 0.5, for a clean bubble (a) Ma = 0, and increasing contamination so that (b) Ma = 0.2, (c) Ma = 2, and (d) Ma = 20. The thick red curve shows the bubble's (sharp) interface (VoF T = 0.5 contour). Flow around the bubble is visualized through the streamlines in the reference frame of the bubble, shown with lines and arrows colored blue-white-red; vorticity is colored blackgreen-white. Surfactant concentration is colored black-white (and localized at the interface). The addition of surfactants shows the immobilization of the interface with increasing Ma number. The streamlines show the formation of a wake. Examining the magnitude of the velocity inside the bubble confirms the reduction of the internal circulation.

FIG. 8 .

 8 FIG. 8. Interfacial velocity in the reference frame of the bubble for Ga = 100, Bo = 0.5 at steady (terminal) state t g/d b = 5 as a function of the polar angle θ measured from the nose of the bubble. (a) Ma = 0 (b) Ma = 0.2 (c) Ma = 2 (d) Ma = 20. The tangential component us/u b (black) and normal un/u b (blue) of the bubble are shown. The normal component is zero as there is no flow across the bubble. The normalized surfactant concentration Γ = Γ/Γ0 is shown in red read with the axis on the right. The tangential velocity vanishes at the stagnant cap angle θ0.The stagnant cap angle corresponds to the vanishing tangential velocity and is θ0 = π for a clean bubble. The stagnant cap angle decreases as the Marangoni number is increased, and also corresponds to the maximum surfactant concentration. The presence of a small concentration of Γ below the stagnant cap angle θ < θ0 is caused by the non-zero surface diffusivity of surfactant Ds > 0.

FIG. 9 .

 9 FIG. 9. (a,b) Evolution of the bubble rise velocity as a function of time for Bo = 0.5 and (a) Ga = 5, (b) Ga = 100 for a grid resolution ∆x = d b /50. The Ma number increases from 0 (clean bubble) to 2 from blue to red. The terminal velocity of the bubble decreases with an increase in Ma, while the reduction of the rise speed saturates at high Ma. (c) Effect of the Marangoni number on the terminal velocity, normalized by the terminal velocity in the clean case u b (∞)/u c b (∞). Marangoni effects start to be important for smaller values of Ma for smaller Ga. All curves saturate for M a > 10. (d) Normalized drag coefficient (CD -C rigid )/(C rigid -C clean ) as a function of the stagnant cap angle, which decreases with Ma. The clean bubble has θ0 = π and large surfactant effects (high Ma) lead to θ0 = 0, which corresponds to the rigid body translation. Numerical data are compared to the theoretical curve from Sadhal and Johnson [92] (black line). Inset shows the variation of the stagnant cap angle with an increase of the Marangoni number. The reduction of rise velocity plateaus as θ0 → 0.

FIG. 10 .

 10 FIG. 10. Helicoidal bubble paths. (a) Trajectory evolution, t/ g/d b = 0 -82.1, for Ga = 100.25, Bo = 10, matching the numerical simulations from Cano-Lozano et al. [101] and probing similar regimes to Tagawa et al. [102]. The clean case (Ma=0) is shown in black while a contaminated case (Ma=1) is shown in red. We find a helical trajectory in the clean case while the contaminated bubble follows a zigzag trajectory. (b,c) Bubble interface, streamlines and surfactant concentration at the interface at t g/d0 = 60 (squares in (a)) for Ga = 100.25, Bo = 10, (b) Ma = 0, (c) Ma = 1 (with surfactant concentration Γ = Γ/Γ0). The blue-colored field in (c) shows the presence of surfactants and forming a stagnant cap. The immobile interface at the cap has modified the flow at the rear of the bubble forming a wake.

  and Trygvasson

	Interface	Dimensions Mesh			Citations	Statement of Open-
	Tracking						Source
		2D	Structured uniform mesh		Stone and Leal [35]
	Boundary	3D	Structured uniform mesh		Li and Pozrikidis
	Element					[36], Feigl et al. [52]
		3D	Unstructured mesh		Yon and Pozrikidis
						[37]
		2D	Structured uniform mesh for veloc-	Muradoglu and Tryg-
			ity/pressure, Unstructured for in-	gvason [39]
			terface tracking
	Front Tracking	3D	Structured uniform mesh for veloc-	Muradoglu and Tryg-
			ity/pressure, Unstructured for in-	gvason [29]
			terface tracking
		3D	Structured Adaptive Mesh Refine-	De Jesus et al. [40]
			ment (AMR) for velocity/pressure,
			Unstructured for interface tracking
	Level Set-Front	3D	Structured uniform for veloc-	Shin et al. [44]
	Tracking		ity/pressure,	Unstructured for
			interface tracking
		3D	Structured uniform mesh		Xu et al. [46]
		2D	Structured uniform		Lakshmanan	and
						Ehrhard [53]
	Level-Set	2D	Structured uniform mesh		Piedfert et al. [47]
		2D	Structured uniform mesh		Atasi et al. [27]
	Arbitrary	3D	Unstructured	Adaptive	Mesh	Dieter-Kissling et al.
	Lagrangian-		Refinement			[54]
	Eulerian (ALE)
	Diffuse-Interface 2D & 3D	Structure	Adaptive	Mesh	Teigen et al. [50, 55]
			Refinement
	Algebraic VoF	2D	Unstructured	Adaptive	Mesh	Antritter et al. [56]
			Refinement
		2D & 3D	Structured uniform mesh		Renardy et al. [57]
	Geometric VoF	2D	Unstructured	Adaptive	Mesh	James and Lowen-
			Refinement(AMR)		grub [58]
	Geometric	2D & 3D	Structured Adaptive Mesh Re-	Present-Study	Open-Source
	VoF-Phase		finement (AMR)
	field

TABLE I

 I

. Summary of representative studies of numerical methods for surfactant-laden flows in interfacial fluid dynamics. All these methods solve for surfactant transport equations coupled with Marangoni stresses and two-phase Navier Stokes equations, with different interface tracking and mesh structure. The present work proposes a fully open-source geometric VoF method coupled to a phase-field representation.

TABLE II .

 II List of simulations for various Ga, Bo, and Ma. C. Effect of the Marangoni number on the terminal axisymmetric rise velocity

ACKNOWLEDGMENTS

This work was supported by NSF grant 2242512 and NSF CAREER 1844932 to LD. The authors are pleased to acknowledge that the work reported on in this paper was substantially performed using the Princeton Research Computing resources at Princeton University which is a consortium of groups led by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's Research Computing.

Declaration of Interests. The authors report no conflict of interest.

Appendix A: Grid Convergence

We demonstrate grid convergence in the bubble rise velocity configuration. Simulations at increasing levels of resolution, corresponding to an increasing number of grid points per bubble diameter are performed. Figure 11(a,b) shows the normalized rise velocities of the bubble for Ga = 5, 100, Bo = 0.5. Three different resolutions were used for simulations with adaptivity over dependent variables, L9≡

The grid convergence is achieved for Ga = 5 in the resolution L9 and L10. As the Ga number is increased the resolution requirement increases (as the boundary layer is thinner). Fig. 11(b) shows the grid convergence for the highest Ga = 100, converging between L10 and L11. The L10 hence is the resolution where all the simulations show convergence, we used this for the analysis in the present study.