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

Samuel D. Tomlinson, Frédéric Gibou, Paolo Luzzatto-Fegiz, Fernando

Temprano-Coleto, Oliver E. Jensen, Julien R. Landel

► To cite this version:

Samuel D. Tomlinson, Frédéric Gibou, Paolo Luzzatto-Fegiz, Fernando Temprano-Coleto, Oliver E. Jensen, et al.. Drag reduction in surfactant-contaminated superhydrophobic channels at high Péclet numbers. Journal of Fluid Mechanics, In press. hal-05034999

HAL Id: hal-05034999

https://hal.science/hal-05034999v1

Submitted on 15 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Drag reduction in surfactant-contaminated superhydrophobic channels at high Péclet numbers

Samuel D. Tomlinson<sup>1, 2†</sup>, Frédéric Gibou<sup>3</sup>, Paolo Luzzatto-Fegiz<sup>3</sup>, Fernando Temprano-Coleto<sup>4, 5</sup>, Oliver E. Jensen<sup>6</sup> and Julien R. Landel<sup>6, 7</sup>

(Received xx; revised xx; accepted xx)

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 to 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.

**Key words:** Marangoni convection, drag reduction, microfluidics

† Email address for correspondence: sdt50@cam.ac.uk

<footnote n="1"><sup>1</sup>Department of Engineering, University of Cambridge, Trumpington Street, Cambridge, CB2 1PZ, UK</footnote>

<footnote n="2"><sup>2</sup>Centre for Climate Repair, Centre for Mathematical Sciences, Wilberforce Road, Cambridge, CB3 0WA</footnote>

<footnote n="3"><sup>3</sup>Department of Mechanical Engineering, University of California, Santa Barbara, CA 93106, USA</footnote>

<footnote n="4"><sup>4</sup>Andlinger Center for Energy and the Environment, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="6"><sup>6</sup>Department of Mathematics, University of Manchester, Oxford Road, Manchester, M13 9PL, UK</footnote>

<footnote n="7"><sup>7</sup>Universite Claude Bernard Lyon 1, Laboratoire de Mécanique des Fluides et d'Acoustique (LMFA), UMR5509, CNRS, Ecole Centrale de Lyon, INSA Lyon, 69622 Villeurbanne, France</footnote>
# 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 to 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 gathered significant attention in the academic literature (Schönecker et al. 2014; Park 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 (Pe), forming concentration boundary layers when 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 (2D) channel flow featuring streamwise-periodic SHSs using COMSOL, employing surfactant properties that are characteristic of sodium dodecyl sulfate (SDS). Temprano-Coleto et al. (2023) extended this methodology to encompass a three-dimensional (3D) channel flow featuring streamwise and spanwise-periodic SHSs, while Sundin & Bagheri (2022) investigated 2D channel flow laden with surfactants over streamwise-periodic liquid-infused surfaces (LISs). 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 to
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 an 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 an 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 will term a quasi-stagnant cap.

Scaling theories have been developed to analyse the slip length over SHSs and LISs 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 characterizes 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 spatially-1D 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 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 will 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 1D long-wave model first breaks down due to 2D 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 SHS. 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 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 2D problem. We non-dimensionalise the fluid and surfactant equations in the long-wave limit and derive 1D and 2D 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 2D model

We explore a steady 2D 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{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
<figure><img src="image_2.png" /><figcaption>FIGURE 1. A 2D laminar pressure-driven channel flow transporting soluble surfactant confined between a streamwise-periodic SHS and 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{D}_1$ and $\hat{D}_2$ respectively, as specified in (2.1).</figcaption></figure>

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,

$$ \hat{D}_1 = \{\hat{x} \in [-\phi \hat{P}, \phi \hat{P}] \} \times \{\hat{y} \in [0, 2\hat{H}]\}, \qquad (2.1a) $$

$$ \hat{D}_2 = \{\hat{x} \in [\phi \hat{P}, (2-\phi) \hat{P}]\} \times \{\hat{y} \in [0, 2\hat{H}]\}, \qquad (2.1b) $$

as illustrated in figure 1.

A comprehensive discussion of the equations governing fluid and surfactant behaviour in a 3D 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 2D. Within domains $\hat{D}_1$ and $\hat{D}_2$, we have Stokes equations with an advection-diffusion equation for bulk surfactant

$$ \hat{\nabla} \cdot \hat{u} = 0, \quad \hat{\mu} \hat{\nabla}^2 \hat{u} - \hat{\nabla} \hat{p} = 0, \quad \hat{D} \hat{\nabla}^2 \hat{c} - \hat{u} \cdot \hat{\nabla} \hat{c} = 0, \qquad (2.2a-c) $$

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

$$ \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, \\ \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, \qquad (2.3a-d) $$

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,

$$ \hat{u}\hat{\Gamma} - \hat{D}_I\hat{\Gamma}_{\hat{x}} = 0. \qquad (2.4) $$

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

$$ \hat{u} = 0, \quad \hat{v} = 0, \quad \hat{c}_{\hat{y}} = 0. \qquad (2.5a-c) $$

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

$$ \hat{q}((\phi \hat{P})^{-}, \hat{y}) = \hat{q}((\phi \hat{P})^{+}, \hat{y}), \quad \hat{q}(-\phi \hat{P}, \hat{y}) = \hat{q}((2 - \phi) \hat{P}, \hat{y}), \qquad (2.6) $$

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

$$ \frac{d}{d\hat{x}} \int_{0}^{2\hat{H}} (\hat{u}\hat{c} - \hat{D}\hat{c}_{\hat{x}}) d\hat{y} - (\hat{K}_{d}\hat{\Gamma} - \hat{K}_{a}\hat{c}(\hat{x}, 0)) = 0, \qquad (2.7) $$

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

$$ \int_{0}^{2\hat{H}} (\hat{u}\hat{c} - \hat{D}\hat{c}_{\hat{x}}) d\hat{y} + (\hat{u}(\hat{x}, 0)\hat{\Gamma} - \hat{D}_{I}\hat{\Gamma}_{\hat{x}}) = \hat{K}. \qquad (2.8) $$

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

$$ \int_{-\phi \hat{P}}^{\phi \hat{P}} (\hat{K}_d \hat{\Gamma} - \hat{K}_a \hat{c}(\hat{x}, 0)) d\hat{x} = 0. \qquad (2.9) $$

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

$$ \int_{0}^{2\hat{H}} (\hat{u}\hat{c} - \hat{D}\hat{c}_{\hat{x}}) d\hat{y} = \hat{K}. \qquad (2.10) $$

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

$$ \hat{Q} = \int_{0}^{2\hat{H}} \hat{u} d\hat{y}. \qquad (2.11) $$

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} = \langle \hat{p} \rangle(-\phi \hat{P}) - \langle \hat{p} \rangle((2 - \phi) \hat{P}) > 0$, where $\langle \cdot \rangle = \int_{0}^{2\hat{H}} \cdot d\hat{y}/(2\hat{H})$. We define the normalised drag reduction

$$ DR = \frac{\Delta \hat{p}_I - \Delta \hat{p}}{\Delta \hat{p}_I - \Delta \hat{p}_U}, \qquad (2.12) $$

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

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

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

$$
\mathcal{D}_1 = \{x \in [-\phi, \phi]\} \times \{y \in [0, 2]\}, \qquad (2.14a)
$$

$$
\mathcal{D}_2 = \{x \in [\phi, 2-\phi]\} \times \{y \in [0, 2]\}. \qquad (2.14b)
$$

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

$$
\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})/Pe - \epsilon u c_x - \epsilon v c_y = 0, \tag{2.15a-d}
\end{align}
$$

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

$$
\begin{align}
& u_y - \epsilon M a \Gamma_x = 0, \quad v = 0, \quad c_y - D a (c - \Gamma) = 0, \nonumber \\
& \epsilon^2 \Gamma_{xx} / Pe_I + B i (c - \Gamma) - \epsilon (u \Gamma)_x = 0, \tag{2.16a-d}
\end{align}
$$

with $Ma = \hat{A}\hat{K}_a\hat{K}\hat{H}/(\hat{\mu}\hat{K}_d\hat{Q}^2)$ the Marangoni number, $Da = \hat{K}_a\hat{H}/\hat{D}$ the Damköhler number, $Pe_I = \hat{Q}/\hat{D}_I$ the interfacial Péclet number and $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

$$
u\Gamma - \epsilon\Gamma_x / Pe_I = 0. \qquad (2.17)
$$

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

$$
u = 0, \quad v = 0, \quad c_y = 0. \tag{2.18a-d}
$$

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

$$
\mathbf{q}(\phi^{-}, y) = \mathbf{q}(\phi^{+}, y), \quad \mathbf{q}(-\phi, y) = \mathbf{q}(2-\phi, y). \qquad (2.19a, b)
$$

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

$$
\frac{d}{dx} \int_{0}^{2} \left( u c - \frac{\epsilon c_x}{Pe} \right) dy - \frac{D a}{\epsilon Pe} \left( \Gamma - c(x, 0) \right) = 0 \quad \text{in } \mathcal{D}_1, \qquad (2.20a)
$$

$$
\int_{0}^{2} \left( u c - \frac{\epsilon c_{x}}{P e} \right) d y + \frac{D a}{B i P e} \left( u(x, 0) \Gamma - \frac{\epsilon \Gamma_{x}}{P e_{I}} \right) = 1 \quad \text{in } \mathcal{D}_{1}, \qquad (2.20b)
$$

$$
\int_{-\phi}^{\phi} (\Gamma - c(x, 0)) dx = 0 \quad \text{in } \mathcal{D}_1, \qquad (2.20c)
$$

$$
\int_0^2 \left( u c - \frac{\epsilon c_x}{P e} \right) d y = 1 \quad \text{in } \mathcal{D}_2, \qquad (2.20d)
$$

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

$$
\int_0^2 u \, dy = 1. \qquad (2.21)
$$
Finally, the drag reduction (2.12) becomes

$$DR = \frac{\Delta p_I - \Delta p}{\Delta p_I - \Delta p_U}, \qquad (2.22)$$

where $\Delta p \equiv \langle p \rangle(-\phi) - \langle p \rangle(2 - \phi)$ and $\langle \cdot \rangle \equiv \int_0^2 \cdot dy/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, Pe, Ma, Pe_I, Bi, 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 below.

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

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

$$(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) + \dots \qquad (2.23)$$

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

In $D_1$ and $D_2$, flow is driven by a streamwise pressure gradient, and wall-normal diffusion is comparable to advection. The bulk equations (2.15) become

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

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

$$u_{0y} - \mathcal{M}\Gamma_{0x} = 0, \quad v_0 = 0, \quad c_{0y} - Da(c_0 - \Gamma_0) = 0, \quad \mathcal{B}(c_0 - \Gamma_0) - (u_0\Gamma_0)_x = 0. \quad (2.25a-d)$$

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

$$u_0 = 0, \quad v_0 = 0, \quad c_{0y} = 0. \qquad (2.26a-c)$$

Short 2D Stokes-flow regions arise at the junctions between $D_1$ and $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$,

$$u_0\Gamma_0 = 0. \qquad (2.27)$$

Continuity of bulk concentration (2.19) between $D_1$ and $D_2$ requires

$$c_0(\phi^{-}, y) = c_0(\phi^{+}, y), \quad c_0(-\phi, y) = c_0(2 - \phi, y). \qquad (2.28a, b)$$

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.20) becomes

