Transient bubble rising in the presence of a surfactant at very low concentrations

Contents lists available at ScienceDirect International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmulflow Transient bubble rising in the presence of a surfactant at very low concentrations D. Fernández-Martínez a, M.G. Cabezas a,∗, J.M. López-Herrera b, M.A. Herrada a, J.M. Montanero b aDepartamento de Ingeniería Mecánica, Energética y de los Materiales and Instituto de Computación Científica Avanzada (ICCAEx), Universidad de Extremadura, E-06006 Badajoz, Spain bE.T.S.I., Depto. de Ingeniería Aeroespacial y Mecánica de Fluidos, Universidad de Sevilla, Camino de los Descubrimientos s/n 41092, Spain A R T I C L E I N F O Keywords: Bubble rising Surfactant Dynamic adsorption layer A B S T R A C T We study the formation of the dynamic adsorption layer when a bubble is released in a tank containing water with a tiny amount of surfactant. The influence of the sorption kinetic constants is examined by comparing the experiments with Sodium Dodecyl Sulfate (SDS) and Triton X-100. The experiments allowed us to determine the parameter conditions that lead to a stable bubble rising and to validate the simulation. A simple scaling analysis and the simulation show that the formation of the dynamic adsorption layer can be split into three phases characterized by disparate time scales. The mechanisms controlling those phases are surfactant convection, adsorption–desorption, and diffusion. The amount of surfactant adsorbed onto the interface increases monotonously throughout the three phases. The experiments and the simulation show that the rising velocity reaches a maximum at times of the order of 𝑘−1 𝑑 (𝑘𝑑 is the desorption constant) when the dynamic adsorption layer is practically formed. This occurs even when only traces of surfactant are present in the liquid. The non-monotonous behavior of the maximum surfactant surface concentration is explained in terms of the reverse flow in the rear of the bubble right after the bubble release. This work contributes to the understanding of the complex interplay between hydrodynamics and surfactant transport and kinetics over bubble rising. 1. Introduction Many technological and industrial processes involve the interaction between gases and liquids. Among them, we can mention waste-water treatment, chemical (Jakobsen, 2014) and biochemical (Doran, 2013) reactors, nuclear engineering (Lahey, 1990), and metallurgical bubble column reactors. The heat and mass transfer across the interface of bubbles dominate the phenomena occurring in these processes. This interfacial transfer is considerably influenced by the size, shape, trajectory, and velocity of these bubbles. The processes mentioned above involve large populations of bubbles interacting among them. Single-bubble dynamics serve as a basis for analyzing complex multi-bubble systems. The rising of a bubble in still water is a paradigmatic problem that has been analyzed for several centuries and continues to capture the attention of many researchers today (Saffman, 1956; Sanada et al., 2008; Tripathi et al., 2015; Cano-Lozano et al., 2016a,b; Kure et al., 2021; Bonnefis et al., 2023). ∗Corresponding author. E-mail address: [email protected] (M.G. Cabezas). When a bubble is released in a liquid bath containing surfactant, the surfactant molecules adsorb onto the free surface during the bubble rising. The molecules are advected by the outer current toward the rear of the bubble. Surfactant accumulates in that region, where desorption considerably increases. This process results in an even distribution of surfactant over the interface, producing a Marangoni stress that substantially alters the forces acting on the bubble (Takagi and Matsumoto, 2011). The bubble rising in the presence of a surfactant is a complex phenomenon in which fluid-dynamic and physicochemical processes are coupled, especially for large non-spherical bubbles. The bubble shape change implies variations of the interfacial area, which induces sorption processes counteracting the interface expansion or compression (Levich, 1962). Interface expansion/compression, sorption kinetics, and convection over the bubble surface compete to establish the socalled dynamic adsorption layer (Dukhin et al., 1998, 2015; Ulaganathan et al., 2016; Zawala et al., 2023). https://doi.org/10.1016/j.ijmultiphaseflow.2025.105205 Received 30 September 2024; Received in revised form 28 February 2025; Accepted 28 February 2025 International Journal of Multiphase Flow 188 (2025) 105205 Available online 8 March 2025 0301-9322/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). D. Fernández-Martínez et al. The evolution of the bubble velocity in the presence of a surfactant significantly differs from that in clean water (Saffman, 1956; Sanada et al., 2008; Tripathi et al., 2015; Cano-Lozano et al., 2016a,b; Herrada and Eggers, 2023; Bonnefis et al., 2023). For sufficiently large surfactant concentrations, the bubble velocity reaches its maximum value after initial acceleration. Then, the bubble decelerates until a plateau value (the terminal velocity) is attained. This non-steady rising process corresponds to the major part of the bubble trajectory in many practical situations. The non-monotonous dependence of the bubble velocity on the vertical position (time) is usually referred to as the local velocity profile (Krzan and Malysa, 2002; Krzan et al., 2007; Kracht and Finch, 2010; Dukhin et al., 2015). The local velocity profile of a rising bubble is a fingerprint of the dynamic adsorption layer’s transient behavior. It is strongly affected by the adsorption of surface-active species at the bubble surface. For this reason, the bubble path allows one to obtain information about the surfactant sorption kinetics that no other technique can provide, particularly when comprising proteins at various solvent conditions, such as pH and ionic strength. Adsorption and desorption occur at the same localization in most experimental configurations. Conversely, the surfactant is essentially adsorbed onto the front mobile part of a rising bubble, while desorption occurs mainly within the rear stagnant cap. This feature may favor the study of the adsorption kinetics of large molecules (Dukhin et al., 2015). There is abundant experimental information about surfactants’ effect on the rising single bubble. For instance, it is well known that the height and width of the velocity profile maximum decrease as the surfactant concentration increases (Krzan et al., 2007). Kracht and Finch (2010) found that the instantaneous values of the bubble velocity and aspect ratio were correlated in their experiments. This suggests that surfactants can essentially affect bubble rise velocity through the bubble shape. Fdhila and Duineveld (1996) showed that the bubble slowdown occurs when the surfactant coverage reaches approximately half the bubble surface. The steady-state velocity was independent of the surfactant concentration for sufficiently large concentrations (Sam et al., 1996; Zhang and Finch, 2001). Several studies have examined the effects of specific surfactants. For instance, adding a surfactant to 𝛽-lactoglobulin (Ulaganathan et al., 2014) solutions or varying the pH of bovine serum albumin solutions (Zawala et al., 2010) dramatically alter the velocity profile of a rising bubble. The time for establishing an immobile rigid surface layer at the rising bubble surface becomes shorter with increasing pH (Ulaganathan et al., 2016). The bubble mobility is more influenced than the bubble shape when adding sodium dodecyl sulfate (SDS) to a pseudo-plastic liquid (Tzounakos et al., 2004). Tagawa et al. (2014) studied the bubble path instability in 1-pentanol and Triton X-100 solutions, observing that the drag force monotonically increases with the surfactant concentration, while the lift force showed a non-monotonic behavior. Cetyltrimethyl ammonium bromide (CTAB) dissolved in ultrapure water exhibits two distinct velocity stages, whereas Tween 80 solutions display a three-stage velocity profile (Luo et al., 2022). The information in all the experimental studies is limited to the time-dependent bubble shape and velocity. Interpreting these experimental data requires theoretical models to describe the unsteady dynamic adsorption layer. Approximate analytical methods have been proposed to describe the unsteady dynamic adsorption layer and its influence on bubble rising (Zholkovskij et al., 2000). The results only apply to spherical shapes and low Reynolds numbers (Zholkovskij et al., 2000). In the quasisteady approximation, the instantaneous velocity has been assumed to be the steady velocity corresponding to the amount of the solute accumulated at that instant (Dukhin et al., 2015). Previous studies have used this theory to obtain the parameters of adsorption–desorption kinetics from the local velocity profile measured in the experiments. For non-spherical bubbles and large Reynolds numbers, accurate transient numerical simulations are required to produce reliable results. Most numerical studies focus mainly on the final stationary (stable) (Mestel, 1994; Lakshmanan and Ehrhard, 2010; Li and Mao, 2001; Takemura, 2005; Fukuta et al., 2008; Tasoglu et al., 2008; Kentheswaran et al., 2023; Dani et al., 2022) or oscillating (unstable) (Pesci et al., 2018; Farsoiy et al., 2024) regime reached by the rising bubble. The initial stage following the bubble release has received less attention. Transient axisymmetric simulations show that the bubble attains a maximum velocity before slowing down to its steady-state velocity (Liao and McLaughlin, 2000). A speed comparable to the steady-state speed in pure water is reached even if the bubble is contaminated before releasing it (Liao and McLaughlin, 2000). The elasticity number and the bulk Peclet number significantly affect the initial transient motion as well (Tasoglu et al., 2008). Cuenot et al. (1997) considered a simplified dynamic problem in which the bubble remained spherical and its velocity was constant throughout the contamination process, focusing on the physicochemical processes neglected in many previous approaches. They described the temporal evolution of the relevant interfacial quantities to show the formation of the dynamic adsorption layer. The rising of a bubble covered by Triton X-100 was simulated by decoupling the solution of the fluid flow and mass transfer under different approximations (Zhang et al., 2001). The best results were obtained assuming the stagnant cap model and surfactant transfer to the interface controlled by diffusion. Matsumoto et al. (2006) showed the influence of the adsorption and desorption constants on the local velocity profile and the temporal evolution of the drag coefficient. The height and width of the velocity profile maximum decrease as the adsorption (desorption) constant increases (decreases). Pesci et al. (2018) conducted three-dimensional direct numerical simulations of a bubble rising in a liquid containing a soluble surfactant. The numerical results agreed with their experiments (Pesci et al., 2017). They concluded that the initial transient stage is very sensitive to the initial surface concentration. The quasi-steady state of the rise velocity is reached without adsorption and desorption being necessarily in equilibrium. Rubio et al. (2024) have recently explained how traces of surfactant can significantly change the steady (terminal) regime of bubble rising. The tiny surface tension variation occurs within an extremely thin, diffusive surface boundary layer. This produces a Marangoni stress confined within that layer that immobilizes the free surface and drastically increases the outer viscous stress only in the diffusive layer. The dynamic adsorption layer described by Rubio et al. (2024) qualitatively differs from those described in other works under different conditions, in which the Marangoni and outer viscous stresses are spread over a considerable portion of the bubble surface. In this work, we will study numerically how the singular dynamic adsorption layer appearing at very low surfactant concentrations (Rubio et al., 2024) grows over time. We will analyze the temporal evolution of the relevant interfacial quantities and explain the effect of surfactant convection, kinetics, and diffusion on the growth of the dynamic adsorption layer. The numerical results of Matsumoto et al. (2006) for deformable bubbles are restricted to small Reynolds and Peclet numbers. Our boundary-fitted numerical method is similar to that employed in that work. We use a spectral collocation technique (Khorrami, 1989) to accumulate the grid points next to the interface, which allows us to resolve the extremely thin surfactant boundary layer on the outer side of the bubble surface. In this way, we solve the complete model (Matsumoto et al., 2006) even for realistic values of the Reynolds and Peclet numbers. It is well known that the surfactant monolayer enhances the path instability, reducing the critical radius above which helical and zigzagging trajectories are observed (Tagawa et al., 2014; Rubio et al., 2024). Non-axisymmetric instabilities are not allowed in our axisymmetric simulations. For this reason, we will conduct experiments to determine the parameter conditions leading to a straight path (axisymmetric flow) and to validate our numerical solutions for those conditions. International Journal of Multiphase Flow 188 (2025) 105205 2 D. Fernández-Martínez et al. Fig. 1. Sketch of the numerical domain. The blue and red outer boundaries correspond to the inlet and non-reflecting boundary conditions, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 2. Methods 2.1. Governing equations Consider a bubble of radius 𝑅= [3𝑉∕(4𝜋)]1∕3 (𝑉 is the bubble volume), density 𝜌(𝑖), and viscosity 𝜇(𝑖) rising in a liquid of density 𝜌(𝑜) and viscosity 𝜇(𝑜). The surface tension of the clean interface is 𝜎𝑐, while the gravity acceleration is 𝑔. We dissolve a surfactant in the liquid at the concentration 𝑐∞. The surfactant can be considered soluble because the sorption characteristic times are smaller than the hydrodynamic characteristic time. In the framework of the model considered here, the surfactant properties are the volumetric diffusion coefficient 𝑜, the surface diffusion coefficient 𝑠, the adsorption and desorption constants, 𝑘𝑎 and 𝑘𝑑, and the maximum packing density 𝛤∞. The hydrodynamic equations are solved in a cylindrical system of coordinates (𝑟, 𝑧) whose origin solidly moves with the bubble’s upper point (Fig. 1). The continuity and momentum equations are 𝛁⋅𝐯(𝑗)= 0,(1) 𝜌(𝑗)𝐷𝐯(𝑗) 𝐷𝑡 = −𝜌(𝑗)(𝑔+𝑑2ℎ 𝑑𝑡2)𝐞𝐳+𝛁⋅𝝈(𝑗),(2) where 𝐯(𝑗)=𝑢(𝑗)𝐞𝐫+𝑤(𝑗)𝐞𝐳 is the axisymmetric velocity field, the superscripts 𝑗=𝑖 and 𝑜 refer to the inner and outer phases, respectively, 𝐞𝐫 and 𝐞𝐳 are the unit vectors along the axis 𝑟 and 𝑧, respectively, 𝐷∕𝐷𝑡 is the material derivative, ℎ=𝑍−𝑧 is the vertical position of the bubble’s upper point (Fig. 1), 𝝈(𝑗)= −𝑝(𝑗)𝐈+𝝉(𝑗)(3) is the stress tensor, 𝑝(𝑗) is the hydrostatic pressure, 𝐈 is the identity matrix, and 𝝉(𝑗)=𝜇(𝑗)[𝛁𝐯(𝑗)+(𝛁𝐯(𝑗))𝑇](4) is the viscous stress tensor. The velocity field is continuous (𝐯(𝑖)=𝐯(𝑜)) at the free surface. The interface is parametrized in terms of the meridional arc length 𝑠 (0≤𝑠≤𝑠𝑓) as 𝑟𝑠=𝑓(𝑠) and 𝑧𝑠=𝑔(𝑠), where (𝑟𝑠, 𝑧𝑠) is the interface location. Here, 𝑠𝑓 is the arc length corresponding to the bubble’s rear point (Fig. 1). The kinematic compatibility condition reads 𝜕𝑓 𝜕𝑡 +𝑢(𝑖)𝑔′−𝑤(𝑖)𝑓′= 0,(5) where the apostrophe (′) indicates the derivative with respect to 𝑠. The equilibrium of normal and tangential stresses on the interface leads to the equations 𝐞𝐧⋅(𝝈(𝑜)−𝝈(𝑖))⋅𝐞𝐧=𝜎𝜅, 𝐞𝐭⋅(𝝈(𝑜)−𝝈(𝑖))⋅𝐞𝐧=𝜏Ma,(6) where 𝜎 is the local value of the interfacial tension, 𝜅= ∇⋅𝐞𝐧 is (twice) the mean curvature, and 𝐞𝐧=𝑔′𝐞𝑟−𝑓′𝐞𝑧 (𝑓′2 +𝑔′2)1∕2 and 𝐞𝐭=𝑓′𝐞𝑟+𝑔′𝐞𝑧 (𝑓′2 +𝑔′2)1∕2 (7) are the unit vectors normal and tangential to the interface, respectively. In addition, 𝜏Ma =𝜎′ is the Marangoni stress. We neglect the viscous surface stresses in Eqs. (6) because the shear and dilatational viscosities of the surfactant monolayer are very small (Ponce-Torres et al., 2020). We restrict ourselves to low surfactant concentrations, which implies that the surfactant is present as monomers. The monomer volumetric concentration 𝑐(𝑜)(𝐫, 𝑡) in the outer phase is calculated from the conservation equation (Craster et al., 2009; Kalogirou and Blyth, 2019) 𝜕𝑐(𝑜) 𝜕𝑡 +𝐯(𝑗)⋅𝛁𝑐(𝑜)=𝑜𝛁2𝑐(𝑜).(8) Now, we model the transfer of monomers between the bulk and the bubble surface. The net sorption flux =𝑎−𝑑 is calculated as the difference between the adsorption 𝑎 and desorption 𝑑 flux. The kinetic model 𝑎=𝑘𝑎𝑐(𝑜) 𝑠(1 − 𝛤 𝛤∞),𝑑=𝑘𝑑𝛤(9) is adopted to calculate these fluxes. As mentioned above, 𝑘𝑎 and 𝑘𝑑 are the adsorption and desorption constants, respectively, 𝑐(𝑜) 𝑠 is the bulk surfactant concentration evaluated at the interface, 𝛤 is the surfactant surface concentration (the surface coverage, measured in mols per unit area), and 𝛤∞ is the maximum packing density. The surfactant surface concentration 𝛤 verifies the advectiondiffusion equation (Craster et al., 2009) 𝜕𝛤 𝜕𝑡 +𝛁𝑠⋅(𝛤𝐯𝑠) + 𝛤(𝛁𝑠⋅𝐧)(𝐯⋅𝐞𝐧) = 𝑠∇2 𝑠𝛤+,(10) where 𝛁𝑠 is the tangential intrinsic gradient along the free surface, 𝐯𝑠=𝗜𝑠𝐯 is the (two-dimensional) surface velocity, 𝗜𝑠=𝗜−𝐞𝐧𝐞𝐧 is the tensor that projects any vector on that surface, 𝗜 is the identity tensor, and 𝑠 is the surface diffusion coefficient. The dependence of the surface tension 𝜎 on the surface concentration 𝛤 is given by the Langmuir equation of state (Tricot, 1997) 𝜎=𝜎𝑐+𝛤∞𝑅𝑔𝑇ln (1 − 𝛤 𝛤∞),(11) where 𝜎𝑐 is the surface tension of the clean interface, 𝑅𝑔 is the gas constant, and 𝑇 is the temperature. The kinetic model (9) yields the Langmuir isotherm at equilibrium (𝑎=𝑑). Combining the Langmuir and Gibbs isotherms leads to the Langmuir equation of state (11) (Tricot, 1997). This means the model (9) and (11) are consistent. For this reason, we selected Eq. (11) to model the SDS effect, even though it may be less accurate than fitting the experimental data (Ponce-Torres et al., 2020). We assume that the liquid bath is not perturbed by the bubble in the upstream region 𝑟=𝑅𝑜 and 𝑧 > 0 (Fig. 1). Therefore, 𝑤(𝑜)= − 𝑑ℎ 𝑑𝑡 , 𝑢(𝑜)= 0, 𝑝(𝑜)+𝜌(𝑜)𝑔𝑧 =const. (12) in that boundary. The non-reflecting boundary conditions 𝜕𝑢(𝑜) 𝜕𝑧 =𝜕𝑤(𝑜) 𝜕𝑧 = 0 (13) are applied in the downstream region far from the bubble (𝑟=𝑅𝑜 and 𝑧 < 0) to capture the wake (Fig. 1). The surfactant concentration at 𝑟=𝑅𝑜 is 𝑐∞. We consider the regularity conditions 𝑢(𝑗)=𝜕𝑤(𝑗) 𝜕𝑟 =𝜕𝑝(𝑗) 𝜕𝑟 =𝜕𝑐(𝑜) 𝜕𝑟 = 0 (14) International Journal of Multiphase Flow 188 (2025) 105205 3 D. Fernández-Martínez et al. Table 1 Properties for SDS and Triton X-100. The values for SDS were those considered by Rubio et al. (2024). The values for Triton X-100 were taken from Refs. Tagawa et al. (2014) and Lin et al. (1990). CMC (mol/m3)𝑘𝑎 (m/s) 𝑘𝑑 (s−1)𝛤∞ (mol/m2)𝑜 (m2/s) SDS 8.0 2.34 × 10−4 6.4 3.9 × 10−6 8.0 × 10−10 Triton X-100 0.23 1.45 × 10−4 3.3 × 10−2 2.9 × 10−6 2.6 × 10−10 at the symmetry axis 𝑟= 0. The condition 𝑤(𝑜)= 0 at the interface upper point allows us to calculate the bubble’s vertical position ℎ(𝑡). Finally, we specify the bubble’s volume through the equation 𝑉=𝜋∫𝑠𝑓 0 𝑓2𝑔′𝑑𝑠. (15) We start the simulation from the hydrostatic solution corresponding to a spherical bubble. Consider a bubble ‘‘instantaneously’’ inflated in a tank of still water with surfactant. The Damkohler number Da=𝜏𝐷∕𝜏𝑘 relates the time scales 𝜏𝐷 and 𝜏𝑘 corresponding to the surfactant transfer limited by diffusion and sorption, respectively (see the Appendix A) (Manikantan and Squires, 2020). Considering the values shown in the next section, the Damkohler number takes values of the order of 10 and 102 for SDS and Triton X-100, respectively, indicating that the surfactant transfer to the bubble is limited by diffusion in both cases. In the experiment, the bubble remained attached to the capillary for approximately 5 s. The surface concentrations of SDS and Triton X-100 at that instant were similar to the equilibrium concentration for SDS (see Fig. A.13 in the Appendix A). For this reason, we consider this value as the initial surface concentration in all the simulations. The above equations are numerically integrated with a variant of the boundary-fitted spectral method proposed by Herrada and Montanero (2016). The major difficulty associated with a soluble surfactant is the existence of a very thin diffusive boundary layer next to the interface for the small diffusion coefficients of most surfactants (large Peclet numbers). We use Chebyshev spectral collocation points to accumulate the grid points next to the interface (Herrada and Montanero, 2016), facilitating the resolution of this layer (Rubio et al., 2024). A grid sensitivity analysis is shown in the Supplemental Material. 2.2. Experimental method In an experiment, nitrogen was injected through a needle to form the bubble in the center of the tank bottom. The bubble detached from the needle and rose across the tank until it reached the free surface. We used a virtual binocular stereo vision system (Luo et al., 2022) to image two perpendicular views of the rising bubble. The images were processed at the pixel level (Canny, 1986) to determine the bubble shape and velocity. Details of the experimental setup and the image processing analysis can be found in the Supplemental Material. This experimental method was validated by Rubio et al. (2024) from comparison with the results for clean water obtained by Duineveld (1995). We consider SDS and Triton X-100 in our experiments. SDS is an anionic surfactant with a molecular weight of 288.4 g/mol, while Triton X-100 is a nonionic surfactant with a molecular weight of 647 g/mol. Table 1 shows the relevant properties of these surfactants. The diffusion coefficients of SDS and Triton X-100 are commensurate with each other. SDS is slightly more active than Triton X-100, as indicated by the value of 𝛤∞. The SDS adsorption constant is also similar to that of Triton X-100. Interestingly, the desorption rate 𝑘𝑑 of Triton X-100 is two orders of magnitude smaller than that of SDS. For small surfactant concentrations, 𝛤=𝑘𝑎𝑐∕𝑘𝑑 at equilibrium [Eqs. (9)]. This means that one needs to dissolve much fewer moles of Triton X-100 in water to achieve the same surfactant surface concentration at equilibrium. This is reflected in the critical micelle concentration, around 35 times smaller in the Triton X-100 case. 3. Experimental results As mentioned in the Introduction, the goal of our experimental study is twofold: (i) to determine the parameter conditions that lead to an axisymmetric, straight (stable) bubble rising and (ii) to validate the numerical method. Fig. 2 shows the three types of trajectories observed in our experiments with Triton X-100. Similar results were obtained for SDS (Rubio et al., 2024). Fig. 2a corresponds to a bubble following a quasi-straight trajectory. The small amplitude oscillations can be attributed to the camera vibration while moving (we did not observe those oscillations when the camera was at rest). The path is slightly tilted (the tilt angle is smaller than 0.1◦), probably due to a small asymmetry in the bubble formation process. For 𝑐∞∕𝑐cmc = 5 × 10−4 [case (b)], the bubble motion remains stable until 𝑧≃ 350 mm. Helical instability develops at larger distances from the ejector. This experiment highlights the importance of conducting experiments with large tanks and may explain the discrepancies among the critical conditions found in previous works. The distance from the ejector at which the helical instability develops decreases as the surfactant concentration increases (Fig. 2c). For 𝑐∞∕𝑐cmc = 3.3 × 10−3 [case (d)], the helical instability evolves toward a zig-zag motion, as shown by the projection of the bubble path onto the 𝑥𝑦 plane (the red line corresponds to the zig-zag motion). The zig-zag instability arises close to the bubble ejector for the highest surfactant concentration considered in our analysis [case (e)]. Figs. 3–5 compare the effect of SDS and Triton X-100 on the bubble velocity and aspect ratio. We monitor the bubble velocity and the aspect ratio to safely determine whether the bubble has reached a steady motion. The results in Figs. 3and 4 were obtained for the same value of the relative concentration 𝑐∞∕𝑐cmc of SDS and Triton X-100. For small concentrations, Eqs. (9) yield 𝛤eq∕𝛤∞= [𝑘𝑎𝑐cmc∕(𝑘𝑑𝛤∞)] 𝑐∕𝑐cmc at equilibrium, where 𝑘𝑎𝑐cmc∕(𝑘𝑑𝛤∞) = 75 and 348 for SDS and Triton X-100, respectively. This means that the same value of 𝑐∕𝑐cmc leads to a much higher equilibrium surface coverage 𝛤eq∕𝛤∞ and density 𝛤eq in the Triton X-100 case. However, the Triton X-100 equilibrium surface coverage is reached at much longer times due to the much smaller value of the desorption constant. This suggests that Triton X-100 must produce a smaller effect on a relatively short time scale, but this effect must eventually exceed that produced by SDS at sufficiently large times. The results in Figs. 3and 4 are consistent with the above prediction. The bubble velocity and aspect ratio for SDS are smaller than those for Triton X-100. Nevertheless, a crossover is expected to occur beyond the maximum height analyzed in the experiment. The results in Figs. 5and 6 were obtained for practically the same value of the absolute concentration 𝑐∞ of SDS and Triton X-100. In this case, the adsorption flux 𝑎≃𝑘𝑎𝑐∞ is expected to take similar values for the two surfactants. This explains why the bubbles covered with SDS and Triton X-100 exhibit similar velocities and aspect ratios in the first stage of the bubble rising. As explained in Section 2.2, the same volumetric concentration leads to a much higher Triton X-100 surface concentration at equilibrium. This explains why the Triton X100 effects are larger than those caused by SDS at the same absolute concentration (Figs. 5and 6). The experimental local velocity profile exhibits the so-called ‘‘overshooting’’ phenomenon: the velocity reaches a maximum and subsequently decreases to its terminal value. Consider the following balance of forces acting on the bubble: =𝐷0(𝑣𝑧(𝑡)) + (𝑡),(16) where  is the approximately constant buoyancy force, 𝐷0 represents the drag force exerted on the same bubble but in the absence of surfactant, and  is the extra drag force associated with the surfactant monolayer. Interestingly, the trajectories with and without SDS are practically the same for 𝑧 ≲ 50 mm (Fig. 5). This implies that ≃ 0 during this part of the bubble rising probably because the increase International Journal of Multiphase Flow 188 (2025) 105205 4 D. Fernández-Martínez et al. Fig. 2. Bubble trajectory for 𝑅= 0.76 mm and 𝑐∞∕𝑐cmc = 10−4 (a), 5 × 10−4 (b), 10−3 (c), 3.6 × 10−3 (d), and 10−2 (e). The graphs also show the projections of the trajectories onto the planes (𝑥, 𝑦), (𝑥, 𝑧), and (𝑦, 𝑧). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 3. Bubble velocity 𝑣𝑧 and aspect ratio 𝜒 as a function of the vertical position 𝑧 of the center of gravity for 𝑅= 0.76 mm and 𝑐∞∕𝑐cmc = 5×10−4 . The solid and dashed lines correspond to the stable and oscillatory parts of the bubble trajectory, respectively. in the viscous friction is balanced by the decrease in the bubble deformation (the bubble adopts a more aerodynamic shape) due to the surfactant. For 𝑧 ≳ 50 mm, the bubble shape hardly changes while (i) the viscous friction increases due to the growth of the Marangoni stress Fig. 4. Bubble velocity 𝑣𝑧 and aspect ratio 𝜒 as a function of the vertical position 𝑧 of the center of gravity for 𝑅= 0.76 mm and 𝑐∞∕𝑐cmc = 10−3 . The solid and dashed lines correspond to the stable and oscillatory parts of the bubble trajectory, respectively. within the surfactant surface boundary layer and (ii) the separation point slightly displaces upstream due to the displacement of that layer, as explained in Section 4.2. This results in an increase in the drag force and a decrease in the bubble velocity (the overshooting phenomenon). International Journal of Multiphase Flow 188 (2025) 105205 5 D. Fernández-Martínez et al. Fig. 5. Bubble velocity 𝑣𝑧 and aspect ratio 𝜒 as a function of the vertical position 𝑧 of the center of gravity for 𝑅= 0.66 mm and the concentrations 𝑐∞= 8.3 × 10−4 mol/m3 of Triton X-100 and 𝑐∞= 8 × 10−4 mol/m3 of SDS. The solid and dashed lines correspond to the stable and oscillatory parts of the bubble trajectory, respectively. The red circles correspond to the experimental data of Zhang and Finch (2001) for 𝑅= 0.70 mm and Triton X-100 at the concentration 𝑐∞= 7.5 × 10−4 mol/m3. The arrow indicates the maximum in the velocity for SDS due to the overshooting phenomenon. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 6. Bubble velocity 𝑣𝑧 and aspect ratio 𝜒 as a function of the vertical position 𝑧 of the center of gravity for 𝑅= 0.76 mm and the concentrations 𝑐∞= 8.3×10−4 mol/m3 of Triton X-100 and 𝑐∞= 8 × 10−4 mol/m3 of SDS. The solid and dashed lines correspond to the stable and oscillatory parts of the bubble trajectory, respectively. Figs. 5and 6 allow us to analyze the influence of the bubble radius. The surfactant effect increases as the bubble radius decreases. Specifically, the terminal velocity reduction caused by SDS is larger for the 𝑅= 0.66 mm case, and Triton X-100 destabilizes the bubble path at smaller values of 𝑧 in that case. This result can be expected because the surfactant produces interfacial effects, and the surface-to-volume ratio increases as 𝑅 decreases. All the trajectories analyzed in Figs. 3–5 become unstable except for the case SDS at 𝑐∞= 8 × 10−4 mol/m3, where the bubble reaches a steady state. In the next section, we numerically analyze the temporal evolution of the dynamic adsorption layer for that case. The analysis of this axisymmetric flow is realistic because the flow does not suffer from non-axisymmetric instabilities over the entire bubble trajectory. 4. Numerical results 4.1. Validation of the numerical method We calculated the numerical solution for 𝜌(𝑖)= 1.2 kg∕m3, 𝜌(𝑜)= 997 kg∕m3, 𝜇(𝑖)= 1.83 × 10−5 kg/(m⋅s), 𝜇(𝑜)= 9.3 × 10−4 kg/(m s), and 𝜎𝑐= 72 mN/m. The values of 𝛤∞ were taken from Table 1. The diffusion coefficient 𝑜 was set equal to that of SDS to favor the convergence of the numerical method, which may slightly affect the terminal velocity (Rubio et al., 2024). We considered 𝑠= 4.2 × 10−7 m2/s. The value of 𝑠 was increased with respect to the expected experimental value to favor the convergence of the numerical method. We verified that this does not significantly affect the result. Specifically, an increase in 𝑠 of the order of 10−1 produces an increase in the terminal velocity of the order of 10−5. This occurs because the Peclet number is sufficiently large for the surfactant surface distribution to be dominated by the advection and adsorption–desorption kinetics (Tagawa et al., 2014). The reported values of the adsorption and desorption constants of most surfactants are relatively inconsistent. For this reason, these values are adjusted in the simulation to reproduce the bubble path observed experimentally. In fact, bubble rising has been proposed to determine the values of those constants. We followed different strategies to determine the adsorption and desorption constants of Triton X-100 and SDS. For the small surfactant concentrations considered in our analysis, 𝛤∕𝛤∞≪1, which implies that 𝜎𝑐−𝜎≃𝐿𝑑𝑐∞𝑅𝑔𝑇(17) at equilibrium. This equation shows the critical role played by the depletion length 𝐿𝑑=𝑘𝑎∕𝑘𝑑 (or, equivalently, the Langmuir equilibrium adsorption constant 𝐾𝐿=𝐿𝑑∕𝛤∞) in the bubble dynamics. A reliable value of the depletion length, 𝐿𝑑= 4.4 mm, has been determined for Triton X-100 in previous studies (Prosser and Franses, 2001; Tagawa et al., 2014). For this reason, it is reasonable to fit the numerical bubble path to the experimental one by varying 𝑘𝑎 and 𝑘𝑑 while keeping 𝐿𝑑 equal to that value (Lin et al., 1990; Fdhila and Duineveld, 1996). We have done this with Triton X-100. The surface tension measurement for SDS does not lead to a reliable value of 𝐿𝑑 (Fdhila and Duineveld, 1996; Prosser and Franses, 2001). Alternatively, Rubio et al. (2024) took the desorption constant 𝑘𝑑 from accurate experimental measurements and used 𝑘𝑎 as the fitting parameter, whose experimental value exhibits more uncertainty. The value shown in Table 1 led to the best agreement between the numerical terminal velocity and the experimental one for the case analyzed in this work. We take that value for the transient simulations conducted here. The depletion length 𝐿𝑑 obtained from this procedure is consistent with the value given by Fdhila and Duineveld (1996). Fig. 7 compares the numerical and experimental local velocity profiles for Triton X-100 and SDS. The numerical results for Triton X100 were calculated for 𝑘𝑎= 7.25 × 10−4 m∕s and 14.5 × 10−4 m∕s while keeping 𝐿𝑑 constant (𝐿𝑑= 4.4 mm), as explained above. The simulation with 𝑘𝑎= 14.5 × 10−4 m∕s approximately captures the dependency of the bubble velocity on the vertical coordinate. Significant deviations arise in the oscillatory part of the experimental bubble trajectory, which suggests that the oscillatory instability reduces the average vertical velocity. As observed, the bubble trajectory considerably depends on 𝑘𝑎 for a fixed value of the depletion length 𝐿𝑑. This contrasts with the time evolution of the surface coverage when an initially clean spherical bubble remains at rest in a liquid bath, which depends on 𝐿𝑑=𝑘𝑎∕𝑘𝑑 (not on 𝑘𝑎 and 𝑘𝑑 separately) in the diffusion-controlled limit applicable to our case (see the Appendix A). Therefore, the sorption constants enter our problem separately due to surfactant convection caused by the bubble motion. The results for SDS show a good agreement between the simulation and the experiment. For 𝑧 ≳ 80 mm, both the bubble velocity and the aspect ratio take approximately constant values (Fig. 5), indicating that a quasi-steady regime has been reached. In this regime, the bubble velocity is significantly smaller than that in pure water. This difference is accurately captured by the simulation. The overshooting phenomenon is more noticeable in the experiment. This suggests that the formation of the surfactant surface boundary layer is slightly delayed with respect to its numerical counterpart. This delay may occur because the surface diffusion coefficient is smaller than the value considered in the simulation. International Journal of Multiphase Flow 188 (2025) 105205 6 D. Fernández-Martínez et al. Fig. 7. Bubble velocity 𝑣𝑧 as a function of the vertical position 𝑧 of the center of gravity for Triton X-100 (upper graph) and SDS (lower graph). The experiments with Triton X-100 were conducted for 𝑅= 0.76 mm and 𝑐∞= 8.3 × 10−4 mol/m3. The experiments with SDS were obtained for 𝑅= 0.66 mm and 𝑐∞= 8 × 10−4 mol/m3. The solid and dashed lines correspond to the stable and oscillatory parts of the experimental bubble trajectory, respectively. The symbols are the simulation results. The simulation results for Triton X-100 were calculated with 𝑘𝑎= 7.25 × 10−4 m∕s and 14.5 × 10−4 m∕s. In the two cases, 𝐿𝑑= 4.39 mm. The results for SDS were obtained with the values of 𝑘𝑎 and 𝑘𝑑 shown in Table 1. Table 2 Dimensionless numbers defined based on the characteristic quantities without surfactant. 𝜌=𝜌(𝑖)∕𝜌(𝑜)𝜇=𝜇(𝑖)∕𝜇(𝑜)𝐵=𝜌(𝑜)𝑔𝑅2∕𝜎𝑐Ar=𝜌(𝑜)2𝑔𝑅3∕𝜇(𝑜)2 1.20 × 10−3 1.96 × 10−2 5.92 × 10−2 3.24 × 103 4.2. Growth of the dynamic adsorption layer This section analyzes numerically the growth of the dynamic adsorption layer formed when a bubble is released in water containing a surfactant at a very low concentration. Specifically, we consider the same configuration as that studied by Rubio et al. (2024) in the steady regime, a bubble with 𝑅= 0.66 mm in water containing SDS at the concentration 𝑐∞= 8 × 10−4 mol/m3 (see Figs. 5and 7). Tables 2 and 3 show the values of the dimensionless numbers characterizing the problem. In this section, length, time, velocity, flux, and stresses are measured in terms of 𝑅, 𝑡𝑐=𝑅∕𝑣𝑡, 𝑣𝑡, 𝑣𝑡𝑐∞, and 𝜌(𝑜)𝑣2 𝑡, respectively. Fig. 8 shows the surfactant monolayer evolution during the bubble rising. Since the bubble shape changes over time, the polar angle measured from the bubble’s center of gravity also changes, and its use as the independent variable is confusing. For this reason, we have selected instead the surface area 𝑎 measured from the north pole. The surfactant surface concentration 𝛤 sharply increases in the rear of the bubble. The maximum of 𝛤(𝑎) displaces toward the bubble equator until it reaches the approximately constant location 𝑎∕(4𝜋𝑅2)≃0.9, where 𝛤 becomes almost ten times the initial concentration 𝛤0. This large increase in 𝛤 occurs within a very thin diffusive layer, producing a significant surface tension gradient even though this quantity changes in less than 1%. The surface tension gradient gives rise to a Marangoni stress 𝜏Ma three orders of magnitude larger than the tangential viscous stress in a surfactant-free bubble (Rubio et al., 2024). Although the Marangoni stress is confined within the surface boundary layer, it immobilizes part of the bubble’s south hemisphere and significantly reduces the terminal velocity (Fig. 5). The results for the net sorption flux =𝑎−𝑑 show that surfactant adsorbs onto (desorbs from) the interface in front of (behind) the diffusive layer over the entire bubble motion. The hydrostatic pressure distributions over the bubble surface reflect the bubble acceleration. Pressure is built up in front of the bubble as the bubble velocity increases, while a negative gauge pressure arises around the equator, flattening the bubble (Fig. 9). Fig. 8. Surfactant surface concentration 𝛤, surface tension 𝜎, surface velocity 𝑣, net sorption flux =𝑎−𝑑, Marangoni stress 𝜏Ma, and hydrostatic pressure 𝑝 as a function of the area 𝑎. Fig. 9. Evolution of the bubble shape. Fig. 10 shows the evolution of the main quantities describing the bubble dynamics. The symbols indicate the instants corresponding to the profiles in Fig. 8. In clean water, the bubble velocity monotonously increases until it asymptotically reaches its terminal value. However, the velocity in our simulation reaches a maximum at 𝑡∕𝑡𝑐∼ 50 (𝑡∼ 0.1 s) and then slightly decreases (the overshooting phenomenon) (Fig. International Journal of Multiphase Flow 188 (2025) 105205 7 D. Fernández-Martínez et al. Table 3 Dimensionless numbers defined involving the physical properties of the surfactant monolayer. Pe =𝑅𝑣𝑡∕𝑜Pe𝑠=𝑅𝑣𝑡∕𝑠𝐾𝑎=𝑘𝑎𝑐cmc𝑅∕(𝛤∞𝑣𝑡)𝐾𝑑=𝑘𝑑𝑅∕𝑣𝑡Ma=𝛤∞𝑅𝑔𝑇∕𝜎𝑐𝑐∞∕𝑐cmc 2.57 × 1054.89 × 1021.02 1.36 × 10−2 0.132 10−4 Fig. 10. Bubble vertical position 𝑧, velocity 𝑣𝑧, mean surfactant concentration 𝛤, maximum surfactant concentration 𝛤max, normalized net sorption rate 𝛷, bubble surface area 𝐴, and aspect ratio 𝜒 as a function of time 𝑡. The symbols indicate the instants corresponding to the profiles in Fig. 8. The dashed horizontal lines indicate the terminal values. The arrow indicates the maximum velocity due to the overshooting phenomenon. 10b). This decrease is produced by the extra drag resulting from the late growth of the dynamic surfactant layer. For 𝑡 ≳ 0.1 s, there is still net surfactant adsorption at the interface (Fig. 10e). The concentration jump at the surfactant surface boundary layer grows. The Marangoni stress peak becomes higher and wider (Fig. 8e), increasing the normal outer viscous stress that contributes to the drag (Rubio et al., 2024). Additionally, the surfactant surface boundary layer slightly displaces upstream (Fig. 8a), shifting the separation point in the same direction and enlarging the low-pressure area in the wake that also contributes to the drag (Fig. 8f). These two effects are responsible for the bubble deceleration for 𝑡∕𝑡𝑐≳50 (𝑡 ≳ 0.1 s). The mean surfactant concentration 𝛤 monotonously increases (Fig. 10c) due to the continuous net adsorption of surfactant. As explained below, the sharp increase in 𝛤max for 𝑡∕𝑡𝑐≲5 (𝑡 ≲ 0.01 s) (Fig. 10d) reflects the intense surfactant surface convection right after the bubble release. To analyze this aspect of the problem, we calculate the normalized net sorption rate 𝛷=1 4𝜋𝑅2𝛤0 𝑑 𝑑𝑡 (𝐴𝛤 ),(18) where 𝐴 is the bubble surface area. 𝛷 is the time derivative of the surfactant mass trapped in the bubble, 𝐴𝛤 , relative to its initial value. As explained below, it reaches its maximum value at 𝑡∕𝑡𝑐≃ 7.64 (𝑡≃ 0.016 s) (Fig. 10e), when the bubble front (clean) surface has grown. Fig. 10f shows that the bubble surface slightly expands (𝐴 increases) during the bubble acceleration and compresses during the deceleration. This expansion/compression is accompanied by a significant change in the aspect ratio 𝜒 (Fig. 10g), as also observed in Fig. 9. The instantaneous values of the bubble velocity and aspect ratio are correlated, as Kracht and Finch (2010) observed in their experiments. Four time scales describe the formation of the dynamic adsorption layer. The fastest process is convection to the rear part of the bubble of the surfactant adsorbed before the bubble is released. The time scale of this mechanism is 𝑅∕𝑣𝑧𝑐 ∼ 10−2 s, where 𝑣𝑧𝑐 ≃ 0.1 m∕s is the characteristic rising velocity during the bubble acceleration. This first process is followed by that influenced by adsorption–desorption kinetics. The adsorption time scale is 𝑡𝑎= (𝑘𝑎𝑐𝑠∕𝛤)−1. The initial surfactant concentration is approximately that at equilibrium. Then, 𝑐𝑠∕𝛤≃𝑘𝑑∕𝑘𝑎, and, therefore, 𝑡𝑎=𝑘−1 𝑑. We conclude that the adsorption time scale is commensurate with the desorption characteristic time 𝑡𝑑= 𝑘−1 𝑑∼ 0.1 s, which suggests that both mechanisms compete with each other during the same phase of the dynamic adsorption layer formation. The surfactant diffusion governs the last stage of this process. The surface and volumetric diffusion times are 𝑡𝐷𝑠 =𝑅2∕𝑠∼ 1 s and 𝑡𝐷𝑜 =𝑅2∕𝑜∼ 500 s, respectively. The time scales mentioned above can be observed when analyzing the maximum surfactant concentration 𝛤max and mean surfactant concentration 𝛤. The surfactant adsorbed onto the interface at 𝑡= 0 s is swept toward the rear of the bubble very fast. At 𝑡∕𝑡𝑐= 3.64 (𝑡= 0.008 s), when the bubble has traveled a distance smaller than its radius (Fig. 10a), the maximum surfactant concentration 𝛤max has increased over five times (Fig. 10d). However, the mean surfactant concentration 𝛤 has remained almost constant (Fig. 10c). This confirms that surfactant convection dominates the process on the time scale 𝑅∕𝑣𝑧𝑐 ∼ 10−2 s. Now, we analyze the evolution of the surfactant surface distribution during the convection-dominated stage of the process. For this purpose, we have re-plotted in Fig. 11 the profiles shown in Fig. 8, separating them in several graphs and zooming into the rear of the bubble. The maximum surfactant concentration is initially located at the bubble bottom (see Fig. 11-a1 for 𝑡∕𝑡𝑐<6.44 (𝑡 < 14 ms)). The surface tension gradient produces a Marangoni stress 𝜏Ma that significantly slows the surfactant-loaded part of the interface. This stress becomes large enough to reverse the interface motion in the bubble rear at 𝑡∕𝑡𝑐≃ 5.24 (𝑡≃ 11 ms) (Fig. 11-a3). The reverse flow transports the surfactant away from the bubble bottom, making 𝛤max decrease even though there is net adsorption of surfactant, as indicated by the increase in 𝛤 (Fig. 10c). At 𝑡∕𝑡𝑐≃ 8.83 (𝑡≃ 19 ms), the surfactant concentration profile exhibits the shape characteristic of the steady solution. Specifically, more than 80% of the bubble surface is almost clean, and a thin diffusive boundary layer separates this region from the surfactant-loaded rear (Fig. 8a). The Marangoni stress peak in the International Journal of Multiphase Flow 188 (2025) 105205 8 D. Fernández-Martínez et al. Fig. 11. Surfactant surface concentration 𝛤, surface tension 𝜎, surface velocity 𝑣 and net sorption flux =𝑎−𝑑, Marangoni stress 𝜏Ma, and hydrostatic pressure 𝑝 as a function of 𝑎. surface boundary layer practically immobilizes the surfactant-loaded part of the interface (Fig. 8c). The adsorption–desorption of surfactant controls the second phase of the dynamic adsorption layer growth. The flow drags the surfactant toward the rear of the bubble. The uneven surfactant concentration enhances adsorption in the almost-clean region and favors desorption in the loaded rear (Fig. 11-b4). This results in a positive net sorption flux until the bubble reaches its steady velocity, as indicated by the continuous increase in 𝛤 (Fig. 10c). The net sorption rate 𝛷 [Eq. (18)] initially grows (Fig. 10e) as the clean part of the bubble surface enlarges. This occurs until the Marangoni stress peak almost immobilizes the surfactant-loaded region at 𝑡∕𝑡𝑐≃ 6.44 (𝑡≃ 14 ms) (Fig. 11-a3). Then, the bubble deformation increases the surface area (Fig. 10f). This allows for the growth of the surfactant-loaded region. This effect and the increase in the surfactant concentration favor desorption, reducing the net sorption rate. After the maximum velocity and deformation are reached (𝑡∕𝑡𝑐>50, 𝑡 > 0.1 ms), the slight recovery of the spherical shape (Fig. 9) contributes to enlarging the loaded area at the expense of the clean region (Fig. 11c1), which also reduces the net adsorption rate. The process continues until reaching the delicate balance between the surfactant convection, diffusion, and sorption corresponding to the steady regime (Rubio et al., 2024). Fig. 12 shows the volumetric surfactant concentration in the stationary state as a function of the distance 𝑛 normal to the interface at the interface points indicated in the figure. In Fig. 12a, these points are located upstream of the momentum boundary layer separation. In this case, a thin surfactant boundary layer forms with a thickness 𝛿 verifying the scaling law 𝛿∕𝑅∼Pe−1∕2 (Rubio et al., 2024). As can be observed, the numerical results are consistent with the prediction 𝛿∕𝑅∼ 0.002 of that law. The interface points considered in Fig. 12b are within the surfactant wake. Diffusion transports the surfactant molecules away from the interface. 5. Conclusions We analyzed the transient bubble rising in the presence of a surfactant at very low concentrations. The experiments for 𝑅= 0.76 mm and 𝑐∞∕𝑐cmc ≥5×10−4 of Triton X-100 showed that a helical instability eventually develops. The helical instability evolved toward a zig-zag motion as the Triton X-100 concentration increased. For the same relative concentrations 𝑐∞∕𝑐cmc of Triton X-100 and SDS, Triton X-100 produces a smaller effect on a relatively short time scale, but this effect eventually exceeds that produced by SDS at sufficiently large times. The Triton X-100 effects are larger than those caused by SDS at the same absolute concentration. The magnitude of the surfactant effect increases as the bubble radius decreases for both Triton X-100 and SDS. The experiments allowed us to determine the parameter conditions that lead to an axisymmetric, straight (stable) bubble rising. The bubble motion remained stable for 𝑅= 0.66 mm and the concentration 𝑐∞= 8 × 10−4 mol/m3 of SDS. The experimental results satisfactorily agreed with the transient numerical solution for this case. We numerically analyzed the growth of the dynamic adsorption layer of SDS for 𝑅= 0.66 mm and 𝑐∞= 8 × 10−4 mol/m3. The vertical velocity of a bubble in clean water increases monotonously. Conversely, the bubble velocity in our simulation attains a maximum at around 𝑡∕𝑡𝑐≃ 65 (𝑡≃ 0.1 s). This occurs despite the tiny surfactant concentration considered in our analysis, which suggests that the so-called ‘‘overshooting’’ phenomenon occurs for any surfactant concentration. Both the Marangoni stress confined in the boundary layer and the increase in the low-pressure region area are responsible for the bubble deceleration. The bubble surface slightly expands during the bubble acceleration and compresses during the deceleration. A non-monotonous variation of the aspect ratio accompanies this expansion/compression. The amount of surfactant in the monolayer continuously increases over the bubble motion. One can distinguish four consecutive phases in the formation of the dynamic adsorption layer: the almost instantaneous convection of International Journal of Multiphase Flow 188 (2025) 105205 9