Drag reduction in surfactant-contaminated superhydrophobic channels at high Péclet numbers

Abstract

Motivated by microfluidic applications, we investigate drag reduction in laminar pressure-driven flows in channels with streamwise-periodic superhydrophobic surfaces (SHSs) contaminated with soluble surfactant. We develop a model in the long-wave and weak-diffusion limit, where the streamwise SHS period is large compared with the channel height and the Péclet number is large. Using asymptotic and numerical techniques, we determine the influence of surfactant on drag reduction in terms of the relative strength of advection, diffusion, Marangoni effects and bulk–surface exchange. In scenarios with strong exchange, the drag reduction exhibits a complex dependence on the thickness of the bulk-concentration boundary layer and surfactant strength. Strong Marangoni effects immobilise the interface through a linear surfactant distribution, whereas weak Marangoni effects yield a quasi-stagnant cap. The quasi-stagnant cap has an intricate structure with an upstream slip region, followed by intermediate inner regions and a quasi-stagnant region that is mediated by weak bulk diffusion. The quasi-stagnant region differs from the immobile region of a classical stagnant cap, observed for instance in surfactant-laden air bubbles in water, by displaying weak slip. As exchange weakens, the bulk and interface decouple: the surfactant distribution is linear when the surfactant is strong, whilst it forms a classical stagnant cap when the surfactant is weak. The asymptotic solutions offer closed-form predictions of drag reduction across much of the parameter space, providing practical utility and enhancing understanding of surfactant dynamics in flows over SHSs.

1. Introduction

Superhydrophobic surfaces (SHSs) offer a promising avenue for decreasing drag in laminar and turbulent flows (Lee et al. 2016; Park et al. 2021). Hydrophobic chemistry and microscopic topography on the SHSs entrap gas bubbles that lubricate the flow and decrease the wall-average shear stress compared with solid walls. This design presents numerous opportunities for industrial and environmental applications, including enhanced cooling (Lam et al. 2015) and reduced emissions (Xu et al. 2020). Consequently, SHSs have attracted significant attention in the academic literature (Park et al. 2014; Schönecker et al. 2014; Türk et al. 2014; Cheng et al. 2015; Seo & Mani 2018; Rastegari & Akhavan 2019; Kirk et al. 2020; Landel et al. 2020; Tomlinson et al. 2023b ). However, recent investigations have highlighted potential challenges. Interfacial displacement and deflection (Ng & Wang 2009; Teo & Khoo 2010), gas-phase flow dynamics (Game et al. 2017) and heat- and surfactant-induced Marangoni stresses (Peaudecerf et al. 2017; Kirk et al. 2020) have been shown to impede the anticipated drag reduction. Consequently, in certain applications, SHSs may not outperform solid walls. Motivated by these questions, this study investigates drag reduction in a surfactant-contaminated flow confined between a SHS and a solid wall. Specifically, we focus on the computationally challenging regime of weak bulk molecular diffusion, a regime which has remained largely unexplored in prior literature.

Surfactants have been measured in both artificial (Hourlier-Fargette et al. 2018; Temprano-Coleto et al. 2023) and natural (Pereira et al. 2018; Frossard et al. 2019) environments. When a liquid flow becomes contaminated with soluble surfactant, these chemical compounds are transported throughout the flow, adsorbing onto interfaces and desorbing into the bulk. The surfactant distribution within the bulk is regulated by the bulk Péclet number ( $\textit {Pe}$ ), forming concentration boundary layers when $\textit {Pe}$ is large. As the flow carries surfactant towards the downstream stagnation point of an interface, it accumulates, forming a concentration gradient. Depending on factors such as the flow properties, surfactant characteristics and geometry, the resulting adverse Marangoni force can cause the theoretically shear-free interface to behave like a no-slip (or partially no-slip) wall (Peaudecerf et al. 2017). Experimental support for this mechanism can be found in studies by Kim & Hidrovo (2012), Bolognesi et al. (2014), Schäffel et al. (2016), Peaudecerf et al. (2017), Song et al. (2018) and Temprano-Coleto et al. (2023).

Peaudecerf et al. (2017) and Landel et al. (2020) conducted numerical simulations of a two-dimensional (2-D) channel flow featuring streamwise-periodic SHSs using COMSOL, employing surfactant properties that are characteristic of sodium dodecylsulphate (SDS). Temprano-Coleto et al. (2023) extended this methodology to encompass a three-dimensional channel flow featuring streamwise- and spanwise-periodic SHSs, while Sundin & Bagheri (2022) investigated 2-D channel flow laden with surfactants over streamwise-periodic liquid-infused surfaces. In all of these configurations, numerical simulations of the Navier–Stokes equations and advection–diffusion equations for bulk and interfacial surfactant incurred substantial computational costs, particularly in the high-Péclet-number regime, where bulk-concentration boundary layers can be extremely thin. To address these challenges of high Péclet number values, we develop a numerical method based on Chebyshev collocation. This numerical method clusters nodes around the bulk-concentration boundary layers, reducing computational demand compared with numerical simulations. We complement numerical simulations with asymptotic analysis that reveals the underlying flow structures and offers direct predictions of flow properties.

The distinction between strictly insoluble and very weakly soluble surfactant can be subtle. In the former case, there is no exchange between bulk and interface and the amount of surfactant present on a SHS plastron must be prescribed as an input parameter. In the latter case, adsorption and desorption processes control bulk and interface exchange and the amount of surfactant present on a SHS plastron is determined by coupled transport processes on the interface and in the bulk. In many practical applications, diffusive effects are weak, and might not be expected to influence the manner in which surfactant compromises drag reduction. However, the surfactant flux between bulk and interface is mediated by bulk diffusion, enabling it to play a central role in determining drag reduction, even at very high Péclet numbers. In such circumstances, weak bulk diffusion can be the cause of weak apparent solubility, even when adsorption and desorption processes are very fast. Furthermore, while bulk diffusion cannot be neglected, interfacial diffusion can be, but then adsorption and desorption must be accommodated by weak stretching or compression of the interface. Taking these factors together, we will see how structures such as the classical stagnant cap, of some prescribed size when surfactant is strictly insoluble, must be replaced by a more mobile structure that we term a quasi-stagnant cap.

Scaling theories have been developed to analyse the slip length over SHSs and liquid-infused surfaces with streamwise-periodic gratings (Landel et al. 2020; Sundin & Bagheri 2022) and SHSs with streamwise- and spanwise-periodic gratings (Temprano-Coleto et al. 2023). These theories typically assume small surfactant concentrations and uniform shear stresses at the interface. For common surfactants such as SDS, Temprano-Coleto et al. (2023) demonstrated a significant slip length when the grating length exceeds a modified depletion length and a mobilisation length. The modified depletion length depends on the depletion length (Manikantan & Squires 2020), the height of the channel and the ratio of interfacial and bulk diffusivities. It characterises the interfacial length above which interfacial diffusion is weak compared to bulk–interface exchanges. The mobilisation length depends on the surfactant concentration and surfactant properties through the Marangoni, Damköhler and Biot numbers, and the SHS geometry (Temprano-Coleto et al. 2023). The mobilisation length being generally larger than the modified depletion length, it is the key length scale above which an interface can display significant slip. Landel et al. (2020), Temprano-Coleto et al. (2023) and Sundin & Bagheri (2022) used numerical simulations to calibrate the empirical coefficients linked to the bulk-concentration boundary-layer thickness. Unlike these prior studies, our asymptotic approach to the problem bypasses the need for empirical coefficient fitting, extending the slip-length scalings identified in previous studies (Landel et al. 2020; Temprano-Coleto et al. 2023) and uncovering new drag reduction regimes for surfactant-contaminated SHSs.

Tomlinson et al. (2023a ) developed an asymptotic theory for laminar channel flow featuring streamwise- and spanwise-periodic grooves and then investigated the impact of spatio–temporal fluctuations in surfactant concentration over the SHS (Tomlinson et al. 2024). Their theory assumed that the channel is long and bulk diffusion is strong enough to eliminate cross-channel concentration gradients. They derived a spatiallyone-dimensional (1-D) long-wave model that accommodated non-uniform shear stresses at the interface and enabled an analysis of the soluble stagnant-cap regime at low bulk Péclet numbers. The numerical simulations performed at moderate bulk Péclet number by Sundin & Bagheri (2022) show adsorption at the upstream end of the stagnant cap, but they did not offer a scaling theory for this regime. Similar stagnant-cap scenarios, where there is no adsorption at the upstream end of the interface, were studied by Baier & Hardt (2021) and Mayer & Crowdy (2022) in the context of insoluble surfactant, with a linear and nonlinear equation of state, respectively. These studies were extended to include protrusion at the liquid–gas interface: Baier & Hardt (2022) found recirculating interfacial flows for weak protrusions, which Rodriguez-Broadbent & Crowdy (2023) extended to larger protrusions, deriving an effective slip length. Crowdy et al. (2023) investigated a linear extensional flow between two fluids. Assuming that $\textit {Pe}$ is zero, they showed that solubility and strong exchange with the bulk can reduce the length of the stagnant cap and make the interface less immobile. We show how ‘remobilisation’ due to strong bulk–surface exchange can also arise when $Pe$ is large, revealing the coupling between the bulk-concentration boundary layer and the underlying interface, which is particularly intricate when Marangoni effects are weak.

Tomlinson et al. (2023a ) delineated three key asymptotic regions of parameter space: one dominated by Marangoni effects, with minimal drag reduction; and others dominated by advection and diffusion, with significant drag reduction. They also identified when their 1-D long-wave model first breaks down due to 2-D effects. In this paper, we employ a combination of asymptotic and numerical methods to investigate the impact of soluble surfactant on a laminar pressure-driven channel flow confined between a streamwise-periodic SHS and a solid wall. We assume that the period of the SHS is significantly longer than the height of the channel and that the bulk and surface Péclet numbers are large. To address the numerical challenges of these high-Péclet-number flows, we develop asymptotic solutions for weak diffusion, which may offer insights into a broader range of surfactant-contaminated flows than flows over SHSs. For example, numerical simulations of flows involving soluble-surfactant-contaminated drops and bubbles exhibit velocity and surfactant distributions reminiscent of stagnant caps at large but finite $\textit {Pe}$ (Oguz & Sadhal 1988; Tasoglu et al. 2008). However, the structure of these stagnant caps has only been described in the insoluble limit (Harper 2004; Palaparthi et al. 2006).

Additionally, we construct asymptotic expressions for practical quantities such as the slip length and drag reduction. These asymptotic expressions not only aid in understanding the physics of surfactant-contaminated flows over SHSs but also facilitate the prediction of the slip length and drag reduction in experiments.

In § 2, we formulate the 2-D problem. We non-dimensionalise the fluid and surfactant equations in the long-wave limit and derive 1-D and 2-D long-wave models for bulk and interfacial surfactant behaviour. The calculations underpinning drag reduction predictions are laid out in appendices. In § 3, we present our main asymptotic and numerical findings. We analyse the drag reduction and compare the predicted flow and surfactant fields with results obtained from COMSOL simulations and experiments. In § 4, we illustrate the uses of our study and discuss the key asymptotic results that may be generalisable to broader surfactant-contaminated flows.

2. Formulation

2.1. Governing equations of the 2-D model

