The motion of a surfactant-laden bubble in a channel or a Hele-Shaw cell

Abstract

We investigate how the addition of surfactant affects the governing equations for a bubble in a two-dimensional channel in the small-capillary-number limit. In the limit where the time scale for absorption of surfactant is much shorter than the time scales for transport of surfactant along the surface of the bubble, we derive a set of idealised free-surface boundary conditions that capture the effects of surfactant in a single dimensionless ‘elasticity parameter’, and apply them to the front and rear of the bubble separately. At the front of the bubble, there are three regions of interest: the front cap, the thin film region and a transition region that smoothly connects the other two regions. Through matched asymptotic expansions, we derive predictions for the thin film height and the pressure drop across the front meniscus. We find that the viscous pressure drop across the front meniscus is always larger for a surfactant-laden bubble than for a surfactant-free bubble, by an order-one factor of up to $4^{2/3}$. Using a similar analysis at the rear of the bubble, we find that the viscous pressure drop across the rear meniscus is also always larger in magnitude for a surfactant-laden bubble than for a surfactant-free bubble, again up to a maximum factor of $4^{2/3}$. Finally, we use these results to show that, for the same flow conditions, an isolated surfactant-laden bubble in a Hele-Shaw cell will travel more slowly than an isolated surfactant-free bubble.

1. Introduction

Bretherton (1961) first analysed the motion of a surfactant-free bubble through a viscous liquid in a capillary tube. His analysis was modernised by Park & Homsy (1984) into the language of matched asymptotic expansions. We follow a similar methodology in this paper for a surfactant-laden bubble. We then use the results of the analysis to show how the equation of motion derived by Booth, Griffiths & Howell (2023) for an approximately circular bubble in a Hele-Shaw cell is modified by the presence of surfactant, resulting in a model that systematically accounts for the effects of surfactant-laden thin liquid films above and below the bubble, without the inclusion of any ad hoc fitting parameters.

There is a plethora of literature studying the propagation of a bubble or a finger of air in a viscous liquid containing surfactants (see, for example, Ratulowski & Chang 1990; Park 1992; Stebe & Barthes-Biesel 1995; Ghadiali & Gaver 2003; Halpern & Gaver 2012, and references therein). In particular, Ratulowski & Chang (1990) extended the results of Bretherton (1961) to a finger of air propagating into a viscous liquid containing soluble surfactants. They proposed five different distinguished limits based on the convective, diffusive and kinetic time scales. Park (1992) built on the work of Ratulowski & Chang to consider the flow of a long bubble in a capillary tube and found that surfactants can rigidify the bubble’s surface. Maruvada & Park (1996) used these results to describe a rigid elliptical surfactant-laden bubble in a Hele-Shaw cell, generalising the Taylor & Saffman (1959) solution for the motion of a bubble with constant Laplace pressure to include the rigidifying effects of surfactants that can occur in the `convective equilibrium’ regime proposed by Ratulowski & Chang (1990).

In this paper, we focus on what Ratulowski & Chang term the `bulk equilibrium’ model, in which there is an abundance of surfactant in solution, such that the bulk concentration remains approximately constant. Physically, this limit could correspond to the bulk surfactant concentration being significantly above the critical micelle concentration. Furthermore, we consider the regime in which the surfactant is highly soluble, so the surface concentration is close to equilibrium with the bulk. Such a system has been recently studied experimentally by Baué et al. (2025). While the work we present in this paper is focused on the motion of a bubble, the methodology can be adapted to study the propagation of a liquid plug or film coating. This problem thus has applications also to the transport of fluid in the lung (see Waters & Grotberg 2002; Grotberg 2011; Halpern & Gaver 2012; Shemilt et al. 2023), and in fibre coating (see Shen et al. 2002; Delacotte et al. 2012).

Ginley & Radke (1989) also study a surfactant-laden bubble in a capillary tube, while Waters & Grotberg (2002) study the similar problem of the propagation of a surfactant-laden plug. Both of these papers investigate a regime where the surfactant concentration is close to equilibrium, with the resulting Marangoni stress along the bubble interface being negligible at leading order. Both then provide asymptotic results for the pressure drop across the propagating bubble or plug and the height of the thin film deposited behind. We consider a similar regime in this paper; however, we study the distinguished limit in which Marangoni effects enter the problem at leading order. In the limit of small Marangoni stress we recover the asymptotic predictions of Ginley & Radke (1989) and Waters & Grotberg (2002); however, we find that their prediction that the film height decreases with increasing Marangoni stress quickly fails. The thin film height actually starts to increase as the strength of the Marangoni effect is increased, to a maximum factor of $4^{2/3}$ times the Bretherton prediction for a surfactant-free bubble. The factor $4^{2/3}$ has been seen in multiple studies of large Marangoni stress (see, for example, Ratulowski & Chang 1990; Park, Maruvada & Yoon 1994; Stebe & Barthes-Biesel 1995; Shen et al. 2002), and occurs due to the bubble’s surface being stationary in the laboratory frame of reference, rather than satisfying the zero-stress condition used in the original Bretherton problem.

Figure 1. Schematic of a two-dimensional surfactant-laden bubble propagating at speed $\hat {U}_b$ along a channel of height $2\hat {H}$ . We take the origin to be at the start of the thin film region, whose length and height are denoted by $\hat {l}$ and $\hat {h}_\infty$ , respectively. The pressures in the bubble and in the fluid outside are denoted by $\hat {p}_b$ and $\hat {p}$ , respectively.

The basic problem studied in this paper consists of a surfactant-laden bubble propagating along a two-dimensional channel, as sketched in figure 1. For the surfactant-free case, Bretherton (1961) obtained approximate formulae for the height $\hat {h}_\infty$ of the thin film deposited by the front meniscus and the pressure drop $\hat {p}_b-\hat {p}$ across the bubble meniscus, in the limit as the capillary number ${\textit{Ca}}$ tends to zero, namely

(1.1) \begin{equation} \hat {h}_\infty \sim a\hat {H}{\textit{Ca}}^{2/3}, \end{equation}

and

(1.2) \begin{equation} \hat {p}_b - \hat {p} \sim \frac {\hat {\gamma }_0 }{\hat {H}} +\beta \,\frac {\hat {\gamma }_0 }{\hat {H}}\,{\textit{Ca}}^{2/3}. \end{equation}

The dimensionless prefactors $a$ and $\beta$ are determined numerically: $a\approx 1.337$ , while the value of $\beta$ depends on whether the meniscus is advancing or retreating, with

(1.3) \begin{equation} \beta =\begin{cases} \beta _1\approx 3.88\quad &\text{at front (advancing) meniscus}, \\ \beta _2\approx -1.13\quad &\text{at rear (retreating) meniscus}. \end{cases} \end{equation}

Our aim in this paper is to determine how the values of these constants $\{a,\beta _1,\beta _2\}$ are modified for a surfactant-laden bubble, using systematic matched asymptotic expansions. We then apply these results to the propagation of a surfactant-laden bubble through a Hele-Shaw cell and thus derive a drag law that contains no ad hoc fitting parameters, in a similar fashion to Booth et al. (2023).

This paper is structured as follows. We begin in § 2 by writing down the governing equations for the flow around a bubble propagating through a two-dimensional channel containing a surfactant-laden viscous liquid. In § 3 we describe the asymptotic structure of the front region of the bubble, following the methodology of Park & Homsy (1984). We find that there are three regions of interest: the front cap, the thin film region and the transition region, which allows for a smooth transition from the front cap into the thin film region. The equations in the transition region are then analysed in § 4, whereby we find the surfactant-laden analogues of $a$ in (1.1) and $\beta _1$ in (1.3). In § 5 we describe the asymptotic structure of the rear of the bubble, which has an additional matching condition that the height $\hat {h}_\infty$ of the thin film flowing towards the rear of the bubble must equal the height of the thin film deposited by the front meniscus. Then, in § 6, we analyse the rear transition region to find the surfactant-laden analogue of $\beta _2$ in (1.3). Combining the results of §§ 4 and 6 for $\beta _1$ and $\beta _2$ , in § 7 we investigate how the inclusion of surfactant affects the velocity of an isolated bubble in a Hele-Shaw cell. Finally, in § 8 we summarise our key findings.

2. Governing equations

2.1. Dimensional modelling

We consider the steady propagation of a two-dimensional bubble inside a channel of height $2\hat {H}$ (see figure 1). We orient the $\hat {x}$ -axis and the $\hat {z}$ -axis along and perpendicular to the lower wall, respectively. We define $\boldsymbol{i}$ and $\boldsymbol{k}$ as the unit vectors in the $\hat {x}$ - and $\hat {z}$ -directions, respectively. We assume that buoyancy effects are negligible, so the flow is symmetric across the centreline of the channel and we can restrict our attention to the region $0\leqslant \hat {z} \leqslant \hat {H}$ . We define the bubble’s surface in this region to be at $\hat {z}=\hat {h}(\hat {x})$ . The normal to the liquid surface pointing into the bubble and the tangent to the bubble’s surface are denoted by $\boldsymbol{n}$ , and $\boldsymbol{t}$ , respectively, and are given by