$$ \frac{d}{dx} \int_{0}^{2} u_{0}c_{0} \, dy - \frac{Da}{\mathcal{D}} (\Gamma_{0} - c_{0}(x, 0)) = 0 \quad \text{in } \mathcal{D}_{1}, \qquad (2.29a) $$

$$ \int_{0}^{2} u_{0}c_{0} \, dy + \frac{Da}{\mathcal{B}\mathcal{D}} (u_{0}(x, 0)\Gamma_{0}) = 1 \quad \text{in } \mathcal{D}_{1}, \qquad (2.29b) $$

$$ \int_{-\phi}^{\phi} (\Gamma_{0} - c_{0}(x, 0)) \, dx = 0 \quad \text{in } \mathcal{D}_{1}, \qquad (2.29c) $$

$$ \int_{0}^{2} u_{0}c_{0} \, dy = 1 \quad \text{in } \mathcal{D}_{2}. \qquad (2.29d) $$

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

$$ \int_{0}^{2} u_{0} \, dy = 1. \qquad (2.30) $$

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

$$ DR_{0} = \frac{\Delta p_{I} - \Delta p_{0}}{\Delta p_{I} - \Delta p_{U}}. \qquad (2.31) $$

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

$$ u_0 = \tilde{U} p_{0x} + \mathcal{M} \bar{U} \Gamma_{0x} \quad \text{in } \mathcal{D}_1, \quad u_0 = \check{U} p_{0x} \quad \text{in } \mathcal{D}_2, \qquad (2.32a, b) $$

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

$$ \tilde{U}_{yy} = 1, \quad \text{subject to} \quad \tilde{U}_{y}(0) = 0, \quad \tilde{U}(2) = 0; \qquad (2.33a-c) $$

$$ \bar{U}_{yy} = 0, \quad \text{subject to} \quad \bar{U}_{y}(0) = 1, \quad \bar{U}(2) = 0; \qquad (2.34a-c) $$

$$ \check{U}_{yy} = 1, \quad \text{subject to} \quad \check{U}(0) = 0, \quad \check{U}(2) = 0. \qquad (2.35a-c) $$

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

$$ \tilde{Q}p_{0x} + \mathcal{M}\bar{Q}\Gamma_{0x} = 1 \quad \text{in } \mathcal{D}_1, \quad \check{Q}p_{0x} = 1 \quad \text{in } \mathcal{D}_2, \qquad (2.36a, b) $$

where $\tilde{Q} \equiv -8/3$, $\bar{Q} \equiv -2$ and $\check{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),

$$ v_0 = \tilde{V} p_{0xx} + \mathcal{M} \bar{V} \Gamma_{0xx} \quad \text{in } \mathcal{D}_1, \quad v_0 = 0 \quad \text{in } \mathcal{D}_2, \qquad (2.37a, b) $$

where $\tilde{V} \equiv -y^3/6 + 2y$ and $\bar{V} \equiv -y^2/2 + 2y$ (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} + \mathcal{M}\bar{V}(2)\Gamma_{0xx} = -\tilde{Q}p_{0xx} - \mathcal{M}\bar{Q}\Gamma_{0xx} = 0$, using the volume-flux condition (2.36a).

We substitute (2.32, 2.37) into the bulk and interfacial surfactant equations (2.24)–(2.29) to derive the 2D 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
<table><caption>TABLE 1. A summary of the transport coefficients in the 2D long-wave model, (2.38)–(2.42), with their definition and physical interpretation. The transport coefficients α, β, γ, δ and ν are written in terms of Pe, Da, Bi, Ma, Pe<sub>I</sub> and ε defined in §2.2. Compared to the 3D transport coefficients in Tomlinson et al. (2023a), α<sub>3D</sub> ∝ ε<sup>2</sup>α, β<sub>3D</sub> ∝ β, γ<sub>3D</sub> ∝ γ, δ<sub>3D</sub> ∝ ε<sup>2</sup>δ and ν<sub>3D</sub> ∝ ε<sup>2</sup>ν.</caption><thead><tr><th>Coefficient</th><th>Name</th><th>Definition</th><th>Physical interpretation</th></tr></thead><tbody><tr><td>α</td><td>Bulk diffusion coefficient</td><td><sup>1</sup>/<sub>εPe</sub> ≡ <sup>ΔP̂</sup>/<sub>Q̂H</sub></td><td><u>Bulk diffusion</u><br><u>Bulk advection</u></td></tr><tr><td>β</td><td>Partition coefficient</td><td><sup>Da</sup>/<sub>BiPe</sub> ≡ <sup>K̂<sub>a</sub></sup>/<sub>ΔĤK̂<sub>d</sub></sub></td><td><u>Adsorption</u><br><u>Desorption</u></td></tr><tr><td>γ</td><td>Surfactant strength</td><td><sup>εMaDa</sup>/<sub>BiPe</sub> ≡ <sup>ΔĤK̂K̂<sub>a</sub><sup>2</sup></sup>/<sub>μK̂<sub>d</sub><sup>2</sup>P̂Q<sup>2</sup></sub></td><td><u>Marangoni effects</u><br><u>Interfacial advection</u></td></tr><tr><td>δ</td><td>Interfacial diffusion coefficient</td><td><sup>Da</sup>/<sub>εBiPePe<sub>I</sub></sub> ≡ <sup>ΔÎK̂<sub>a</sub>P̂</sup>/<sub>ΔH<sup>2</sup>K̂<sub>d</sub>Q̂</sub></td><td><u>Interfacial diffusion</u><br><u>Interfacial advection</u></td></tr><tr><td>ν</td><td>Exchange strength</td><td><sup>Da</sup>/<sub>εPe</sub> ≡ <sup>K̂<sub>a</sub>P̂</sup>/<sub>Q̂</sub></td><td><u>Adsorption</u><br><u>Bulk advection</u></td></tr></tbody></table>

equation (2.24c) gives

$$ \alpha(\epsilon^2 c_{0xx} + c_{0yy}) - \left(\frac{3}{4} - \frac{3y^2}{16}\right)c_{0x} - \frac{\gamma}{\beta}\left(-\frac{1}{2} + y - \frac{3y^2}{8}\right)\Gamma_{0x}c_{0x} - \frac{\gamma}{\beta}\left(\frac{y}{2} - \frac{y^2}{2} + \frac{y^3}{8}\right)\Gamma_{0xx}c_{0y} = 0 \quad \text{in } \mathcal{D}_1, \quad (2.38a) $$

$$ \alpha(\epsilon^2 c_{0xx} + c_{0yy}) - \left(\frac{3y}{2} - \frac{3y^2}{4}\right)c_{0x} = 0 \quad \text{in } \mathcal{D}_2, \quad (2.38b) $$

where α, β and γ 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

$$ 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, \quad (2.39a, b) $$

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

$$ \frac{3}{4}\beta\Gamma_0 - \frac{1}{2}\gamma\Gamma_{0x}\Gamma_0 - \epsilon^2\delta\Gamma_{0x} = 0. \quad (2.40) $$

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

$$ c_{0y} = 0. \quad (2.41) $$
The surfactant-flux conditions (2.29a, b, d) become

$$ \frac{d}{dx} \int_{0}^{2} \left( \left( \frac{3}{4} - \frac{3y^2}{16} \right) c_0 + \frac{\gamma}{\beta} \left( -\frac{1}{2} + y - \frac{3y^2}{8} \right) \Gamma_{0x} c_0 - \epsilon^2 \alpha c_{0x} \right) dy - \nu(\Gamma_0 - c_0(x, 0)) = 0 \quad \text{in } \mathcal{D}_1, \quad (2.42a) $$

$$ \int_{0}^{2} \left( \left( \frac{3}{4} - \frac{3y^2}{16} \right) c_0 + \frac{\gamma}{\beta} \left( -\frac{1}{2} + y - \frac{3y^2}{8} \right) \Gamma_{0x} c_0 - \epsilon^2 \alpha c_{0x} \right) dy + \frac{3}{4} \beta \Gamma_0 - \frac{1}{2} \gamma \Gamma_{0x} \Gamma_0 - \epsilon^2 \delta \Gamma_{0x} = 1 \quad \text{in } \mathcal{D}_1, \quad (2.42b) $$

$$ \int_{0}^{2} \left( \left( \frac{3y}{2} - \frac{3y^2}{4} \right) c_0 - \epsilon^2 \alpha c_{0x} \right) dy = 1 \quad \text{in } \mathcal{D}_2, \quad (2.42c) $$

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

$$ \tilde{Q}p_{0x} + (\gamma/\beta)\tilde{Q}\Gamma_{0x} = 1 \quad \text{in } \mathcal{D}_1, \quad \tilde{Q}p_{0x} = 1 \quad \text{in } \mathcal{D}_2. \quad (2.43a, b) $$

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)/\tilde{Q}$. When the interface is immobilised by surfactant, $\Delta p_I = -2/\tilde{Q}$. Between these two limits, $\Delta p_0 = -2\phi/\tilde{Q} - 2(1 - \phi)/\tilde{Q} + \gamma\tilde{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

$$ DR_0 = 1 - \frac{\gamma \Delta \Gamma_0}{3 \phi \beta}. \quad (2.44) $$

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)

$$ \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)]}. \quad (2.45) $$

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

## 2.4. The 1D long-wave model

Figure 2 shows maps of $(\alpha, \gamma)$-parameter space, illustrating how the present 2D model relates to existing 1D 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

$$c_0 = \bar{c}_0(x) + \bar{c}_1(x, y)/\alpha + \dots, \quad \Gamma_0 = \bar{\Gamma}_0(x) + \bar{\Gamma}_1(x)/\alpha + \dots, \quad (2.46)$$

(2.38)–(2.42) reduce to the “1D long-wave model”, given by

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

$$\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}, \quad (2.47b)$$

$$\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, \quad (2.47c)$$

$$\bar{c}_0 - 2\epsilon^2 \alpha \bar{c}_{0x} = 1 \quad \text{in } \mathcal{D}_2. \quad (2.47d)$$

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

$$\bar{c}_0(\phi^-) = \bar{c}_0(\phi^+), \quad (2.48a)$$

$$\bar{c}_0(-\phi) = \bar{c}_0(2-\phi), \quad (2.48b)$$

$$\bar{c}_0(\phi^-) - 2\epsilon^2 \alpha \bar{c}_{0x}(\phi^-) = \bar{c}_0(\phi^+) - 2\epsilon^2 \alpha \bar{c}_{0x}(\phi^+), \quad (2.48c)$$

$$\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), \quad (2.48d)$$

$$\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. \quad (2.48e)$$