Figure 1. A 2-D laminar pressure-driven channel flow transporting soluble surfactant confined between a streamwise-periodic SHS and a solid wall. We position the origin of the Cartesian coordinate system, $(\hat {x}, \, \hat {y})$ , at the centre of the liquid–gas interface. Each periodic cell has channel height $2\hat {H}$ and period length $2\hat {P}$ . The SHS, characterised by gas fraction $\phi$ , has an interface region of length $2\phi \hat {P}$ and a solid region of length $2(1 - \phi )\hat {P}$ ; the area above these regions defines subdomains $\hat {\mathcal{D}}_1$ and $\hat {\mathcal{D}}_2$ , respectively, as specified in (2.1).

We explore a steady 2-D laminar pressure-driven channel flow contaminated with soluble surfactant, confined between a streamwise-periodic SHS and a solid wall (see figure 1). The streamwise and wall-normal directions are represented by $\hat {x}$ and $\hat {y}$ coordinates respectively, where hats indicate dimensional quantities. Assuming the fluid to be incompressible and Newtonian, we define the velocity field $\hat {\boldsymbol {u}}=(\hat {u}(\hat {x},\,\hat {y}),\,\hat {v}(\hat {x},\,\hat {y}))$ , pressure field $\hat {p}=\hat {p}(\hat {x},\,\hat {y})$ , bulk surfactant concentration $\hat {c}=\hat {c}(\hat {x},\,\hat {y})$ and interfacial surfactant concentration $\hat {\Gamma }=\hat {\Gamma }(\hat {x})$ . Due to the periodicity of the SHS in the streamwise direction, our analysis focuses on a single periodic cell with dimensions $2\hat {P}$ in length and $2\hat {H}$ in height. At $\hat {y}=0$ , there is a solid ridge of length $2(1-\phi )\hat {P}$ and a liquid–gas interface, or plastron, of length $2\phi \hat {P}$ , where $\phi$ is the gas fraction. The interface is assumed to be flat. There is a solid boundary at $\hat {y}= 2\hat {H}$ . The periodic domain is partitioned into two subdomains,

(2.1a) \begin{align} \hat {\mathcal{D}}_1 &= \{\hat {x}\in [-\phi \hat {P},\, \phi \hat {P}]\} \times \{\hat {y}\in [0, \,2 \hat {H}]\}, \end{align}

(2.1b) \begin{align} \hat {\mathcal{D}}_2 &= \{\hat {x}\in [\phi \hat {P}, \, (2 - \phi )\hat {P}]\} \times \{\hat {y}\in [0,\, 2 \hat {H}]\}, \end{align}

as illustrated in figure 1.

A comprehensive discussion of the equations governing fluid and surfactant behaviour in a three-dimensional geometry featuring streamwise- and spanwise-periodic SHSs was presented in Tomlinson et al. (2023a ). Here, we provide a succinct overview of the model for a streamwise-periodic SHS in two dimensions. Within domains $\hat {\mathcal{D}}_1$ and $\hat {\mathcal{D}}_2$ , we have Stokes equations with an advection–diffusion equation for bulk surfactant:

(2.2a--c) \begin{align} \hat {\boldsymbol {\nabla }} \cdot \hat {\boldsymbol {u}} = 0, \quad \hat {\mu } \hat {\nabla }^2 \hat {\boldsymbol {u}}- \hat {\boldsymbol {\nabla }} \hat {p} = \boldsymbol {0}, \quad \hat {D} \hat {\nabla }^2 \hat {c} - \hat {\boldsymbol {u}}\cdot \hat {\boldsymbol {\nabla }} \hat {c} = 0, \end{align}

where $\hat {\mu }$ is the dynamic viscosity and $\hat {D}$ is the bulk diffusivity. Linearising the equation of state, we assume $\hat {\sigma }_{\hat {x}} = - \hat {A}\hat {\Gamma }_{\hat {x}}$ , where $\hat {\sigma }$ is the surface tension and $\hat {A}$ is the surface activity. For $\hat {x} \in [-\phi \hat {P}, \, \phi \hat {P}]$ and $\hat {y}=0$ , we impose the tangential stress balance, no penetration, linear adsorption–desorption kinetics and an advection–diffusion equation for interfacial surfactant:

(2.3a--d) \begin{align} \hat {\mu } \hat {u}_{\hat {y}} - \hat {A} \hat {\Gamma }_{\hat {x}} = 0, \quad \hat {v}=0, \quad \hat {D} \hat {c}_{\hat {y}} - \hat {K}_a \hat {c} + \hat {K}_d \hat {\Gamma } =0,\nonumber \\ \quad \hat {D}_I \hat {\Gamma }_{\hat {x}\hat {x}} + \hat {K}_a \hat {c} - \hat {K}_d \hat {\Gamma } -(\hat {u} \hat {\Gamma })_{\hat {x}} = 0, \end{align}

where $\hat {D}_I$ is the interfacial diffusivity, $\hat {K}_a$ is the adsorption rate and $\hat {K}_d$ is the desorption rate. At $\hat {x} = \pm \phi \hat {P}$ and $\hat {y}=0$ , there is no interfacial flux of surfactant:

(2.4) \begin{align} \hat {u}\hat {\Gamma } - \hat {D}_I \hat {\Gamma }_{\hat {x}} = 0. \end{align}

For $\hat {x} \in [\phi \hat {P}, \, (2-\phi )\hat {P}]$ and $\hat {y}=0$ , and likewise for $\hat {x} \in [-\phi \hat {P}, \, (2 - \phi )\hat {P}]$ and $\hat {y}=2\hat {H}$ , we impose no slip, no penetration and no bulk flux of surfactant:

(2.5a--c) \begin{align} \hat {u} = 0, \quad \hat {v} = 0, \quad \hat {c}_{\hat {y}} =0. \end{align}

For $\hat {\boldsymbol {q}} \equiv (\hat {\boldsymbol {u}}, \,\hat {p}_{\hat {x}},\,\hat {c})$ , continuity and periodicity conditions on bulk quantities are given by

(2.6) \begin{align} \hat {\boldsymbol {q}}((\phi \hat {P})^-,\,\hat {y}) = \hat {\boldsymbol {q}}((\phi \hat {P})^+,\,\hat {y}), \quad \hat {\boldsymbol {q}}(-\phi \hat {P},\,\hat {y}) = \hat {\boldsymbol {q}}((2- \phi )\hat {P},\,\hat {y}), \end{align}

where the superscript $-$ ( $+$ ) indicates evaluation in $\hat {\mathcal{D}}_1$ ( $\hat {\mathcal{D}}_2$ ). The pressure field is continuous across $x=\pm \phi \hat {P}$ but is not periodic. Within $\hat {\mathcal{D}}_1$ , the bulk surfactant equations and boundary conditions, (2.2c ), (2.3c ) and (2.5c ), can be combined to give

(2.7) \begin{align} \frac {\text{d}}{\text{d}\hat {x}} \int _{0}^{2 \hat {H}} (\hat {u}\hat {c} - \hat {D} \hat {c}_{\hat {x}}) \, \text{d}\hat {y} - (\hat {K}_d \hat {\Gamma } - \hat {K}_a \hat {c}(\hat {x},\, 0)) = 0, \end{align}

which, combined with the interfacial surfactant equation (2.3d ), leads to an expression for the total flux of surfactant, $\hat {K}$ , which must be uniform:

(2.8) \begin{align} \int _{0}^{2 \hat {H}} (\hat {u}\hat {c} - \hat {D} \hat {c}_{\hat {x}}) \, \text{d}\hat {y} + (\hat {u}(\hat {x},\, 0)\hat {\Gamma } - \hat {D}_I \hat {\Gamma }_{\hat {x}}) = \hat {K}. \end{align}

Integrating the interfacial surfactant equation (2.3d ) over the plastron and utilising the no-flux condition (2.4) yields a constraint on net adsorption and desorption:

(2.9) \begin{align} \int _{ - \phi \hat {P}}^{\phi \hat {P}} (\hat {K}_d \hat {\Gamma } - \hat {K}_a \hat {c}(\hat {x},\, 0)) \, \text{d}\hat {x} = 0. \end{align}

Similarly, within $\hat {\mathcal{D}}_2$ , (2.3c ) and (2.5c ) can be combined to give

(2.10) \begin{align} \int _{0}^{2 \hat {H}} (\hat {u}\hat {c} - \hat {D} \hat {c}_{\hat {x}}) \, \text{d}\hat {y} = \hat {K}. \end{align}

From (2.2a ), (2.3b ) and (2.5b ), it follows that the volume flux of fluid in $\hat {\mathcal{D}}_1$ and $\hat {\mathcal{D}}_2$ , $\hat {Q}$ , is also uniform, where

(2.11) \begin{align} \hat {Q} = \int _{0}^{2 \hat {H}} \hat {u} \, \text{d}\hat {y}. \end{align}

The fluxes $\hat {Q}$ and $\hat {K}$ are prescribed in this model.

The flow is driven in the streamwise direction by a cross-channel-averaged pressure drop per period, $\Delta \hat {p} \equiv \langle \hat {p}\rangle (-\phi \hat {P}) - \langle \hat {p}\rangle ((2-\phi )\hat {P}) \gt 0$ , where $\langle \cdot \rangle \equiv \int _{0}^{2 \hat {H}} \cdot \, \text{d}\hat {y} / (2\hat {H})$ . We define the normalised drag reduction

(2.12) \begin{align} {DR} = \frac {\Delta \hat {p}_I - \Delta \hat {p}}{\Delta \hat {p}_I - \Delta \hat {p}_U}, \end{align}

where $\Delta \hat {p} = \Delta \hat {p}_I$ when the interface is immobilised due to surfactant ( ${DR}=0$ ) and $\Delta \hat {p} = \Delta \hat {p}_U$ when the interface is unaffected by surfactant ( ${DR}=1$ ). We aim to determine the dependence of $DR$ on the dimensional parameters of the problem ( $\phi$ , $\hat {P}$ , $\hat {H}$ , $\hat {\mu }$ , $\hat {D}$ , $\hat {A}$ , $\hat {D}_I$ , $\hat {K}_a$ , $\hat {K}_d$ , $\hat {K}$ and $\hat {Q}$ ) in the singular high-Péclet-number regime.

2.2. Non-dimensionalisation

We non-dimensionalise the governing equations (2.1)–(2.12) using $\hat {Q}/\hat {H}$ for the velocity scale, $\hat {K}/\hat {Q}$ for the bulk concentration scale and $\hat {K}_a \hat {K}/(\hat {K}_d\hat {Q})$ for the interfacial concentration scale. We write

(2.13a--g) \begin{align} x = \frac {\hat {x}}{\hat {P}}, \quad y = \frac {\hat {y}}{\epsilon \hat {P}}, \quad u = \frac {\hat {u}}{\hat {Q}/\hat {H}}, \quad v = \frac {\hat {v}}{\epsilon \hat {Q}/\hat {H}}, \nonumber \\ p = \frac {\hat {p}}{\hat {\mu } \hat {Q}/(\epsilon \hat {H}^2)}, \quad c = \frac {\hat {c}}{\hat {K}/\hat {Q}}, \quad \Gamma = \frac {\hat {\Gamma }}{\hat {K}_a\hat {K}/(\hat {K}_d\hat {Q})}, \end{align}