(2.1) \begin{align} \boldsymbol{n} & = \dfrac {-\hat {h}_{\hat {x}} \boldsymbol{i}+\boldsymbol{k}}{\sqrt {1+\hat {h}_{\hat {x}}^2}}, \quad \quad \boldsymbol{t} = \dfrac {\boldsymbol{i}+\hat {h}_{\hat {x}}\boldsymbol{k}}{\sqrt {1+\hat {h}_{\hat {x}}^2}}, \end{align}

where the subscript variable means partial differentiation with respect to that variable.

We assume the motion of the bubble is sufficiently slow that the flow is in the Stokes regime and we move into a frame of reference in which the bubble is stationary and the walls travel at a velocity $-\hat {U}_b\boldsymbol{i}$ , where $\hat {U}_b$ is the bubble propagation speed. We denote the liquid velocity as $\hat {\boldsymbol{u}}=(\hat {u},\hat {v})$ and the pressure as $\hat {p}$ . The motion of the liquid is then governed by the Stokes equations

(2.2a) \begin{align} \hat {\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol{u}} &=0, \quad \quad\quad\;\;\; \text{in}\quad \hat {\varOmega }, \end{align}

(2.2b) \begin{align} \hat {\boldsymbol{\nabla \!}} \hat {p} &=\hat {\mu }\hat {{\nabla} }^2\hat {\boldsymbol{u}} \quad\quad\, \text{in}\quad \hat {\varOmega }, \end{align}

where $\hat {\mu }$ is the constant liquid viscosity, $\hat {\varOmega }$ denotes the liquid region and $\hat {\boldsymbol{\nabla }}=(\partial /\partial \hat {x},\partial /\partial \hat {z} )$ is the two-dimensional gradient operator.

The surfactant concentration on the surface of the bubble, $\hat {\varGamma }$ , is governed by the advection–diffusion–reaction equation (Stone 1990)

(2.3) \begin{align} \frac {{\rm d}}{\textrm{d} \hat {s}}\big (\hat {\varGamma \!} u_{\!{S}} \big ) &= \hat {D}\frac {{d}^2\hat {\varGamma }}{\textrm{d}\hat {s}^2} +\hat {k}_1 \hat {C}-\hat {k}_2 \hat {\varGamma }, \quad \quad \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

where $\hat {s}$ is arclength, $\hat {D}$ is the surface diffusion coefficient, $\hat {C}$ is the bulk concentration and

(2.4) \begin{equation} {u \!}_{{S}}=\bigl .\hat {\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{t}\bigr |_{\hat {z}=\hat {h}} \end{equation}

is the tangential surface velocity. We have supplemented (2.3) with a linear reaction term, with rate constants $\hat {k}_1$ and $\hat {k}_2$ , because the surfactant concentrations in the bulk and on the surface are assumed to be close to equilibrium. There are numerous numerical studies that include nonlinear reaction kinetics, such as Fujioka & Grotberg (2005) and Muradoglu et al. (2019). In general one would solve a coupled advection–diffusion equation for $\hat {C}$ but, as mentioned above, we focus on the ‘bulk equilibrium’ limit in which the bulk concentration does not vary significantly, so we take $\hat {C}$ to be a known constant. We further assume that the surfactant quickly adsorbs or desorbs onto the surface of the bubble and thus the surfactant concentration is close to equilibrium. This implies that surfactant cannot accumulate and thus the bubble cannot rigidify, a phenomenon seen in surfactant systems without fast reactions (see e.g. Park 1992).

On the Hele-Shaw cell wall, we supply the no-slip boundary condition

(2.5a) \begin{align} \hat {\boldsymbol{u}} &= -\hat {U}_b\boldsymbol{i} \quad \quad \text{on} \quad \hat {z}=0. \end{align}

On the bubble’s surface, we supply a kinematic condition and normal and tangential stress balances as follows:

(2.5b) \begin{align} \hat {\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{n} &= 0 \quad \quad\quad \quad\quad \quad \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

(2.5c) \begin{align} \boldsymbol{n}\boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}\boldsymbol{\cdot }\boldsymbol{n} &= -\hat {p}_b+\hat {\gamma } \hat {\kappa } \quad \quad\; \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

(2.5d) \begin{align} \boldsymbol{t}\boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}\boldsymbol{\cdot }\boldsymbol{n} &= \frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {s}} \quad \quad\quad \quad\quad \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

where $\hat {\gamma }$ is the (no longer constant) surface tension, $\hat {\kappa }$ is the curvature, $\hat {p}_b$ is the constant pressure inside the bubble and $\hat {\boldsymbol{\sigma }}$ is the viscous stress tensor, given by

(2.5e) \begin{align} \hat {\boldsymbol{\sigma }}=-\hat {p}\boldsymbol{I}+\hat {\mu }\big (\hat {\boldsymbol{\nabla }}\hat {\boldsymbol{u}} +{\hat {\boldsymbol{\nabla }}\hat {\boldsymbol{u}}}^{\text{T}}\big ), \end{align}

where $\boldsymbol{I}$ denotes the two-dimensional identity tensor and the superscript $\text{T}$ denotes the transpose.

We close our model with an equation of state, which relates the surface tension, $\hat {\gamma }$ , to the surfactant surface concentration, $\hat {\varGamma }$ . Since we are assuming that the surfactant concentrations in the bulk and on the surface are close to equilibrium, we also supply a linear equation of state, i.e.

(2.6) \begin{align} \hat {\gamma } = \hat {\gamma }_0 + \frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {\varGamma }}\bigg |_{\hat {\varGamma }=\hat {\varGamma }_0 }(\hat {\varGamma }-\hat {\varGamma }_0), \end{align}

where $\hat {\varGamma }_0 = \hat {k}_1\hat {C}/\hat {k}_2$ is the equilibrium concentration of surfactant and $\hat {\gamma }_0$ is the surface tension at equilibrium.

2.2. Non-dimensionalisation

We non-dimensionalise the system (2.2)–(2.6) as follows (in which dimensionless variables are denoted without hats):

(2.7) \begin{align} (\hat {x},\hat {z},\hat {s},\hat {h}(\hat {x}))&=\hat {H}(x,z,s,h(x)), \quad \quad \hat {\boldsymbol{u}}= \hat {U}_b\boldsymbol{u}, \quad \quad \hat {\gamma }=\hat {\gamma }_0 \gamma , \nonumber \\ (\hat {p},\hat {p}_b)&=\frac {\hat {\gamma }_0}{\hat {H}}(p,p_b), \quad \quad\quad \quad\;\, \hat {\kappa }=\frac {\kappa }{\hat {H}}, \quad \quad \hat {\varGamma } = \hat {\varGamma }_0 \varGamma . \end{align}

The dimensionless governing equations are given by

(2.8) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{u} &= 0, \quad\quad \boldsymbol{\nabla \!}p \quad \quad= {\textit{Ca}} {\nabla} ^2\boldsymbol{u} \quad\quad \text{in }0\lt z\lt h(x). \end{align}

The dimensionless versions of the boundary conditions (2.5) are given by

(2.9a) \begin{align} \boldsymbol{u} &= -\boldsymbol{i} \quad\quad \text{on} \quad z=0, \end{align}

and, on the bubble’s surface $z=h(x)$ ,

(2.9b) \begin{align} \frac {{\rm d}}{\textrm{d} s}\big (\varGamma\! u_{\!{S}}\big ) &= \frac {1}{{\textit{Pe}}}\frac {\rm{d}^2\varGamma }{\textrm{d} s^2} +k (1-\varGamma ),\\[-12pt]\nonumber \end{align}

(2.9c) \begin{align} v&=uh_x,\\[-12pt]\nonumber \end{align}

(2.9d) \begin{align} p_b-p-\frac {2{\textit{Ca}}}{1+h_x^2}\big (\big (1-h_x^2\big )u_x +h_x(u_z +v_x)\big ) &= \frac {(1-{{\textit{Ma}}}(\varGamma -1) )h_{xx}} {\big (1+h_x^2\big )^{3/2}}, \end{align}

(2.9e) \begin{align} \frac {1}{\big (1+h_x^2\big )^{1/2}}\big (-4h_x u_x +\big (1-h_x^2\big )(u_z+v_x)\big ) &= -{\frac {{\textit{Ma}}}{{\textit{Ca}}}}\, \varGamma _x, \end{align}