The system of ODEs 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 1D 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- ($\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 2D 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 2D 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 2D 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.
<figure><img src="image_3.png" /><figcaption>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 1D long-wave model of Tomlinson et al. (2023a) are combined with predictions derived in Appendix C from the 2D long-wave model for (a) strong exchange and (b) weak exchange. The Marangoni-dominated region (M) incorporates sub-regions $M^{1D}$, $M_E^{1D}$, $M^{2D}$ and $M_E^{2D}$, for which $DR_0 \ll 1$; the advection-dominated region (A) incorporates sub-regions $A^{1D}$, $A_E^{1D}$, $A^{2D}$ and $A_E^{2D}$, for which $1 - DR_0 \ll 1$; and the diffusion-dominated region (D) incorporates sub-regions $D^{1D}$ and $D_E^{1D}$, for which $1 - DR_0 \ll 1$. Here it is assumed that $\alpha \sim \delta$ and $\beta \sim 1$.</figcaption></figure>

### 3.1. Strong exchange

Figure 3 summarises the surfactant concentration profiles and drag reduction in the strong-exchange limit, computed using the 2D long-wave model, (2.38)–(2.42). For large values of $\gamma$, with strong Marangoni effect, the interface is almost immobilised (figure 3a–b), the $\Gamma_0$ profile is almost linear, and the bulk surfactant concentration transitions from a 2D (figure 3a) to a 1D (figure 3b) 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 3f).

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 2D. 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 sub-regions $M^{1D}$ and $A^{1D}$ for $1 \ll \alpha \ll 1/\epsilon^2$ in figure 2(a), which formally provide a bridge between the 1D and 2D long-wavelength 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}$.

Region $M^{2D}$ (figures 2a and 3a) 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
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

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–(A) and diffusion-dominated (D) regions, the interface is shear-free ($DR_0 \approx 1$). $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 to asymptotic predictions (B 1b) in $M^{1D}$, (B 2b) in $D^{1D}$, (C 21b) 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$. $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 (C 54).

concentration is close to unity, and $\Gamma_0 < 1$ ($\Gamma_0 > 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
<figure><img src="image_11.png" /><figcaption>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. (a) The bulk-concentration boundary layer when Marangoni effects are strong (region $M^{2D}$) and the surfactant distribution almost immobilises the interface. (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.</figcaption></figure>

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

$$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}, \qquad (3.1a, b)$$

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 (C21) below, where the extension must be evaluated numerically. Therefore, we give the simpler formula above, but compare (C21) with numerical results from the 2D 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 3e). 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 2a and 3f), 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 below; 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

$$ c_0(x, 0) \approx \Gamma_0 \approx \begin{cases} \frac{4(\alpha(x+\phi))^{1/2}}{\sqrt{3\pi}\beta} & \text{for } x \in [-\phi, x_0), \\ \frac{3\beta}{2\gamma}(x-x_0) & \text{for } x \in (x_0, \phi], \end{cases} \quad DR_0 = 1 - \frac{a_1\gamma^{3/5}}{\alpha^{1/10}}, \\ u_0(x, 0) \approx \begin{cases} 3/4 & \text{for } x \in [-\phi, x_0), \\ \frac{3}{4}\left[G\left(\frac{x-x_0}{\gamma\alpha^{1/2}}\right)\right]^{-1} & \text{for } x \approx x_0, \\ \frac{3\alpha^{2/3}C'(0)}{5\beta}\left((x-x_0)^{2/3} - \frac{(\phi-x_0)^{5/3}}{x-x_0}\right) & \text{for } x \in (x_0, \phi]. \end{cases} \quad (3.2a-c) $$

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

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

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 2D 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 (both in terms of the required resolution and run time) of the 2D long-wave model and COMSOL simulations (below).

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 2D long-wave model, (2.38)-(2.42), to 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
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

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 (C.54).

slight deviation from no-slip at the liquid–gas interface (figure 5c), generating a weak vertical flow towards the interface. The concentration profiles in the two simulations match well (figure 5a), with the leading-order drag reduction predicted from the 2D 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.2c) 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 (C.40) 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.2c) across the quasi-stagnant cap region. The predictions for the drag reduction and concentration field from the 2D 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} > 0$ for $x \in [-\phi, x_0]$ in figure 5d) and stronger compression (allowing desorption) in the quasi-stagnant region ($u_{0x} < 0$ for $x \in [x_0, \phi]$ in figure 5d). There are discrepancies in $u_0$ at both the upstream and downstream contact lines. At the upstream contact line, the 2D 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 2D 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 2D long-wave model as streamwise diffusion vanishes, the two solutions align in the limit $\epsilon \to 0$, where the 2D 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

When bulk-surface exchange is weak ($\nu = 0.01$), 2D 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 6a-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 6f), which is smoothed by streamwise diffusion as $\alpha = \delta$ increase (figure 6g).

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 (B 4, B 5)) for large $\alpha$ and introducing new limits $M_E^{2D}$ and $D_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}$ (figure 2b and 6f), 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 1D case, is supported by asymptotic calculations in §C.3 below, which yield

$$c_0 \approx 1, \quad \Gamma_0 \approx \begin{cases} 0 & \text{for } x \in [-\phi, x_0], \\ \frac{3\beta}{2\gamma}(x-x_0) & \text{for } x \in [x_0, \phi], \end{cases} \quad DR_0 \approx 1 - \left(\frac{2\gamma}{3\phi\beta}\right)^{1/2}, \quad (3.4a-c)$$

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 2D 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 non-linear 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_E^{2D}$ from region $M_E^{2D}$, using (B 4b), i.e. for decreasing $\gamma$.