where $\epsilon \equiv \hat {H}/ \hat {P}$ is the slenderness parameter. The subdomains (2.1) become

(2.14a) \begin{align} \mathcal{D}_1 &= \{x\in [-\phi, \,\phi ]\} \times \{y\in [0, \,2]\}, \end{align}

(2.14b) \begin{align} \mathcal{D}_2 &= \{x\in [\phi, \,2 - \phi ]\} \times \{y\in [0,\, 2]\}. \end{align}

In $\mathcal{D}_1$ and $\mathcal{D}_2$ , the bulk equations (2.2) become

(2.15a--d) \begin{align} u_x + v_{y} = 0, \quad \epsilon ^2 u_{xx} + u_{yy} - p_x = 0,\nonumber \\ \epsilon ^4 v_{xx} + \epsilon ^2 v_{yy} - p_y = 0, \quad (\epsilon ^2 c_{xx} + c_{yy})/{\textit {Pe}} - \epsilon u c_x - \epsilon v c_y = 0, \end{align}

with ${\textit {Pe}} = \hat {Q}/\hat {D}$ the bulk Péclet number. For $x\in [-\phi, \, \phi ]$ and $y=0$ , the interface conditions (2.3) give

(2.16a--d) \begin{align} u_y - \epsilon {\textit {Ma}} \Gamma _{x} = 0, \quad v =0, \quad c_y - {\textit {Da}} (c -\Gamma ) =0, \nonumber \\ \quad \epsilon ^2 \Gamma _{xx}/{\textit {Pe}}_I + {\textit {Bi}} ( c - \Gamma ) - \epsilon (u \Gamma )_{x} = 0, \end{align}

with ${\textit {Ma}} = \hat {A}\hat {K}_a \hat {K}\hat {H}/(\hat {\mu }\hat {K}_d\hat {Q}^2)$ the Marangoni number, ${\textit {Da}}= \hat {K}_a \hat {H}/\hat {D}$ the Damköhler number, ${\textit {Pe}}_I= \hat {Q}/\hat {D}_I$ the interfacial Péclet number and ${\textit {Bi}} = \hat {K}_d \hat {H}^2/\hat {Q}$ the Biot number. At $x = \pm \phi$ and $y=0$ , the no-flux condition (2.4) becomes

(2.17) \begin{align} u \Gamma - \epsilon \Gamma _{x}/{\textit {Pe}}_I = 0. \end{align}

For $x\in [\phi, \, 2 - \phi ]$ and $y=0$ (and for $x\in [-\phi, \, 2-\phi ]$ and $y=2$ ), the solid boundary conditions (2.5) give

(2.18a--d) \begin{align} u = 0, \quad v = 0, \quad \ c_{y} =0. \end{align}

For $\boldsymbol {q} = (\boldsymbol {u}, \, p_x, \,c)$ , the matching conditions (2.6) become

(2.19a,b) \begin{align} \boldsymbol {q}(\phi ^-,\,y) = \boldsymbol {q}(\phi ^+,\,y), \quad \boldsymbol {q}(-\phi, \,y) = \boldsymbol {q}(2- \phi, \,y). \end{align}

For $x\in [-\phi, \, 2 - \phi ]$ , the surfactant-flux conditions (2.7)–(2.10) yield

(2.20a) \begin{align} \frac {\text{d}}{\text{d}x} \int _{0}^{2} \left ( u c - \frac {\epsilon c_x}{{\textit {Pe}}} \right ) \, \text{d} y - \frac {{\textit {Da}}}{\epsilon {\textit {Pe}}}(\Gamma - c(x,\, 0)) &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.20b) \begin{align} \int _{0}^{2} \left ( u c - \frac {\epsilon c_x}{{\textit {Pe}}}\right ) \, \text{d}y + \frac {{\textit {Da}}}{{\textit {Bi}} {\textit {Pe}}} \left ( u(x,\, 0) \Gamma - \frac {\epsilon \Gamma _{x}}{{\textit {Pe}}_I}\right ) &= 1 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.20c) \begin{align} \int _{ - \phi }^{\phi } ( \Gamma - c(x,\, 0) ) \, \text{d}x &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.20d) \begin{align} \int _{0}^{2} \left (u c - \frac {\epsilon c_x}{{\textit {Pe}}} \right ) \, \text{d}y &= 1 \quad \text{in} \quad \mathcal{D}_2, \\[2pt] \nonumber \end{align}

and in $\mathcal{D}_1$ and $\mathcal{D}_2$ , the volume-flux condition (2.11) gives

(2.21) \begin{align} \int _{0}^{2} u \, \text{d}y = 1. \end{align}

Finally, the drag reduction (2.12) becomes

(2.22) \begin{align} {DR} = \frac {\Delta p_I - \Delta p}{\Delta p_I - \Delta p_U}, \end{align}

where $\Delta {p} \equiv \langle p\rangle (-\phi ) - \langle p\rangle (2-\phi )$ and $\langle \cdot \rangle \equiv \int _{0}^{2} \cdot \, \text{d} y / 2$ .

In the limit $\epsilon \ll 1$ , we proceed to solve the leading-order boundary-value problem and flux constraints defined in (2.15)–(2.21) for $u$ , $v$ , $p$ , $c$ and $\Gamma$ using asymptotic and numerical methods. Subsequently, we evaluate the drag reduction (2.22) and investigate its dependence on the seven dimensionless groups ( $\epsilon$ , $\textit {Pe}$ , $\textit {Ma}$ , ${\textit {Pe}}_I$ , $\textit {Bi}$ , $\textit {Da}$ and $\phi$ ) that characterise the geometry, flow, liquid and surfactant. Moreover, full solutions of (2.15)–(2.21) obtained as COMSOL simulations are presented in § 3.

2.3. The long-wave limit of the 2-D model at high Péclet numbers

We investigate distinguished limits of (2.15)–(2.21) to uncover different physical balances. We assume that $1/{\textit {Pe}} = O(\epsilon )$ , $1/{\textit {Pe}}_I = O(\epsilon )$ , ${\textit {Bi}} = O(\epsilon )$ , ${\textit {Da}} = O(1)$ and ${\textit {Ma}} = O(1/\epsilon )$ as $\epsilon \rightarrow 0$ , denoting this as the ‘2-D long-wave model’ in which streamwise and cross-channel concentration gradients are comparable at leading order. We rescale $1/{\textit {Pe}} = \epsilon / \mathscr{P}$ , $1/{\textit {Pe}}_I = \epsilon /\mathscr{P}_I$ , ${\textit {Bi}} = \epsilon \mathscr{B}$ and ${\textit {Ma}} = \mathscr{M}/\epsilon$ , assuming that $\mathscr{P}$ , $\mathscr{P}_I$ , $\mathscr{B}$ and $\mathscr{M}$ remain $O(1)$ as $\epsilon \rightarrow 0$ . We substitute the expansions

(2.23) \begin{align} (u,\,v, \,p,\, c,\, \Gamma ) = (u_0,\,v_0, \,p_0,\, c_0, \,\Gamma _0) + \epsilon ^2 (u_1,\,v_1, \,p_1,\, c_1, \,\Gamma _1) + \cdots \end{align}

into the governing equations (2.15)–(2.22) and take the leading-order approximation.

In $\mathcal{D}_1$ and $\mathcal{D}_2$ , flow is driven by a streamwise pressure gradient, and wall-normal diffusion is comparable to advection. The bulk equations (2.15) become

(2.24a--d) \begin{align} u_{0x} + v_{0y} = 0, \quad u_{0yy} - p_{0x} = 0, \quad p_{0y} = 0, \quad c_{0yy}/\mathscr{P} - u_0 c_{0x} - v_0 c_{0y} = 0. \end{align}

For $x\in [-\phi, \, \phi ]$ and $y=0$ , the interfacial flow is inhibited by the surfactant gradient and exchange is comparable to advection. The interface conditions (2.16) yield

(2.25a--d) \begin{align} u_{0y} - \mathscr{M} \Gamma _{0x} = 0, \quad v_0 =0, \quad c_{0y} - {\textit {Da}} (c_0 -\Gamma _0) =0, \nonumber \\ \mathscr{B} ( c_0 - \Gamma _0) - (u_0 \Gamma _0)_{x} = 0. \end{align}

For $x\in [\phi, \, 2 - \phi ]$ and $y=0$ (and for $x\in [-\phi, \, 2-\phi ]$ and $y=2$ ), the solid wall boundary conditions (2.18) become

(2.26a--c) \begin{align} u_0 = 0, \quad v_0 = 0, \quad \ c_{0y} =0. \end{align}

Short 2-D Stokes-flow regions arise at the junctions between $\mathcal{D}_1$ and $\mathcal{D}_2$ . In the present long-wave theory, in which surface diffusion appears at the next order, we impose the no-flux condition (2.17) at $x = \pm \phi$ and $y=0$ :

(2.27) \begin{align} u_0 \Gamma _0 = 0. \end{align}

Continuity of bulk concentration (2.19) between $\mathcal{D}_1$ and $\mathcal{D}_2$ requires

(2.28a,b) \begin{align} c_0(\phi ^-,\,y) = c_0(\phi ^+,\,y), \quad c_0(-\phi, \,y) = c_0(2- \phi, \,y). \end{align}

In § 3, comparison with COMSOL simulations will allow us to investigate the effect of these Stokes-flow regions. Streamwise surfactant transport in the bulk and interface is dominated by advection and exchange, so (2.20a ) becomes

(2.29a) \begin{align} \frac {\text{d}}{\text{d}x} \int _{0}^{2} u_0 c_0 \, \text{d} y - \frac {Da}{\mathscr{P}}(\Gamma _0 - c_0(x,\, 0)) &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.29b) \begin{align} \int _{0}^{2} u_0 c_0 \, \text{d}y + \frac {{\textit {Da}}}{\mathscr{B}\mathscr{P}} \left ( u_0(x,\, 0) \Gamma _0 \right ) &= 1 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.29c) \begin{align} \int _{ - \phi }^{\phi } ( \Gamma _0 - c_0(x,\, 0) ) \, \text{d}x &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.29d) \begin{align} \int _{0}^{2} u_0 c_0 \, \text{d}y &= 1 \quad \text{in} \quad \mathcal{D}_2. \end{align}

In $\mathcal{D}_1$ and $\mathcal{D}_2$ , the volume-flux condition (2.21) gives

(2.30) \begin{align} \int _{0}^{2} u_0 \, \text{d}y = 1. \end{align}

The leading-order drag reduction (2.22) is given by

(2.31) \begin{align} {DR}_0 = \frac {\Delta p_I - \Delta p_0}{\Delta p_I - \Delta p_U}. \end{align}

The leading-order pressure field simplifies to $p_0 = p_0(x)$ using the wall-normal equation (2.24c ). The streamwise velocity field is driven by the pressure gradient $p_{0x}$ and is inhibited by the surfactant gradient $\Gamma _{0x}$ . Using linear superposition, we can write

(2.32a,b) \begin{align} u_0 = \tilde {U} p_{0x} + \mathscr{M} \bar {U}\Gamma _{0x} \quad \text{in} \quad \mathcal{D}_1, \quad u_0 = \breve {U} p_{0x} \quad \text{in} \quad \mathcal{D}_2, \end{align}