where

(2.10) \begin{align} &\frac {\rm{d}}{\textrm{d} s} =\frac {1}{\sqrt {1+h_x^2}}\frac {\rm{d}}{\textrm{d} x} \quad\quad \text{and} \quad\quad u_{\!{S}} = \left .\frac {u+h_x\! v}{\sqrt {1+h_x^2}}\right |_{z=h(x)}. \end{align}

Figure 2. Schematic of the front and rear of a bubble, showing the three regions of interest in each: regions 1, front- and rear-cap regions, regions 2, thin film regions and regions 3 transition regions.

The model (2.8)–(2.9) contains the dimensionless parameters

(2.11) \begin{align} {\textit{Ca}} &= \frac {\hat {\mu}\hat {U}_b}{\hat {\gamma }_0}, &{\textit{Pe}} &= \frac {\hat {H}\hat {U}_b}{\hat {D}}, & {{\textit{Ma}}}& {= -\frac {\hat {\varGamma }_0}{\hat {\gamma }_0}\,\frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {\varGamma }},} & k& = \frac {\hat {k}_2 \hat {H}}{\hat {U}_b}, \end{align}

namely the capillary number, the Péclet number, the Marangoni number and the dimensionless reaction constant, respectively. We consider slow flow in which ${\textit{Ca}}\ll 1$ . Furthermore, we assume that ${\textit{Pe}}\gg 1$ , and hence we take the limit ${\textit{Pe}}\rightarrow \infty$ , so (2.9b ) reduces to an advection–reaction equation for $\varGamma$ . Finally, we also consider the regime where $k\gg 1$ so that the surface and bulk surfactants are approximately in equilibrium. We observe from (2.9b ) that $\varGamma \rightarrow 1$ as $k\rightarrow \infty$ and therefore define $\tau =k{\textit{Ca}}^{1/3}(\varGamma -1)$ . Here, the scaling of $\tau$ with ${\textit{Ca}}^{1/3}$ is designed to achieve dominant balances in the transition region between the capillary-static meniscus and the thin film on the channel wall (marked Region 3 in figure 2), as we will show below. We then identify a distinguished limit in which, although the surface concentration of surfactant is almost constant, the Marangoni stress at the free surface is retained in the model at leading order. To this end, we define the dimensionless elasticity parameter

(2.12) \begin{equation} E=\frac {{\textit{Ma}}}{{\textit{Ca}}k}=-\frac {\hat {\varGamma }_0}{\hat {\mu}\hat {H}\hat {k}_2}\,\frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {\varGamma }}, \end{equation}

and take $E=O(1)$ while letting $k\rightarrow \infty$ . The surfactant conservation (2.9b ) and boundary conditions (2.9d )–(2.9e ) thus reduce to

(2.13a) \begin{align} {\textit{Ca}}^{1/3}\frac {\textrm{d}u_{\!{S}}}{\textrm{d}s} &= \tau, \end{align}

(2.13b) \begin{align} p_b-p-\frac {2{\textit{Ca}}}{1+h_x^2}\big (\big (1-h_x^2\big )u_x +h_x(u_z +v_x)\big ) &= \frac {(1-{\textit{ECa}}^{2/3}\tau )h_{xx}}{\big (1+h_x^2\big )^{3/2}}, \end{align}

(2.13c) \begin{align} \frac {{\textit{Ca}}^{1/3}}{\big (1+h_x^2\big )^{1/2}}\big (-4h_x u_x +\big (1-h_x^2\big )(u_z+v_x)\big ) &= -E\tau _x, \end{align}

all on $z=h(x)$ .

Following the above simplifying assumptions, the model now contains just two dimensionless parameters: ${\textit{Ca}}$ , which is assumed to be small, and $E$ , which we take to be $O(1)$ . Waters & Grotberg (2002) study a similar system in which $E\ll 1$ and the Marangoni stress is therefore negligible at leading order. The boundary conditions (2.13) represent the simplest model that contains both the Marangoni stress due to surface gradients in surfactant concentration (2.13c ) and the depletion of adsorbed surfactant by surface dilatation (2.13a ). The resulting system can also be interpreted as a model for surface viscosity (Gounley et al. 2016).

2.3. Summary of assumptions

Here, we summarise the assumptions that have been made in this section.

1. (i) The bulk surfactant concentration $\hat {C}$ is sufficiently large that it does not vary significantly due to interaction with the surface, and we can treat it as effectively constant.

2. (ii) The surfactant is highly soluble, so the surface and bulk surfactant concentrations are approximately in equilibrium, i.e. $k\gg 1$ .

3. (iii) Surface diffusion of surfactant is negligible, so ${\textit{Pe}}\ll 1$ .

Following these assumptions, our model for a highly soluble surfactant is given by (2.8) with the boundary conditions (2.9a ), (2.9c ) and (2.13). In addition, we will assume shortly that the bubble propagates sufficiently slowly for the capillary number ${\textit{Ca}}$ to be small.

3. The front of the bubble

3.1. Regions and scalings

We consider the small- ${\textit{Ca}}$ limit and perform a perturbation expansion in powers of ${\textit{Ca}}^{1/3}$ , following the analysis presented by Park & Homsy (1984), for the similar problem of a surfactant-free bubble. As shown by Park & Homsy, in the small- ${\textit{Ca}}$ limit the problem splits into three regions of interest (see figure 2).

1. (i) The front-cap region. Here, the free surface is capillary static and hence, to leading order, is a circular cap.

2. (ii) The thin film region. Through a lubrication analysis, we find there is a thin liquid film of constant thickness between the channel wall and the bubble.

3. (iii) The transition region. Here, both viscous and capillary forces are important, which allows us to smoothly transition from the circular cap to the thin film region.

We include for completeness the analysis of the front-cap and thin film regions, which follows directly from Park & Homsy. It is in the transition region where the effect of the surfactant is included, and thus the equations deviate from those found by Park & Homsy.

3.2. Region 1: front-cap region

In this region we expand our variables in powers of ${\textit{Ca}}^{1/3}$ : $p(x,z)\sim p_0(x,z) +{\textit{Ca}}^{1/3} p_1(x,z) +{\cdots}$ and so forth. At leading order, the equations of motion (2.8) become

(3.1a,b) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_0 &= 0, \quad\quad \boldsymbol{\nabla \!}p_0 =\boldsymbol{0}, \end{align}

while the boundary conditions (2.13a ) and (2.13c ) are both satisfied identically by $\tau _0=0$ . Equation (3.1b ) implies that $p_0=\text{constant}$ , and the normal stress balances at the bubble’s surface (2.13b ) at leading order reads

(3.2) \begin{align} \Delta p_0 &= \frac {h_0^{\prime\prime}}{(1+h_0^{\prime 2})^{3/2}} \quad\quad \text{on} \quad z=h_0(x), \end{align}

where $\Delta p_0 = p_b - p_0$ is the leading-order difference between the constant pressure inside the bubble, $p_b$ , and the fluid pressure.

To find the leading-order shape of the bubble, we impose the conditions $h_0(x)\rightarrow 1$ and $h_0^{\prime}(x)\rightarrow \infty$ at the front tip of the bubble (which we define to be at $x=0$ ), and $h_0(x_c)=h_0^{\prime}{}(x_c)=0$ , where $x_c\lt 0$ is the a priori unknown location of the point where the leading-order meniscus encounters the cell wall. We thus find that $\Delta p_0=1$ and $x_c=-1$ , and the leading-order shape of the meniscus a circular cap of radius 1 (Park & Homsy 1984), given by

(3.3) \begin{align} h_0(x) &= 1-\sqrt {1-(x+1)^2}, \quad\quad \text{for}\quad x\in (-1,0). \end{align}

In (1.2) we require knowledge of the pressure drop across the meniscus up to $O({\textit{Ca}}^{2/3})$ . Following Park & Homsy (1984), we find that $p_1,\,h_1\equiv 0$ ; then at $O({\textit{Ca}}^{2/3})$ we find that $p_2$ is a constant, denoted by $-\beta _1$ , and

(3.4) \begin{align} \beta _1 = \frac {h_2^{\prime\prime}}{(1+h_0^{\prime 2})^{3/2}}, \end{align}

which can be solved to give

(3.5) \begin{align} h_2(x)=\beta _1\sqrt {1-(x+1)^2}, \end{align}

where $\beta _1$ is the a priori unknown $O({\textit{Ca}}^{2/3})$ pressure correction.

3.3. Region 2: thin film region

In the thin film region, we rescale the variables as