A comparison between the solutions evaluated using the 2D long-wave model, (2.38)-
<figure><img src="image_16.png" /><figcaption>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- ($A_E$) and diffusion-dominated ($D_E$) regions, the interface is mostly shear-free ($DR_0 \approx 1$). $c_0$, $\Gamma_0$, $\langle c_0 \rangle$ and $c_0(x, 0)$ are plotted in (a) $M_E^{2D}$ and (b) $M_E^{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 to (B 4b) in $M_E^{1D}$, (B 5b) in $D_E^{1D}$, (B 4b) in $M_E^{2D}$ and (3.4b) in $A_E^{2D}$. The dashed line in (e) represents the largest $\gamma$ for which $\Gamma_0(-\phi) = 0$ when $\alpha \ll 1$. $c_0$, $\Gamma_0$, $\langle c_0 \rangle$ and $c_0(x, 0)$ are plotted in (f) $A_E^{2D}$ (notice the classical stagnant cap profile) and (g) $D_E^{1D}$ (where streamwise diffusion has smoothed the stagnant cap), where $x = x_0$ is plotted using (C 61).</figcaption></figure>

(2.42), and COMSOL simulations is given in figure 7, focussing 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 2D long-wave model predicts $DR_0 = 0.027$,
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>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_E^{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_E^{2D}$) and $x = x_0$ is plotted using (C 61).</figcaption></figure>

while the numerical simulations yield $DR_{NS} = 0.056$ (due to the 4% difference in $\Delta\Gamma_0$ shown in figure 7a). 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 2D 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\tilde{P}/\tilde{H}$ is half the length of the liquid–gas interface divided by the channel height, $L_m^2 = (\tilde{n}\tilde{R}\tilde{T}\tilde{K}_a^2\tilde{c}_0)/(\tilde{D}\tilde{\mu}\tilde{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
<table><caption>TABLE 2. Summary of the leading-order drag reduction (DR₀) scaling in regions M²D, M¹D and M_E¹D and the corresponding leading-order slip length (λ_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₀ ≪ 1 and λ_e ~ (φ/ε)λ₀, where L̂₂ = (n̂R̂T̂Ĥ²K̂ₐ²c₀)/(D̂μK̂ₐ²), φ/ε = φP̂/Ĥ, L_m = L̂_m/Ĥ and b = Q̂/D̂.</caption><thead><tr><th></th><th>M<sup>2D</sup> (α &ll; 1)</th><th>M<sup>1D</sup> (1 &ll; α &ll; 1/ε<sup>2</sup>)</th><th>M<sup>1D</sup> and M<sub>E</sub><sup>1D</sup> (α &gg; 1/ε<sup>2</sup>)</th></tr></thead><tbody><tr><td>Drag reduction scaling</td><td><sup>φ<sup>5/3</sup>α<sup>2/3</sup></sup>&frasl;<sub>γ</sub> = <sup>φ<sup>5/3</sup>P̂<sup>5/3</sup>Ĥ<sup>1/3</sup>Q̂<sup>1/3</sup></sup>&frasl;<sub>L<sub>m</sub><sup>2</sup>D<sup>1/3</sup></sub></td><td><sup>φ</sup>&frasl;<sub>L<sub>m</sub><sup>2</sup>D̂</sub> = <sup>φP̂ĤQ̂</sup>&frasl;<sub>L<sub>m</sub><sup>2</sup>D̂</sub></td><td><sup>ε<sup>2</sup>α</sup>&frasl;<sub>γ</sub> = <sup>Ĥ<sup>2</sup></sup>&frasl;<sub>L<sub>m</sub><sup>2</sup></sub></td></tr><tr><td>Slip length scaling</td><td><sup>b<sup>1/3</sup>(φ/ε)<sup>8/3</sup></sup>&frasl;<sub>L<sub>m</sub><sup>2</sup></sub></td><td><sup>b(φ/ε)<sup>2</sup></sup>&frasl;<sub>L<sub>m</sub><sup>2</sup></sub></td><td><sup>φ/ε</sup>&frasl;<sub>L<sub>m</sub><sup>2</sup></sub></td></tr></tbody></table>

simplifies to λ_e = (Δp_I - Δp_U)/(εΔp_IΔp_U) (A¹D and A²D). Hence, outside of the specific microchannel applications considered in Temprano-Coleto et al. (2023), whose predictions agree with those in region M¹D (as discussed below), we have found new parameters in regions M²D that influence surfactant impairment of SHS drag reduction. As diffusion weakens, the deviation between the drag-reduction predictions in M¹D and M²D will increase.

For laboratory experimental studies documented in the literature, we can use our 2D 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 P̂, Q̂, Ĥ and K̂ vary across experiments. Peaudecerf et al. (2017) use α ≈ 1 and 4.7 × 10¹ ≲ γ ≲ 2.9 × 10⁴ (at the M¹D/M²D boundary, figure 2); and Temprano-Coleto et al. (2023) use α ≈ 25 and 10⁻¹ ≲ γ ≲ 3.7 × 10² (at the A¹D/M¹D 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 K̂ₐ and K̂₉ given in table 3. However, as K̂ₐ and K̂₉ are merely estimates based on fitting to microchannel experiments, weak exchange could be achieved, for instance, when K̂ₐ = 3.4 × 10⁻⁸m²s⁻¹ and K̂₉ = 7.5 × 10⁻⁵s⁻¹. 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 α ≈ 25 and ε ≈ 0.002, we confirm 1 ≪ α ≪ 1/ε² and they are within the regime M¹D. From (B 1b), we calculate DR₀, convert to λ_e using table 2, and then to u_Ic/u_Ic clean following Temprano-Coleto et al. (2023). Agreement is observed across φ/ε = 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 φ/ε = 762 could be attributed to experimental error, which is consistent with the expected behaviour of increasing slip with grating length.

## 4. Discussion

In this study, we develop a long-wave theory to analyse the behaviour of a 2D 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-
<table><caption>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).</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>Adsorption rate</td><td>K̂<sub>a</sub></td><td>m<sup>2</sup> s<sup>-1</sup></td><td>3.4 × 10<sup>-4</sup></td></tr><tr><td>Desorption rate</td><td>K̂<sub>d</sub></td><td>s<sup>-1</sup></td><td>0.75</td></tr><tr><td>Surface activity</td><td>Å</td><td>J mol<sup>-1</sup></td><td>4913.6</td></tr><tr><td>Dynamic viscosity</td><td>μ̂</td><td>kg m<sup>-1</sup> s<sup>-1</sup></td><td>8.9 × 10<sup>-4</sup></td></tr><tr><td>Surface diffusivity</td><td>D̂</td><td>m<sup>2</sup> s<sup>-1</sup></td><td>7 × 10<sup>-10</sup></td></tr><tr><td>Bulk diffusivity</td><td>D̂<sub>I</sub></td><td>m<sup>2</sup> s<sup>-1</sup></td><td>7 × 10<sup>-10</sup></td></tr></tbody></table>

<table><caption>TABLE 4. A comparison between the slip velocity normalised by the clean (surfactant-free) value, $u_{Ic}/u_{Ic}^{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 $φ/ε$ is half the length of the liquid–gas interface divided by the channel height.</caption><thead><tr><th></th><th>φ/ε = 253</th><th>φ/ε = 418</th><th>φ/ε = 587</th><th>φ/ε = 762</th></tr></thead><tbody><tr><td>Peaudecerf et al. (2017)</td><td>-</td><td>-</td><td>0.1410</td><td>-</td></tr><tr><td>Temprano-Coleto et al. (2023)</td><td>0.0661</td><td>0.1187</td><td>0.1903</td><td>0.1450</td></tr><tr><td>Current</td><td>0.0610</td><td>0.1164</td><td>0.1859</td><td>0.2716</td></tr></tbody></table>

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 2D long-wave model, we delineate asymptotic regions within the parameter space. Under conditions of strong cross-channel diffusion, our findings align with the 1D results in Tomlinson et al. (2023a). Conversely, when cross-channel diffusion is weak, we unveil 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 2D 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 ($I_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_E^{1D}$), advection-dominated ($A^{1D}$ and $A_E^{1D}$) and diffusion-dominated ($D^{1D}$ and $D_E^{1D}$) regions identified in Tomlinson et al. (2023a) (see figure 2). However, under weak cross-channel diffusion, we unveil new subregions dominated by Marangoni effects ($M^{2D}$)
and $M_E^{2D}$) and advection ($A^{2D}$ and $A_E^{2D}$). 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) above. 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 $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 4b). 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 $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 2D 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.
## Acknowledgements

We acknowledge support from CBET-EP SRC (EPSRC Ref. EP/T030739/1, NSF #2054894), as well as partial support from ARO MURI W911NF-17-1-0306. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Numerical solutions for weak cross-channel diffusion

The numerical solution to (2.38)–(2.42) can be divided into two subroutines, employing Chebyshev collocation techniques as outlined by Trefethen (2000). First, given $\Gamma_0^{old}$, we evaluate $c_0^{new}$ in $D_1$ and $D_2$. The subdomains of the periodic cell are mapped to the 2D canonical Chebyshev collocation domain, $D_n = \{\xi \in [-1, 1]\} \times \{\eta \in [-1, 1]\}$, using the transformations

$$ (\xi_1, \xi_2, \eta) = ((x+\phi)/\phi - 1, (x-\phi)/(1-\phi) - 1, y-1). \quad (A1) $$

We then discretise using $N = (N_\xi + 1) \times (N_\eta + 1)$ nodes in each subdomain, at points

$$ (\xi_1, i, \xi_2, i, \eta_j) = (\cos(i\pi/N_\xi), \cos(i\pi/N_\eta), \cos(j\pi/N_\eta)), \quad (A2) $$

where $i = 0, 1, ..., N_\xi$ and $j = 0, 1, ..., N_\eta$. Continuity boundary conditions (2.28) are enforced on $c_0^{new}$ at $i = 0$ and $N_\xi$, supplemented with continuity of bulk diffusive surfactant flux between subdomains

$$ c_{0,1}^{\text{new}} = c_{0,2}^{\text{new}}, \quad \frac{\partial c_{0,1}^{\text{new}}}{\partial \xi_1} = \frac{\partial c_{0,2}^{\text{new}}}{\partial \xi_2}. \quad (A3a, b) $$

Boundary conditions (2.39)–(2.41) are enforced at $j = 0$ and $N_\eta$, ensuring continuity of surfactant flux at the liquid-gas interface, no-flux of surfactant at the solid wall in $D_1$ and no-flux of surfactant at the solid walls in $D_2$. Discretising the bulk surfactant equation (2.38) in conservative form in domains $D_1$ and $D_2$ and incorporating the discretised matching conditions (represented by $M_{12}$ and $M_{21}$) gives

$$ \begin{pmatrix} \mathbf{A}_1 & \mathbf{M}_{12} \\ \mathbf{M}_{21} & \mathbf{A}_2 \end{pmatrix} \begin{pmatrix} c_{0,1}^{\text{new}} \\ c_{0,2}^{\text{new}} \end{pmatrix} = \begin{pmatrix} \mathbf{f}_1 \\ \mathbf{0} \end{pmatrix}, \quad (A4) $$

where $\mathbf{f}_1$ is the forcing due to the interfacial surfactant concentration $\Gamma_0^{old}$ in (2.39a).

Second, given a $c_0^{new}$, we evaluate $\Gamma_0^{new}$ for $x \in [-\phi, \phi]$ and $y = 0$. The liquid-gas interface is mapped to the 1D canonical Chebyshev collocation domain, $D_i = \{\xi \in [-1, 1]\}$, using $\xi_1$ in (A1), discretised using $N_\xi + 1$ nodes at points $\xi_1, i$ in (A2). Discretising the interfacial surfactant equation (2.39b) for $x \in [-\phi, \phi]$ and $y = 0$ and incorporating the discretised no-flux conditions (2.40), results in

$$ \mathbf{B}(\Gamma_0^{\text{old}})\Gamma_0^{\text{new}} = \mathbf{g}, \quad (A5) $$

where $\mathbf{g}$ is the forcing due to the bulk surfactant concentration $c_0^{new}$ in (2.39b).

Therefore, to evaluate $DR_0$ for a given $\alpha, \beta, \gamma, \delta, \epsilon, \nu$ and $\phi$, we choose an initial guess based on the 1D long-wave model discussed in §2.4. The initial guess can be substituted into the forcing $\mathbf{f}_1$ in the bulk concentration subroutine and the linear system inverted to calculate $c_0^{new}$. The bulk concentration can then be substituted into the forcing $\mathbf{g}$ in
the interfacial concentration subroutine and the linear system inverted to calculate $\Gamma_0^{\text{new}}$. We then set $\Gamma_0^{\text{old}} = \Gamma_0^{\text{new}}$ and repeat the above steps. This procedure of updating $c_0$ and $\Gamma_0$ is then continued until convergence is achieved for both the bulk and interfacial concentration fields.

## Appendix B. Key results of the 1D long-wave model

Here we outline key limits of the 1D long-wave model (2.47, 2.48), exploiting analogous findings as in Tomlinson et al. (2023a). For strong exchange, with $\alpha \gg 1$ and $\gamma \gg \max(1, \epsilon^2\alpha)$, a region dominated by Marangoni effects ($M^{1D}$, figure 2a) is characterised by the immobilisation of the liquid-gas interface. The surfactant distribution and drag reduction are given by

$$ \bar{c}_0 \approx \bar{\Gamma}_0 \approx 1 + \frac{3\beta}{2\gamma} \left( x - \frac{\phi(E+1)}{(E-1)} \right), \quad DR_0 \approx \frac{2}{\gamma} \left( \epsilon^2(2\alpha + \delta) + \frac{\phi(E+1)}{(E-1)} \right), \quad (\text{B 1a, b}) $$

where $E \equiv \exp((1 - \phi)/(\epsilon^2\alpha))$. Thus, $DR_0 \approx 2\phi/\gamma$ for $1 \ll \alpha \ll 1/\epsilon^2$ and $DR_0 \approx 2\epsilon^2(2\alpha + \delta)/\gamma$ for $\alpha \gg 1/\epsilon^2$. Region $M^{1D}$ transitions into a diffusion-dominated region ($D^{1D}$, figure 2a) when $\alpha \gg 1/\epsilon^2$ and $\gamma \sim \epsilon^2\alpha$, and into an advection-dominated region ($A^{1D}$, figure 2a) when $1 \ll \alpha \ll 1/\epsilon^2$ and $\gamma \sim 1$. For $\gamma \ll \epsilon^2\alpha$ and $\alpha \gg 1/\epsilon^2$ (region $D^{1D}$),

$$ \bar{c}_0 \approx \bar{\Gamma}_0 \approx \frac{4\alpha + 2\delta(1 - \phi)}{\alpha(4 + 3\beta\phi) + 2\delta(1 - \phi)}, \quad DR_0 \approx 1 - \frac{\gamma(1 - \phi)}{\epsilon^2(\alpha(4 + 3\beta\phi) + 2\delta(1 - \phi))}, \quad (\text{B 2a, b}) $$

and for $\gamma \ll 1$ and $1 \ll \alpha \ll 1/\epsilon^2$ (region $A^{1D}$),

$$ \bar{c}_0 \approx \bar{\Gamma}_0 \approx \frac{4}{3\beta + 4} + \frac{3\beta}{3\beta + 4} \exp \left( \frac{(4 + 3\beta)(x - \phi)}{4\epsilon^2(2\alpha + \delta)} \right), \quad DR_0 \approx 1 - \frac{\gamma}{\phi(3\beta + 4)}. \quad (\text{B 3a, b}) $$

The onset of 2D effects in the 1D long-wave model was estimated in Tomlinson et al. (2023a) to arise via the emergence of shear dispersion, where $\alpha \sim 1$ and $\gamma \ll 1$ or $\alpha \sim 1/\gamma$ and $\gamma \gg 1$. We refine such estimates by solving the 2D long-wave model, which reveals regions $M^{2D}$ and $A^{2D}$ shown in figure 2(a).

A similar picture emerges when bulk-surface exchange is weak. A region dominated by Marangoni effects ($M^{1D}$, figure 2b) arises for $\delta \gg 1$ and $\gamma \gg \max(1, \epsilon^2\alpha)$. Here, $\bar{c}_0$ decouples from $\bar{\Gamma}_0$ and is approximately equal to its background value (i.e. $\bar{c}_0 \approx 1$). $\bar{\Gamma}_0$ varies around this value to fulfil the net adsorption-desorption condition, (2.29), decreasing the drag reduction compared to the strong-exchange problem

$$ \bar{c}_0 \approx 1, \quad \bar{\Gamma}_0 \approx 1 + \frac{3\beta x}{2\gamma} + O(\epsilon^2), \quad DR_0 \approx \frac{2\epsilon^2\delta}{\gamma}. \quad (\text{B 4a-c}) $$

The region $M^{1D}$ transitions to a diffusion-dominated region ($D^{1D}$, figure 2b) when $\gamma \sim \epsilon^2\delta$ and $\delta \gg 1/\epsilon^2$ and an advection-dominated region ($A^{1D}$, figure 2b) when $\gamma \sim 1$ and $1 \ll \delta \ll 1/\epsilon^2$. For $\gamma \ll \epsilon^2\delta$ and $\delta \gg 1/\epsilon^2$ (region $D^{1D}$),

$$ \bar{c}_0 \approx 1, \quad \bar{\Gamma}_0 \approx 1 + \frac{3\beta x}{4\epsilon^2\delta}, \quad DR_0 \approx 1 - \frac{\gamma}{2\epsilon^2\delta}, \quad (\text{B 5a-c}) $$

and for $\gamma \ll 1$ and $1 \ll \delta \ll 1/\epsilon^2$ (region $A^{1D}$),

$$ \bar{c}_0 \approx 1, \quad \bar{\Gamma}_0 \approx \begin{cases} 0 & \text{for } x \in [-\phi, x_0], \\ \frac{3\beta}{2\gamma}(x - x_0) & \text{for } x \in [x_0, \phi], \end{cases} \quad DR_0 \approx 1 - \left(\frac{2\gamma}{3\phi\beta}\right)^{1/2}, \quad (\text{B 6a-c}) $$
where $x_0 = \phi - (8\phi\gamma/(3\beta))^{1/2}$ gives the length of the stagnant cap.

## Appendix C. Asymptotic solutions for weak cross-channel diffusion

We derive asymptotic solutions to the 2D long-wave model, (2.38)–(2.42), when bulk diffusive effects are weak and confined to a boundary layer along $y = 0$. In this limit, streamwise diffusion terms do not appear at leading order, which means that we do not resolve inner regions around $x = -\phi$, $\phi$ and $2 - \phi$. First, we consider the strong-exchange regime, where Marangoni effects (Appendix C.1, figure 4a) and advection (Appendix C.2, figure 4b) dominate. Second, we consider the weak-exchange regime (Appendix C.3). Throughout, we assume $\beta = O(1)$ and $\phi = O(1)$.

Assuming that $c_0 \approx 1$ in the core of the channel, we anticipate that the transport equations (2.38) for the bulk surfactant involve a diffusive boundary layer near $y = 0$. The diffusive boundary layer is made up of components over the interface and ridge, which are governed by

$$ \alpha c_{0yy} - \left( \frac{3}{4} + \frac{\gamma}{\beta} \left( -\frac{1}{2} + y \right) \Gamma_{0x} \right) c_{0x} - \frac{\gamma}{\beta} \left( \frac{y}{2} - \frac{y^2}{2} \right) \Gamma_{0xx} c_{0y} = 0 \quad \text{in } \mathcal{D}_1, \quad (C1a) \\ \alpha c_{0yy} - 3yc_{0x}/2 = 0 \quad \text{in } \mathcal{D}_2, \quad (C1b) $$

where we have assumed $y \ll 1$ and have neglected second (third) order terms in $O(y^2)$ ($O(y^3)$) in the streamwise (wall-normal) velocity. For $x \in [-\phi, \phi]$ and $y = 0$, the bulk-interface flux condition and the interfacial transport equation for the surfactant in (2.39) give

$$ \alpha c_{0y} = \nu(c_0 - \Gamma_0) = \beta \left[ \left( \frac{3}{4} - \frac{\gamma}{2\beta} \Gamma_{0x} \right) \Gamma_0 \right]_x. \quad (C2) $$

At $x = \pm \phi$ and $y = 0$, the no-surfactant-flux conditions at the stagnation points (2.40) become

$$ \left( \frac{3}{4} - \frac{\gamma}{2\beta} \Gamma_{0x} \right) \Gamma_0 = 0. \quad (C3) $$

For $x \in [\phi, 2 - \phi]$ and $y = 0$, the no-flux condition at the solid wall (2.41) gives

$$ c_{0y} = 0. \quad (C4) $$

For $x \in [-\phi, 2 - \phi]$ and $y \to 1$ (shorthand for $y$ values in the core of the channel, far outside the boundary layer along $y = 0$), the core condition is

$$ c_0 \to 1. \quad (C5) $$

The boundary-layer equations (C1)–(C5) can be integrated across the boundary layer to derive excess-surfactant-flux conditions. The total flux constraint (2.42b, c) is satisfied at leading order in the core of the channel, implying that

$$ \beta \left( \frac{3}{4} - \frac{\gamma}{2\beta} \Gamma_{0x} \right) \Gamma_0 + \int_0^1 \left( \frac{3}{4} + \frac{\gamma}{\beta} \left( -\frac{1}{2} + y \right) \Gamma_{0x} \right) (c_0 - 1) \, dy = 0 \quad \text{in } \mathcal{D}_1, \quad (C6a) $$

$$ \int_{-\phi}^{\phi} c_{0y} \, dx = \int_{-\phi}^{\phi} (c_0 - \Gamma_0) \, dx = 0 \quad \text{in } \mathcal{D}_1, \quad (C6b) $$

$$ \int_0^1 \frac{3}{2} y (c_0 - 1) \, dy = 0 \quad \text{in } \mathcal{D}_2, \quad (C6c) $$

provided $c_0 \to 1$ sufficiently quickly where the boundary layer meets the core flow. Equation (C6b) is derived by integrating (C2) and using (C3). The condition in (C6a)
indicates that the extra surfactant flux carried by the flow at the interface balances the flux deficit in the boundary layer (where $c_0 \le 1$), which happens over the upstream part of the interface. Conversely, over the downstream part of the interface, the reduced interfacial surfactant flux is compensated by a surplus of bulk surfactant in the boundary layer (where $c_0 \ge 1$). This can be seen in the boundary-layer profiles at small $\alpha$ values, as depicted in figures 3(a, f) and 6(a, f). The condition in (C6b) specifies that the total amounts of surfactant adsorbed and desorbed along the interface must match, with continuity between the bulk diffusive flux (first integral) and the adsorption-desorption flux (second integral). Motivated by the simulations and scaling arguments provided in §3.1, we propose the asymptotic structures for regions $M^{2D}$ and $A^{2D}$ outlined in figure 4 and derive asymptotic solutions of the surfactant distributions in the following subsections for the different regions in the parameter space where boundary layers exist.

### C.1. Strong Marangoni effect and strong bulk-surface exchange: region $M^{2D}$

We first consider the boundary-layer equations, (C1)–(C6), in the limit where bulk diffusion is weak and bulk-surface exchange is strong, i.e. $\alpha \ll 1$ and $\nu$ large. Suppose the interfacial surfactant gradient leads to a significant adverse Marangoni effect, immobilising the interface. Surfactant adsorbs onto the interface along its upstream part and desorbs along the downstream part. Motivated by the scaling arguments given in §3.1, we set

$$ \Gamma_0 = M_1 + \frac{3\beta x}{2\gamma} - \alpha^{2/3}G(x) + \dots, \quad c_0 = C(x, Y), \quad (C 7a, b) $$

where $y = \alpha^{1/3}Y$ balances advection and diffusion in the bulk-surfactant equation (C1) and $M_1$ is a constant to be determined. Substituting (C7) into (C1)–(C6), (C1) becomes, at leading order,

$$ C_{YY} - 3YC_x/2 = 0. \quad (C 8) $$

For $x \in [-\phi, \phi]$ and $Y = 0$, we enforce $\Gamma_0 = c_0$ as bulk-surface exchange is strong, so that (C2) becomes

$$ C_Y - \frac{\gamma}{2} \left[ G_x \left( \frac{3\beta x}{2\gamma} + M_1 \right) \right]_x = 0. \quad (C 9) $$

At $x = \pm\phi$ and $Y = 0$, the no-flux conditions (C3) at the stagnation points are

$$ G_x = 0, \quad (C 10) $$

and for $x \in [\phi, 2-\phi]$ and $Y = 0$, the no-flux condition (C4) at the solid wall is

$$ C_Y = 0. \quad (C 11) $$

In the bulk of the periodic domain, for $x \in [-\phi, 2-\phi]$ and $Y \to \infty$, the core condition (C5) becomes

$$ C \to 1. \quad (C 12) $$

The excess-surfactant-flux conditions (C6) reduce to

$$ \int_0^\infty \frac{3Y}{2} (C-1) dY + \frac{G_x}{2} \left( \frac{3\beta x}{2} + \gamma M_1 \right) = 0 \quad \text{in } \mathcal{D}_1, \quad (C 13a) $$

$$ \int_0^\infty \frac{3Y}{2} (C-1) dY = 0 \quad \text{in } \mathcal{D}_2. \quad (C 13b) $$

We construct a similarity solution to (C8)–(C13), assuming that the length of the interface is long compared to the boundary-layer thickness, and that there is no interaction
between adjacent plastrons. The bulk concentration is split into two contributions that grow from the leading edge of the plastron,

$$C = 1 - \left(1 - M_1 + \frac{3\beta\phi}{2\gamma}\right) F(\eta) + \frac{3\beta}{2\gamma}(x + \phi)H(\eta) \quad (C 14)$$

for some functions F and H depending on the similarity variable $\eta = Y/(x + \phi)^{1/3}$. The function F captures the response of the bulk field to the sudden drop in concentration along $Y = 0$ at $x = -\phi$; H captures the response to the linear profile of $\Gamma_0$ in $x > -\phi$. Substituting (C 14) into (C 8)-(C 12) leads to the boundary-value problems

$$F_{\eta\eta} + \eta^2 F_{\eta}/2 = 0, \quad F(0) = 1, \quad F(\infty) = 0, \quad (C 15a-c)$$

$$H_{\eta\eta} + \eta^2 H_{\eta}/2 - 3\eta H/2 = 0, \quad H(0) = 1, \quad H(\infty) = 0, \quad (C 16a-c)$$

which have solutions

$$F = \frac{\Gamma(\frac{1}{3}, \eta^3/6)}{\Gamma(\frac{1}{3})}, \quad H = \frac{(4 + \eta^3)\Gamma(\frac{1}{3}, \eta^3/6) - 6^{2/3}\eta \exp(-\eta^3/6)}{4\Gamma(\frac{1}{3})}. \quad (C 17a, b)$$

Here $\Gamma(z)$ is the gamma function and $\Gamma(s, z)$ is the upper incomplete gamma function. We then use the surfactant-flux conditions (C 13) to calculate the response of the interfacial concentration field to the bulk field (C 14),

$$\frac{\gamma G_x}{2} \left( \frac{3\beta x}{2\gamma} + M_1 \right) = \frac{3}{2} (x + \phi)^{2/3} \left( 1 - M_1 + \frac{3\beta\phi}{2\gamma} \right) \int_0^\infty \eta F \, d\eta - \frac{9\beta}{4\gamma} (x + \phi)^{5/3} \int_0^\infty \eta H \, d\eta. \quad (C 18)$$

The total flux of surfactant must be continuous between domains. The no-flux condition (C 10) at $x = \phi$ requires that

$$M_1 = 1 + \frac{3\beta\phi}{2\gamma} \left( \int_0^\infty \eta F \, d\eta - 2 \int_0^\infty \eta H \, d\eta \right) / \int_0^\infty \eta F \, d\eta = 1 - \frac{3\beta\phi}{10\gamma}. \quad (C 19)$$

Integrating the total-flux condition (C 18) across the plastron reveals that

$$\Delta G = \int_{-\phi}^{\phi} \left\{ \frac{27\beta\phi}{10\gamma} (x + \phi)^{2/3} \int_{0}^{\infty} \eta F \, d\eta - \frac{9\beta}{4\gamma} (x + \phi)^{5/3} \int_{0}^{\infty} \eta H \, d\eta \right\} / \left[ \frac{\gamma}{2} \left( 1 + \frac{3\beta}{2\gamma} \left( x - \frac{\phi}{5} \right) \right) \right] \, dx. \quad (C 20)$$

Hence, in region $M^{2D}$, substituting (C 7) into (2.44), the leading-order interfacial concentration and drag reduction are given by

$$\Gamma_0 = 1 + \frac{3\beta}{2\gamma} \left(x - \frac{\phi}{5}\right) + \dots, \quad DR_0 = \frac{\gamma\alpha^{2/3}\Delta G}{3\phi\beta}. \quad (C 21a, b)$$

The leading-order interfacial surfactant distribution $\Gamma_0$ in region $M^{2D}$ is a linear profile with mean value close to 1, as found numerically (see figure 5a). This approximation requires $\Gamma_0 > 0$ at $x = -\phi$, so that $\gamma > 9\beta\phi/5$. For $\gamma \gg 1$, (C 21) reduces to

$$\Delta G = \frac{81\beta(2\phi)^{8/3}}{50\gamma^2} \int_0^\infty \eta F \, d\eta - \frac{27\beta(2\phi)^{8/3}}{16\gamma^2} \int_0^\infty \eta H \, d\eta \equiv \frac{\Delta J}{\gamma^2}, \quad (C 22)$$

$$\Gamma_0 = 1 + \frac{3\beta}{2\gamma} \left(x - \frac{\phi}{5}\right) + \dots, \quad DR_0 = \frac{\alpha^{2/3} \Delta J(\beta, \phi)}{3\phi\beta\gamma} \equiv \frac{m_1 \alpha^{2/3} \phi^{5/3}}{\gamma}, \quad (C 23a, b)$$
where $m_1 = 81(3/2)^{2/3}/(50\Gamma(\frac{1}{3})) \approx 0.79$, as in (3.1). This expression for the leading-order drag reduction (C23b) assumes strong Marangoni effects ($\gamma \gg 1$), weak diffusion ($\alpha \ll 1$) and strong bulk-surface exchange, and shows how drag reduction in this limit arises from differences between the two self-similar components F and H of the bulk concentration field. As mentioned previously, this asymptotic prediction agrees with the numerical simulations in the $M^{2D}$ region shown in figure 3(a, c, d, e).

### C.2. Weak Marangoni effect and strong bulk-surface exchange: region A$^{2D}$

We now consider the boundary-layer equations (C1)-(C6) when $\gamma \ll 1$ and $\alpha \ll 1$. Weak Marangoni effects imply that the interface cannot support any stress over most of its length. Consequently, there exists a slip flow along most of the plastron, characterised by a slip velocity in (C1a) close to 3/4. A short transition region, mediated by surface diffusion, exists between a no-slip and a slip flow at the upstream contact line, but this is disregarded here. Surfactant adsorbs onto the interface in the slip region (in $-\phi < x < x_0$, say) and accumulates (and desorbs) in a region at the downstream end of the plastron (in $x_0 < x < \phi$). Across this “quasi-stagnant” region, the surface velocity smoothly decreases to zero at the downstream contact line. Motivated by simulations and the scaling arguments given in §3.1, we propose the asymptotic structure illustrated in figure 4(b), in which “deceleration” and “transition” regions are nested at the tip of the quasi-stagnant region over length scales indicated. The formal requirement that the deceleration, transition and quasi-stagnant layers are nested ($\gamma\alpha^{1/2} \ll \gamma^{3/4}\alpha^{1/8} \ll \gamma^{3/5}/\alpha^{1/10} \lesssim 1$) is $\gamma \ll \alpha^{-3/2}$. We now discuss these regions in turn.

#### C.2.1. The slip region ($-\phi < x < x_0$)

In the slip region, surfactant adsorbs onto the interface, causing the interface to carry a spatially growing component of the overall flux. Consequently, the bulk advective flux falls with $x$. Here, we set

$$ \Gamma_0 = \alpha^{1/2}G(x), \quad c_0 = C(x, Y), \qquad (C24a, b) $$

where $y = \alpha^{1/2}Y$ balances advection and diffusion in the bulk-transport equation (C1). When we substitute (C24) into (C1)-(C6), the bulk-transport equation (C1) in $D_1$ becomes

$$ C_{YY} - 3C_x/4 = 0. \qquad (C25) $$

At the interface $\alpha^{1/2}G(x) = C(x, 0)$ in the strong-exchange limit and the flux condition (C2) reduces to

$$ C_Y(x, 0) = 3\beta G_x/4. \qquad (C26) $$

The no-flux condition at the upstream stagnation point (C3) is

$$ G(-\phi) = 0. \qquad (C27) $$

For $x \in [-\phi, x_0]$ and $Y \to \infty$, the core condition (C5) is

$$ C \to 1. \qquad (C28) $$

We construct a similarity solution to (C25)-(C28) by expanding the surfactant fields for $\alpha \ll 1$ as

$$ C(x, Y) = C_0(\eta) + (\alpha(x + \phi))^{1/2}C_1(\eta) + \dots, \quad G(x) = (x + \phi)^{1/2}C_1(0), \qquad (C29a, b) $$

with the similarity variable given by $\eta = Y/(\phi + x)^{1/2}$. Substituting (C29) into (C25)-(C28) yields the boundary-value problems

$$ C_0'' + 3\eta C_0'/8 = 0, \quad C_0(0) = 0, \quad C_0(\infty) = 1, \qquad (C30a-c) $$
$$C_1'' + 3\eta C_1'/8 - 3C_1/8 = 0, \quad 3\beta C_1(0) = 8C_0'(0), \quad C_1(\infty) = 0. \quad (C31a-c)$$

These are solved by

$$C_0 = \operatorname{erf}(\sqrt{3}\eta/4), \quad C_1 = \frac{4\operatorname{ierfc}(\sqrt{3}\eta/4)}{(\sqrt{3}\beta)}, \quad (C32a, b)$$

where erf(z) is the error function, erfc(z) is the complementary error function and ierfc(z) ≡ e⁻z²/√π - z erfc(z) is the integrated complementary error function. We find that the interfacial concentration field in the slip region is given by

$$\Gamma_0 = \frac{4}{\beta} \left( \frac{\alpha(x + \phi)}{3\pi} \right)^{1/2} + \dots, \quad (C33)$$

contributing to (3.2a). The interfacial surfactant distribution Γ₀ in region A²D is a non-linear boundary layer profile scaling as α¹/²x¹/² in the slip region of the interface (up to x₀), as found numerically (see figure 5b). There is a large vertical concentration gradient across the boundary layer. Adsorption of surfactant onto the interface is accommodated by the horizontal gradient of the advective flux βu₀Γ₀ ~ 3βΓ₀/4. Marangoni effects remain subdominant, so the slip flow is uniform to this order, but the concentration rises with x to accommodate the adsorbed material. The total flux adsorbed across the slip region, up to the location x = x₀ (to be determined), is given by

$$\int_{-\phi}^{x_0} c_{0y}(x, 0) \, dx = \left( \frac{9\alpha(x_0 + \phi)}{3\pi} \right)^{1/2} + \dots, \quad (C34)$$

which is subdominant to the flux in the core. The slip region terminates abruptly, meeting a very short deceleration region near x = x₀ nested at the start of the quasi-stagnant region (figure 4b).

### C.2.2. Deceleration region ($x - x_0 = O(\gamma\alpha^{1/2})$)

In the deceleration region (figure 4b), u₀(x, 0) varies by O(1) over a short horizontal scale Δx and Γ₀ ~ α¹/² (where ~ denotes ‘scales like’), based on continuity with the slip region. Thus, Δx ~ γα¹/² for γΓ₀ₓ to vary by O(1) in (C2). Bulk diffusion balances advection by u₀ over a vertical lengthscale y ~ (αΔx)¹/² ~ γ¹/²α³/⁴ using (C1). Diffusion therefore influences the concentration field near the interface, while advection dominates transport in the upper part of the boundary layer. At the downstream limit of this region, u₀ falls towards zero, with Γ₀ₓ ≈ 3β/(2γ). The deceleration region is too short for there to be appreciable desorption so flux conservation requires that u₀ ~ α¹/²/Γ₀ ~ α¹/²γ/(x - x₀) at the outlet of the region. We write the interfacial flux as

$$u_0(x, 0)\Gamma_0 = \frac{3}{\beta} \left( \frac{\alpha(x_0 + \phi)}{3\pi} \right)^{1/2} \equiv \frac{3\alpha^{1/2}}{4}\lambda, \quad \text{where} \quad \lambda = \frac{4}{\beta} \left( \frac{x_0 + \phi}{3\pi} \right)^{1/2}, \quad (C35a, b)$$

as u₀ = 3/4 + ... and Γ₀ is (C33) in the slip region. We then rescale, using

$$x = x_0 + \gamma\alpha^{1/2}X, \quad y = \gamma^{1/2}\alpha^{3/4}Y, \quad \Gamma_0 = \alpha^{1/2}\lambda G(X), \quad c_0 = \alpha^{1/2}\lambda C(X, Y). \quad (C36a-d)$$

Substituting (C36) into (C1)-(C6), the bulk-surfactant equation (C1) in D₁ is given by

$$C_{YY} - \left( \frac{3}{4} - \frac{\lambda}{2\beta} G_X \right) C_X - \frac{\lambda Y}{2\beta} G_{XX} C_Y = 0, \quad (C37)$$
neglecting terms of $O(\alpha^{3/4})$. We have $G(X) = C(X, 0)$ in the strong-exchange limit and the interfacial-surfactant equation (C 2) simplifies to

$$ \left[ G \left( \frac{3}{4} - \frac{\lambda}{2\beta} G_X \right) \right]_X = 0. \qquad (C 38) $$

Integrating (C 38) and matching it to the slip region, where $G \to 1$ for $X \to -\infty$, gives

$$ \frac{3}{4}G - \frac{\lambda}{2\beta}GG_X = \frac{3}{4}, \qquad (C 39) $$

which can be solved directly to give

$$ G + \log(G - 1) = (3\beta/(2\lambda))(X - X_0) \qquad (C 40) $$

for some $X_0$ (as in Jensen & Halpern (1998), for example). Substituting $G$ into (C 35a), we compare $u_0(x, 0)$ with the COMSOL simulations in figure 5(d), showing how $u_0$ falls steeply from 3/4. The core condition (C 5) becomes

$$ C \to \infty \text{ as } Y \to \infty. \qquad (C 41) $$

Equations (C 37)–(C 41) constitute a nonlinear boundary-layer problem that describes the abrupt ejection of bulk concentration towards the core of the channel in figure 3(f). We only consider the downstream limit of (C 37)–(C 41) here. The upstream condition on $C$ is provided from the near-wall component of the concentration in the slip region, which grows in amplitude with $Y$ in (C 32). The downstream limit has self-similar form

$$ G \approx \frac{3\beta}{2\lambda}(X - X_0) + \dots, \quad C \approx \frac{3\beta}{2\lambda}(X - X_0)\hat{C}(\eta) + \dots \qquad (C 42a, b) $$

where the similarity variable is given by $\eta = Y/(X - X_0)$ for some $X_0$ and $\hat{C}$ satisfies $\hat{C}_{\eta\eta} = (\lambda/2\beta)\hat{C}$, which has solutions $\exp[\pm\sqrt{\lambda/2\beta\eta}]$. Concentration contours are constant along lines $Y = \sqrt{2\beta/\lambda}(X - X_0)$, representing ejection of bulk surfactant away from the interface. Expressed in terms of the original variables using (C 36), the downstream limit (C 42) of this region yields

$$ \Gamma_0 \approx \frac{3\beta}{2\gamma}(x - x_0) \quad u_0(x, 0) \approx \frac{\lambda\alpha^{1/2}\gamma}{2\beta(x - x_0)} \quad \text{for } \gamma\alpha^{1/2} \ll x - x_0 \ll \gamma^{3/4}\alpha^{1/8}, \quad (C 43a, b) $$

which delivers an interfacial advective flux $3\lambda\alpha^{1/2}/4$ into a very short transition region near $x = x_0$. Surface diffusion is a subdominant effect in this region provided $\epsilon^2 \ll \gamma/\alpha^{1/2}$. The leading-order interfacial concentration profile in the deceleration region, (C 43a), contributes to the second part of (3.2a) (for $x_0 \le x$) and (C 43b) explains the steep fall in surface velocity in figure 5(d) at the end of the slip region; the incident boundary layer of thickness $\alpha^{1/2}$ is thickened to $\alpha^{1/2}/u_0$ as $u_0$ diminishes (see the arrow in figure 3f). The surface velocity in this regime, computed using (C 35) with $\Gamma_0$ evaluated using (C 39), is shown in figure 5(d).

### C.2.3. Transition region ($x - x_0 = O(\gamma^{3/4}\alpha^{1/8})$)

In the transition region, there is a balance in the bulk between vertical diffusion, advection by $u_0$ (which is $o(1)$) and advection by shear (which is affected by viscous and Marangoni effects). Hence, over a horizontal lengthscale $\Delta x$, $\alpha/y^2 \sim u_0/\Delta x \sim y/\Delta x$; thus $y \sim u_0 \sim (\alpha\Delta x)^{1/3}$. Matching $u_0$ to the deceleration region requires $u_0 \sim \alpha^{1/2}\gamma/\Delta x$, implying $\Delta x \sim \gamma^{3/4}\alpha^{1/8}$ and $y \sim u_0 \sim \gamma^{1/4}\alpha^{3/8}$. The interfacial concentration is linear at leading order ($\Gamma_0 \approx 3\beta(x - x_0)/(2\gamma)$) with a correction of size $u_0\Delta x/\gamma \sim \alpha^{1/2}$. The bulk concentration is set by the dominant interfacial concentration and is $O(\Delta x/\gamma) \sim$
$\alpha^{1/8}/\gamma^{1/4}$. The region is again too short for there to be appreciable desorption, so the interfacial flux is still $O(\alpha^{1/2})$. To provide a bridge between the deceleration region and the quasi-stagnant region further downstream, we reintroduce shear $3Y/2$ to the bulk transport equation, using the scaling

$$x = x_0 + \alpha^{1/8}\gamma^{3/4}X, \quad y = \alpha^{3/8}\gamma^{1/4}Y,$$

$$\Gamma_0 = \frac{3\beta}{2\gamma}(x-x_0) + \alpha^{1/2}G(X), \quad c_0 = \frac{\alpha^{1/8}}{\gamma^{1/4}}C(X, Y). \quad (C44a-d)$$

Substituting (C 44) into (C 1)–(C 6), the bulk transport equation (C 1) in $D_1$ simplifies to

$$C_{YY} + \left( \frac{G_X}{2\beta} - \frac{3Y}{2} \right) C_X - \frac{G_{XX}}{2\beta} Y C_Y = 0. \quad (C45)$$

Again $C(X, 0) = G(X)$ and the interfacial transport equation (C 2) becomes

$$-\left[\frac{3\beta X}{4} G_X\right]_X = 0. \quad (C46)$$

Using (C 46), there is no exchange of surfactant with the bulk to leading order. Integrating (C 46) in $X$ and matching to the deceleration region, $X G_X \to \lambda$ as $X \to -\infty$. This can be integrated in $X$ to give $G = -\lambda \log X + A_1$ for some constant $A_1$, which we substitute into (C 45), to obtain

$$C_{YY} - \frac{\lambda}{2\beta} \left( \frac{1}{X} C_X - \frac{Y C_Y}{X^2} \right) - \frac{3}{2} Y C_X = 0. \quad (C47)$$

This linear boundary-layer problem bridges the scalings $Y \sim X$ for $X \ll 1$ and $Y \sim X^{1/3}$ for $X \gg 1$, but we do not solve the full problem here. Expressed in terms of original variables using (C 44), the downstream limit of the transition region demonstrates that

$$\Gamma_0 \approx \frac{3\beta(x-x_0)}{2\gamma}, \quad u_0(x, 0) \approx \frac{\lambda\alpha^{1/2}\gamma}{2\beta(x-x_0)} \quad \text{for} \quad \gamma^{3/4}\alpha^{1/8} \ll x-x_0 \ll 1. \quad (C48a, b)$$

Note that the interfacial concentration profile in this sub-region, (C 48a), remains linear and exactly the same as in the deceleration region, see (C 48a), by continuity. Over this part of the interface, the interfacial flux is dominated by advection: $u_0(x, 0)\Gamma_0 \approx 3\lambda\alpha^{1/2}/4$. This flux is delivered to the rest of the interface, in the last region that we call the quasi-stagnant region for $x_0 < x < \phi$. The drag reduction can only be determined once the interfacial concentration profile across the whole interface is determined.

### C.2.4. The quasi-stagnant region ($x_0 < x < \phi$)

In the quasi-stagnant region, compression of the almost-immobile interface promotes desorption, returning surfactant to the bulk. Therefore, the interfacial flux of surfactant decreases and the bulk flux increases with $x$. Following the scaling arguments and numerical simulations outlined in §3.1, we write

$$\Gamma_0 = \frac{3\beta}{2\gamma}(x-x_0) + \dots, \quad c_0 = \frac{3\beta}{2\gamma}(x-x_0)C(\eta) + \dots, \quad (C49a, b)$$

where the similarity variable $\eta = y/(\alpha(x-x_0))^{1/3}$ is obtained by balancing advection and diffusion in the bulk-transport equation (C 1). Substituting (C 49) into (C 1)–(C 6), the bulk-surfactant equation (C 1) in $D_1$ simplifies to

$$C_{\eta\eta} + \eta^2 C_{\eta}/2 - 3\eta C/2 = 0, \quad C(0) = 1, \quad C(\infty) = 0. \quad (C50a-c)$$
Equation (C 50) is the same as (C 16), whose solution is given in (C 17b). At the interface $x \in [x_0, \phi]$ and $Y = 0$, the interfacial-surfactant equation (C 2) becomes

$$ \beta[u_0(x, 0)(x - x_0)]_x = (\alpha(x - x_0))^{2/3}C'(0). \quad (C 51) $$

Integrating (C 51) and applying the no-flux condition $u_0(\phi, 0) = 0$, gives

$$ \beta u_0(x, 0) = -\frac{3\alpha^{2/3}C'(0)}{5} \left( \frac{(\phi - x_0)^{5/3}}{x - x_0} - (x - x_0)^{2/3} \right), \quad (C 52) $$

where $C'(0) = -3^{5/3}/(2^{4/3}\Gamma(\frac{1}{3})) \approx -0.92$. This profile successfully captures the COM-SOL prediction in figure 5(d). In figure 5(d), the small discrepancy between the asymptotic results (C 35, C 52) and the 2D long-wave model arises from the retention of surface diffusion in its numerical approximation. Matching (C 52) to the surface velocity in the transition region (C 48) yields

$$ 5\lambda\gamma = -6\alpha^{1/6}C'(0)(\phi - x_0)^{5/3}. \quad (C 53) $$

The constants $\lambda$ and $x_0$ have to be evaluated numerically for a given $\alpha, \beta, \gamma$ and $\phi$ using (C 35) and (C 53). Thus, in the quasi-stagnant region of $A^{2D}$, the interfacial concentration is given by

$$ \Gamma_0 \approx \frac{3\beta}{2\gamma} (x - x_0) = \frac{3\beta}{2\gamma} \left( x - \phi + \left( \frac{-5\lambda\gamma}{6\alpha^{1/6}C'(0)} \right)^{3/5} \right), \quad (C 54) $$

using (C 53). As anticipated, the leading-order interfacial surfactant concentration (C 54) obtained in the quasi-stagnant cap region is linear and equal to the profiles obtained in the two previous short sub-regions: the deceleration (§C.2.2) and the transition regions (§C.2.3). This linear profile agrees with our numerical results in region $A^{2D}$ for the quasi-stagnant part of the interface, i.e. for $x_0 < x < \phi$, as depicted in figure 5(b). Only now can we substitute (C 54) into (2.44), to obtain the leading-order drag reduction in region $A^{2D}$, which is given by

$$ DR_0 = 1 - \frac{1}{2\phi} \left( \frac{-5\lambda}{6C'(0)} \right)^{3/5} \frac{\gamma^{6/10}}{\alpha^{1/10}} = 1 - \frac{a_1\gamma^{6/10}}{\alpha^{1/10}}, \quad (C 55) $$

as in (3.2b). Equations (C 54) and (C 55) show how bulk diffusion, captured in the non-dimensional coefficient $\alpha^{1/10}$, has a weak but non-zero influence on the interfacial surfactant profile that determines $DR_0$. The solution is valid as long as the layers are nested ($\gamma \lesssim \alpha^{-3/2}$) and $x_0 > -\phi$, such that $1 - DR_0 < 1$. In order to determine $x_0$, $\lambda$ and $a_1$, we solve (C 35) and (C 53) numerically to get $\lambda$ and $x_0$ for a given $\alpha, \beta, \gamma$ and $\phi$, from which, we can evaluate $a_1 = (-5\lambda/(6C'(0)))^{3/5}/(2\phi)$ and therefore $DR_0$.

In summary, to acquire the solution in the quasi-stagnant region, where the surfactant gradient decelerates the flow towards the downstream stagnation point, (C 54)–(C 55), we have resolved a nested asymptotic structure which includes deceleration (§C.2.2) and transition (§C.2.3) regions. The nested asymptotic structure has been matched to a slip region at the upstream end of the plastron (§C.2.1), where the surfactant gradient has a weak effect on the flow and slips from the upstream stagnation point. These four asymptotic regions are detailed in figure 4(b), which we used to describe and explain the numerical solution that we see in figure 3(f). Using the expressions for $u_0(x, 0)$ and $\lambda$ in the quasi-stagnant and slip regions, (C 52) and (C 35b), respectively, the surfactant desorbed across the quasi-stagnant region is given by

$$ \int_{x_0}^{\phi} c_{0y}(x, 0) \, dx = - \left( \frac{9\alpha(x_0 + \phi)}{3\pi} \right)^{1/2} + \dots, \quad (C 56) $$
which matches the total amount of surfactant adsorbed in the slip region in (C 34), as required by the conservation of surfactant.

### C.3. Weak Marangoni effect and weak bulk-surface exchange: region $A_E^{2D}$

In region $A_E^{2D}$, we simplify the boundary-layer equations (C 1)–(C 6), where bulk diffusion and bulk-surface exchange are weak compared to Marangoni effects and advection, i.e. $\alpha \ll 1$ and small $\nu$. In this regime, there is a one-way coupling of the bulk concentration onto the interfacial concentration, with $c_0 \approx 1$ everywhere, and the interfacial surfactant distribution is in the classical stagnant-cap regime (He et al. 1991), with a streamwise-averaged concentration imposed by the bulk concentration that is approximately equal to 1. Here, the surfactant gradient generated by the flow renders the downstream end of the interface no slip and the upstream end of the interface shear-free. We expand the surfactant field as follows

$$c_0 = 1 + \alpha C_1(x, y) + \dots, \quad \Gamma_0 = G_0(x) + \alpha G_1(x) + \dots, \quad (C 57a, b)$$

and substitute (C 57) into (C 1)–(C 6). In the interfacial surfactant equation (C 2), there is no flux of bulk surfactant onto the interface, such that $c_{0y} = 0$. For $x \in [-\phi, \phi]$ and $y = 0$, we have

$$\frac{3\beta}{4} G_0 - \frac{\gamma}{2} G_{0x} G_0 = 0, \quad (C 58)$$

at leading order, and where we used the no-flux of interfacial surfactant condition (C 3). The surfactant flux condition (C 6b) in $D_1$ with $c_0 \sim 1$ reduces to

$$\int_{-\phi}^{\phi} G_0 \, dx = 2\phi. \quad (C 59)$$

A piecewise linear solution to the nonlinear ordinary differential equation (C 58) exists for $2\phi > 4\gamma/3\beta$, such that in $A_E^{2D}$ the interfacial-surfactant distribution is given by

$$G_0 = \begin{cases} 0 & \text{for } -\phi < x < x_0, \\ \frac{3\beta}{2\gamma}(x - x_0) & \text{for } x_0 < x < \phi. \end{cases} \quad (C 60)$$

This piecewise linear stagnant-cap profile for the interfacial surfactant concentration agrees with our numerical results shown in figure 7(b). We evaluate $x_0$ using the no-net-flux condition (C 59), such that

$$x_0 = \phi - \left( \frac{8\phi\gamma}{3\beta} \right)^{1/2}, \quad (C 61)$$

where $x_0 > -\phi$ to ensure existence of the no-surfactant part of the profile in (C 60). Hence, for $2\phi > 4\gamma/3\beta$, the leading-order drag reduction calculated by substituting (C 60) into (2.44) is given by

$$DR_0 = 1 - \left( \frac{2\gamma}{3\phi\beta} \right)^{1/2}. \quad (C 62)$$

This solution for the $A_E^{2D}$ region, already given in (3.4c), is equivalent to that in region $A_E^{1D}$ in the strong cross-channel diffusion problem (Tomlinson et al. 2023a). It agrees with numerical results shown in figure 6.
BAIER, T. & HARDT, S. 2021 Influence of insoluble surfactants on shear flow over a surface in Cassie state at large Péclet numbers. *J. Fluid Mech.* **907**, A3.

BAIER, T. & HARDT, S. 2022 Shear flow over a surface containing a groove covered by an incompressible surfactant phase. *J. Fluid Mech.* **949**, A34.

BOLOGNESI, G., COTTIN-BIZONNE, C. & PIRAT, C. 2014 Evidence of slippage breakdown for a superhydrophobic microchannel. *Phys. Fluids* **26**, 082004.

CHENG, Y., XU, J. & SUI, Y. 2015 Numerical study on drag reduction and heat transfer enhancement in microchannels with superhydrophobic surfaces for electronic cooling. *Appl. Therm. Eng.* **88**, 71–81.

CROWDY, D.G., CURRAN, A.E. & PAPAGEORGIOU, D.T. 2023 Fast reaction of soluble surfactant can remobilize a stagnant cap. *J. Fluid Mech.* **969**, A8.

FROSSARD, A.A., GÉRARD, V., DUPLESSIS, P., KINSEY, J.D., LU, X., ZHU, Y., BISGROVE, J., MABEN, J.R., LONG, M.S., CHANG, R.Y. & OTHERS 2019 Properties of seawater surfactants associated with primary marine aerosol particles produced by bursting bubbles at a model air–sea interface. *Environ. Sci. Technol.* **53** (16), 9407–9417.

GAME, S.E., HODES, M., KEAVENY, E.E. & PAPAGEORGIOU, D.T. 2017 Physical mechanisms relevant to flow resistance in textured microchannels. *Phys. Rev. Fluids* **2**, 094102.

HARPER, J.F. 2004 Stagnant-cap bubbles with both diffusion and adsorption rate-determining. *J. Fluid Mech.* **521**, 115–123.

HE, Z., MALDARELLI, C. & DAGAN, Z. 1991 The size of stagnant caps of bulk soluble surfactant on the interfaces of translating fluid droplets. *J. Colloid Interface Sci.* **146**, 442–451.

HOURLIER-FARGETTE, A., DERVAUX, J., ANTKOWIAK, A. & NEUKIRCH, S. 2018 Extraction of silicone uncrosslinked chains at air–water–polydimethylsiloxane triple lines. *Langmuir* **34**, 12244–12250.

JENSEN, O.E. & HALPERN, D. 1998 The stress singularity in surfactant-driven thin-film flows. part 1. viscous effects. *J. Fluid Mech.* **372**, 273–300.

KIM, T.J. & HIDROVO, C. 2012 Pressure and partial wetting effects on superhydrophobic friction reduction in microchannel flow. *Phys. Fluids* **24**, 112003.

KIRK, T.L., KARAMANIS, G., CROWDY, D.G. & HODES, M. 2020 Thermocapillary stress and meniscus curvature effects on slip lengths in ridged microchannels. *J. Fluid Mech.* **894**, A15.

LAM, L.S., HODES, M. & ENRIGHT, R. 2015 Analysis of galinstan-based microgap cooling enhancement using structured surfaces. *J. Heat Transf.* **137**, 091003.

LANDEL, J.R., PEAUDECERF, F.J., TEMPRANO-COLETO, F., GIBOU, F., GOLDSTEIN, R.E. & LUZZATTO-FEGIZ, P. 2020 A theory for the slip and drag of superhydrophobic surfaces with surfactant. *J. Fluid Mech.* **883**, A18.

LEE, C., CHOI, C.H. & KIM, C.J. 2016 Superhydrophobic drag reduction in laminar flows: a critical review. *Exp. Fluids* **57**, 1–20.

MANIKANTAN, H & SQUIRES, T M 2020 Surfactant dynamics: hidden variables controlling fluid flows. *J. Fluid Mech.* **892**, P1.

MAYER, M.D. & CROWDY, D.G. 2022 Superhydrophobic surface immobilisation by insoluble surfactant. *J. of Fluid Mech.* **949**, A18.

NG, C.O. & WANG, C.Y. 2009 Stokes shear flow over a grating: implications for superhydrophobic slip. *Phys. Fluids* **21**, 087105.

OGUZ, H.N. & SADHAL, S.S. 1988 Effects of soluble and insoluble surfactants on the motion of drops. *J. Fluid Mech.* **194**, 563–579.

PALAPARTHI, R., PAPAGEORGIOU, D.T. & MALDARELLI, C. 2006 Theory and experiments on the stagnant cap regime in the motion of spherical surfactant-laden bubbles. *J. Fluid Mech.* **559**, 1–44.

PARK, H., CHOI, C.H. & KIM, C.J. 2021 Superhydrophobic drag reduction in turbulent flows: a critical review. *Exp. Fluids* **62**, 1–29.

PARK, H., SUN, G. & KIM, C.J. 2014 Superhydrophobic turbulent drag reduction as a function of surface grating parameters. *J. Fluid Mech.* **747**, 722–734.

PEAUDECERF, F.J., LANDEL, J.R., GOLDSTEIN, R.E. & LUZZATTO-FEGIZ, P. 2017 Traces of surfactants can severely limit the drag reduction of superhydrophobic surfaces. *Proc. Natl. Acad. Sci.* **114**, 7254–7259.

PEREIRA, R., ASHTON, I., SABBAGHZADEH, B., SHUTLER, J.D. & UPSTILL-GODDARD, R.C.
2018 Reduced air-sea co 2 exchange in the atlantic ocean due to biological surfactants. *Nat. Geosci.* **11**, 492–496.

RASTEGARI, A. & AKHAVAN, R. 2019 On drag reduction scaling and sustainability bounds of superhydrophobic surfaces in high reynolds number turbulent flows. *J. Fluid Mech.* **864**, 327–347.

RODRIGUEZ-BROADBENT, H. & CROWDY, D.G. 2023 Superhydrophobic surfaces with recirculating interfacial flow due to surfactants are ‘effectively’ immobilized. *J. Fluid Mech.* **956**, R3.

SCHÄFFEL, D., KOYNOV, K., VOLLMER, D., BUTT, H.J. & SCHÖNECKER, C. 2016 Local flow field and slip length of superhydrophobic surfaces. *Phys. Rev. Lett.* **116**, 134501.

SCHÖNECKER, C., BAIER, T. & HARDT, S. 2014 Influence of the enclosed fluid on the flow over a microstructured surface in the Cassie state. *J. Fluid Mech.* **740**, 168–195.

SEO, J. & MANI, A. 2018 Effect of texture randomization on the slip and interfacial robustness in turbulent flows over superhydrophobic surfaces. *Phys. Rev. Fluids* **3**, 044601.

SONG, D., SONG, B., HU, H., DU, X., DU, P., CHOI, C.H. & ROTHSTEIN, J.P. 2018 Effect of a surface tension gradient on the slip flow along a superhydrophobic air-water interface. *Phys. Rev. Fluids* **3**, 033303.

SUNDIN, J. & BAGHERI, S. 2022 Slip of submerged two-dimensional liquid-infused surfaces in the presence of surfactants. *J. Fluid Mech.* **950**, A35.

TASOGLU, S., DEMIRCI, U. & MURADOGLU, M. 2008 The effect of soluble surfactant on the transient motion of a buoyancy-driven bubble. *Phys. Fluids* **20** (4).

TEMPRANO-COLETO, F., SMITH, S.M., PEAUDECERF, F.J., LANDEL, J.R., GIBOU, F. & LUZZATTO-FEGIZ, P. 2023 A single parameter can predict surfactant impairment of superhydrophobic drag reduction. *Proc. Natl. Acad. Sci.* **120** (3), e2211092120.

TEO, C.J. & KHOO, B.C. 2010 Flow past superhydrophobic surfaces containing longitudinal grooves: effects of interface curvature. *Microfluid Nanofluidics* **9** (2), 499–511.

TOMLINSON, S.D., GIBOU, F., LUZZATTO-FEGIZ, P., TEMPRANO-COLETO, F., JENSEN, O.E. & LANDEL, J.R. 2023a Laminar drag reduction in surfactant-contaminated superhydrophobic channels. *J. Fluid Mech.* **963**, A10.

TOMLINSON, S.D., GIBOU, F., LUZZATTO-FEGIZ, P., TEMPRANO-COLETO, F., JENSEN, O.E. & LANDEL, J.R. 2024 Unsteady evolution of slip and drag in surfactant-contaminated superhydrophobic channels. *J. Fluid Mech.* **1000**, A62.

TOMLINSON, S.D., PEAUDECERF, F.J., TEMPRANO-COLETO, F., GIBOU, F., LUZZATTO-FEGIZ, P., JENSEN, O.E. & LANDEL, J.R. 2023b A model for slip and drag in turbulent flows over superhydrophobic surfaces with surfactant. *Int. J. Heat Fluid Flow* **103**, 109171.

TREFETHEN, L.N. 2000 *Spectral methods in MATLAB*. SIAM.

TÜRK, S., DASCHIEL, G., STROH, A., HASEGAWA, Y. & FROHNAPFEL, B. 2014 Turbulent flow over superhydrophobic surfaces with streamwise grooves. *J. Fluid Mech.* **747**, 186–217.

XU, M., GRABOWSKI, A., YU, N., KEREZYTE, G., LEE, J., PFEIFER, B.R. & KIM, C. 2020 Superhydrophobic drag reduction for turbulent flows in open water. *Phys. Rev. Appl.* **13** (3), 034056.