where $\tilde {U} \equiv y^2/2 - 2$ , $\bar {U} \equiv y - 2$ and $\breve {U} \equiv y^2/2 - y$ are velocity contributions satisfying the following boundary-value problems:

(2.33a–c) \begin{align} \tilde {U}_{yy} = 1, \quad \text{subject to} \quad \tilde {U}_y(0)=0, \quad \tilde {U}(2)=0; \\[-9pt] \nonumber \end{align}

(2.34a–c) \begin{align} \bar {U}_{yy} = 0, \quad \text{subject to} \quad \bar {U}_y(0)= 1, \quad \bar {U}(2) = 0; \\[-9pt] \nonumber \end{align}

(2.35a–c) \begin{align} \breve {U}_{yy} = 1, \quad \text{subject to} \quad \breve {U}(0)= 0, \ \quad \breve {U}(2)=0. \\[-9pt] \nonumber \end{align}

Substituting (2.32) into (2.30), the volume flux constraint gives

(2.36a,b) \begin{align} \tilde {Q} p_{0x} + \mathscr{M} \bar {Q} \Gamma _{0x} =1 \quad \text{in} \quad \mathcal{D}_1, \quad \breve {Q} p_{0x} = 1 \quad \text{in} \quad \mathcal{D}_2, \end{align}

where $\tilde {Q} \equiv -8/3$ , $\bar {Q} \equiv -2$ and $\breve {Q} \equiv -2/3$ . Using the continuity equation (2.24a ), the wall-normal velocity field is forced by $p_{0xx}$ and $\Gamma _{0xx}$ , with velocity contributions that complement those given in (2.32)–(2.34):

(2.37a,b) \begin{align} v_0 = \tilde {V} p_{0xx} + \mathscr{M} \bar {V} \Gamma _{0xx} \quad \text{in} \quad \mathcal{D}_1, \quad v_0 = 0 \quad \text{in} \quad \mathcal{D}_2, \end{align}

where $\tilde {V} \equiv -y^3/6 + 2 y$ and $\bar {V} \equiv -y^2/2 + 2 y$ (so that $\tilde {V}_{y} = -\tilde {U}$ , $\tilde {V}(0)=0$ , $\bar {V}_{y} = -\bar {U}$ , $\bar {V}(0)= 0$ ). The no-penetration condition (2.26b ) is satisfied, as $v_0(2) = \tilde {V}(2) p_{0xx} + \mathscr{M}\bar {V}(2)\Gamma _{0xx} = -\tilde {Q} p_{0xx} - \mathscr{M}\bar {Q}\Gamma _{0xx} = 0$ , using the volume-flux condition (2.36a ).

Table 1. A summary of the transport coefficients in the 2-D long-wave model, (2.38)–(2.42), with their definition and physical interpretation. The transport coefficients $\alpha$ , $\beta$ , $\gamma$ , $\delta$ and $\nu$ are written in terms of $\textit {Pe}$ , $\textit {Da}$ , $\textit {Bi}$ , $\textit {Ma}$ , ${\textit {Pe}}_I$ and $\epsilon$ defined in § 2.2. Compared with the three-dimensional transport coefficients in Tomlinson et al. (2023a ), $\alpha _{{3D}}\propto \epsilon ^2\alpha$ , $\beta _{{3D}}\propto \beta$ , $\gamma _{{3D}}\propto \gamma$ , $\delta _{{3D}}\propto \epsilon ^2\delta$ and $\nu _{{3D}}\propto \epsilon ^2\nu$ .

We substitute (2.32) and (2.37) into the bulk and interfacial surfactant equations (2.24)–(2.29) to derive the 2-D long-wave model. To facilitate its numerical computation, we retain $O(\epsilon ^2)$ bulk and interfacial streamwise diffusion operators to make the problem more elliptic, simplifying the continuity conditions and smoothing the flow between $\mathcal{D}_1$ and $\mathcal{D}_2$ , later seeking the limit in which these effects are weak. The bulk surfactant equation (2.24c ) gives

(2.38a) \begin{align} \alpha (\epsilon ^2 c_{0xx} + c_{0yy}) - \bigg (\frac {3}{4} - \frac {3y^2}{16}\bigg ) & c_{0x} \nonumber \\ - \frac {\gamma }{\beta } \bigg (-\frac {1}{2} + y - \frac {3y^2}{8}\bigg )& \Gamma _{0x} c_{0x} - \frac {\gamma }{\beta } \bigg (\frac {y}{2} - \frac {y^2}{2} + \frac {y^3}{8}\bigg ) \Gamma _{0xx} c_{0y} = 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.38b) \begin{align} &\alpha (\epsilon ^2 c_{0xx} + c_{0yy}) - \bigg (\frac {3y}{2} - \frac {3y^2}{4}\bigg ) c_{0x} = 0 \quad \text{in} \quad \mathcal{D}_2, \\[6pt] \nonumber \end{align}

where $\alpha$ , $\beta$ and $\gamma$ are given in table 1. This table defines the five primary parameter combinations that we use below and gives their physical interpretation. For $x\in [-\phi, \, \phi ]$ and $y=0$ , the boundary condition specifying continuity of surfactant flux and the interfacial surfactant equation (2.25c , d ) become

(2.39a,b) \begin{align} c_{0y} - \frac {\nu }{\alpha }(c_{0} - \Gamma _0) = 0, \quad \nu ( c_{0} - \Gamma _0) - \frac {3}{4}\beta \Gamma _{0x} + \frac {1}{2}\gamma (\Gamma _{0x}\Gamma _0)_x + \epsilon ^2 \delta \Gamma _{0xx} = 0, \end{align}

where $\nu$ and $\delta$ are given in table 1. At $x= \pm \phi$ and $y=0$ , no flux of interfacial surfactant (2.27) gives

(2.40) \begin{align} \frac {3}{4}\beta \Gamma _{0} - \frac {1}{2}\gamma \Gamma _{0x}\Gamma _0 - \epsilon ^2 \delta \Gamma _{0x} = 0. \end{align}

For $x\in [\phi, \, 2 - \phi ]$ and $y=0$ (and for $x\in [-\phi, \, 2-\phi ]$ and $y=2$ ), no flux of bulk surfactant (2.26c ) becomes

(2.41) \begin{align} c_{0y} =0. \end{align}

The surfactant-flux conditions (2.29a , b , d ) become

(2.42a) \begin{align} \frac {\text{d}}{\text{d}x} \int _{0}^{2} \left (\bigg (\frac {3}{4} - \frac {3y^2}{16}\bigg ) c_{0} + \frac {\gamma }{\beta } \bigg (-\frac {1}{2} + y - \frac {3y^2}{8}\bigg ) \Gamma _{0x} c_{0} - \epsilon ^2 \alpha c_{0x} \right )\text{d} y \nonumber \\ - \nu (\Gamma _0 - c_0(x,\, 0)) &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.42b) \begin{align} \int _{0}^{2} \left (\bigg (\frac {3}{4} - \frac {3y^2}{16}\bigg ) c_{0} + \frac {\gamma }{\beta } \bigg (-\frac {1}{2} + y - \frac {3y^2}{8}\bigg ) \Gamma _{0x} c_{0} - \epsilon ^2 \alpha c_{0x} \right )\text{d}y \nonumber \\ + \frac {3}{4}\beta \Gamma _{0} - \frac {1}{2}\gamma \Gamma _{0x}\Gamma _0 - \epsilon ^2 \delta \Gamma _{0x} &= 1 \quad \text{in} \quad \mathcal{D}_1, \end{align}

(2.42c) \begin{align} \int _{0}^{2} \left ( \bigg (\frac {3y}{2} - \frac {3y^2}{4}\bigg ) c_0 - \epsilon ^2 \alpha c_{0x} \right ) \, \text{d}y &= 1 \quad \text{in} \quad \mathcal{D}_2, \\[6pt] \nonumber \end{align}

and in $\mathcal{D}_1 \cup \mathcal{D}_2$ , the volume-flux condition (2.30) becomes

(2.43a,b) \begin{align} \tilde {Q} p_{0x} + (\gamma / \beta ) \bar {Q} \Gamma _{0x} = 1 \quad \text{in} \quad \mathcal{D}_1, \quad \breve {Q} p_{0x} = 1 \quad \text{in} \quad \mathcal{D}_2. \end{align}

To summarise, the seven-parameter Stokes-flow system (2.15)–(2.21) involving five variables $(\hat {u}, \hat {v}, \, \hat {p}, \, \hat {c}, \, \hat {\Gamma })$ has been reduced to the seven-parameter long-wave system (2.38)–(2.42) involving two variables $(\hat {c}, \, \hat {\Gamma })$ . Further, by neglecting streamwise diffusion in the bulk, $-\epsilon ^2 \alpha c_{0x}$ , and at the interface, $-\epsilon ^2 \delta \Gamma _{0x}$ , the system reduces to a five-parameter system ( $\alpha$ , $\beta$ , $\gamma$ , $\nu$ , $\phi$ ). This long-wave model disregards short (but passive) inner regions associated with the small parameters $\epsilon$ and $\delta$ .

The numerical solution of (2.38)–(2.42) (Appendix A) is generally less expensive (in terms of memory and runtime) than COMSOL simulations. The numerical method iteratively solves linearised bulk- and interfacial-concentration problems until convergence is achieved. Convergence in the nonlinear problem can be slow, necessitating small tolerances. Accurate resolution of singularities at $x=\pm \phi$ is essential to achieving convergence; accordingly, we employ domain decomposition and Chebyshev collocation techniques that cluster nodes around $y=0,\,2$ and $x = -\phi, \, \phi, \,2-\phi$ . Once $c_0$ and $\Gamma _0$ are determined, $u_0$ and $v_0$ can be calculated using (2.32) and (2.37), respectively.

Following Tomlinson et al. (2023a ), we express the leading-order drag reduction ( ${DR}_0$ ) in terms of the transport parameters and geometry ( $\beta$ , $\gamma$ and $\phi$ ) as follows. When the interface is unaffected by surfactant, $\Delta p_U = -2\phi /\tilde {Q} - 2(1-\phi )/\breve {Q}$ . When the interface is immobilised by surfactant, $\Delta p_I = -2/\breve {Q}$ . Between these two limits, $\Delta p_0 = -2\phi /\tilde {Q} - 2(1-\phi )/\breve {Q} + \gamma \bar {Q} \Delta \Gamma _0 / (\beta \tilde {Q})$ , where $\Delta \Gamma _0 = \Gamma _0(\phi ) - \Gamma _0(-\phi )$ . Substituting these expressions into (2.31) gives the leading-order drag reduction as

(2.44) \begin{align} {DR}_0 = 1 - \frac {\gamma \Delta \Gamma _0}{3\phi \beta }. \end{align}

When evaluated using COMSOL, we write ${DR}_0 = {DR}_{NS}$ . An alternative measure used commonly in the SHS literature is the effective slip length, $\lambda _0$ , which can be evaluated from ${DR}_0$ using (Tomlinson et al. 2023a )

(2.45) \begin{align} \lambda _0 = \frac {{DR}_0(\Delta p_I - \Delta p_U)}{\Delta p_I[\Delta p_U {DR}_0 + \Delta p_I(1 - {DR}_0)]}. \end{align}

Before presenting the results of the 2-D long-wave model, it is helpful to recall key features of the problem at low Péclet numbers, when a 1-D description applies.