(3.6) \begin{align} x& = \tilde {x}, \quad\quad z= {\textit{Ca}}^{2/3}\tilde {z}, \quad\quad h= {\textit{Ca}}^{2/3}\tilde {h}, \quad\quad p-p_b={\textit{Ca}}^{2/3}\tilde {p},\nonumber \\ u&=\tilde {u}, \quad\quad v= {\textit{Ca}}^{2/3}\tilde {v} \quad\quad\, \tau ={\textit{Ca}}^{1/3} \tilde {\tau }. \end{align}

Again we expand the variables in powers of ${\textit{Ca}}^{1/3}$ . At leading order, the equations of motion (2.8) become

(3.7a–c) \begin{align} \tilde {u}_{0\tilde {x}}+\tilde {v}_{0\tilde {z}} &= 0, \quad\quad \tilde {u}_{0\tilde {z}\tilde {z}} =0, \quad\quad \tilde {p}_{0\tilde {z}} =0. \end{align}

The kinematic boundary conditions (2.9a ) and (2.9c ) become

(3.8a) \begin{align} \tilde {u}_0=-1,\quad \tilde {v}_0&=0 \quad \quad\quad \quad \text{on} \quad \tilde {z}=0, \end{align}

(3.8b) \begin{align} \tilde {v}_0&=\tilde {u}_0\tilde {h}_{0\tilde {x}} \quad \quad\, \text{on} \quad \tilde {z}=\tilde {h}_0(\tilde {x}), \end{align}

and the remaining boundary conditions (2.13a )–(2.13c ) reduce to

(3.9a–c) \begin{align} \tilde {\tau }_0&=-\tilde {u}_{0x}, & \tilde {p}_0&=\tilde {h}_{0xx}, & \tilde {u}_{0\tilde {z}}&=0 \quad \quad \text{on} \quad \tilde {z}=\tilde {h}_0(\tilde {x}). \end{align}

Hence, we find that $\tilde {\tau }_0\equiv 0$ , $\tilde {u}_0\equiv -1$ , $\tilde {p}_0\equiv 0$ and $\tilde {h}_0(\tilde {x}) = \text{constant}$ . Thus, region 2 has a constant (a priori unknown) film thickness

(3.10) \begin{equation} \tilde {h}_0=\frac {\hat {h}_\infty }{\hat {H}\!{\textit{Ca}}^{2/3}}, \end{equation}

where $\hat {h}_\infty$ is the dimensional film thickness indicated in figure 1.

Next, we examine the transition region, which allows us to smoothly transition from the constant film region to the circular cap at the front of the bubble.

3.4. Region 3: transition region

In the transition region, we shift the origin to $x=-1$ and rescale the variables as

(3.11) \begin{align} \begin{aligned} x+1 &= {\textit{Ca}}^{1/3}X, & z&= {\textit{Ca}}^{2/3}Z,& h &= {\textit{Ca}}^{2/3}H, & p=P,\\ u&=U,& v&= {\textit{Ca}}^{1/3}V, &\tau &=G. \end{aligned} \end{align}

We again expand the variables in powers of ${\textit{Ca}}^{1/3}$ .

The leading-order equations of motion (2.8) become

(3.12a–c) \begin{align} U_{0X}+V_{0Z}&=0, & P_{0X}&=U_{0ZZ}, & P_{0Z}&=0 \quad \quad \text{in }\quad 0\lt Z\lt H_0(X). \end{align}

The boundary conditions (2.9a ) and (2.9c ) become

(3.13a) \begin{align} U_0=-1,\quad V_0&=0 \quad \quad\quad \quad\;\; \text{on} \quad Z=0, \end{align}

(3.13b) \begin{align} V_0 &=U_0 H_{0X} \quad \quad \text{on} \quad Z=H_0(X), \end{align}

and the remaining boundary conditions (2.13) become

(3.14a–c) \begin{align} G_0&=-U_{S0}^{\prime}, & -P_0&=H_0^{\prime\prime}, & U_{0Z}&= -{\textit{EG}}_0^{\prime} \quad \quad \text{on} \quad Z=H_0(X), \end{align}

where $U_{S0}(X) = U_0(X,H_0(X))$ .

Using (3.12)–(3.14), we find that $U_0$ is given by

(3.15) \begin{align} U_0(X,Z) =-1 +\frac {1}{2}H_0^{\prime\prime\prime}(X) (2H_0(X) Z-Z^2) - {\textit{EG}}_0^{\prime}(X) Z. \end{align}

By integrating (3.15) across the liquid layer (from $Z=0$ to $Z=H_0(X)$ ) we find that the flux of liquid in the $X$ -direction is given by

(3.16) \begin{align} Q= \frac {1}{3}H_0^{\prime\prime\prime} H_0^3 - \frac {1}{2}{\textit{EG}}_0^{\prime} H_0^2 - H_0. \end{align}

By conservation of mass this quantity must equal the flux of liquid in the thin film region, where $Q=-\tilde {h}_0$ . We thus obtain an equation for the height of the bubble’s surface in the transition region , namely

(3.17) \begin{align} H_0^{\prime\prime\prime} = \frac {3(H_0-\tilde {h}_0)}{H_0^3} +\frac {3{\textit{EG}}_0^{\prime}}{2H_0}. \end{align}

Next, using (3.15) and (3.17), we find that (3.14a ) becomes

(3.18) \begin{align} \frac {{d}}{\textrm{d}X}\left (\frac {1}{4}{\textit{EG}}_0^{\prime}H_0 +\frac {3\tilde {h}_0}{2H_0}\right )&=G_0. \end{align}

Equations (3.17) and (3.18) form a closed system for the film profile, $H_0(X)$ , and the perturbation to the surface concentration of surfactant, $G_0(X)$ , in the transition region. In addition, we enforce the matching conditions

(3.19a) \begin{align} H_0(X) &\rightarrow \tilde {h}_0, \quad \quad\quad\quad\;\;\;\, G_0(X) \rightarrow 0 \quad \quad \text{as}\quad X\rightarrow -\infty , \end{align}

(3.19b) \begin{align} H_0(X)&\sim \frac {1}{2}X^2 +\beta _1, \quad \quad G_0(X) \rightarrow 0 \quad \quad \text{as}\quad X\rightarrow \infty . \end{align}

In principle, the solution of the system (3.17)–(3.18) subject to the far-field behaviour (3.19) determines the a priori unknown constants $\tilde {h}_0$ and $\beta _1$ along with $H_0$ and $G_0$ .

In the next section, we analyse the problem (3.17)–(3.19).

4. Analysis of the transition region equations

4.1. Normalisation

We begin by normalising the equations (3.17) and (3.18) by scaling the variables as

(4.1) \begin{align} &\xi = \frac {X+\mathcal{S}}{\tilde {h}_0}, \quad \quad \eta (\xi ) = \frac {H_0(X)}{\tilde {h}_0}, \quad \quad g(\xi ) = \tilde {h}_0G_0(X). \end{align}

The (3.17) and (3.18) are translation invariant, so we introduce an arbitrary shift $\mathcal{S}$ to simplify the forthcoming analysis. Under these scalings, (3.17) and (3.18) become

(4.2a) \begin{align} \eta ''' &= \frac {3(\eta -1)}{\eta ^3} +\frac {3\mathcal{E}\!g'}{2\eta }, \end{align}

(4.2b) \begin{align} \frac {\mathcal{E}}{4} (\eta g')'&=g+\frac {3\eta '}{2\eta ^2}, \end{align}

where

(4.3) \begin{equation} \mathcal{E}=\frac {E}{\tilde {h}_0}. \end{equation}

The boundary conditions (3.19) imply that

(4.4a) \begin{align} \eta (\xi ) \rightarrow 1,\quad\quad &g(\xi ) \rightarrow 0 \quad\quad \text{as}\quad \xi \rightarrow -\infty , \end{align}

(4.4b) \begin{align} &g(\xi ) \rightarrow 0 \quad\quad \text{as}\quad \xi \rightarrow \infty . \end{align}

The solution of (4.2a ) can be shown to behave quadratically as $\xi \rightarrow \infty$ so

(4.5) \begin{align} \eta (\xi ) \sim \frac {1}{2}a \xi ^2 +b\xi +c \quad\quad \text{as}\quad \xi \rightarrow \infty , \end{align}

where $a$ , $b$ and $c$ are constants. Notice that the coefficients are not uniquely determined due to the arbitrary choice of origin for $\xi$ . However, the translation-invariant groups $a$ and $ac -( {1}/{2})b^2$ are uniquely determined. By comparison with (3.19), we see that they are related to the a priori unknown constants $\tilde {h}_0$ and $\beta _1$ by

(4.6) \begin{align} a&=\tilde {h}_0, \quad\quad ac-\frac {1}{2}b^2 = \beta _1. \end{align}

The solution strategy for the problem (4.2)–(4.4) is explained in the following subsection. Once $\eta$ and $g$ have been computed for a given value of $E$ , the surface velocity of the thin film is calculated from (3.15), giving

(4.7) \begin{align} U_{\!S} = \frac {1}{2}-\frac {\mathcal{E}}{4}\eta g' -\frac {3}{2\eta }. \end{align}

4.2. Solution

We solve (4.2) numerically by shooting from $\xi \rightarrow -\infty$ . Linearising (4.2) about the far-field behaviour (4.4a ), we find that

(4.8) \begin{align} \eta (\xi ) &\sim 1 + \sum _{n=1}^5 A_n \textrm{e}^{\lambda _n \xi },& g(\xi ) &\sim \sum _{n=1}^5 B_n \textrm{e}^{\lambda _n \xi } \quad \quad \text{as} \quad \xi \rightarrow -\infty , \end{align}

where the $\lambda _n$ are roots of the quintic polynomial

(4.9) \begin{align} \lambda ^5 - \frac {4}{\mathcal{E}}\lambda ^3 - 12\lambda ^2 +\frac {12}{\mathcal{E}}=0. \end{align}

This equation has two real and positive roots (which we label $\lambda _1$ and $\lambda _2$ ), one real and negative (labelled $\lambda _3$ ), and a complex conjugate pair with negative real part (labelled $\lambda _c$ and $\overline {\lambda }_c$ ). We require the solution to decay as $\xi \rightarrow -\infty$ , so only the positive eigenvalues are permitted. Hence, the decaying linearised far-field behaviour is given by

(4.10a) \begin{align} \eta (\xi ) &\sim 1 + A_1 \textrm{e}^{\lambda _1 \xi } + A_2 \textrm{e}^{\lambda _2 \xi } \quad \quad \text{as} \quad \xi \rightarrow -\infty , \end{align}

(4.10b) \begin{align} g(\xi ) &\sim \frac {\lambda _1\big (\lambda _1^3-12\big )A_1}{6}\,\textrm{e}^{\lambda _1\xi } +\frac {\lambda _2\big (\lambda _2^3-12\big )A_2}{6}\,\textrm{e}^{\lambda _2\xi } \quad \quad \text{as} \quad \xi \rightarrow -\infty , \end{align}

where $A_1$ and $A_2$ are a priori unknown constants. Due to the translation invariance we may (e.g.) set $A_1 = \pm 1$ by choice of $\mathcal{S}$ in (4.1). (Although translation allows us to set the coefficient of one exponential to have magnitude 1, we do not know its sign in advance; however, we always find that in the front transition region $A_1\gt 0$ .) We then determine $A_2$ via the shooting method to ensure our solution satisfies $g(\xi )\rightarrow 0$ as $\xi \rightarrow \infty$ .

For each value of $\mathcal{E}$ , we use the above shooting method to solve for $\eta$ and $g$ , then read off the coefficients $\{a,b,c\}$ in the quadratic behaviour (4.5) as $\xi \rightarrow \infty$ . We then use (4.3) and (4.6) to determine ${a=}\tilde {h}_0$ and $\beta _1$ parametrically as functions of $E$ . However, the shooting problem can become delicate for small or large values of $\mathcal{E}$ . In the next two subsections, we present asymptotic results for these two limits.

4.3. Small- $\mathcal{E}$ limit

In the limit where $\mathcal{E}$ is small we expand

(4.11) \begin{align} \eta &\sim \eta _{0}+\mathcal{E}\eta _{1}+{\cdots} , \quad \quad g\sim g_{0}+\mathcal{E}g_{1} +{\cdots} . \end{align}

This regime is similar to that studied by Waters & Grotberg (2002) for a surfactant-laden liquid plug and by Ginley & Radke (1989) who considered a bubble in a capillary tube.

At $O(1)$ in (4.2), we find that

(4.12a) \begin{align} \eta _0^{\prime\prime\prime} &= \frac {3(\eta _{0}-1)}{\eta _{0}^3}, \end{align}

(4.12b) \begin{align} g_{0} &=-\frac {3\eta _{0}^{\prime}}{2\eta _{0}^2}. \end{align}

The decoupled (4.12a ) for $\eta _{0}$ is the same Landau–Levich equation used by Bretherton (1961) to determine the shape of a surfactant-free bubble in the transition region. The solution for $\eta _0$ is uniquely determined, up to an arbitrary translation, and the corresponding leading-order surfactant concentration profile $g_0(\xi )$ , given by (4.12b ), is plotted in figure 3. Although the limit $\mathcal{E}\rightarrow 0$ is singular, removing the highest derivative in (4.2b ), we see that $g_0$ tends to zero in the far field, as required, and no boundary-layer behaviour is produced, as also found by Ginley & Radke (1989) and Waters & Grotberg (2002). The coefficients in the quadratic behaviour

(4.13a) \begin{align} \eta _{0}(\xi )&\sim \frac {1}{2}a_0 \xi ^2 +b_0 \xi +c_0 \quad \quad \text{as}\quad \xi \rightarrow \infty , \end{align}

are also determined uniquely and, in particular, we have $a_0\approx 1.337$ and $a_0c_0-({1}/{2})b_0^2\approx 3.88$ , as found by Bretherton (1961).

Figure 3. The surfactant concentration in the front transition region, $g(\xi )$ , with $\mathcal{E}\rightarrow 0$ (black), $\mathcal{E}=1$ (blue), $\mathcal{E}=4$ (red) and $\mathcal{E}\rightarrow \infty$ (purple).

To find the correction to $\tilde {h}_0$ and $\beta _1$ due to the effect of surfactants, we proceed to first order in (4.2a ) to obtain the equation

(4.14) \begin{align} \eta _{1}^{\prime\prime\prime} = \frac {12(3-2\eta _{0})\eta _{1} + 9(2\eta _{0}^{\prime 2}-\eta _{0}\eta _0^{\prime\prime})}{4\eta _{0}^4}, \end{align}

for the correction to the thin film height. We solve (4.14) in the same fashion as (4.12a ) by shooting from $\xi \rightarrow -\infty$ . Again the solution to (4.14) behaves quadratically for large positive $\xi$

(4.15a) \begin{align} \eta _{01} &\sim \frac {1}{2}a_1 \xi ^2 +b_1 \xi +c_1 \quad \quad \text{as}\quad \xi \rightarrow \infty , \end{align}

where the constants $a_1,\,b_1,\,c_1$ are in principle determined (up to an arbitrary translation) by the solution of (4.14). In particular we find that $a_1\approx -0.0146$ and $a_0c_1+c_0a_1-b_0b_1\approx 0.58$ .

Finally, we obtain the small- $E$ expansions for ${a=}\tilde {h}_0$ and $\beta _1$ from (4.6). Note that the definition (4.3) of $\mathcal{E}$ involves $\tilde {h}_0$ , so we have to manipulate the expansions to remove the dependence on $\tilde {h}_0$ to get

(4.16a) \begin{align} {a}&\sim a_0+\frac {a_1}{a_0}\,E\approx 1.337-0.011E&&\text{as}\quad E\rightarrow 0, \end{align}

(4.16b) \begin{align} \beta _1&\sim a_0c_0-\frac {1}{2}b_0^2+\frac {a_0c_1+c_0a_1-b_0b_1}{a_0}\,E\approx 3.88+0.43E&&\text{as}\quad E\rightarrow 0. \end{align}

We note that (4.16a ) and (4.16b ) differ from (26) in Waters & Grotberg (2002) because they include a factor of 3 in their ${\textit{Ca}}$ , which induces a factor of $3^{1/3}$ in the definition of $E$ . We also note that Waters & Grotberg’s expression for the pressure drop is twice ours, due to the cylindrical instead of two-dimensional geometry that they study.

4.4. Large- $\mathcal{E}$ limit

In the other extreme where $\mathcal{E}$ is large, we expand

(4.17) \begin{align} \eta &\sim \eta _{0}+\frac {1}{\mathcal{E}}\,\eta _{1}+{\cdots} , \quad \quad g\sim \frac {1}{\mathcal{E}}\,g_{1} +{\cdots} . \end{align}

At $O(1)$ in (4.2) we find

(4.18a) \begin{align} \eta _0^{\prime\prime\prime} &= \frac {3(\eta _{0}-1)}{\eta _{0}^3}+\frac {3g_{1}^{\prime}}{2\eta _{0}}, \end{align}

(4.18b) \begin{align} \frac {1}{4}(g_{1}^{\prime}\eta _{0})' &= \frac {3\eta _{0}^{\prime}}{2\eta _{0}^2}. \end{align}

We can integrate (4.18b ) and substitute into (4.18a ) to obtain

(4.19a) \begin{align} \eta _0^{\prime\prime\prime}&=\frac {12(\eta _{0}-1)}{\eta _{0}^3}, \end{align}

(4.19b) \begin{align} g_{1}&=\frac {1}{2}\,\eta _{0}\eta _0^{\prime\prime}-\frac {1}{4} (\eta _{0}^{\prime})^2-\frac {1}{2}\,\beta _{10}, \end{align}

where $\beta _{10}$ is the leading-order approximation for the coefficient $\beta _1$ .

Once again, (4.19a ) is the Landau–Levich equation and it is similar to the surfactant-free equation (4.12a ) found by Bretherton (1961) except with an additional factor of 4 in the numerator. This additional factor of 4 induces an increase in the thin film height and correction to the pressure drop by a factor of $4^{2/3}$ , i.e.

(4.20a) \begin{align} {a} &\rightarrow 4^{2/3}\boldsymbol{\cdot }1.337\approx 3.369 \quad \quad\quad \;\;\, \text{as} \quad E\rightarrow \infty , \end{align}

(4.20b) \begin{align} \beta _1 &\rightarrow \beta _{10}\approx 4^{2/3}\boldsymbol{\cdot }3.88\approx 9.78 \quad \quad \text{as} \quad E\rightarrow \infty . \end{align}

In this limit, we find that the surface velocity (4.7) is given by $U_{\!S}\equiv -1$ , which corresponds to the bubble interface travelling at the same velocity as the walls of the cell.

These results reproduce the large-Marangoni-number limit reported in previous studies (see, for example, Ratulowski & Chang 1990; Park 1992; Stebe & Barthes-Biesel 1995; Shen et al. 2002). However, we also evaluate the correction to the surfactant concentration, given by (4.19b ) and plotted in figure 5, where it is evident that $g_1$ does not satisfy the far-field condition $g_1(\xi )\rightarrow 0$ as $\xi \rightarrow -\infty$ . This apparent inconsistency can be resolved by examining an outer region in which

(4.21) \begin{align} \xi &=\mathcal{E}^{1/2}\varXi , & g(\xi )&=\mathcal{E}^{-1}\varPsi (\varXi ), & \eta (\xi )&\sim 1+\text{exponentially small terms}, \end{align}

so that (4.2b ) is transformed to

(4.22) \begin{equation} \varPsi ''=4\varPsi , \end{equation}

up to exponentially small corrections. By matching with (4.19b ) we thus obtain the leading-order outer solution

(4.23) \begin{equation} \varPsi _0(\varXi )=-\frac {1}{2}\,\beta _{10}\textrm{e}^{2\varXi }. \end{equation}

4.5. Results

In figure 3, we plot the correction from equilibrium to the surfactant concentration, $g$ , in the front transition region when $\mathcal{E}\rightarrow 0$ , $\mathcal{E}=1$ , $\mathcal{E}=4$ and $\mathcal{E}\rightarrow \infty$ . We use the arbitrary shift $\mathcal{S}$ introduced in § 4.1 to align the peaks of the concentration profiles. In the limit $\mathcal{E}\rightarrow \infty$ , $g$ vanishes across the entire domain, but in all other cases, we observe that $g\lt 0$ and so the surfactant concentration is everywhere below equilibrium in the front transition region. Similar concentration profiles were observed by Stebe & Barthes-Biesel (1995) in a system with an elevated bulk concentration. In figure 4 we plot the film height in the transition region and, for all values of $\mathcal{E}$ , we observe similar profiles to those found by Bretherton (1961) and Park & Homsy (1984) for a surfactant-free bubble.

Figure 4. The free-surface profile $\eta (\xi )$ in the front transition region, with $\mathcal{E}\rightarrow 0$ (black), $\mathcal{E}=1$ (blue), $\mathcal{E}=4$ (red).

Figure 5. The leading-order (a) perturbation to the surfactant concentration, $g_1(\xi )$ , and (b) free-surface profile, $\eta _0(\xi )$ , in the front transition region in the limit $\mathcal{E}\rightarrow \infty$ .

We plot the lowest-order perturbation to the surfactant concentration, $g_1(\xi )$ , and the film height, $\eta _0(\xi )$ , in the limit $\mathcal{E}\rightarrow \infty$ in figure 5. The leading-order solution evidently does not satisfy the downstream boundary condition $g_1(\xi )\rightarrow 0$ as $\xi \rightarrow -\infty$ , implying that there must be a boundary layer at infinity, as explained in § 4.4. In figure 6 we plot the leading-order surface velocity, $U_{\!S}$ , for $\mathcal{E}\rightarrow 0$ , $\mathcal{E}=1$ , $\mathcal{E}=4$ and $\mathcal{E}\rightarrow \infty$ . We observe that, for finite $\mathcal{E}$ , there is a stagnation point (in the frame of the bubble) within the transition region. Its location is close to the minimum point of $g$ in figure 3, because the flow directed outwards advects surfactant away from the stagnation point. The presence of a stagnation point along the front of the bubble is a prevalent feature of gas bubbles in Hele-Shaw cells or capillary tubes, even in systems with more complicated surfactant dynamics and non-zero Reynolds numbers (Fujioka & Grotberg 2005; Zheng, Fujioka & Grotberg 2007).

Figure 6. The surface velocity, $U_S(\xi )$ , in the front transition region with $\mathcal{E}\rightarrow 0$ (black), $\mathcal{E}=1$ (blue), $\mathcal{E}=4$ (red) and $\mathcal{E}\rightarrow \infty$ (purple).

Figure 7. The normalised thin film height $a$ versus elasticity parameter $E$ . The solid curve is from the numerical solution of (4.2), and dashed curves are the asymptotic predictions: (4.16a ) for small $E$ and (4.20a ) for large $E$ . (a) A log–linear plot to show the full range of $E$ . (b) The solution for $0\leqslant E\leqslant 5.77$ .

The normalised height of the thin film, $a$ , is plotted as a function of the elasticity parameter $E$ in figure 7. As $E\rightarrow 0$ , surfactant effects become negligible and the thin film height approaches Bretherton’s result of $1.337$ for a surfactant-free bubble (Bretherton 1961). At the other extreme, when $E$ is large, $a$ approaches 3.369, which is larger by a factor of $4^{2/3}$ , as expected. Interestingly, (4.16a ) predicts a decrease in the thin film height for small $E$ (see figure 7 b), however, figure 7(b) shows that the asymptotic result (4.16a ) (also obtained by Waters & Grotberg 2002) quickly becomes redundant, and the normalised thin film height, $a$ , increases with $E$ thereafter.

The correction to the pressure drop across the front meniscus, $\beta _1$ , is plotted as a function of $E$ in figure 8. Again when $E$ is small we recover the Bretherton (1961) result that $\beta _1\approx 3.88$ . We observe that $\beta _1$ is a monotonic increasing function and when $E$ is large $\beta _1$ approaches 9.78, in agreement with (4.20b ). In the numerical simulation of (4.2), accurate convergence for the value of $\beta _1$ could not be achieved for values of $E\lesssim 0.{5}$ due to the sensitivity of the numerical shooting method, caused by the singular nature of the system (4.2) as $E\rightarrow 0$ . In general, it is harder to compute the value of $\beta _1$ than $a$ at small $E$ because a significantly larger value of $\xi$ is needed to robustly extract the value of $ac-b^2/2$ from the quadratic function (4.5) than to determine $a$ . The numerical approach is thus useful provided $E\gtrsim O(1)$ , while the asymptotic approximation (4.16b ) is useful when $E$ is small, and we are reassured by figure 8 that there is at least a small overlap region where they approximately agree.

Figure 8. The correction to the pressure drop across the front meniscus, $\beta _1$ , versus elasticity parameter $E$ . The solid curve is from the numerical solution of (4.2), and dashed curves are the asymptotic predictions: (4.16b ) for small $E$ and (4.20b ) for large $E$ .

5. Rear of the bubble

5.1. Regions

As for the front meniscus, in the small- ${\textit{Ca}}$ limit the problem at the rear of the bubble splits into three regions of interest (see figure 2). In particular, for the rear-cap region, we can follow the same analysis as in § 3.2 and find that the leading-order shape and surfactant concentration are given by

(5.1a) \begin{align} h_0(x) &= 1-\sqrt {1-(l-1 +x)^2}, \end{align}

(5.1b) \begin{align} h_2(x) &= \beta _2\sqrt {(x+l)(2-l -x)}, \end{align}

(5.1c) \begin{align} \varGamma _0(x)&=1, \end{align}

where $l$ is the dimensionless length of the bubble, and $\beta _2$ is the a priori unknown $O({\textit{Ca}}^{2/3})$ correction to the pressure drop across the rear meniscus.

As in § 3.3, in the thin film region (region 2), the film height $\tilde {h}_0$ is constant, and now in principle known from the solution for the front meniscus (see figure 7). In the next section we analyse the equation in the rear transition region in a similar manner to § 4.

6. Analysis of the rear transition region equations

6.1. Normalisation

We again normalise by scaling the variables as

(6.1) \begin{align} &\xi = \frac {X+\mathcal{S}}{\tilde {h}_0}, \quad\quad \eta (\xi ) = \frac {H_0(X)}{\tilde {h}_0}, \quad \quad g(\xi ) = \tilde {h}_0 G_0(X), \end{align}

where $\mathcal{S}$ is an arbitrary shift of our coordinates. We obtain exactly the same (4.2) as the front transition region, i.e.

(6.2a) \begin{align} \eta '''&= \frac {3(\eta -1)}{\eta ^3} + \frac {3\mathcal{E}g'}{2\eta }, \end{align}

(6.2b) \begin{align} \frac {\mathcal{E}}{4}(g'\eta )'&=g+\frac {3\eta '}{2\eta ^2}. \end{align}

We solve (6.2) numerically now by shooting from $\xi \rightarrow \infty$ . We find that the decaying linearised solution is given by

(6.3a) \begin{align} \eta &\sim 1\pm \textrm{e}^{\lambda _3 \xi } + S \textrm{e}^{\lambda _R\xi }\cos (\lambda _I \xi -q), \end{align}

(6.3b) \begin{align} g&\sim \pm \frac {6 \lambda _3}{\mathcal{E} \lambda _3^2 -4}\textrm{e}^{\lambda _3 \xi } -S \! \varLambda _c(\lambda _R,\lambda _I) \textrm{e}^{\lambda _R\xi }\cos (\lambda _I \xi -q)-S\! \varLambda _s(\lambda _R,\lambda _I) \textrm{e}^{\lambda _R\xi }\sin (\lambda _I \xi -q), \end{align}

as $\xi \rightarrow \infty$ . Here, $S$ and $q$ are a priori unknown shooting parameters, $\lambda _3$ is the real negative solution of (4.9) and $\lambda _c= \lambda _R +\textrm{i}\lambda _I$ is the complex root with negative real part. The coefficients are given by

(6.4a) \begin{align} \varLambda _c(\lambda _R,\lambda _I) = \frac {6 \lambda _R \left(-4+\mathcal{E}\left(\lambda _R^2+\lambda _I^2\right)\right)}{16+8\mathcal{E}\left(\lambda _I^2-\lambda _R^2\right)+\mathcal{E}^2 \left(\lambda _R^2+\lambda _I^2\right)^2}, \end{align}

(6.4b) \begin{align} \varLambda _s(\lambda _R,\lambda _I) = \frac {6 \lambda _I \left(4+\mathcal{E}\left(\lambda _R^2+\lambda _I^2\right)\right)}{16+8\mathcal{E}\left(\lambda _I^2-\lambda _R^2\right)+\mathcal{E}^2 \left(\lambda _R^2+\lambda _I^2\right)^2}. \end{align}

Note again that the $\pm$ occurs in (6.3) because, although translation allows us to set the coefficient of the exponential to be of magnitude 1, we do not know its sign in advance. Finally, we now have that $\eta _0$ behaves quadratically for large negative $\xi$ , i.e.

(6.5) \begin{align} \eta &\sim \frac {1}{2}A \xi ^2 +B\xi +C \quad \quad \text{as}\quad \xi \rightarrow -\infty . \end{align}

We solve (6.2) for each value of $\mathcal{E}=E/\tilde {h}_0(E)$ by applying a shooting method with two unknown parameters $S$ and $q$ , which are fixed by ensuring $g(\xi ) \rightarrow 0$ and $\eta ''(\xi )\rightarrow \tilde {h}_0(E)$ as $\xi \rightarrow -\infty$ , where $\tilde {h}_0(E)$ is as shown in figure 7. The first condition corresponds to matching the surfactant concentration in the thin film to the equilibrium concentration in the rear cap (see § 5.1), and the second ensures that the thin film height at the rear meniscus matches the height of the thin film deposited at the front meniscus. Following the matching procedure laid out in § 4.1, we then obtain the $O({\textit{Ca}}^{2/3})$ correction to the pressure drop across the rear meniscus as

(6.6) \begin{align} \beta _2 = AC-\frac {1}{2}B^2. \end{align}

This two-parameter shooting problem can be extremely sensitive, so we examine the limiting cases using asymptotic analysis.

6.2. Small- $\mathcal{E}$ limit

In the extreme where $\mathcal{E}$ is small we expand

(6.7a) \begin{align} \eta &\sim \eta _{0}+\mathcal{E}\eta _{1}+{\cdots} , \end{align}

(6.7b) \begin{align} g&\sim g_{0}+\mathcal{E}g_{1} +{\cdots} . \end{align}

Then at $O(1)$ in (6.2) we again find that the equations reduce to (4.12). We note that (6.2b ) is singular in the limit $\mathcal{E}\rightarrow 0$ ; however, for the same reasons as presented in § 4.3 there is no boundary-layer behaviour and the solution of (4.12b ) satisfies all the relevant boundary conditions.

Again, the (4.12a ) for $\eta _{0}$ decouples and is just the usual Landau–Levich equation obtained for a surfactant-free bubble. At first order in (6.2) we again obtain (4.14) for the correction to the bubble’s surface. We solve (4.14) by shooting from $\xi \rightarrow \infty$ . Following the matching methodology laid out in § 4.1, we thus find that the $O({\textit{Ca}}^{2/3})$ correction to the pressure drop is given by

(6.8) \begin{align} \beta _2&\sim -1.13-0.73E \quad \quad \text{as}\quad E\rightarrow 0. \end{align}

The leading term in (6.8) is Bretherton’s classical result for the rear meniscus of a surfactant-free bubble (Bretherton 1961), and the second term is the first correction due to the presence of surfactant.

6.3. Large- $E$ limit

In this limit we follow the same methodology as in § 4.4 to obtain

(6.9) \begin{align} \eta _0^{\prime\prime\prime}=\frac {12(\eta _{0}-1)}{\eta _{0}^3}. \end{align}

Matching with the rear-cap solutions (5.1) we find that the correction to the pressure drop is then given by

(6.10) \begin{align} \beta _2 \sim -4^{2/3}\boldsymbol{\cdot }1.13 \approx -2.85 \quad \quad \text{as}\quad E\rightarrow \infty . \end{align}

Again this is a factor of $4^{2/3}$ larger than the original Bretherton (1961) result. This extends the large-Marangoni-number limit reported in many studies (see, for example, Ratulowski & Chang 1990; Park 1992) to the rear meniscus.

6.4. Results

In figure 9 we show example solutions for $E=1.36$ ( $\mathcal{E}=1)$ and $E=5.76$ ( $\mathcal{E}=4)$ . In figure 9(b) we observe that the surfactant concentration can be both above and below the equilibrium concentration in the rear transition region, in contrast to the front transition region, where the concentration is always below equilibrium. For these specific solutions, we find that $\beta _2 \approx -1.50$ for $E=1$ and $\beta _2\approx -2.33$ for $E=4$ , which are greater in magnitude than the pressure drop $\beta _2\approx -1.13$ for a surfactant-free bubble found by Bretherton (1961). In figure 9(c) we plot the corresponding surface velocities, $U_{\!S}$ for the same values of $E$ . Similarly to the front meniscus, we observe a stagnation point (in the frame of the bubble) in the transition region. However, here the flow is directed into the stagnation point, resulting in a local increase in the surfactant concentration.

Figure 9. (a) The surface profile, $\eta$ , (b) the surfactant concentration, $g_1$ , and (c) the surface velocity, $U_{\!S}$ , in the rear transition region with $\mathcal{E}=1$ , (black) $\mathcal{E}=4$ (red).

7. Application to the motion of bubbles in a Hele-Shaw cell

7.1. Force balance

We are now in a position to include the effect of surfactants in the models presented by Booth et al. (2023, 2025a , b ) and Wu et al. (2024) for the motion of an approximately circular bubble in a Hele-Shaw cell moving due to a uniform background flow $\hat {U}_{\!{f}}\boldsymbol{i}$ . Booth et al. (2023) find that the dimensionless velocity $\boldsymbol{U \!}_b$ of such a bubble is determined by the force balance

(7.1) \begin{equation} \frac {\boldsymbol{U \!}_b}{|\boldsymbol{U \!}_b|^{1/3}} =\frac {\delta }{\pi }\oint _{\partial \varOmega } -p\boldsymbol{n}\, \textrm{d}s, \end{equation}

where $\partial \varOmega$ is the bubble’s surface as viewed from above (see figure 10) and the Bretherton parameter is defined by

(7.2) \begin{equation} \delta = \frac {3\sqrt {\pi }\varGamma (11/6)}{(\beta _1 - \beta _2)\varGamma (4/3)}\, \frac {{\textit{Ca}}_{\!{f}}^{1/3}}{\epsilon }. \end{equation}

Here, $\epsilon = \hat {H}/\hat {R}$ , where $\hat {R}$ is the radius of the bubble (measured from above), and ${\textit{Ca}}_{\!{f}} = \hat {\mu }\hat {U_f}/\hat {\gamma }$ is the capillary number based on the background flow speed, $\hat {U}_{\!{f}}$ , both of which are assumed to be small. In the distinguished limit where ${\textit{Ca}}_{\!{f}} = O(\epsilon ^3)$ as $\epsilon \rightarrow 0$ , so the viscous lubrication pressure balances the pressure drop across the menisci, the bubble is circular to leading order (Booth et al. 2023). For a surfactant-free bubble, $\beta _1$ and $\beta _2$ are given by the values $\beta _1(0)\approx 3.88$ and $\beta _2(0)\approx -1.13$ originally calculated by Bretherton. This result is now easily generalised for a surfactant-laden bubble by using the expressions for $\beta _1(E)$ and $\beta _2(E)$ found in §§ 4 and 6, respectively. Crucially, we recall that the elasticity parameter $E$ , given by (2.11), is independent of the capillary number, so the values of $\beta _1$ and $\beta _2$ depend only on the given surfactant properties and concentration. Note that the values of $\beta _{1}$ and $\beta _2$ are both $O(1)$ for the entire range of values of $E\in [0,\infty ).$

Figure 10. Plan view of a surfactant-laden bubble in a Hele-Shaw cell in a uniform background flow.

From (7.1), we find that the velocity of an isolated bubble in a uniform background flow is given by $\boldsymbol{U \!}_b = U_{\!{b}}\boldsymbol{i}$ , where

(7.3) \begin{align} \frac {U_b^{2/3}}{2-U_b} = \delta = \left (\frac {\beta _1(0)-\beta _2(0)}{\beta _1(E)-\beta _2(E)}\right )\delta _B, \end{align}

and we define the surfactant-free Bretherton parameter

(7.4) \begin{equation} \delta _B = \frac {3\sqrt {\pi }\varGamma (11/6)}{(\beta _1(0)-\beta _2(0))\varGamma (4/3)}\, \frac {{\textit{Ca}}_{\!{f}}^{1/3}}{\epsilon } \approx 1.12 \frac {{\textit{Ca}}_{\!{f}}^{1/3}}{\epsilon }. \end{equation}

7.2. Results

In figure 11 we plot $U_b$ versus $\delta _B$ , for a range of values of $E$ . Note that, if we plotted versus the Bretherton parameter, $\delta$ given by (7.2), then all the curves would collapse. Plotting $U_b$ versus $\delta _B$ allows us to analyse the effect of surfactant on the bubble velocity in comparison with a surfactant-free bubble experiencing the same flow conditions. We observe that the velocity of a surfactant-laden bubble ( $E\gt 0)$ at each $\delta _B$ is less than that of a surfactant-free bubble ( $E=0$ ) at the same value of $\delta _B$ . This trend continues as we increase $E$ , up to the limiting case $E\rightarrow \infty$ when $\delta = 4^{2/3}\delta _B$ , the maximum value that $\delta$ can take for a fixed $\delta _B$ . Hence, we always find that a surfactant-laden bubble travels more slowly than a surfactant-free bubble under the same flow conditions.

The form of (7.3) has the same structure as the expression found by Baué et al. (2025) (their (5.11)) for the velocity of a droplet in a highly soluble surfactant solution, in the limit as the droplet viscosity tends to zero. Note that their capillary number is calculated from the droplet velocity, whereas we use the capillary number ${\textit{Ca}}_{\!{f}}$ based on the background flow speed. One can make an analogy between their constant $1/K$ and the prefactor of $\delta _B$ in (7.4). Our model provides the dependence of this prefactor on the surfactant properties, which is missing in their work. Baué et al. (2025) found experimentally that, with highly soluble surfactants, the velocity of a bubble increases with its size, which is consistent with (7.3).

Figure 11. The dimensionless bubble velocity, $U_b$ (7.3) as a function of the surfactant-free Bretherton parameter, $\delta _B$ (7.4) for a range of values of $E=0$ (black), $E=1.36$ ( $\mathcal{E}=1)$ (blue), $E=5.76$ ( $\mathcal{E}=4)$ (red), $E=\infty$ (purple), with $\beta _1$ and $\beta _2$ given by (4.6) and (6.6), respectively.

8. Conclusions

In this paper, we develop a model for the propagation of a two-dimensional surfactant-laden bubble in a channel. We adopt the so-called bulk equilibrium model, in which there is assumed to be an abundance of surfactant in the liquid. We then identify a distinguished asymptotic limit in which the reaction kinetics are so fast that the surface concentration of surfactant remains close to equilibrium, but the Marangoni stress is large enough still to enter the model at leading order. The resulting boundary conditions (2.13) capture the important physical effects of surfactant in a single dimensionless parameter $E$ .

Through the method of matched asymptotic expansions, we derive results for the dimensionless height of the thin films between the bubble and the channel walls and for the corrections to the pressure drop across the front and rear menisci of the bubble. Such an analysis is reliant on the bubble being long, so we can treat the front and rear of the bubble separately. Our bulk equilibrium surfactant model produces results analogous to Bretherton’s, in which the thin film height and the pressure corrections scale with ${\textit{Ca}}^{2/3}$ (Bretherton 1961), but where the prefactors are now numerically determined functions of $E$ , with the surfactant-free case corresponding to $E=0$ . Previous work (Waters & Grotberg 2002) found that the height of the deposited film is a decreasing function of $E$ in the limit $E\rightarrow 0$ . Strikingly, we show that this asymptotic prediction fails for $E$ as small as $0.2$ , and in fact the film height almost always increases with $E$ , up to a maximum value larger than Bretherton’s by a factor of $4^{2/3}$ . Likewise, we find that the net pressure difference across both menisci increases with $E$ , again by a factor of up to $4^{2/3}$ in the limit as $E\rightarrow \infty$ . The factor of $4^{2/3}$ comes from the bubble’s surface being stationary in the laboratory frame of reference, rather than satisfying the zero-stress condition as in the original Bretherton problem.

The key outputs from our analysis are the normalised corrections to the pressure drop $\beta _1$ and $\beta _2$ across the front and rear menisci, respectively. In practice, the computation of these parameters across a range of values of $E$ is very challenging because of the extreme sensitivity of the relevant shooting problems, especially for the rear meniscus, where there are two shooting parameters. To perform an exhaustive parameter sweep, particularly in the singular limit where $E\rightarrow 0$ , it may be necessary to adopt an alternative numerical approach, for example solving the boundary-value problem directly by discretising the whole domain.

We use our results for the modified pressure drop across the bubble to obtain a generalised equation of motion for a bubble in a Hele-Shaw cell that includes the effects of surfactants. As in Booth et al. (2023), the effective viscous drag on the bubble is measured by a dimensionless `Bretherton parameter’, $\delta \propto {\textit{Ca}}^{1/3}/\epsilon$ , with just the prefactor now a function of $E$ (see (7.2)). We find that, for the same flow conditions, an isolated surfactant-laden bubble will travel more slowly than an isolated surfactant-free bubble. Crucially, $E$ depends only on the physical properties of the surfactant and the fluid, as well the cell height, but not on any local flow properties (e.g. the local capillary number). The model thus easily generalises to an arbitrary number of bubbles by modifying the prefactor in $\delta$ in the same way for each bubble.

Our modelling relies on the surfactant being highly soluble, in the sense that the time scale for adsorption is much shorter than that for surfactant transport, i.e. the Damköhler number is large. It also relies on the capillary number being small, i.e. the bubble propagates slowly enough for the free surface to be dominated by capillary effects. The former can be achieved in practice using surfactants such as sodium alkyl sulphates, or alkyl trimethylammonium bromides with fewer than 11 carbons in the alkyl chain (Baué et al. 2025). The latter is almost always satisfied in microfluidic devices (Stone, Stroock & Ajdari 2004).

Our analysis relies on the front and rear menisci being well separated in the flow direction, which is not true near the `poles’ of the bubble (in plane view), where the bubble meniscus is parallel to the background velocity. In these regions, a different asymptotic scaling allows one to explain how the parameter $\beta$ varies smoothly between the constant values $\beta _1$ and $\beta _2$ , as shown by Burgess & Foster (1990) for clean bubbles. However, Booth et al. (2023) showed that these regions provide a correction to the bubble velocity that is $O(\epsilon ^{6/5})$ , and thus negligible to lowest order.