Figure 2. Maps of ( $\alpha$ , $\gamma$ )-parameter space for $\epsilon \ll 1$ , using parameters given in table 1, showing the magnitude of drag reduction ( $DR_0$ ) across different asymptotic regions. Existing predictions from the 1-D long-wave model of Tomlinson et al. (2023a ) are combined with predictions derived in Appendix C from the 2-D long-wave model for (a) strong exchange and (b) weak exchange. The Marangoni-dominated region ( $M$ ) incorporates subregions $M^{{1D}}$ , $M^{{1D}}_E$ , $M^{{2D}}$ and $M^{{2D}}_E$ , for which ${DR}_0 \ll 1$ ; the advection-dominated region ( $A$ ) incorporates subregions $A^{{1D}}$ , $A^{{1D}}_E$ , $A^{{2D}}$ and $A^{{2D}}_E$ , for which $1-{DR}_0 \ll 1$ ; and the diffusion-dominated region ( $D$ ) incorporates subregions $D^{{1D}}$ and $D^{{1D}}_E$ , for which $1-{DR}_0 \ll 1$ . Here it is assumed that $\alpha \sim \delta$ and $\beta \sim 1$ .

2.4. The 1-D long-wave model

Figure 2 shows maps of $(\alpha, \,\gamma )$ -parameter space, illustrating how the present 2-D model relates to existing 1-D results. In the strong-cross-channel-diffusion limit, for which $\alpha \gg 1$ (with $\epsilon ^2\alpha =O(1)$ and $\epsilon ^2\delta =O(1)$ ), cross-channel concentration gradients are small. Expanding the bulk and interfacial concentration fields

(2.46) \begin{align} c_0 = \bar {c}_0(x) + \bar {c}_1(x,\,y)/\alpha + \cdots, \quad \Gamma _0 = \bar {\Gamma }_0(x) + \bar {\Gamma }_1(x)/\alpha + \cdots, \end{align}

equations (2.38)–(2.42) reduce to the ‘1-D long-wave model’, given by

(2.47a) \begin{align} \bar {c}_{0x} - 2 \epsilon ^2\alpha \bar {c}_{0xx} - \nu (\bar {\Gamma }_0 - \bar {c}_0) &= 0 \quad \text{in} \ \mathcal{D}_1, \end{align}

(2.47b) \begin{align} \frac {3}{4} \beta \bar {\Gamma }_{0x} - \frac {1}{2} \gamma (\bar {\Gamma }_{0x} \bar {\Gamma }_0)_x - \epsilon ^2\delta \bar {\Gamma }_{0xx} + \nu (\bar {\Gamma }_0 - \bar {c}_0) &= 0 \quad \text{in} \ \mathcal{D}_1, \end{align}

(2.47c) \begin{align} \bar {c}_0 - 2 \epsilon ^2\alpha \bar {c}_{0x} +\frac {3}{4} \beta \bar {\Gamma }_0 - \frac {1}{2} \gamma \bar {\Gamma }_{0x} \bar {\Gamma }_0 - \epsilon ^2\delta \bar {\Gamma }_{0x} &= 1 \quad \text{in} \ \mathcal{D}_1, \end{align}

(2.47d) \begin{align} \bar {c}_0 - 2 \epsilon ^2\alpha \bar {c}_{0x} &= 1 \quad \text{in} \ \mathcal{D}_2. \\[6pt] \nonumber \end{align}

The system in (2.47) is solved subject to continuity of bulk surfactant, continuity of bulk-surfactant flux and no flux of interfacial surfactant:

(2.48a) \begin{align} \bar {c}_0(\phi ^-) = \bar {c}_0(\phi ^+), \\[-33pt] \nonumber \end{align}

(2.48b) \begin{align} \bar {c}_0(-\phi ) = \bar {c}_0(2-\phi ), \\[-33pt] \nonumber \end{align}

(2.48c) \begin{align} \bar {c}_0(\phi ^-) - 2\epsilon ^2\alpha \bar {c}_{0x}(\phi ^-) = \bar {c}_0(\phi ^+) - 2\epsilon ^2\alpha \bar {c}_{0x}(\phi ^+), \\[-33pt] \nonumber \end{align}

(2.48d) \begin{align} \bar {c}_0(-\phi ) - 2\epsilon ^2\alpha \bar {c}_{0x}(-\phi ) = \bar {c}_0(2-\phi ) - 2\epsilon ^2\alpha \bar {c}_{0x}(2-\phi ), \\[-33pt] \nonumber \end{align}

(2.48e) \begin{align} \frac {3}{4} \beta \bar {\Gamma }_0(\pm \phi ) - \frac {1}{2} \gamma \bar {\Gamma }_0(\pm \phi ) \bar {\Gamma }_{0x}(\pm \phi ) - \epsilon ^2\delta \bar {\Gamma }_{0x}(\pm \phi ) = 0. \end{align}

The system of ordinary differential equations in (2.47)–(2.48) can be solved numerically using the method outlined in Appendix A of Tomlinson et al. (2023). The drag reduction can be evaluated using (2.44). As illustrated in figure 2, the 1-D long-wave model (2.47)–(2.48) exhibits multiple asymptotic regimes in the $(\alpha, \,\gamma )$ -parameter space, which provides useful context for the results to follow. We highlight regions dominated by Marangoni effects ( $M$ ), interfacial advection ( $A$ ) and interfacial diffusion ( $D$ ); a subscript $E$ denotes cases in which $\Gamma _0(x)$ and $c_0(x,0)$ can be out of equilibrium. We summarise key asymptotic results from Tomlinson et al. (2023a ) that apply in the strong-exchange ( $\bar {c}_0 \approx \bar {\Gamma }_0$ , large $\nu$ ) and weak-exchange ( $\bar {c}_0 \neq \bar {\Gamma }_0$ , small $\nu$ ) limits in Appendix B, in regions labelled with superscript $1D$ . We now proceed to investigate the emergence of 2-D flow structures as bulk diffusion weakens, i.e. for $\alpha$ small ( $\ll 1$ ) or at intermediate values $\sim 1$ .

3. Results

In § 3.1, we investigate the leading-order drag reduction ( ${DR}_0$ ), interfacial and bulk surfactant concentrations ( $\Gamma _0$ and $c_0$ ), as well as the streamwise and wall-normal velocity components ( $u_0$ and $v_0$ ), for strong bulk–surface exchange ( $\nu = 100$ ). We vary the bulk diffusion $(\alpha )$ and surfactant strength ( $\gamma$ ), exploring both strong- and weak-cross-channel-diffusion regimes by solving the 2-D long-wave model (2.38)–(2.42). We focus primarily on the weak-cross-channel-diffusion regime ( $\alpha = \delta \ll 1$ ), as the strong-cross-channel-diffusion regime (§ 2.4) was addressed in Tomlinson et al. (2023a ). For simplicity, we constrain the partition coefficient ( $\beta$ ) and interfacial diffusion ( $\delta$ ), assuming that $\alpha = \delta$ , $\beta = 1$ and $\phi =0.5$ throughout. However, asymptotic solutions are derived in Appendix C for generic variables $\alpha$ , $\beta$ , $\delta$ and $\phi$ , which we use to validate numerical simulations and identify dominant physical mechanisms. In § 3.2, we consider the weak-exchange limit ( $\nu = 0.01$ ). As illustrated in figure 2, in both strong- and weak-exchange limits, we identify regions of parameter space dominated by Marangoni effects ( $M$ ), advection ( $A$ ) and diffusion ( $D$ ). Additionally, we compare numerical solutions of the 2-D long-wave model (2.38)–(2.42) with numerical solutions of the Stokes and advection–diffusion equations (2.15)–(2.21) in the weak-cross-channel-diffusion limit.

3.1. Strong exchange

Figure 3 summarises the surfactant concentration profiles and drag reduction in the strong-exchange limit, computed using the 2-D long-wave model, (2.38)–(2.42). For large values of $\gamma$ , with strong Marangoni effect, the interface is almost immobilised (figure 3 a,b), the $\Gamma _0$ profile is almost linear and the bulk surfactant concentration transitions from a 2-D (figure 3 a) to a 1-D (figure 3 b) distribution as $\alpha$ increases. Surfactant adsorbs onto the interface across the upper half of the plastron and desorbs from it in the lower half. In contrast, for small $\gamma$ with weaker Marangoni effects, interfacial surfactant accumulates near the downstream end of the plastron, lengthening the adsorption region and compressing the desorption region (figure 3 f).

Figure 3. The leading-order drag reduction ( ${DR}_0$ ), bulk ( $c_0$ ) and interfacial ( $\Gamma _0$ ) surfactant concentration fields for $\beta = 1$ , $\epsilon = 0.1$ , $\nu = 100$ and $\phi = 0.5$ , computed using (2.38)–(2.42) when bulk–surface exchange is strong. In the Marangoni-dominated ( $M$ ) region, the SHS is no-slip ( ${DR}_0\ll 1$ ), and in the advection-dominated ( $A$ ) and diffusion-dominated ( $D$ ) regions, the interface is shear-free ( ${DR}_0\approx 1$ ). Fields $c_0$ , $\Gamma _0$ , $\langle c_0 \rangle$ and $c_0(x, \, 0)$ are plotted in (a) $M^{{2D}}$ and (b) $M^{{1D}}$ . (c) Contours of ${DR}_0$ , (d) plots of ${DR}_0$ for different $\gamma$ and (e) plots of ${DR}_0$ for different $\alpha$ , where (c,d,e) are compared with asymptotic predictions (B1b ) in $M^{{1D}}$ , (B2b ) in $D^{{1D}}$ , (C22b ) in $M^{{2D}}$ and (3.2b ) in $A^{{2D}}$ . The dashed line in (e) represents the largest $\gamma$ for which $\Gamma _0(-\phi )=0$ when $\alpha \ll 1$ . Fields $c_0$ , $\Gamma _0$ , $\langle c_0 \rangle$ and $c_0(x, \, 0)$ are plotted in (f) $A^{{2D}}$ (notice the quasi-stagnant-cap profile) and (g) $D^{{1D}}$ , where $x=x_0$ is plotted using (C55).

Figure 3(c–e) illustrates the variation of ${DR}_0$ with $\alpha$ and $\gamma$ . In the present strong-exchange limit, drag reduction decreases as $\alpha$ decreases and the bulk concentration becomes more 2-D. Reduced bulk diffusion promotes the formation of bulk-concentration boundary layers, while slowing fluxes between the bulk and interface. We characterise this transition by comparing the established large- $\alpha$ asymptotes (in regions $M^{{1D}}$ and $D^{{1D}}$ of the $(\alpha, \gamma )$ plane; see (B1), (B2)) with new small- $\alpha$ limits (derived below) in regions $M^{{2D}}$ and $A^{{2D}}$ . We use these limits to characterise the immobilisation of the interface as $\gamma$ increases, highlighting in particular the transition from region $A^{{2D}}$ to region $M^{{2D}}$ for small $\alpha$ . The central asymptotic subregions $M^{{1D}}$ and $A^{{1D}}$ for $1\ll \alpha \ll 1/\epsilon ^2$ in figure 2(a), which formally provide a bridge between the 1-D and 2-D long-wave models, are smoothed out by transitional effects associated with the nominally small geometric parameter $\epsilon$ taking the value 0.1 in figure 3. We now discuss these new regions, starting with $M^{{2D}}$ .

Figure 4. Schematics of the asymptotic structure of the bulk-concentration boundary layer. Weak diffusion ensures that $c_0$ is approximately uniform in the core of the channel, varying primarily in a thin concentration boundary layer near the SHS. Blue (pink) regions illustrate regions where surfactant is drawn from (released into) the bulk onto (from) the interface. $(\textit {a})$ The bulk-concentration boundary layer when Marangoni effects are strong (region $M^{{2D}}$ ) and the surfactant distribution almost immobilises the entire interface. $(\textit {b})$ The bulk-concentration boundary layer when Marangoni effects are weak (region $A^{{2D}}$ ), creating a slip region with low surfactant concentration upstream of a quasi-stagnant region in which interfacial surfactant accumulates. Details of each asymptotic region are provided in Appendix C.

Region $M^{{2D}}$ (figures 2 a and 3 a) emerges when $\alpha \ll 1$ and $\gamma$ is large. In this region, Marangoni effects dominate, with weak bulk diffusion resulting in a bulk-concentration boundary layer near $y=0$ . We provide a scaling argument for ${DR}_0$ as follows, with a detailed derivation given in Appendix C.1 and a schematic of the asymptotic structure in figure 4(a). At the interface, $\Gamma _0$ is in equilibrium with $c_0$ , and the interface is nearly immobile, exhibiting a linear surfactant distribution with a slope of size $1/\gamma$ . The bulk concentration is close to unity, and $\Gamma _0\lt 1$ ( $\Gamma _0\gt 1$ ) in the upstream (downstream) half of the plastron, leading to adsorption (desorption). Perturbations to the bulk-concentration boundary layer are driven by the linear surface concentration distribution, giving the boundary layer a self-similar structure. The thickness of the boundary layer is $\alpha ^{1/3} \ll 1$ (characteristic of a shear flow), resulting in fluxes of size $\alpha c_{0y} \sim \alpha ^{2/3}$ onto and off the interface. Adsorption and desorption are accommodated by weak stretching and compression of the interface, generating smaller contributions to $\Gamma _0$ of size $\alpha ^{2/3}/\gamma ^2$ and to ${DR}_0$ of size $\alpha ^{2/3}/\gamma$ , using (2.44). The constraint $\alpha \ll 1$ ensures that the boundary layer is thin and we require $\gamma \gg 1$ and $\alpha ^{2/3} \ll \gamma$ for the correction to $\Gamma _0$ to be small (although a more precise condition will emerge below). The leading-order surfactant distribution and drag reduction are given by

(3.1a,b) \begin{align} c_0(x, \, 0) \approx \Gamma _0 \approx 1 + \frac {3\beta }{2\gamma }\left (x - \frac {\phi }{5}\right ), \quad {DR}_0 \approx \frac {m_1\alpha ^{2/3}\phi ^{5/3}}{\gamma }, \end{align}

where the coefficient $m_1 \approx 0.79$ is given explicitly in Appendix C.1. The range of validity of the leading-order solution in (3.1) is extended from $\gamma \gg 1$ to $\gamma \sim 1$ in (C22), where the extension must be evaluated numerically. Therefore, we give the simpler formula above, but compare (C22) with numerical results from the 2-D long-wave model at small $\alpha$ in figure 3(c–e), which capture numerical results successfully. In particular, decreasing $\gamma$ leads to a steepening of the interfacial surfactant gradient until $\Gamma _0$ approaches zero at the upstream contact line. This transition provides a lower bound on region $M^{{2D}}$ , namely $\Gamma _0(-\phi )=0$ when $\gamma = 9\phi \beta /5$ (see vertical dashed lines in figure 3 e). Hence, we plot asymptotic predictions up to this limit in figure 3(e).

When $\alpha \ll 1$ and $\gamma$ is sufficiently small, we transition into region $A^{{2D}}$ (figures 2 a and 3 f), where advection dominates Marangoni effects at the interface, while interacting with a strongly coupled bulk-concentration boundary layer. We distinguish two primary regions along the interface: (i) a slip region, $x\in [-\phi, \,x_0)$ (for some $x_0$ ), at the upstream end of the interface, where Marangoni effects are weak and surfactant adsorbs from the bulk; and (ii) a quasi-stagnant region, $x\in (x_0,\,\phi ]$ , at the downstream end of the interface, where surfactant gradients decelerate the interface through Marangoni effects without completely immobilising the interface (in contrast with a classic stagnant cap), while desorbing the accumulated surfactant. We now provide a scaling argument for ${DR}_0$ that is supported by a detailed derivation in Appendix C.2; a schematic of the asymptotic structure is given in figure 4(b). In the slip region, of length $x_0+\phi =O(1)$ , $c_0$ varies by $O(1)$ across a boundary-layer thickness $y\sim \alpha ^{1/2}$ , yielding a flux per unit length of size $\alpha c_{0y}\sim \alpha ^{1/2}$ onto the interface. This allows $\Gamma _0$ to grow along the slip region from zero to $O(\alpha ^{1/2})$ . Although $c_0$ is effectively zero at the interface at leading order, it is coupled to $\Gamma _0$ through an $O(\alpha ^{1/2})$ correction. The slip region therefore delivers an interfacial flux of size $O(\alpha ^{1/2})$ at $x\approx x_0$ , which is carried through short regions across which the interface decelerates rapidly at the leading edge of the quasi-stagnant region (i.e. at $x\approx x_0$ ). The sizes of the short regions are indicated in figure 4(b). Across the quasi-stagnant region, of length $L=\phi -x_0$ (to be determined), $\Gamma _0$ is approximately linear with a slope of size $1/\gamma$ , setting $c_0$ and $\Gamma _0$ to be of size $L/\gamma$ . In the bulk, diffusion balances shear-flow transport over a vertical length $y \sim (\alpha L)^{1/3}$ and therefore $\alpha c_{0y} \sim (\alpha L)^{2/3}/\gamma$ . Integrated over $L$ , this accommodates the $O(\alpha ^{1/2})$ interfacial flux, resulting in $L \sim \gamma ^{3/5}/\alpha ^{1/10}$ . Consequently, the drag reduction arising from the surfactant adsorbed and desorbed in the slip and quasi-stagnant regions is proportional to ${DR}_0 \sim \gamma \Delta \Gamma \sim L$ .

In Appendix C.2, we present a more detailed analysis showing that

(3.2a--c) \begin{align} c_0(x, \, 0) \approx \Gamma _0 \approx \begin{cases} \dfrac {4(\alpha (x+\phi ))^{1/2}}{\sqrt {3\pi }\beta } \quad \text{for} \quad x\in [-\phi, \, x_0), \\ \dfrac {3\beta }{2\gamma }(x-x_0) \quad \text{for} \quad x\in (x_0,\,\phi ], \end{cases} \quad {DR}_0 = 1 - \dfrac {a_1 \gamma ^{3/5}}{\alpha ^{1/10}},\nonumber \\ { u_0(x, \, 0) \approx \begin{cases}3/4 \quad \text{for} \quad x\in [-\phi, x_0), \\ \dfrac {3}{4} \left [G\left (\dfrac {x-x_0}{\gamma \alpha ^{1/2}}\right )\right ]^{-1} \text{for} \quad x\approx x_0, \\[6pt] \dfrac {3\alpha ^{2/3} C^{\prime}(0)}{5\beta } \Bigg ((x-x_0)^{2/3} - \dfrac {(\phi -x_0)^{5/3}}{x-x_0}\Bigg ) \quad \text{for} \ x\in (x_0,\,\phi ]. \end{cases} } \\[-6pt] \nonumber \end{align}

The coefficient $a_1 = (-5\lambda /(6C^{\prime}(0)))^{3/5}/(2\phi )$ in (3.2b ) where $\lambda = (16(x_0+\phi )/ (3\pi \beta ^2))^{1/2}$ ; $C^{\prime}(0) \approx -0.92$ in (3.2c ) is given explicitly in Appendix C.2; the function $G$ in (3.2 $c$ ) is given implicitly in (C41). The length of the quasi-stagnant cap is

(3.3) \begin{align} L=\phi -x_0= \frac {2 \phi a_1 \gamma ^{3/5}}{\alpha ^{1/10}}. \end{align}

In (3.2), $x\in [-\phi, x_0)$ , $x\approx x_0$ and $x\in (x_0,\phi ]$ are shorthand for the slip, deceleration and quasi-stagnant regions shown in figure 4(b); the leading-order approximations in (3.2) match together as explained in Appendix C. Equation (3.2c ) highlights the role of interfacial compression in the quasi-stagnant region, accommodating desorption.

Asymptotes for ${DR}_0$ , $c_0(x,\,0)$ and $\Gamma _0$ in figure 3(c–f) are evaluated numerically for a given $\alpha$ and $\gamma$ and capture the behaviour exhibited by the numerical simulations of the 2-D long-wave model. Increasing $\gamma$ from low values causes lengthening of the quasi-stagnant region, pushing $x_0$ towards the upstream contact line and reducing $DR_0$ . Figure 3(e) shows how (3.2b ) applies almost until $x_0$ approaches $-\phi$ . For small $\alpha$ , $c_0(-\phi, \, 0) \approx \Gamma _0(-\phi ) \approx 0$ for $\gamma \ll 1$ , causing the $c_0$ field to develop a singular first derivative at the upstream contact point. This concentration-gradient singularity significantly impacts the numerical solution (in terms of both the required resolution and run time) of the 2-D long-wave model and COMSOL simulations (below).

Figure 5. The leading-order streamwise velocity field ( $u_0$ ) and surfactant concentration fields ( $c_0$ and $\Gamma _0$ ) for $\beta = 1$ , $\nu = 100$ , $\epsilon = 0.1$ and $\phi = 0.5$ evaluated using (2.38)–(2.42) (lines) and COMSOL simulations (2.15)–(2.21) (symbols), when bulk–surface exchange is strong. (a,c) Plots of $c_0$ and $u_0$ , respectively, for $\alpha = 0.1$ and $\gamma = 1$ , when the flow is in the Marangoni-dominated region with weak cross-channel diffusion ( $M^{{2D}}$ ). (b,d) Plots of $c_0$ and $u_0$ , respectively, for $\alpha = 0.1$ and $\gamma = 0.1$ , when the flow is in the advection-dominated region with weak cross-channel diffusion ( $A^{{2D}}$ ), where $x=x_0$ is plotted using (C55).

Further insight into the $M^{{2D}}$ and $A^{{2D}}$ solutions is provided in figure 5, which compares $u_0$ , $c_0$ and $\Gamma _0$ computed using the 2-D long-wave model, (2.38)–(2.42), with the Stokes and advection–diffusion equations in COMSOL, (2.15)–(2.21), with the same parameters as the examples shown in figure 3(a, f). For the solution in the $M^{{2D}}$ region, $u_0$ exhibits a slight deviation from no slip at the liquid–gas interface (figure 5 c), generating a weak vertical flow towards the interface. The concentration profiles in the two simulations match well (figure 5 a), with the leading-order drag reduction predicted from the 2-D long-wave model, ${DR}_0 = 0.106$ , in close agreement with the COMSOL simulations, ${DR}_{NS} = 0.109$ . For the solution in the region $A^{{2D}}$ , depicted in figure 5(b,d), the liquid–gas interface is nearly shear-free and $u_0\approx 3/4$ in the slip region, where the interfacial surfactant concentration is low, before falling towards zero across the quasi-stagnant region, where the gradient of $\Gamma _0$ is large. The asymptotic predictions (3.2 c) capture respectively the slight decrease of the slip velocity in the deceleration region and then the rapid fall towards zero across the quasi-stagnant region; the shift parameter in (C41) was chosen as $X_0 = 0.5$ to minimise the error against the COMSOL simulations.

Figure 5 illustrates the distinction between a classical stagnant cap of a strictly insoluble surfactant, for which $u_0$ vanishes abruptly in the cap region where $\Gamma _0$ has a large gradient, and the present quasi-stagnant-cap structure, where desorption of the soluble surfactant in the bulk boundary layer is mediated by molecular diffusion, and which leads to a small but non-zero slip velocity (3.2 $c$ ) across the quasi-stagnant-cap region. The predictions for the drag reduction and concentration field from the 2-D long-wave model and numerical simulations in COMSOL agree, with ${DR}_0 = 0.745$ and ${DR}_{NS} = 0.747$ . Both methods capture weak interfacial stretching (allowing adsorption) in the slip region ( $u_{0x}\gt 0$ for $x\in [-\phi, \, x_0]$ in figure 5 d) and stronger compression (allowing desorption) in the quasi-stagnant region ( $u_{0x}\lt 0$ for $x\in [x_0,\,\phi ]$ in figure 5 d). There are discrepancies in $u_0$ at both the upstream and downstream contact lines. At the upstream contact line, the 2-D long-wave theory does not capture the Stokes-flow region, where $u_0$ transitions sharply from zero to the slip velocity. At the downstream contact line, the retention of weak streamwise diffusion in the 2-D long-wave model leads to a small but non-zero value of $u_0$ , ensuring zero surfactant flux. In contrast, the small- $\alpha$ asymptotic solution, which neglects streamwise diffusion, correctly predicts that the slip velocity vanishes at the contact line. Since the small- $\alpha$ solution represents a limiting case of the 2-D long-wave model as streamwise diffusion vanishes, the two solutions align in the limit $\epsilon \to 0$ , where the 2-D long-wave model’s prediction for $u_0(\phi, 0)$ also tends to zero at the downstream edge. Further features of the surface velocity profile at the tip of the quasi-stagnant region, which explain the dramatic thickening of the boundary layer in figures 3(f) and 4(b), are discussed in Appendix C.2.

3.2. Weak exchange

Figure 6. The leading-order drag reduction ( ${DR}_0$ ), bulk ( $c_0$ ) and interfacial ( $\Gamma _0$ ) surfactant concentration fields for $\beta = 1$ , $\epsilon = 0.1$ , $\nu = 0.01$ and $\phi = 0.5$ , computed using (2.38)–(2.42) when bulk–surface exchange is weak. In the Marangoni-dominated ( $M_E$ ) region, the SHS is mostly no slip ( ${DR}_0\ll 1$ ), and in the advection-dominated ( $A_E$ ) and diffusion-dominated ( $D_E$ ) regions, the interface is mostly shear-free ( ${DR}_0\approx 1$ ). Fields $c_0$ , $\Gamma _0$ , $\langle c_0 \rangle$ and $c_0(x, \, 0)$ are plotted in (a) $M^{{2D}}_E$ and (b) $M^{{1D}}_E$ . (c) Contours of ${DR}_0$ , (d) plots of ${DR}_0$ for different $\gamma$ and (e) plots of ${DR}_0$ for different $\alpha$ , where (c,d,e) are compared to (B4b ) in $M^{{1D}}_E$ , (B5b ) in $D^{{1D}}_E$ , (B4b ) in $M^{{2D}}_E$ and (3.4b ) in $A^{{2D}}_E$ . The dashed line in (e) represents the largest $\gamma$ for which $\Gamma _0(-\phi )=0$ when $\alpha \ll 1$ . Fields $c_0$ , $\Gamma _0$ , $\langle c_0 \rangle$ and $c_0(x, \, 0)$ are plotted in (f) $A^{{2D}}_E$ (notice the classical stagnant-cap profile) and (g) $D^{{1D}}_E$ (where streamwise diffusion has smoothed the stagnant cap), where $x=x_0$ is plotted using (C62).

When bulk–surface exchange is weak ( $\nu = 0.01$ ), 2-D effects are again prevalent in bulk-concentration profiles at small $\alpha$ (figure 6). However, the impact of bulk boundary layers on the drag reduction is more modest than in the strong-exchange limit, in that they have less of an effect on the leading-order drag reduction. Solutions at large $\gamma$ again have approximately linear interfacial surfactant profiles that immobilise the interface (figure 6 a,b). These profiles are almost decoupled from the bulk. They act as a weak sink/source combination driving adsorption/desorption. Solutions with weak Marangoni effects (small $\gamma$ ) show the formation of a more classical stagnant cap at small $\alpha$ (figure 6 f), which is smoothed by streamwise diffusion as $\alpha =\delta$ increases (figure 6 g).

Figure 6(c–e) illustrates the dependency of ${DR}_0$ on $\alpha$ and $\gamma$ , evaluated against existing limits $M_E^{{1D}}$ and $D_E^{{1D}}$ (see (B4), (B5)) for large $\alpha$ and introducing new limits $M_E^{{2D}}$ and $A_E^{{2D}}$ at small $\alpha$ . Importantly, weak coupling between the bulk and the interface suppresses the dependence of $DR_0$ on $\alpha$ in the limit of very small bulk diffusivity, making the drag-reduction contours horizontal in figure 6(c) for $\alpha \to 0$ .

The following scaling argument applies to region $A_E^{2D}$ (figures 2 b and 6 f), where Marangoni effects and bulk diffusion are weak. Here, the interfacial concentration is approximated by the classical distribution for nearly insoluble surfactant, which is close to zero along the interface except at the downstream end where it forms a stagnant cap of length $L$ and slope $1/\gamma$ . Thus, $\Gamma _0$ is of magnitude $L/\gamma$ , and its integral over the whole interface is $O(L^2/\gamma )$ , while the bulk concentration remains close to unity. The surfactant flux-balance constraint (2.29c ), which matches net absorption to net desorption over the interface, therefore requires $L^2/\gamma =O(1)$ , implying that $\Gamma _0$ is $O(\gamma ^{-1/2})$ and $DR_0$ is $O(\gamma ^{1/2})$ . This scaling argument, which mirrors the 1-D case, is supported by asymptotic calculations in Appendix C.3, which yield

(3.4a--c) \begin{align} c_0 \approx 1, \quad \Gamma _0 \approx \begin{cases} \ 0 \quad \text{for} \quad x\in [-\phi, \, x_0], \\ \ \dfrac {3\beta }{2\gamma }(x-x_0) \quad \text{for} \quad x\in [x_0,\,\phi ], \end{cases} \quad{DR}_0 \approx 1 - \left (\dfrac {2\gamma }{3\phi \beta }\right )^{1/2}, \end{align}

where $\phi - x_0 = (8 \phi \gamma /(3\beta ))^{1/2}$ is the length of the stagnant cap. Bulk diffusion has a higher-order effect, influencing exchange via adsorption and desorption, and hence weak stretching and compression of the interface. In figure 6(c–f), dashed black asymptotes for $c_0$ , $\Gamma _0$ and ${DR}_0$ match those results from the 2-D long-wave model, almost up to $\gamma = 3\beta \phi /2$ in figure 6(e) (see vertical dashed lines), which denotes the conditions where the stagnant cap reaches $x=-\phi$ and the slip region of this nonlinear distribution vanishes. The disappearance of the stagnant-cap distribution with increasing $\gamma$ is linked to a rapid decrease in $DR_0$ , owing to the loss of the slip region. The same rapid transition is recovered if we approach region $A^{2D}$ from region $M_E^{2D}$ , using (B4 b), i.e. for decreasing $\gamma$ .

Figure 7. The leading-order streamwise velocity field ( $u_0$ ) and surfactant concentration fields ( $c_0$ and $\Gamma _0$ ) for $\beta = 1$ , $\nu = 0.01$ , $\epsilon = 0.1$ and $\phi = 0.5$ evaluated using (2.38)–(2.42) (lines) and COMSOL simulations (2.15)–(2.21) (symbols) when bulk–surface exchange is weak. (a,c) Plots of $c_0$ and $u_0$ , respectively, for $\alpha = 0.1$ and $\gamma = 1$ , when the flow is in the Marangoni-dominated region with weak cross-channel diffusion ( $M^{{2D}}_E$ ). (b,d) Plots of $c_0$ and $u_0$ , respectively, for $\alpha = 0.1$ and $\gamma = 0.1$ , when the flow is in the advection-dominated region with weak cross-channel diffusion ( $A^{{2D}}_E$ ) and $x = x_0$ is plotted using (C62).

A comparison between the solutions evaluated using the 2-D long-wave model, (2.38)–(2.42), and COMSOL simulations is given in figure 7, focusing on the examples discussed in figure 6(a,f). Figure 7(a,c) shows $c_0$ , $\Gamma _0$ and $u_0$ in region $M_E^{{2D}}$ . The slip velocity is approximately zero for both symbols and lines, as the interfacial surfactant gradient effectively immobilises the interface. The 2-D long-wave model predicts ${DR}_0 = 0.027$ , while the numerical simulations yield ${DR}_{NS} = 0.056$ (due to the 4 % difference in $\Delta \Gamma _0$ shown in figure 7 a). Figure 7(b,d) depicts stagnant-cap solutions in region $A_E^{{2D}}$ . The upstream end of the liquid–gas interface ( $x\in [-\phi, \, x_0]$ ) is shear-free, resulting in a slip velocity $u_0 = 3/4$ (substituting $\Gamma _{0x} = 0$ into (2.32a ), (2.36a ), $u_0 = \tilde {U} p_{0x}$ , where $\tilde {U}(0) = - 2$ and $p_{0x} = - 3 / 8$ ). Conversely, the downstream end of the liquid–gas interface ( $x\in [x_0,\, \phi ]$ ) does not exhibit slip, leading to a near-zero streamwise velocity at the downstream stagnation point. The transitional and downstream regions induce a strong wall-normal flow that drives bulk surfactant into the core of the channel. We find ${DR}_0 = 0.650$ and ${DR}_{NS} = 0.662$ , demonstrating agreement between the predictions of the 2-D long-wave model and COMSOL simulations, validating the model’s accuracy under a wide range of conditions.

3.3. Comparison with experiments

To facilitate comparison with experiments, in table 2, we convert ${DR}_0$ to slip lengths ( $\lambda _e$ ) using the non-dimensionalisation in Temprano-Coleto et al. (2023). Temprano-Coleto et al. (2023) proposed that for microchannel applications, $\lambda _e \sim (\phi /\epsilon )^2 / L_m^2$ (using our notation), where $\phi /\epsilon = \phi \hat {P}/\hat {H}$ is half the length of the liquid–gas interface divided by the channel height, $L_m^2 = (\hat {n}\hat {R}\hat {T} \hat {K}_a^2 \hat {c}_0)/(\hat {D}\hat {\mu }\hat {K}_d^2)$ is the squared mobilisation length and $\lambda _e \sim (\phi / \epsilon ) \lambda _0$ in (2.45). The slip-length formulae in Landel et al. (2020) and Sundin & Bagheri (2022) have the same dependencies in this limit. We find that, in all regimes, the slip length shares the same dependence on $1/L_m^{2}$ , with quadratic dependence on $\phi /\epsilon$ only when $1 \ll \alpha \ll 1/\epsilon ^{2}$ and $\gamma \gg 1$ (region $M^{{1D}}$ ). However, when $\alpha \ll 1$ and $\gamma \gg \alpha ^{1/6}$ , $\lambda _e$ has a stronger (8/3 power) dependence on $\phi /\epsilon$ (region $M^{{2D}}$ ), and when $\alpha \gg 1/\epsilon ^2$ and $\gamma \gg \epsilon ^2 \alpha$ , we find a weaker (linear) dependence (region $M^{{1D}}$ ). The remaining asymptotic solutions for ${DR}_0$ are summarised in figure 2. When $1-{DR}_0 \ll 1$ , the slip length simplifies to $\lambda _e = (\Delta p_I - \Delta p_U) / (\epsilon \Delta p_I \Delta p_U)$ ( $A^{{1D}}$ and $A^{{2D}}$ ). Hence, outside of the specific microchannel applications considered in Temprano-Coleto et al. (2023), whose predictions agree with those in region $M^{{1D}}$ (as discussed below), we have found new parameters in regions $M^{{2D}}$ that influence surfactant impairment of SHS drag reduction. As diffusion weakens, the deviation between the drag-reduction predictions in $M^{{1D}}$ and $M^{{2D}}$ will increase.

Table 2. Summary of the leading-order drag reduction $({DR}_0)$ scaling in regions $M^{{2D}}$ , $M^{{1D}}$ and $M^{{1D}}_E$ and the corresponding leading-order slip length ( $\lambda _e$ ) scaling using the non-dimensionalisation in Temprano-Coleto et al. (2023), in terms of the mobilisation length $L_m$ . Hats denote dimensional quantities. The drag reduction is converted to the slip length using (2.45) for ${DR}_0 \ll 1$ and $\lambda _e \sim (\phi /\epsilon ) \lambda _0$ , where $\hat {L}_m^2 = (\hat {n}\hat {R}\hat {T} \hat {H}^2 \hat {K}_a^2 \hat {c}_0)/(\hat {D}\hat {\mu }\hat {K}_d^2)$ , $\phi /\epsilon = \phi \hat {P}/\hat {H}$ , $L_m = \hat {L}_m/\hat {H}$ and $b = \hat {Q}/\hat {D}$ .

For laboratory experimental studies documented in the literature, we can use our 2-D long-wave model to estimate surfactant effects via the drag reduction, employing parameters typical of microchannel applications in laminar flows. Estimates for these parameters, based on the surfactant SDS, using models and experimental data from Temprano-Coleto et al. (2023), are given in table 3. The remaining parameters $\hat {P}$ , $\hat {Q}$ , $\hat {H}$ and $\hat {K}$ vary across experiments. Peaudecerf et al. (2017) use $\alpha \approx 1$ and $4.7\times 10^{1} \lessapprox \gamma \lessapprox 2.9\times 10^4$ (at the $M^{{1D}}/M^{{2D}}$ boundary; figure 2) and Temprano-Coleto et al. (2023) use $\alpha \approx 25$ and $10^{-1} \lessapprox \gamma \lessapprox 3.7\times 10^{2}$ (at the $A^{{1D}}/M^{{1D}}$ boundary; figure 2). Qualitatively, all of these studies report surfactant effects, which is consistent with their location in the non-dimensional parameter space. All of these studies operate within the strong-exchange regime for the values of $\hat {K}_a$ and $\hat {K}_d$ given in table 3. However, as $\hat {K}_a$ and $\hat {K}_d$ are merely estimates based on fitting to microchannel experiments, weak exchange could be achieved, for instance, when $\hat{K}_a = 3.4 \times 10^{-8}\,{\textrm{m}\, \textrm{s}}^{-1}$ and $\hat{K}_d = 7.5\times 10^{-5}\,{\textrm{s}}^{-1}$ . Quantitatively, we compare predicted values of slip velocity using our theoretical model with experimental results in table 4. Using the experimental parameters from Temprano-Coleto et al. (2023), where $\alpha \approx 25$ and $\epsilon \approx 0.002$ , we confirm $1\ll \alpha \ll 1/\epsilon ^2$ and they are within the regime $M^{{1D}}$ . From (B1 b), we calculate ${DR}_0$ , convert to $\lambda _e$ using table 2 and then to $u_{Ic}/u_{Ic}^{\textrm {clean}}$ following Temprano-Coleto et al. (2023). Agreement is observed across $\phi /\epsilon = 253$ , 418 and 587, which is remarkable as our theory does not include any fitting parameters. The deviation from Temprano-Coleto et al. (2023) at $\phi /\epsilon = 762$ could be attributed to experimental error, which is consistent with the expected behaviour of increasing slip with grating length.

Table 3. A summary of the parameters in the dimensional problem, (2.1)–(2.12), with their values based on the surfactant SDS, models and experimental data from Temprano-Coleto et al. (2023).

Table 4. A comparison between the slip velocity normalised by the clean (surfactant-free) value, $u_{Ic}/u_{Ic}^{\textrm {clean}}$ , predicted using the current model (2.38)–(2.42), and experimental data from Peaudecerf et al. (2017) and Temprano-Coleto et al. (2023).The parameter $\phi /\epsilon$ is half the length of the liquid–gas interface divided by the channel height.

4. Discussion

In this study, we develop a long-wave theory to analyse the behaviour of a 2-D laminar pressure-driven channel flow contaminated with soluble surfactant. The channel is bounded by a SHS with streamwise-periodic grooves and a solid wall. We linearise the equation of state and adsorption–desorption kinetics, solving Stokes and advection–diffusion equations in the long-wave limit. Our investigation focuses primarily on regimes where concentration gradients in the streamwise and cross-channel directions are comparable. By numerically solving the 2-D long-wave model, we delineate asymptotic regions within the parameter space. Under conditions of strong cross-channel diffusion, our findings align with the 1-D results in Tomlinson et al. (2023a ). Conversely, when cross-channel diffusion is weak, we reveal new regions of the parameter space, where the drag reduction exhibits a non-trivial dependence on the thickness of the bulk-concentration boundary layer and surfactant strength (discussed below). We derive asymptotic solutions for the boundary-layer problem, validating them against both the numerical solution of the 2-D long-wave model (2.38)–(2.42) and the full Stokes and advection–diffusion equations (2.15)–(2.21) in COMSOL. These complementary methodologies, ranging from asymptotic solutions to long-wave numerical simulations, give physical insight for applications where numerical simulations can be computationally prohibitive.

We have explored how the interfacial concentration ( $\Gamma _0$ ), bulk concentration ( $c_0$ ) and leading-order drag reduction ( ${DR}_0$ ) are influenced by bulk diffusion ( $\alpha$ ), surface advection ( $\beta$ ), surfactant strength ( $\gamma$ ), interfacial diffusion ( $\delta$ ) and bulk–surface exchange ( $\nu$ ) (figures 3 and 6); the dimensionless parameters $\alpha$ , $\beta$ , $\gamma$ , $\delta$ and $\nu$ are given in table 1. In cases where cross-channel diffusion dominates, we recover the Marangoni-dominated ( $M^{{1D}}$ and $M^{{1D}}_E$ ), advection-dominated ( $A^{{1D}}$ and $A^{{1D}}_E$ ) and diffusion-dominated ( $D^{{1D}}$ and $D^{{1D}}_E$ ) regions identified in Tomlinson et al. (2023a ) (see figure 2). However, under weak cross-channel diffusion, we reveal new subregions dominated by Marangoni effects ( $M^{{2D}}$ and $M^{{2D}}_E$ ) and advection ( $A^{{2D}}$ and $A^{{2D}}_E$ ). In these regions, a concentration boundary layer forms in the bulk, influencing the drag reduction when bulk–surface exchange is strong. When Marangoni effects dominate, the interface is immobilised and the surfactant distribution along it is approximately linear. The asymptotic solutions for ${DR}_0$ are summarised in figure 2.

When both bulk diffusion and Marangoni effects are weak, part of the interface is shear-free and the surfactant distribution forms either what we call a quasi-stagnant cap (when bulk–interface exchange is strong) or a classical stagnant cap (when exchange is weak). The classical stagnant cap for an insoluble surfactant-contaminated SHS has been described in Baier & Hardt (2021) and Mayer & Crowdy (2022). The quasi-stagnant cap has received much less attention and has a particularly intricate structure, illustrated in figure 4(b). Surface stretching and compression accommodate weak adsorption from, and desorption to, the bulk across thin concentration boundary layers, ‘remobilising’ the interface (analogous to a mechanism described by Crowdy et al. (2023) at low $\textit {Pe}$ – and so without the involvement of bulk boundary layers – for a linear extensional flow). The upstream ‘slip’ region has a bulk-concentration boundary layer of thickness $O(\alpha ^{1/2})$ , which is to be expected above an almost fully mobile interface. Weak stretching draws surfactant from the bulk to the interface via diffusive adsorption. The slip region transitions abruptly to a quasi-stagnant region, which grows to a thickness $O(\gamma ^{1/5}\alpha ^{3/10})$ . This differs from a classical stagnant cap (of a fully insoluble surfactant) in having weak surface compression at a rate that balances diffusive desorption. The transition from the mobile to the quasi-immobile interface takes place across two nested regions at the tip of the quasi-stagnant cap: a short deceleration region, where the surface velocity falls abruptly, displacing the bulk boundary layer upwards towards the core; and a slightly longer transition region, across which shear in the boundary layer balances weakening surface advection (figure 4 $b$ ). Shear then dominates in the bulk-concentration boundary layer along the quasi-stagnant region.

It is likely that the physical balances arising in these regions may emerge in other flow configurations involving soluble surfactant transport near confined interfaces at high Péclet numbers, such as the cap forming at the rear of a rising drop or bubble. For example, computations by Oguz & Sadhal (1988) and Tasoglu et al. (2008) revealed quasi-stagnant caps at large but finite $\textit {Pe}$ . Interfacial flux balances between diffusion-limited adsorption (slip) and desorption (stagnant region) were given by Harper (2004) and Palaparthi et al. (2006), treating the size of the cap as a parameter. Here, we determine the size of the cap by solving the bulk-concentration boundary layer, which leads (for example) to the interfacial surfactant concentration being of size $\beta /(\gamma ^{2/5}\alpha ^{1/10})$ in the $A^{{2D}}$ regime.

In summary, our study compares asymptotic and numerical solutions for a 2-D laminar pressure-driven channel flow, confined by a streamwise-periodic SHS and a solid wall, and contaminated with soluble surfactant. While numerical solutions demand significant computational resources, our asymptotic solutions provide a cost-effective alternative, particularly in regimes at high Péclet numbers where accurate numerical solutions of very thin boundary layers are computationally very demanding. These asymptotic solutions offer accurate predictions devoid of empirical fitting coefficients and provide physical insights into the mechanisms governing drag reduction across a large part of the parameter space, including in particular high-Péclet-number flow regimes.