# Laminar drag reduction in surfactant-contaminated superhydrophobic channels

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.. Laminar drag reduction in surfactant-contaminated superhydrophobic channels. Journal of Fluid Mechanics, 2023, 963, pp.A10. 10.1017/jfm.2023.264. hal-04651750

HAL Id: hal-04651750

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

Submitted on 17 Jul 2024

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
<figure><img src="image_2.png" /></figure>

# Laminar drag reduction in surfactant-contaminated superhydrophobic channels

Samuel D. Tomlinson<fnref n="1" />,†, Frédéric Gibou<fnref n="2" />, Paolo Luzzatto-Fegiz<fnref n="2" />, Fernando Temprano-Coleto<fnref n="3" />, Oliver E. Jensen<fnref n="1" /> and Julien R. Landel<fnref n="1" />

(Received 26 October 2022; revised 12 February 2023; accepted 27 March 2023)

Although superhydrophobic surfaces (SHSs) show promise for drag reduction applications, their performance can be compromised by traces of surfactant, which generate Marangoni stresses that increase drag. This question is addressed for soluble surfactant in a three-dimensional laminar channel flow, with periodic SHSs made of long finite-length longitudinal grooves located on both walls. We assume that bulk diffusion is sufficiently strong for cross-channel concentration gradients to be small. Exploiting long-wave theory and accounting for the difference between the rapid transverse and slower longitudinal Marangoni flows, we derive a one-dimensional model for surfactant transport from the full three-dimensional transport equations. Our one-dimensional model allows us to predict the drag reduction and surfactant distribution across the parameter space. The system exhibits multiple regimes, involving competition between Marangoni effects, bulk and interfacial diffusion, bulk and interfacial advection, shear dispersion and surfactant exchange between the bulk and the interface. We map out asymptotic regions in the high-dimensional parameter space, and derive explicit closed-form approximations of the drag reduction, without any fitting or empirical parameters. The physics underpinning the drag reduction effect and the negative effect of surfactant is discussed through analysis of the velocity field and surfactant concentrations, which show both uniform and non-uniform stress distributions. Our theoretical predictions of the drag reduction compare well with results from the literature solving numerically the full three-dimensional transport problem. Our atlas of maps provides a comprehensive analytical guide for designing surfactant-contaminated channels with SHSs, to maximise the drag reduction in applications.

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

† Email address for correspondence: samuel.tomlinson@manchester.ac.uk

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

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

<footnote n="3">³Andlinger Center for Energy and the Environment, Princeton University, Princeton, NJ 08544, USA</footnote>
# 1. Introduction

When a fluid flows over a superhydrophobic surface (SHS), surface chemistry and microscopic roughness combine to entrap an array of microscopic gas pockets at the SHS. In the non-wetted or Cassie-Baxter state (see, e.g., the reviews of Rothstein 2010; Lee, Choi & Kim 2016; Park, Choi & Kim 2021), the reduced liquid-solid contact area at the SHS reduces viscous drag at the surface when compared with a solid wall, which leads to an overall decrease in drag. Owing to these drag-reducing capabilities, SHSs have been considered in a number of applications, such as low-Reynolds-number laminar flows (e.g. Ou & Rothstein 2005; Sbragaglia & Prosperetti 2007; Rothstein 2010; Schönecker, Baier & Hardt 2014; Lee et al. 2016; Landel et al. 2020), high-Reynolds-number turbulent flows (e.g. Park, Sun & Kim 2014; Türk et al. 2014; Golovin et al. 2016; Seo & Mani 2018; Rastegari & Akhavan 2019; Park et al. 2021) and the thermal management of microelectronics (e.g. Baier, Steffes & Hardt 2010; Cheng, Xu & Sui 2015; Lam, Hodes & Enright 2015; Kirk et al. 2020).

Although the physical mechanism behind laminar drag reduction is well understood, SHSs seldom achieve the high drag-reduction performance predicted by theory (Lee et al. 2016; Park et al. 2021). A number of practical difficulties can explain the lower performance measured experimentally, such as displacement of the liquid-gas interfaces (Biben & Joly 2008; Ng & Wang 2009), contact-angle effects (Sbragaglia & Prosperetti 2007; Teo & Khoo 2010), viscous drag from the gas phase (Schönecker et al. 2014; Game et al. 2017), as well as surfactant-induced Marangoni stresses (Peaudecerf et al. 2017; Landel et al. 2020). In this study, we investigate in detail the effect of soluble surfactant on the drag reduction of SHSs in laminar channel flows. Surfactants can be transported by a liquid, adsorb at liquid-gas interfaces and lower the surface tension of these interfaces (Manikantan & Squires 2020). Transported by the flow along the liquid-gas interface, surfactants can accumulate at stagnation points (liquid-gas-solid contact lines), inducing an adverse Marangoni stress at the interface which increases the drag, thereby negating the drag-reducing effect of the nominally shear-free interface.

The effect of surfactants on SHS drag reduction was suggested by experimental studies in the last decade, which reported only a modest drag reduction when compared with solid walls. For example, Kim & Hidrovo (2012) used microscopy to analyse the location of the liquid-gas interface in a channel flow over transverse ridges. Combined with flow-rate measurements, they found that the frictional properties of the channel were insensitive to the degree of microtexture wetting, and hence closely resembled solid walls. Likewise, using fluorescence microscopy and passive tracers, Bolognesi, Cottin-Bizonne & Pirat (2014) studied the flow field and shape of the liquid-gas interfaces for a channel flow over longitudinal ridges. They found non-zero shear stresses at the interfaces, whereas most theoretical and numerical studies normally employ shear-free conditions to model SHSs (Rothstein 2010).

Following these observations, several studies have explored the effects of surfactants on SHS drag reduction. Schäffel et al. (2016) combined fluorescence correlation spectroscopy with numerical simulations to examine the effective slip length in a channel flow with cylindrical pillars. They found that the effective slip length reduced when compared with simulations of a surfactant-free channel. Peaudecerf et al. (2017) performed experiments and numerical simulations in a channel with longitudinal ridges, demonstrating that low levels of surfactant could yield large changes in interfacial conditions. They showed that liquid-gas interfaces could be rendered no-slip at surfactant concentrations well below typical environmental values. They used channel-flow experiments where the imposed pressure gradient was removed abruptly to show that the observed reverse flow was
only possible due to adverse out-of-equilibrium surfactant gradients. Song *et al.* (2018) performed experiments for finite longitudinal ridges and effectively infinite concentric, annular ridges. For finite longitudinal ridges, surface tension gradients arose due to the presence of downstream transverse contact lines, increasing the drag when compared with infinite concentric ridges, where there were no stagnation points.

The experiments of Peaudecerf *et al.* (2017) and Song *et al.* (2018) were conducted in nominally ‘clean’ channels without added surfactant. Schäffel *et al.* (2016) also performed experiments with added surfactant, finding barely measurable increases in drag relative to their nominally ‘clean’ experiments. As noted by Peaudecerf *et al.* (2017), the counterintuitive result reported by Schäffel *et al.* (2016) was most likely due to the fact that traces of surfactants were already present in their experiment. Surfactant traces are inherently present in most engineered systems, due to general manufacturing conditions, materials, and the high surface-to-volume ratio of the SHS textures. For instance, microfluidic experimental devices are often made of polydimethylsiloxane (PDMS), which is known to lead to surfactant effects (Hourlier-Fargette *et al.* 2018). In nature, surfactant traces have been measured in sea water (Pereira *et al.* 2018), rivers, estuaries and fog (Lewis 1991; Facchini *et al.* 2000).

Due to the potentially strong adverse effect of surfactants on the drag reduction performance of SHSs, it is crucial to model their effect. Theoretical models can explain experimental observations, and also provide predictions for the design of SHSs where surfactant effects are mitigated. However, as discussed by Landel *et al.* (2020), the theoretical modelling of flows inclusive of surfactant over SHSs is complex. To make progress, Landel *et al.* (2020) assumed that the surfactant concentration is small. This assumption is consistent with normal environmental conditions where surfactants are found in trace amounts (Peaudecerf *et al.* 2017). Landel *et al.* (2020) constructed a scaling theory to model the slip and drag in steady two-dimensional (2D) pressure-driven channel flows bounded by SHSs, made of long transverse gratings, in the low-Reynolds-number flow regime. They also performed finite-element 2D numerical simulations to compute the constants in the scaling theory, and validated their scaling predictions. Their results showed that the slip length and drag reduction are affected by surfactants across a broad range of the parameter space. They focused on parameter combinations associated with microfluidic applications; it remains of interest to perform a comprehensive asymptotic analysis of the parameter space. Moreover, their theory did not consider non-uniform interfacial surfactant distributions associated with the ‘stagnant cap’ regime, where the upstream region of the interface has a weak surfactant gradient and is almost shear-free, whereas the downstream region of the interface has a strong surfactant gradient and is effectively no-slip (see, e.g., Bond & Newton 1928; Frumkin & Levich 1947; Levich 1962; He, Maldarelli & Dagan 1991, which describe similar behaviour for air bubbles rising in surfactant-contaminated water). In this regime, the advection of surfactant at the interface dominates relative to surface diffusion and bulk-surface exchange.

The scaling theory of Landel *et al.* (2020) was extended to three-dimensional (3D) low-Reynolds-number channel flows with long but finite longitudinal gratings by Temprano-Coleto *et al.* (2023). They used a long-wave limit and assumed low surfactant concentrations to show that the slip velocity and the slip length scale in a similar fashion to 2D flows. In the limit of large Damköhler numbers, as found for common surfactants, and sufficiently small bulk Péclet numbers, they predicted that significant slip can be achieved provided that the grating length is longer than both a modified depletion length and a mobilisation length. The modified depletion length depends on surfactant properties and the height of the channel, and is generally small, of the order of 1 mm, for small-scale applications. The mobilisation length depends on the normalised surfactant concentration,
Marangoni number, Damköhler number and Biot number. As the mobilisation length can be much longer for typical small-scale applications, they concluded that the mobilisation length alone controls the slip in most small-scale applications. If the grating length is smaller than the mobilisation length, they showed that the slip velocity increases with the square of the grating length, which is in agreement with their numerical simulations and experiments, as well as with the experimental results from Peaudecerf *et al.* (2017) and Song *et al.* (2018).

The contamination of SHS channels with surfactant may be simplified if one considers only the effect of insoluble surfactant at the liquid–gas interface, neglecting the exchange with bulk surfactant which remains at some background concentration. Baier & Hardt (2021) found analytical expressions for the velocity field and interfacial surfactant distribution in a 2D channel flow over transverse grooves in an advection-dominated regime (i.e. large interfacial Péclet numbers), under the assumption that the surface tension depends linearly on the surface concentration. Under the assumption of insoluble surfactant, drag is not reduced by increasing streamwise groove length, in contrast to findings from experiments and soluble-surfactant models, where slip increases with the square of the groove length (Peaudecerf *et al.* 2017; Landel *et al.* 2020; Temprano-Coleto *et al.* 2023). For periodic transverse grooves, Baier & Hardt (2021) determined the effective slip length, which is strongly dependent on the Marangoni number for large gas fractions, decreasing rapidly as the Marangoni number increases from zero. Mayer & Crowdy (2022) also considered surface immobilisation due to insoluble surfactant in a shear flow over periodic transverse SHSs; however, unlike Baier & Hardt (2021), the authors varied the surfactant load, interfacial Péclet and Marangoni number, using a nonlinear equation of state. They combined asymptotic theory with numerical solutions to identify two distinct mechanisms behind surface immobilisation. The first is the previously discussed stagnant cap mechanism. The second immobilisation mechanism emerges when there is a region of near-maximal surfactant concentration close to the downstream stagnation point. In contrast to the stagnant cap mechanism, surfactant gradients need not be large for appreciable Marangoni stresses to develop. The near-maximal surfactant concentration mechanism is captured using a nonlinear equation of state and thus can appear outside of the advection-dominated region, i.e. for small Péclet and Marangoni numbers, provided there is sufficient surfactant present at the interface.

Liquid-infused surfaces (LISs) offer an alternative to SHSs where the SHS pockets are filled with a lubricating immiscible fluid instead of gas (Wong *et al.* 2011; Wexler, Jacobi & Stone 2015). They can self-repair and are more robust than SHSs if designed properly. However, the drag reduction of LISs decreases as the viscosity ratio between internal and external fluids increases (Schönecker *et al.* 2014). Furthermore, recent numerical simulations by Sundin & Bagheri (2022) have indicated that LISs may be more susceptible to surfactant effects than SHSs. Sundin & Bagheri (2022) extended the theory introduced in Landel *et al.* (2020) to account for surfactants in a 2D shear flow, predicting the critical surfactant concentration for the slip to be reduced appreciably: $\hat{C}_c = 4 \times 10^{-4} \text{ mol m}^{-3}$ for water–air SHSs and $\hat{C}_c = 5 \times 10^{-5} \text{ mol m}^{-3}$ for water–dodecane LISs. For low applied shear stresses, Sundin & Bagheri (2022) found that the distribution of surfactant at the interface is approximately uniform and a scaling theory was used to derive an expression for the slip length. For high applied shear stresses, surfactant accumulates at the downstream stagnation point and forms a stagnant cap. Sundin & Bagheri (2022) employed numerical simulations to find that the stagnant cap regime only exists below a particular bulk concentration (above which the scaling theory once again becomes valid). Considering a typical surfactant, e.g. sodium dodecyl sulfate (SDS), the critical bulk concentration is shown to be proportional to the inverse Marangoni number.
Surfactant-contaminated superhydrophobic channels

In this paper, we show that a one-dimensional (1D) asymptotic theory, derived as coupled nonlinear ordinary differential equations (ODEs) from the full 3D transport problem, can capture the impact of surfactant in 3D channels bounded by periodic SHSs made of long finite-length longitudinal grooves. The main assumption behind our theory is that bulk diffusion is strong enough to suppress cross-channel concentration gradients. We thereby sidestep the need for introducing empirical or fitting coefficients associated with scaling analyses, as done in previous studies (Landel et al. 2020; Sundin, Zaleski & Bagheri 2021; Temprano-Coleto et al. 2023). This allows us to address non-uniform shear stresses at the liquid-gas interface. By mapping the drag reduction across a large part of the high-dimensional parameter space, we identify a multitude of asymptotic regions and their boundaries, unexplored by previous studies. Explicit closed-form asymptotic solutions predict the drag reduction in all the regions of the parameter space studied, offering analytical predictions for practical use where numerical simulation of the fluid and surfactant equations are computationally expensive. By addressing the role of shear dispersion, the theory developed here constitutes also a stepping stone towards a wider class of laminar flows with weaker cross-channel diffusion, where the bulk concentration field varies in three dimensions and must be resolved numerically.

The paper is arranged as follows. In § 2, the full 3D transport problem is formulated in terms of nine dimensionless parameters. In § 3, an asymptotic model is derived for the flow and surfactant transport; cross-channel integration reduces the number of independent dimensionless parameters to six. In § 4, key results are presented for the drag reduction using 2D maps that illustrate the structure of the parameter space; the underlying physics is described using the surfactant distribution and the 3D velocity field of the channel flow. In § 5, the implications and extensions of this study are discussed; a table of five dimensionless groups that control drag and the negative effect of surfactant is presented, expressing these dimensionless groups in terms of the dimensional parameters of the problem.

## 2. Formulation

### 2.1. Governing equations

Consider a steady 3D laminar channel flow contaminated with a soluble surfactant and bounded between two SHSs that are separated by a distance $2\hat{H}$, as illustrated in figure 1. Hats indicate dimensional quantities. The $\hat{x}$-, $\hat{y}$- and $\hat{z}$-coordinates are oriented in the streamwise, wall-normal and transverse directions, with $\hat{x} = (\hat{x}, \hat{y}, \hat{z})$. The liquid is assumed to be incompressible and Newtonian with dynamic viscosity $\hat{\mu}$, velocity $\hat{u} = (\hat{u}(\hat{x}), \hat{v}(\hat{x}), \hat{w}(\hat{x}))$, pressure $\hat{p}(\hat{x})$, bulk surfactant distribution $\hat{c}(\hat{x})$ and interfacial surfactant distribution $\hat{\Gamma}(\hat{x}, \hat{z})$. Owing to the periodicity of the geometry, we restrict attention to a single periodic cell with streamwise (transverse) period length $2\hat{P}_x$ ($2\hat{P}_z$), liquid-gas interface length (width) $2\phi_x\hat{P}_x$ ($2\phi_z\hat{P}_z$) and gas fraction $\phi_x$ ($\phi_z$). Liquid-gas interfaces, or plastrons, are assumed to be flat. The domain is partitioned into two subdomains per period, that are bounded by the plastron and solid ridge, namely

$$ \hat{D}_1 = \{\hat{x} \in [-\phi_x \hat{P}_x, \phi_x \hat{P}_x]\} \times \{\hat{y} \in [0, 2\hat{H}]\} \times \{\hat{z} \in [-\hat{P}_z, \hat{P}_z]\}, \quad (2.1a) $$

$$ \hat{D}_2 = \{\hat{x} \in [\phi_x \hat{P}_x, (2 - \phi_x) \hat{P}_x]\} \times \{\hat{y} \in [0, 2\hat{H}]\} \times \{\hat{z} \in [-\hat{P}_z, \hat{P}_z]\}, \quad (2.1b) $$
<figure><img src="image_3.png" /><figcaption>Figure 1. (a) Schematic depicting a plane periodic streamwise channel flow (illustrated by the array of arrows) of a liquid transporting a soluble surfactant. The origin of the Cartesian coordinate system, $\hat{x} = 0$, is located in the middle of the bottom interface. The channel height in the wall-normal direction is $2\hat{H}$, the transverse gas fraction, $\phi_z$, and the transverse period, $2\hat{P}_z$. On the top and bottom SHSs are no-slip ridges, $\hat{R}$ (outlined in red), solid surfaces, $\hat{S}$ (yellow), and liquid-gas interfaces, $\hat{I}$ (green) onto which surfactants can adsorb and desorb, modifying the interfacial stress through the Marangoni effect. (b) Cross-sectional view of the periodic domain at $\hat{z} = 0$, showing the streamwise gas fraction, $\phi_x$, streamwise period, $2\hat{P}_x$, and highlighting domains $\hat{D}_1$ and $\hat{D}_2$ (magenta).</figcaption></figure>

as outlined in figure 1(b). At the SHSs, $\hat{y} = 0$ and $\hat{y} = 2\hat{H}$, we define the liquid-gas interfaces, the ridge surfaces and the solid surfaces, respectively, as

$$ \hat{\mathcal{I}} = \{\hat{x} \in [-\phi_x \hat{P}_x, \phi_x \hat{P}_x]\} \times \{\hat{y} \in \{0, 2\hat{H}\}\} \times \{\hat{z} \in [-\phi_z \hat{P}_z, \phi_z \hat{P}_z]\}, \quad (2.2a) $$

$$ \hat{\mathcal{R}} = \{\hat{x} \in [-\phi_x \hat{P}_x, \phi_x \hat{P}_x]\} \times \{\hat{y} \in \{0, 2\hat{H}\}\} \times \{\hat{z} \in [-\hat{P}_z, -\phi_z \hat{P}_z] \cup [\phi_z \hat{P}_z, \hat{P}_z]\}, \quad (2.2b) $$

$$ \hat{S} = \{\hat{x} \in [\phi_x \hat{P}_x, (2 - \phi_x) \hat{P}_x]\} \times \{\hat{y} \in \{0, 2\hat{H}\}\} \times \{\hat{z} \in [-\hat{P}_z, \hat{P}_z]\}, \quad (2.2c) $$

as outlined in figure 1(a). Note that placing the origin of the coordinate system at the centre of the domain (common practice in the literature for flows over SHSs) would increase the number of matching conditions between subdomains from two to three.

To model the fluid we use the steady Stokes equations, neglecting inertia and any body forces. The bulk surfactant is coupled to the flow field by a steady advection-diffusion equation. In $\hat{D}_1$ and $\hat{D}_2$, we therefore have

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

where $\hat{D}$ is the surfactant bulk diffusivity. The interfacial surfactant is coupled to the flow by a steady advection-diffusion equation and an equation of state. The bulk concentration is coupled to the interfacial surfactant by continuity of flux, where exchange at the interface is modelled using a source-sink term consistent with the Henry isotherm (Chang & Franses 1995). The equation of state and adsorption-desorption kinetics are linearised, which is valid for small deviations in the concentration of surfactant away from some reference value (Manikantan & Squires 2020). On $\hat{I}$, we balance the tangential components of the stress, $\hat{T} \cdot \hat{n}$ where $\hat{T} = -\hat{p}\hat{I} + \hat{\mu}(\hat{\nabla}\hat{u} + (\hat{\nabla}\hat{u})^T)$ and $\hat{n}$ is the unit
normal to $\hat{\mathcal{I}}$ (pointing into the channel), with tangential gradients of surface tension $\hat{\sigma}$. We assume that (i) the effects arising from the gas trapped in the SHSs are negligible, (ii) $\hat{\sigma} = \hat{\sigma}_0 - \hat{A}(\hat{\Gamma} - \hat{\Gamma}_0)$ where $\hat{\sigma}_0$ is the reference surface tension, $\hat{A}$ is the surface activity and $\hat{\Gamma}_0$ is the reference surface concentration of surfactant and (iii) the surface tension remains large enough to suppress deflections of the interface from its assumed flat state. Thus, along $\hat{\mathcal{I}}$, we impose the tangential stress balance in the streamwise and transverse directions, no-penetration of velocity, continuity of surfactant flux and the transport equation for interfacial surfactant

$$ \left. \begin{array}{l} \hat{\mu} \mathbf{n} \cdot \hat{\nabla} \hat{u} - \hat{A} \hat{\Gamma}_{\hat{x}} = 0, \quad \hat{\mu} \mathbf{n} \cdot \hat{\nabla} \hat{w} - \hat{A} \hat{\Gamma}_{\hat{z}} = 0, \quad \hat{v} = 0, \quad \hat{D} \mathbf{n} \cdot \hat{\nabla} \hat{c} - \hat{K}_a \hat{c} + \hat{K}_d \hat{\Gamma} = 0, \\ \hat{D}_I(\hat{\Gamma}_{\hat{x}\hat{x}} + \hat{\Gamma}_{\hat{z}\hat{z}}) + \hat{K}_a \hat{c} - \hat{K}_d \hat{\Gamma} - (\hat{u}\hat{\Gamma})_{\hat{x}} - (\hat{w}\hat{\Gamma})_{\hat{z}} = 0, \end{array} \right\} \quad (2.4a-e) $$

with $\hat{D}_I$ the surfactant interfacial diffusivity, $\hat{K}_a$ the adsorption rate and $\hat{K}_d$ the desorption rate. On $\partial\hat{\mathcal{I}}$ (the contact line bounding the interfaces), no-flux of surfactant requires that

$$ \hat{u}\hat{\Gamma} - \hat{D}_I\hat{\Gamma}_{\hat{x}} = 0 \quad \text{at} \quad \hat{x} = \pm\phi_x\hat{P}_x, \quad \hat{w}\hat{\Gamma} - \hat{D}_I\hat{\Gamma}_{\hat{z}} = 0 \quad \text{at} \quad \hat{z} = \pm\phi_z\hat{P}_z. \quad (2.5a,b) $$

Along $\hat{\mathcal{R}}$ and $\hat{\mathcal{S}}$, we impose no-slip, no-penetration of velocity and no-flux of bulk surfactant

$$ \hat{u} = 0, \quad \hat{w} = 0, \quad \hat{v} = 0, \quad \hat{c}_{\hat{y}} = 0. \quad (2.6a-d) $$

Defining $\hat{q} = (\hat{u}, \hat{p}_{\hat{x}}, \hat{c})$, periodicity across the unit cell $\hat{D}_1 \cup \hat{D}_2$ means that

$$ \hat{q}(-\phi_x \hat{P}_x, \hat{y}, \hat{z}) = \hat{q}((2 - \phi_x) \hat{P}_x, \hat{y}, \hat{z}), \quad \hat{q}(\hat{x}, \hat{y}, -\hat{P}_z) = \hat{q}(\hat{x}, \hat{y}, \hat{P}_z). \quad (2.7a,b) $$

The bulk flow and concentration are continuous at the boundary between $\hat{D}_1$ and $\hat{D}_2$

$$ \hat{q}((\phi_x \hat{P}_x)^-, \hat{y}, \hat{z}) = \hat{q}((\phi_x \hat{P}_x)^+, \hat{y}, \hat{z}), \quad (2.8) $$

where the superscripts - and + mean that the boundary condition is evaluated in $\hat{D}_1$ and $\hat{D}_2$, respectively.

Owing to the symmetry about $\hat{y} = \hat{H}$, the top and bottom interfaces are assumed to have the same distribution of surfactant. We can integrate (2.3a-c)-(2.8) across the channel to show how the streamwise bulk surfactant flux changes as surfactants adsorb and desorb

$$ \frac{d}{d\hat{x}} \int_{\hat{z}=-\hat{P}_z}^{\hat{P}_z} \int_{\hat{y}=0}^{2\hat{H}} (\hat{u}\hat{c} - \hat{D}\hat{c}_{\hat{x}}) d\hat{y} d\hat{z} = 2 \int_{\hat{z}=-\phi_z \hat{P}_z, \hat{y}=0}^{\phi_z \hat{P}_z} (\hat{K}_d \hat{\Gamma} - \hat{K}_a \hat{c}) d\hat{z}; \quad (2.9) $$

and likewise how the streamwise interfacial surfactant flux changes along the plastron

$$ \frac{d}{d\hat{x}} \int_{\hat{z}=-\phi_z \hat{P}_z, \hat{y}=0}^{\phi_z \hat{P}_z} (\hat{u}\hat{\Gamma} - \hat{D}_I\hat{\Gamma}_{\hat{x}}) d\hat{z} = - \int_{\hat{z}=-\phi_z \hat{P}_z, \hat{y}=0}^{\phi_z \hat{P}_z} (\hat{K}_d \hat{\Gamma} - \hat{K}_a \hat{c}) d\hat{z}. \quad (2.10) $$

The boundary-value problem (2.3a-c)-(2.8) can also be integrated to show that the total flux of liquid $\hat{Q}$ is uniform along the streamwise length of the channel

$$ \hat{Q} = \int_{\hat{z}=-\hat{P}_z}^{\hat{P}_z} \int_{\hat{y}=0}^{2\hat{H}} \hat{u} \, d\hat{y} \, d\hat{z}, \quad (2.11) $$
the total flux of surfactant $\hat{K}$ is uniform along the streamwise length of the channel

$$ \hat{K} = \int_{\hat{z}=-\hat{P}_z}^{\hat{P}_z} \int_{\hat{y}=0}^{2\hat{H}} (\hat{u}\hat{c} - \hat{D}\hat{c}_{\hat{x}}) \, d\hat{y} \, d\hat{z} + 2 \int_{\hat{z}=-\phi_z \hat{P}_z, \hat{y}=0}^{\phi_z \hat{P}_z} (\hat{u}\hat{\Gamma} - \hat{D}_I \hat{\Gamma}_{\hat{x}}) \, d\hat{z}, \quad (2.12) $$

where the factor of 2 in front of the second integral accounts for both top and bottom symmetrical SHSs, and the adsorption–desorption flux of surfactant between the interface and the bulk integrates to zero over the whole interface owing to conservation of mass,

$$ \int_{\hat{x}=-\phi_x \hat{P}_x}^{\phi_x \hat{P}_x} \int_{\hat{z}=-\phi_z \hat{P}_z, \hat{y}=0}^{\phi_z \hat{P}_z} (\hat{K}_d \hat{\Gamma} - \hat{K}_a \hat{c}) \, d\hat{z} \, d\hat{x} = 0. \quad (2.13) $$

Integrating (2.9) plus twice (2.10) with respect to $\hat{x}$ recovers (2.12).

A key quantity of interest in the present study is drag reduction. The flow is driven in the streamwise $\hat{x}$-direction by a cross-channel-averaged pressure drop per period given by $\Delta\hat{p} \equiv \langle\hat{p}\rangle(-\phi_x \hat{P}_x) - \langle\hat{p}\rangle((2 - \phi_x)\hat{P}_x) > 0$, where $\langle\cdot\rangle \equiv \int_{\hat{z}=-\hat{P}_z}^{\hat{P}_z} \int_{\hat{y}=0}^{2\hat{H}} \cdot d\hat{y} \, d\hat{z}/(4\hat{P}_z\hat{H})$ is the cross-channel average. In the limit when $\hat{\mathcal{I}}$ is immobilised (yielding an effective no-slip boundary condition for the velocity), we have $\Delta\hat{p} = \Delta\hat{p}_R$, say, and in the uncontaminated (or surfactant-free) limit $\Delta\hat{p} = \Delta\hat{p}_U$ (when $\hat{\mathcal{I}}$ is a shear-free surface). We define the normalised drag reduction as

$$ DR = \frac{\Delta \hat{p}_R - \Delta \hat{p}}{\Delta \hat{p}_R - \Delta \hat{p}_U}, \quad (2.14) $$

which varies from $DR = 0$ to 1 in the cases of minimum and maximum drag reduction, depending on the surfactant-induced Marangoni stresses. Equation (2.14) differs from other definitions of the drag reduction which compare the flow over a SHS with a solid wall, such as $1 - \Delta\hat{p}/\Delta\hat{p}_R$ (e.g. Lee et al. 2016).

## 2.2. Non-dimensionalisation

Non-dimensionalising the governing equations (2.1)–(2.14) using $\epsilon \hat{U} = \hat{Q}/(\hat{H}^2)$ for the velocity scale, with $\epsilon = \hat{H}/\hat{P}_x$ for the slenderness parameter, $\hat{P} = \hat{\mu}\hat{U}/\hat{H}$ for the pressure scale, $\hat{C} = \hat{K}/\hat{Q}$ for the bulk concentration scale and $\hat{G} = \hat{K}_a\hat{C}/\hat{K}_d$ for the interface concentration scale, we write

$$ x = \frac{\hat{x}}{\hat{P}_x}, \quad \mathbf{x}_{\perp} = \frac{\hat{\mathbf{x}}_{\perp}}{\epsilon \hat{P}_x}, \quad u = \frac{\hat{u}}{\epsilon \hat{U}}, \quad \mathbf{u}_{\perp} = \frac{\hat{\mathbf{u}}_{\perp}}{\hat{U}}, \quad p = \frac{\hat{p}}{\hat{P}}, \quad c = \frac{\hat{c}}{\hat{C}}, \quad \Gamma = \frac{\hat{\Gamma}}{\hat{G}}, \quad (2.15a-g) $$

where $\hat{\mathbf{x}}_{\perp} \equiv (\hat{y}, \hat{z})$ and $\hat{\mathbf{u}}_{\perp} \equiv (\hat{v}, \hat{w})$. Assuming that $\epsilon \ll 1$, this normalisation yields a long-wave theory for steady flow in the streamwise direction and introduces a velocity scaling which captures the rapid cross-channel transport that acts to eliminate cross-channel gradients of surfactant. The transverse flow decays exponentially quickly (Mcnair, Jensen & Landel 2022), but should formally be retained to develop consistent expansions. Here we only consider channels with an order-one aspect ratio (such that $\hat{H} \sim \hat{P}_z$); other asymptotic scalings are left for future work. The longitudinal subdomains
(2.1) become

$$ \mathcal{D}_1 = \{x \in [-\phi_x, \phi_x]\} \times \{y \in [0, 2]\} \times \{z \in [-P_z, P_z]\}, \quad (2.16a) $$

$$ \mathcal{D}_2 = \{x \in [\phi_x, 2 - \phi_x]\} \times \{y \in [0, 2]\} \times \{z \in [-P_z, P_z]\}, \quad (2.16b) $$

where $P_z = \hat{P}_z/\hat{H}$ is the non-dimensional pitch, and interfaces (2.2) are

$$ \mathcal{I} = \{x \in [-\phi_x, \phi_x]\} \times \{y \in [0, 2]\} \times \{z \in [-\phi_z P_z, \phi_z P_z]\}, \quad (2.17a) $$

$$ \mathcal{R} = \{x \in [-\phi_x, \phi_x]\} \times \{y \in [0, 2]\} \times \{z \in [-P_z, -\phi_z P_z] \cup [\phi_z P_z, P_z]\}, \quad (2.17b) $$

$$ \mathcal{S} = \{x \in [\phi_x, 2 - \phi_x]\} \times \{y \in [0, 2]\} \times \{z \in [-P_z, P_z]\}. \quad (2.17c) $$

We substitute the non-dimensionalisation (2.15a-g) into the governing equations (2.3a-c)-(2.14) to acquire rescaled governing equations in terms of the non-dimensional variables given in (2.15a-g) and $\epsilon^2$. In $\mathcal{D}_1$ and $\mathcal{D}_2$,

$$ \left. \begin{array}{l} \epsilon^2 u_x + \nabla_{\perp} \cdot \mathbf{u}_{\perp} = 0, \quad \epsilon^2 \mathbf{u}_{xx} + \nabla_{\perp}^2 \mathbf{u} - \nabla p = \mathbf{0}, \\ Pe^{-1}(\epsilon^2 c_{xx} + \nabla_{\perp}^2 c) - \epsilon^2 u c_x - \mathbf{u}_{\perp} \cdot \nabla_{\perp} c = 0, \end{array} \right\} \quad (2.18a-c) $$

with $Pe = \hat{U}\hat{H}/\hat{D}$ the bulk Péclet number; $\nabla_{\perp} \equiv (\partial_y, \partial_z)$ and $\nabla_{\perp}^2 \equiv \partial_{yy} + \partial_{zz}$ are cross-channel differential operators. Along $\mathcal{I}$,

$$ \left. \begin{array}{l} \mathbf{n} \cdot \nabla u - Ma \Gamma_x = 0, \quad \mathbf{n} \cdot \nabla w - Ma \Gamma_z = 0, \quad v = 0, \quad \mathbf{n} \cdot \nabla c - Da(c - \Gamma) = 0, \\ Pe_I^{-1}(\epsilon^2 \Gamma_{xx} + \Gamma_{zz}) + Bi(c - \Gamma) - \epsilon^2(u\Gamma)_x - (w\Gamma)_z = 0, \end{array} \right\} \quad (2.19a-e) $$

with $Ma = \hat{A}\hat{G}/\hat{\mu}\hat{U}$ the Marangoni number, $Da = \hat{K}_a\hat{H}/\hat{D}$ the Damköhler number, $Pe_I = \hat{H}\hat{U}/\hat{D}_I$ the interfacial Péclet number and $Bi = \hat{K}_d\hat{H}/\hat{U}$ the Biot number. On $\partial\mathcal{I}$,

$$ u\Gamma - Pe_I^{-1}\Gamma_x = 0 \quad \text{at } x = \pm\phi_x, \quad w\Gamma - Pe_I^{-1}\Gamma_z = 0 \quad \text{at } z = \pm\phi_z P_z. \quad (2.20a,b) $$

On $\mathcal{R}$ and $\mathcal{S}$,

$$ u = 0, \quad w = 0, \quad v = 0, \quad c_y = 0. \quad (2.21a-d) $$

Periodicity and continuity between subdomains require that, for $\mathbf{q} = (\mathbf{u}, p_x, c)$,

$$ \mathbf{q}(-\phi_x, y, z) = \mathbf{q}(2 - \phi_x, y, z), \quad \mathbf{q}(x, y, -P_z) = \mathbf{q}(x, y, P_z), \quad (2.22a,b) $$

and

$$ \mathbf{q}(\phi_x^-, y, z) = \mathbf{q}(\phi_x^+, y, z). \quad (2.23) $$

The total liquid flux, total surfactant flux and net flux of surfactant from interface to bulk are given respectively by

$$ \int_{z=-P_z}^{P_z} \int_{y=0}^{2} u \, dy \, dz = 1, \quad (2.24) $$

$$ \int_{z=-P_z}^{P_z} \int_{y=0}^{2} \left( u\epsilon - \frac{c_x}{Pe} \right) dy \, dz + \frac{2Da}{BiPe} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} \left( u\Gamma - \frac{\Gamma_x}{Pe_I} \right) dz = 1, \quad (2.25) $$

$$ \int_{x=-\phi_x}^{\phi_x} \int_{z=-\phi_z P_z}^{\phi_z P_z} (\Gamma - c) \, dz \, dx = 0, \quad (2.26) $$
where bulk and surface fluxes are related to adsorption-desorption at the interface via

$$ \frac{d}{dx} \int_{z=-P_z}^{P_z} \int_{y=0}^{2} \left( u c - \frac{c_x}{Pe} \right) dy dz = \frac{2Da}{\epsilon^2 Pe} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} (\Gamma - c) dz, \quad (2.27a) $$

$$ \frac{d}{dx} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} \left( u \Gamma - \frac{\Gamma_x}{Pe_I} \right) dz = - \frac{Bi}{\epsilon^2} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} (\Gamma - c) dz. \quad (2.27b) $$

The group $Da/(BiPe) = \hat{K}_a/(\hat{K}_d\hat{H})$ in (2.25) appears frequently in this problem. The depletion length $\hat{L}_d = \hat{K}_a/\hat{K}_d$ (Manikantan & Squires 2020) corresponds to the depth into the liquid necessary to balance the adsorption and desorption fluxes between the interface and the bulk, at equilibrium. Hence, the group $Da/(BiPe)$ compares $\hat{L}_d$ to the channel height $\hat{H}$. The drag reduction (2.14) becomes

$$ DR = \frac{\Delta p_R - \Delta p}{\Delta p_R - \Delta p_U}, \quad (2.28) $$

where $\Delta p \equiv \langle p \rangle(-\phi_x) - \langle p \rangle(2 - \phi_x)$ and $\langle \cdot \rangle \equiv \int_{z=-P_z}^{P_z} \int_{y=0}^{2} dy dz/(4P_z)$.

The non-dimensional governing equations (2.18a-c), boundary conditions (2.19a-e)-(2.23) and flux constraints (2.24)-(2.26), define a 3D boundary-value problem for $u, p, c$ and $\Gamma$. The solution depends on nine dimensionless groups ($\epsilon, Pe, Ma, Pe_I, Bi, Da, \phi_x, \phi_z$ and $P_z$) that characterise the geometry, flow, liquid and surfactant. The number of dimensionless groups differs from Landel et al. (2020) (where there are eight), as the authors considered 2D geometries, and Temprano-Coleto et al. (2023) (where there are 10), as we have absorbed the non-dimensional background concentration into $Ma$. Our aim is to construct a reduced model in the limit $\epsilon \ll 1$ to predict the drag reduction (2.28).

### 3. Model

With nine non-dimensional parameters in the problem, we choose distinguished limits to reveal different dominant physical balances. We take $\epsilon \to 0$ with some parameters held fixed and others varying proportionally to $\epsilon^2$. Specifically, in § 3.2, we take $Da = Bi = O(1)$ as $\epsilon \to 0$, allowing the leading-order bulk ($c_0$) and interfacial ($\Gamma_0$) surfactant concentrations to remain in equilibrium, such that $c_0 = \Gamma_0$; we call this the ‘strong exchange limit’. In § 3.3, we treat smaller $Da$ and $Bi$, so that $\Gamma_0$ decouples from $c_0$; we call this the ‘moderate exchange limit’. First, however (in § 3.1), we introduce the resulting 1D surfactant transport equations coupling the bulk and the interface, which depend on five non-dimensional parameters that characterise the strength of bulk and interfacial advection, Marangoni effects, diffusion and bulk-interface exchange and one geometrical parameter.

#### 3.1. Surfactant transport equations

To derive the 1D surfactant transport equations from (2.27), we assume that $Pe = Pe_I = Ma = O(1)$ in the limit $\epsilon \ll 1$, which implies that cross-channel concentration gradients are small, so that $c \approx c_0(x)$ and $\Gamma \approx \Gamma_0(x)$. A detailed derivation is provided in § 3.3; however, in this subsection, we introduce the 1D surfactant transport equations that we solve to generate the results in § 4. We briefly describe their physical meaning to provide the reader with an overview of the different physical processes at play in the problem studied. The 1D surfactant transport equations relate bulk and surface fluxes of surfactant
through adsorption–desorption fluxes at the interface. As shown in § 3.3, they reduce to a simpler coupled nonlinear system of second-order ODEs:

$$ (c_0 - \alpha c_{0x})_x - \frac{\nu}{\epsilon^2}(\Gamma_0 - c_0) = 0 \quad \text{in } \mathcal{D}_1, \qquad (3.1a) $$

$$ (\beta\Gamma_0 - \gamma\Gamma_0\Gamma_{0x} - \delta\Gamma_{0x})_x - \frac{\nu}{\epsilon^2}(c_0 - \Gamma_0) = 0 \quad \text{in } \mathcal{D}_1, \qquad (3.1b) $$

$$ (c_0 - \alpha c_{0x})_x = 0 \quad \text{in } \mathcal{D}_2. \qquad (3.1c) $$

The total flux of surfactant (2.25) becomes

$$ c_0 - \alpha c_{0x} + \beta \Gamma_0 - \gamma \Gamma_0 \Gamma_{0x} - \delta \Gamma_{0x} = 1 \quad \text{in } \mathcal{D}_1, \qquad (3.2a) $$

$$ c_0 - \alpha c_{0x} = 1 \quad \text{in } \mathcal{D}_2. \qquad (3.2b) $$

Continuity of bulk surfactant (2.23) and bulk surfactant flux (2.25) between unit cells, and no flux of interfacial surfactant (2.20a,b) through contact lines, are together given by

$$ \left. \begin{array}{l} c_0(\phi_x^-) = c_0(\phi_x^+), \quad c_0(-\phi_x) = c_0(2-\phi_x), \\ c_0 - \alpha c_{0x} = 1, \quad \beta\Gamma_0 - \gamma\Gamma_0\Gamma_{0x} - \delta\Gamma_{0x} = 0 \quad \text{at } x = \pm\phi_x. \end{array} \right\} (3.3a-d) $$

Equations (3.1)–(3.3a–d) are the steady leading-order surfactant transport equations, flux constraints and boundary conditions in $\mathcal{D}_1$ and $\mathcal{D}_2$, accounting for bulk advection by the pressure gradient ($c_{0x}$ in (3.1)), surface advection by the pressure gradient ($\beta\Gamma_{0x}$), Marangoni advection ($\gamma(\Gamma_0\Gamma_{0x})_x$), bulk diffusion ($\alpha c_{0xx}$), surface diffusion ($\delta\Gamma_{0xx}$) and bulk-surface exchange ($\pm\nu(\Gamma_0 - c_0)/\epsilon^2$). In § 3.2, we derive and discuss further the non-dimensional parameters $\alpha, \beta, \gamma$ and $\delta$ associated with the strengths of the above physical processes. In § 3.2, we focus on the limit where bulk-surface exchange of surfactant is strong relative to advection and diffusion, such that $\nu/\epsilon^2 \gg O(1, \alpha, \beta, \delta)$ in (3.1)–(3.3a–d) and, therefore, the $c_0$ and $\Gamma_0$ fields are in equilibrium. In § 3.3, we study the limit where bulk-surface exchange of surfactant is comparable to advection and diffusion, such that $\nu/\epsilon^2 \sim O(1, \alpha, \beta, \delta)$ in (3.1)–(3.3a–d) and, therefore, the $c_0$ and $\Gamma_0$ fields are distinct.

## 3.2. Strong exchange

### 3.2.1. Strong cross-channel diffusion

We begin by assuming that $Pe = Pe_I = Bi = Da = Ma = O(1)$ in the limit $\epsilon \ll 1$. We call this the ‘strong cross-channel diffusion limit’. In this subsection, we derive the $O(1)$ coefficients $\alpha, \beta, \gamma$ and $\delta$, such that $\nu/\epsilon^2 \gg O(1)$ for $\epsilon \ll 1$ and $\nu = O(1)$. We substitute the expansions

$$ (u, v, w, p, c, \Gamma) = (u_0, v_0, w_0, p_0, c_0, \Gamma_0) + \epsilon^2(u_1, v_1, w_1, p_1, c_1, \Gamma_1) + \dots, \qquad (3.4) $$

in the governing equations (2.18a–c)–(2.28) and take the $O(1)$ approximation. In domains $\mathcal{D}_1$ and $\mathcal{D}_2$,

$$ \nabla_{\perp} \cdot \mathbf{u}_{\perp 0} = 0, \quad \nabla_{\perp}^2 \mathbf{u}_0 - \nabla p_0 = 0, \quad Pe^{-1} \nabla_{\perp}^2 c_0 - \mathbf{u}_{\perp 0} \cdot \nabla_{\perp} c_0 = 0. \qquad (3.5a-c) $$

Along the interface $\mathcal{I}$,

$$ \left. \begin{array}{l} \mathbf{n} \cdot \nabla u_0 - Ma\Gamma_{0x} = 0, \quad \mathbf{n} \cdot \nabla w_0 - Ma\Gamma_{0z} = 0, \quad v_0 = 0, \\ \mathbf{n} \cdot \nabla c_0 - Da(c_0 - \Gamma_0) = 0, \quad Pe_I^{-1}\Gamma_{0zz} + Bi(c_0 - \Gamma_0) - (w_0\Gamma_0)_z = 0. \end{array} \right\} (3.6a-e) $$
and on the interface contour $\partial\mathcal{I}$,

$$w_0\Gamma_0 - Pe_I^{-1}\Gamma_{0z} = 0 \quad \text{at } z = \pm\phi_z P_z. \qquad (3.7)$$

There are inner regions near $x = \pm\phi_x$, within which $u_0\Gamma_0 - Pe_I^{-1}\Gamma_{0x} = 0$ is imposed on $\partial\mathcal{I}$, and the boundary conditions $u_0(\phi_x^-, y, z) = u_0(\phi_x^+, y, z)$ and $u_0(-\phi_x, y, z) = u_0(2 - \phi_x, y, z)$ are imposed between $D_1$ and $D_2$. Within the inner regions, the flow and surfactant field are governed by the coupled 3D Stokes and surfactant transport equations. In the present long-wave theory it is sufficient to impose continuity of $c_0$ between $D_1$ and $D_2$, such that

$$c_0(\phi_x^-, y, z) = c_0(\phi_x^+, y, z), \quad c_0(-\phi_x, y, z) = c_0(2 - \phi_x, y, z), \qquad (3.8a,b)$$

and continuity of volume and surfactant flux between $D_1$ and $D_2$, such that

$$\int_{z=-P_z}^{P_z} \int_{y=0}^{2} u_0 \, dy \, dz = 1, \qquad (3.9)$$

$$\int_{z=-P_z}^{P_z} \int_{y=0}^{2} \left( u_0 c_0 - \frac{c_{0x}}{Pe} \right) dy \, dz + \frac{2Da}{BiPe} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} \left( u_0 \Gamma_0 - \frac{\Gamma_{0x}}{Pe_I} \right) dz = 1, \quad (3.10)$$

for all $-\phi_x \le x \le 2 - \phi_x$. Along the ridge $\mathcal{R}$ and solid $\mathcal{S}$,

$$u_0 = 0, \quad v_0 = 0, \quad w_0 = 0, \quad c_{0y} = 0. \qquad (3.11a-d)$$

Transverse periodicity can be rewritten as the symmetry conditions

$$u_{0z} = 0, \quad v_{0z} = 0, \quad w_0 = 0, \quad c_{0z} = 0 \quad \text{at } z = \pm P_z. \qquad (3.12a-d)$$

The drag reduction becomes

$$DR_0 = \frac{\Delta p_R - \Delta p_0}{\Delta p_R - \Delta p_U}. \qquad (3.13)$$

From (3.5a-c)-(3.12a-d), the leading-order solution simplifies to $\Gamma_0 = \Gamma_0(x)$, $c_0 = c_0(x)$ where $c_0 = \Gamma_0$, $p_0 = p_0(x)$ and $v_0 = w_0 = 0$. The leading-order velocities in the cross-section, $v_0$ and $w_0$, vanish due to transverse surfactant gradients decaying exponentially fast in a time-dependent setting, as discussed by Mcnair et al. (2022). Therefore, in our steady problem, the concentration field does not vary in the transverse direction and there are no concentration gradients to generate velocities in the cross-plane. The streamwise flow is driven by $p_{0x}$, from (3.5b), and $\Gamma_{0x}$, from (3.6a). Using linear superposition, we can write

$$u_0 = \tilde{U} p_{0x} + Ma \tilde{U} \Gamma_{0x} \quad \text{in } D_1 \quad \text{and} \quad u_0 = \tilde{U} p_{0x} \quad \text{in } D_2, \qquad (3.14a,b)$$

where $\tilde{U}(y, z)$ and $\tilde{U}(y, z)$ are velocity contributions related to the bulk pressure and surfactant gradient, respectively, in $D_1$, and $\tilde{U}(y)$ is the velocity contribution related to the bulk pressure gradient in $D_2$. Hence, we must solve the following boundary-value
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Contour plot of $\tilde{U}$ defined by (3.15a-g), the leading-order contribution to the streamwise flow due to the pressure gradient, $p_{0x}$. (b) Contour plot of $\tilde{U}$ defined by (3.16a-g), the leading-order contribution to the streamwise flow due to the surfactant gradient, $\Gamma_{0x}$. As $p_{0x} < 0$, $\tilde{U}$ contributes positively to the leading-order streamwise velocity component $u_0$, whilst $\tilde{U}$ contributes negatively, because $\Gamma_{0x} > 0$, following (3.14a,b). The thick black lines represent the solid regions of the SHS with the transverse gas fraction $\phi_z = 0.5$ and width $P_z = 1$.</figcaption></figure>

problems: streamwise flow driven by a pressure gradient over $D_1$,

$$ \left. \begin{array}{l} \nabla_{\perp}^2 \tilde{U} = 1, \quad \text{subject to} \quad \tilde{U}_y(0, z_s) = 0, \quad \tilde{U}(0, z_{ns}) = 0, \quad \tilde{U}_y(2, z_s) = 0, \\ \qquad \tilde{U}(2, z_{ns}) = 0, \quad \tilde{U}_z(y, -P_z) = 0, \quad \tilde{U}_z(y, P_z) = 0; \end{array} \right\} (3.15a-g) $$

streamwise flow driven by a surfactant-induced Marangoni shear stress over $D_1$,

$$ \left. \begin{array}{l} \nabla_{\perp}^2 \tilde{U} = 0, \quad \text{subject to} \quad \tilde{U}_y(0, z_s) = 1, \quad \tilde{U}(0, z_{ns}) = 0, \quad \tilde{U}_y(2, z_s) = -1, \\ \qquad \tilde{U}(2, z_{ns}) = 0, \quad \tilde{U}_z(y, -P_z) = 0, \quad \tilde{U}_z(y, P_z) = 0; \end{array} \right\} (3.16a-g) $$

and streamwise flow driven by a pressure gradient over $D_2$,

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

with $z_s \equiv \{z \in [-\phi_z P_z, \phi_z P_z]\}$ and $z_{ns} \equiv \{z \in [-P_z, -\phi_z P_z]\} \cup \{z \in [\phi_z P_z, P_z]\}$.

Numerical solutions to (3.15a-g), (3.16a-g) can be seen in figure 2 for $\phi_z = 0.5$ and $P_z = 1$, computed using the method outlined in Appendix A. Equation (3.17a-c) can be integrated to give $\tilde{U} = y(y-2)/2$ and $\tilde{Q} = \int_{z=-P_z}^{P_z} \int_{y=0}^{2} \tilde{U} \, dy \, dz = -4P_z/3$. We define the volume and surface fluxes

$$ \tilde{Q} = \int_{z=-P_z}^{P_z} \int_{y=0}^{2} \tilde{U} \, dy \, dz, \quad \tilde{Q} = \int_{z=-P_z}^{P_z} \int_{y=0}^{2} \tilde{U} \, dy \, dz, \qquad (3.18a,b) $$

$$ q = \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} u_0 \, dz, \quad \tilde{q} = \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} \tilde{U} \, dz, \quad \tilde{q} = \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} \tilde{U} \, dz. \qquad (3.19a-c) $$

Bulk $(\tilde{Q}, \tilde{Q})$ and surface fluxes $(\tilde{q}, \tilde{q})$ are plotted as functions of $\phi_z$ and $P_z$ in figure 3.

Next, we substitute the expansions (3.4) into the governing equations (2.18a-c)-(2.22a,b) to evaluate the $O(\epsilon^2)$ equations. The $O(\epsilon^2)$ system is required to calculate the first-order cross-channel flow field $(v_1, w_1)$, pressure field $p_1$, bulk $c_1$ and surface
<figure><img src="image_6.png" /><figcaption>Figure 3. Plot of (a) $\tilde{Q}$ defined in (3.18a) and (b) $\tilde{q}$ defined in (3.19b), the contributions to the bulk and surface flux due to the pressure gradient $p_{0x} < 0$, for varying $\phi_z$ and $P_z$. (c) Plot of $\tilde{Q}$ defined in (3.18b) and (d) plot of $\tilde{q}$ defined in (3.19c), the contributions to the bulk and surface flux due to the surfactant gradient $\Gamma_{0x} > 0$, for varying $\phi_z$ and $P_z$. As $p_{0x} < 0$, $\tilde{Q}$ and $\tilde{q}$ contribute positively to the leading-order velocity flux, whereas $\tilde{Q}$ and $\tilde{q}$ contribute negatively, because $\Gamma_{0x} > 0$, following (3.25a).</figcaption></figure>

concentration field $\Gamma_1$ driven by streamwise gradients of the leading-order quantities ($u_0 = (u_0(x, y, z), 0, 0)$, $c_0(x)$ and $\Gamma_0(x)$), and to close the leading-order problem. In $D_1$ and $D_2$,

$$ \nabla_{\perp} \cdot \mathbf{u}_{\perp 1} = -u_{0x}, \quad \nabla_{\perp}^2 \mathbf{u}_{\perp 1} - \nabla_{\perp} p_1 = \mathbf{0}, \quad Pe^{-1} \nabla_{\perp}^2 c_1 = -Pe^{-1} c_{0xx} + u_0 c_{0x}. \qquad (3.20a-c) $$

Along $\mathcal{I}$,

$$ \left. \begin{array}{l} \mathbf{n} \cdot \nabla w_1 - Ma\Gamma_{1z} = 0, \quad v_1 = 0, \quad \mathbf{n} \cdot \nabla c_1 - Da(c_1 - \Gamma_1) = 0, \\ Pe_I^{-1}\Gamma_{1zz} + Bi(c_1 - \Gamma_1) - \Gamma_0 w_{1z} = -Pe_I^{-1}\Gamma_{0xx} + (u_0\Gamma_0)_x, \end{array} \right\} (3.21a-d) $$

and on $\partial\mathcal{I}$,

$$ w_1\Gamma_0 - Pe_I^{-1}\Gamma_{1z} = 0 \quad \text{at } z = \pm\phi_z P_z. \qquad (3.22) $$

Along $\mathcal{R}$ and $\mathcal{S}$,

$$ v_1 = 0, \quad w_1 = 0, \quad c_{1y} = 0. \qquad (3.23a-c) $$

Transverse periodicity can be rewritten as the symmetry conditions

$$ v_{1z} = 0, \quad w_1 = 0, \quad c_{1z} = 0 \quad \text{at } z = \pm P_z. \qquad (3.24a-c) $$

Solvability conditions imposed on (3.20a-c)-(3.24a-c) constrain $u_0$, $c_0$ and $\Gamma_0$ at $O(1)$. The forcing on the right-hand sides of (3.20a-c), (3.21a-d) must be orthogonal to each vector in the null space of the linear operator that is adjoint to the left-hand sides of (3.20a-c) and (3.21a-d). Conveniently, these conditions are provided by the conservation arguments that result in the leading-order velocity and surfactant flux conditions (3.9), (3.10). Substituting the leading-order streamwise velocity (3.14a,b) into the leading-order velocity flux condition (3.9), we obtain

$$ \tilde{Q}p_{0x} + Ma\tilde{Q}\Gamma_{0x} = 1, \quad q = \tilde{q}p_{0x} + Ma\tilde{q}\Gamma_{0x} \quad \text{in } D_1, \quad \tilde{Q}p_{0x} = 1 \quad \text{in } D_2. \qquad (3.25a-c) $$

Substituting (3.25a-c) into the leading-order surfactant flux condition (3.10) and remembering that $c_0(x) = \Gamma_0(x)$, we obtain the strong cross-channel diffusion problem

$$ (\beta + 1) c_0 - \gamma c_0 c_{0x} - (\alpha + \delta) c_{0x} = 1 \quad \text{in } D_1, \qquad (3.26a) $$

$$ c_0 - \alpha c_{0x} = 1 \quad \text{in } D_2, \qquad (3.26b) $$
subject to continuity of bulk surfactant (3.3a,b). In (3.26), we have introduced

$$ \alpha = \frac{4P_z}{Pe} = \frac{4\hat{D}\hat{P}_z}{\hat{H}^2\hat{U}} \quad (\text{bulk diffusion}), \quad (3.27a) $$

$$ \beta = \frac{2Da\tilde{q}}{BiPe\tilde{Q}} = \frac{2\hat{K}_a\tilde{q}}{\hat{H}\hat{K}_d\tilde{Q}} \quad (\text{partition coefficient}), \quad (3.27b) $$

$$ \gamma = \frac{2MaDa}{BiPe} \left( \frac{\tilde{q}\tilde{Q}}{\tilde{Q}} - \bar{q} \right) = \frac{2\hat{K}_a\hat{A}\hat{G}}{\hat{H}\hat{K}_d\hat{\mu}\hat{U}} \left( \frac{\tilde{q}\tilde{Q}}{\tilde{Q}} - \bar{q} \right) \quad (\text{surfactant strength}), \quad (3.27c) $$

$$ \delta = \frac{4\phi_z P_z Da}{BiPePe_I} = \frac{4\hat{D}_I \hat{K}_a \phi_z \hat{P}_z}{\hat{H}^3 \hat{K}_d \hat{U}} \quad (\text{surface diffusion}). \quad (3.27d) $$

The bulk diffusion coefficient $\alpha > 0$ is a rescaled (by $4\hat{P}_z/\hat{H}$) inverse Péclet number $1/Pe = \hat{D}/(\hat{H}\hat{U})$, nominally assumed $O(1)$. The coefficient $\alpha$ characterises the ratio of streamwise diffusion to streamwise advection in the bulk for a given normalised channel width, assumed $O(1)$. The partition coefficient $\beta > 0$ is the rescaled (by $2\tilde{q}/\tilde{Q}$) normalised surfactant depletion depth $L_d = \hat{K}_a/(\hat{H}\hat{K}_d)$. The coefficient $\beta$ in (3.26a) determines the portion of the surfactant flux transported by the pressure-driven flow along the interface in comparison with the surfactant flux transported by the pressure-driven flow in the bulk. For $\beta \gg 1$ ($\beta \ll 1$), the surfactant advection flux is strongest at the interface (bulk). The surfactant strength coefficient $\gamma > 0$ is the rescaled (by $2(\tilde{q}\tilde{Q}/\tilde{Q} - \bar{q})$) product of $L_d$ and the Marangoni number $Ma = \hat{A}\hat{G}/(\hat{\mu}\hat{U})$, the ratio of surface tension changes due to interfacial surfactant to viscous forces. The coefficient $\gamma$ characterises the nonlinear effect of streamwise surfactant-induced Marangoni stresses on the streamwise surfactant flux at the interface and in the bulk. The surface diffusion parameter $\delta > 0$ is a rescaled (by $4\phi_z P_z$) product of $L_d$ and the inverse surface Péclet number $1/Pe_I = \hat{D}_I/(\hat{H}\hat{U})$, nominally assumed $O(1)$. The coefficient $\delta$ measures the ratio of interfacial streamwise surface diffusion to streamwise advection for a given normalised interface width.

In summary, the leading-order problem (3.5a-c)-(3.12a-d) has been reduced to the five-parameter problem (3.26) subject to (3.3a,b), a special case (where $\nu/\epsilon^2 \gg O(1, \alpha, \beta, \delta)$) of the six-parameter problem (3.1)-(3.3a-d), requiring the solution of a first-order ODE to determine $c_0$ and $DR_0$ (explained in § 3.4). Then, the leading-order streamwise velocity (3.14a,b) and first-order problem (3.20a-c)-(3.24a-c) can be used to construct the 3D flow field.

### 3.2.2. Moderate cross-channel diffusion

In general, when cross-channel diffusion is weak and cross-channel gradients in the concentration field are comparable to the cross-channel average, a numerical technique must be employed to evaluate $c$. We call this the ‘generalised limit’. In Appendix B.1, we identify the regions of parameter space where cross-channel gradients first emerge by perturbing around the cross-channel-averaged concentration. Writing $c_0 = \langle c_0 \rangle(x) + c'_0(x, y, z)$, we assume that $\|c'_0\| = \sqrt{\langle c_0^2 \rangle} \ll \langle c_0 \rangle$ for all $x$, seeking to characterise the leading-order effects of shear dispersion in this problem (Taylor 1953). We expect shear-dispersion effects to arise in surfactant-contaminated SHS channels because the pressure-driven flow advects surfactant at different speeds across the channel. The advected surfactant is then mixed by cross-channel diffusion, which, in turn, results
in longitudinal dispersion of the surfactant distribution. We call this the ‘moderate cross-channel diffusion limit’. The validity of this assumption is evaluated a posteriori in Appendix B.2.

To develop this limit, we assume that $1/Pe = 1/Pe_I = Bi = O(\epsilon^2)$ and $Da = Ma = O(1)$ for $\epsilon \ll 1$, substituting the asymptotic expansion (3.4) into the governing equations (2.18a-c)-(2.27) (as in § 3.2.1). We then use the first-order system to eliminate $c'_0$ and the leading-order surfactant flux constraint to derive the moderate cross-channel diffusion problem, which is effectively a modification to the strong cross-channel diffusion problem (3.26) accounting for shear dispersion,

$$ (\beta + 1) c_0 - \gamma c_0 c_{0x} - \frac{\epsilon^2}{\alpha} \left( s_1 c_{0x} + s_2 \frac{\gamma}{\beta} c_{0x}^2 + s_3 \frac{\gamma^2}{\beta^2} c_{0x}^3 \right) = 1 \quad \text{in } \mathcal{D}_1, \quad (3.28a) $$

$$ c_0 - \frac{\epsilon^2}{\alpha} (s_4 c_{0x}) = 1 \quad \text{in } \mathcal{D}_2, \quad (3.28b) $$

subject to continuity of bulk surfactant (3.3a,b). We have defined the $O(1)$ coefficients

$$ \left. \begin{aligned} s_1 &= -\frac{16P_z^2 \langle \tilde{U}\tilde{C} \rangle}{\tilde{Q}^2}, & s_2 &= -\frac{16P_z^2 \tilde{q}(\tilde{Q}(\langle \tilde{U}\tilde{C} \rangle + \langle \tilde{U}\tilde{C} \rangle) - 2\tilde{Q}\langle \tilde{U}\tilde{C} \rangle)}{\tilde{Q}^2(\tilde{q}\tilde{Q} - \tilde{q}\tilde{Q})}, \\ s_3 &= -\frac{16P_z^2 \tilde{q}^2 (\tilde{Q}^2 \langle \tilde{U}\tilde{C} \rangle - \tilde{Q}\tilde{Q}(\langle \tilde{U}\tilde{C} \rangle + \langle \tilde{U}\tilde{C} \rangle) + \tilde{Q}^2 \langle \tilde{U}\tilde{C} \rangle)}{\tilde{Q}^2(\tilde{q}\tilde{Q} - \tilde{q}\tilde{Q})^2}, & s_4 &= -\frac{16P_z^2 \langle \tilde{U}\tilde{C} \rangle}{\tilde{Q}^2}, \end{aligned} \right\} \quad (3.29a-d) $$

where $\tilde{C}$, $\tilde{C}$ and $\tilde{C}$ are solutions of the boundary-value problems (B12a-g)-(B14a-c), respectively (see Appendix B.1), which depend on geometrical factors and $Da$. Here, $(\tilde{U}, \tilde{U})$ are given by (3.15a-g), (3.16a-g), $(\tilde{Q}, \tilde{Q})$ by (3.18a,b), and $(\tilde{q}, \tilde{q})$ by (3.19a-c). The coefficients $s_1 > 0$, $s_2$ and $s_3$ are plotted for different $\phi_z$, for $P_z = 1$ and as a function of $Da$ in figure 11 (Appendix B.1). The coefficient $s_4 = 3/35$, because $\tilde{U} = y(y-2)/2$, $\tilde{C} = 1/5 - y^3/6 + y^4/24$ and $\tilde{Q} = -4P_z/3$ in $\mathcal{D}_2$. Comparing the strong (3.26) and moderate cross-channel diffusion problems (3.28), streamwise diffusion terms proportional to $\alpha$ and $\delta$ in (3.26) have been replaced by shear dispersion terms proportional to $\epsilon^2/\alpha$ in (3.28), some of which are nonlinear due to Marangoni effects.

### 3.2.3. Composite equation

The ODEs in the strong cross-channel diffusion limit (3.26) can be combined with the ODEs in the moderate cross-channel diffusion limit (3.28) to construct the composite approximation

$$ (\beta + 1) c_0 - \gamma c_0 c_{0x} - (\alpha + \delta) c_{0x} - \frac{\epsilon^2}{\alpha} \left( s_1 c_{0x} + s_2 \frac{\gamma}{\beta} c_{0x}^2 + s_3 \frac{\gamma^2}{\beta^2} c_{0x}^3 \right) = 1 \quad \text{in } \mathcal{D}_1, \quad (3.30a) $$

$$ c_0 - \alpha c_{0x} - \frac{\epsilon^2}{\alpha} (s_4 c_{0x}) = 1 \quad \text{in } \mathcal{D}_2, \quad (3.30b) $$

subject to continuity of bulk surfactant (3.3a,b). The dependence of $\alpha$, $\beta$, $\gamma$ and $\delta$ on $Ma$, $Pe$, $Pe_I$, $Bi$, $Da$, $\phi_z$ and $P_z$ is given in (3.27). The dependence of $s_1$, $s_2$, $s_3$ and $s_4$ on $Da$, $\phi_z$ and $P_z$ is given in (3.29a-d) and discussed in Appendix B.1. Together, the solution to (3.30) depends on nine dimensionless parameters. The numerical
method used to solve (3.3a), (3.3b), (3.30) is outlined in Appendix A.1. Incorporation of shear dispersion effects in (3.30), through perturbation of the cross-channel-averaged concentration ($c_0 = \langle c_0 \rangle(x) + c'_0(x, y, z)$), enables the strong-exchange model to explore regimes where cross-channel concentration gradients first become significant ($\|c'_0\| = O(\langle c_0 \rangle)$ for some $x$), enabling us to estimate approximately the boundary of validity of the 1D model (3.1)–(3.3a–d) (Appendix B.2). When bulk diffusion is strong enough for shear dispersion to become negligible, (3.30) reduces to (3.26), or (3.2) with $c_0 = \Gamma_0$.

### 3.3. Moderate exchange

In § 3.2, we assumed that $Da = O(1)$ in order to have strong exchange at the interface, such that bulk and interfacial surfactants were in equilibrium at leading order, i.e. $c_0 = \Gamma_0$. To study the case where $c_0$ and $\Gamma_0$ are not in equilibrium, we rescale $Bi = \epsilon^2 \mathcal{B}$ and $Da = \epsilon^2 \mathcal{D}$ with $\mathcal{B} = \mathcal{D} = O(1)$, whilst retaining $Pe = Pe_I = Ma = O(1)$. Substituting the expansion (3.4) into the governing equations (2.18a–c)–(2.27), then the leading-order strong cross-channel diffusion equations (3.5a–c), (3.6a–c), (3.7)–(3.9) and (3.11a–d)–(3.12a–d) are recovered. Along the interface $\mathcal{I}$

$$c_{0y} = 0, \quad Pe_I^{-1} \Gamma_{0zz} - (w_0 \Gamma_0)_z = 0. \qquad (3.31a,b)$$

The total flux of surfactant (2.25) is given by

$$\int_{z=-P_z}^{P_z} \int_{y=0}^{2} \left( u_0 c_0 - \frac{c_{0x}}{Pe} \right) dy dz + \frac{2\mathcal{D}}{\mathcal{B}Pe} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} \left( u_0 \Gamma_0 - \frac{\Gamma_{0x}}{Pe_I} \right) dz = 1. \quad (3.32)$$

The bulk and surface fluxes (2.27) are related to adsorption–desorption fluxes at the interface via

$$\frac{d}{dx} \left( \int_{z=-P_z}^{P_z} \int_{y=0}^{2} \left( u_0 c_0 - \frac{c_{0x}}{Pe} \right) dy dz \right) = \frac{2\mathcal{D}}{Pe} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} (\Gamma_0 - c_0) dz, \quad (3.33a)$$

$$\frac{d}{dx} \left( \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} \left( u_0 \Gamma_0 - \frac{\Gamma_{0x}}{Pe_I} \right) dz \right) = -\mathcal{B} \int_{z=-\phi_z P_z, y=0}^{\phi_z P_z} (\Gamma_0 - c_0) dz. \quad (3.33b)$$

The cross-channel velocity field again decays rapidly in time to $v_0 = w_0 = 0$. However, the bulk concentration $c_0 = c_0(x)$ and surface concentration $\Gamma_0 = \Gamma_0(x)$ are no longer equal at leading-order. The streamwise velocity field is given by $u_0 = \tilde{U}p_{0x} + Ma\tilde{U}\Gamma_{0x}$ in $D_1$ and $u_0 = \tilde{U}p_{0x}$ in $D_2$, where $\tilde{U}$, $\tilde{U}$ and $\tilde{U}$ are given by (3.15a–g)–(3.17a–c).

At $O(\epsilon^2)$, the first-order strong cross-channel diffusion equations (3.20a–c), (3.21a,b) and (3.22)–(3.24a–c) are recovered. On $\mathcal{I}$,

$$c_{1y} = \mathcal{D}(c_0 - \Gamma_0), \quad Pe_I^{-1} \Gamma_{1zz} - \Gamma_0 w_{1z} = -Pe_I^{-1} \Gamma_{0xx} + (u_0 \Gamma_0)_x - \mathcal{B}(c_0 - \Gamma_0). \qquad (3.34a,b)$$

This system is required to evaluate the first-order cross-channel flow, ($v_1, w_1$), driven by streamwise gradients of $u_0, c_0$ and $\Gamma_0$, and to close the leading-order problem.

Substituting the streamwise velocity (3.14a,b) into the velocity flux condition (3.9), we recover $\tilde{Q}p_{0x} + Ma\tilde{Q}\Gamma_{0x} = 1$ in $D_1$ and $\tilde{Q}p_{0x} = 1$ in $D_2$, where $\tilde{Q}$ is given by (3.18a), $\tilde{Q}$ is given by (3.18b) and $\tilde{Q} = -4P_z/3$. Substituting the streamwise velocity (3.14a,b) into (3.33) we derive the steady leading-order surfactant transport equations (3.1). Substituting (3.14a,b) into the surfactant transport equations (3.32) (or adding (3.1a) and (3.1b) and
integrating) we recover integral constraints on the surfactant transport equations (3.1), which are given by surfactant flux constraints (3.2). The coefficients $\alpha, \beta, \gamma$ and $\delta$ are defined in (3.27) and

$$ \nu = \frac{4\phi_z P_z Da}{Pe} = \frac{4\hat{K}_a \phi_z \hat{P}_z}{\hat{H}\hat{U}} \quad (\text{surface exchange}). \quad (3.35) $$

The surface exchange coefficient $\nu$ is a rescaled (by $4\phi_z \hat{P}_z$) product of the inverse bulk Péclet and Damköhler numbers, comparing the rate of adsorption of surfactant contained in the layer of fluid $\hat{H}$ with the rate of advection. At the domain boundaries, we impose continuity of bulk surfactant (3.3a,b), continuity of bulk surfactant flux between $D_1$ and $D_2$ (3.3c) and no flux of interfacial surfactant (3.3d). For $\nu/\epsilon^2 \gg O(1, \alpha, \beta, \delta)$, then $c_0 = \Gamma_0$ at leading order, recovering the transport equations (3.26) in the strong cross-channel diffusion and strong exchange limit. The numerical method used to solve the 1D surfactant transport equations (3.1)–(3.3a–d) is outlined in Appendix A.1.

### 3.4. Drag reduction

As mentioned in §2, we study the transition from an immobilised interface where $DR_0 = 0$, to a shear-free interface where $DR_0 = 1$. When the interface is shear-free, $\tilde{Q}p_{0x} = 1$ in $D_1$ and $\tilde{Q}p_{0x} = 1$ in $D_2$, such that integrating $-p_{0x}$ across the period gives

$$ \Delta p_U = - \int_{x=-\phi_x}^{\phi_x} p_{0x} \, dx - \int_{x=\phi_x}^{2-\phi_x} p_{0x} \, dx = - \frac{2\phi_x}{\tilde{Q}} - \frac{2(1-\phi_x)}{\tilde{Q}}. \quad (3.36) $$

In the limit where the interface is immobilised, $\tilde{Q}p_{0x} = 1$ in $D_1$ and $D_2$, giving $\Delta p_R = -2/\tilde{Q}$. Note that when $q=0$, the surface velocity flux condition (3.25b) implies that $\tilde{q}p_{0x} + Ma\tilde{q}\Gamma_{0x} = 0$. We can substitute this into the bulk velocity flux condition (3.25a), to find $p_{0x}(\tilde{Q}\tilde{q} - \tilde{Q}\tilde{q}) = \tilde{q}$ in $D_1$. However, we also have $p_{0x} = 1/\tilde{Q}$ in $D_1$ for $q=0$. Combining these two expressions for $p_{0x}$ in $D_1$, gives $\tilde{Q}\tilde{q} - \tilde{Q}\tilde{q} = \tilde{q}\tilde{Q}$ and $\tilde{q}(\tilde{Q} - \tilde{Q}) = \tilde{Q}\tilde{q}$. It follows that $\tilde{Q}\tilde{Q}\tilde{q}/((\tilde{q}\tilde{Q} - \tilde{q}\tilde{Q})(\tilde{Q} - \tilde{Q})) = 1$, which is valid for any $q$. Between these no-slip and shear-free limits, $\tilde{Q}p_{0x} + Ma\tilde{Q}\Gamma_{0x} = 1$ in $D_1$ and $\tilde{Q}p_{0x} = 1$ in $D_2$, such that integrating $-p_{0x}$ across the period gives

$$ \Delta p_0 = - \int_{x=-\phi_x}^{\phi_x} p_{0x} \, dx - \int_{x=\phi_x}^{2-\phi_x} p_{0x} \, dx = - \frac{2\phi_x}{\tilde{Q}} - \frac{2(1-\phi_x)}{\tilde{Q}} + \frac{Ma\tilde{Q}\Delta\Gamma_0}{\tilde{Q}}, \quad (3.37) $$

where we have defined $\Delta\Gamma_0 = \Gamma_0(\phi_x) - \Gamma_0(-\phi_x) > 0$ ($\Delta\Gamma_0 = \Delta c_0$ in the strong exchange limit). Using the definition of $(\beta, \gamma)$ in (3.27b), (3.27c) and substituting $\Delta p_U$, $\Delta p_R$ and $\Delta p_0$ into the leading-order drag reduction (3.13), we have

$$ DR_0 = 1 - \frac{\gamma \Delta \Gamma_0}{2 \phi_x \beta}. \quad (3.38) $$

In (3.38), the first term gives the drag reduction when the interface is shear-free and the second term measures the effect of surfactant. In order to derive the expression in (3.38), we used the fact that $\tilde{Q}\tilde{Q}\tilde{q}/((\tilde{q}\tilde{Q} - \tilde{q}\tilde{Q})(\tilde{Q} - \tilde{Q})) = 1$.

Drag reduction is one possible measure of the performance of surfactant-contaminated SHSs. An alternative measure, commonly used for laminar flows over SHSs, is the effective slip length $\lambda_e$ (Landel et al. 2020; Temprano-Coleto et al. 2023), defined as the uniform slip length applied to the top and bottom boundaries of an equivalent channel of
the same height as the SHS channel and which has the same flow rate as the SHS channel under the same pressure gradient (Lauga & Stone 2003). To evaluate $\lambda_e$, we integrate the leading-order streamwise momentum equation (3.5b) for an equivalent channel with the mixed boundary conditions (3.6a–e) and (3.11a–d) replaced by $\lambda_e u_{0y} - u_0 = 0$ on $\mathcal{I}$, $\mathcal{R}$ and $\mathcal{S}$. We obtain $u_0 = \tilde{U} \Delta p_{0x}$ where $\tilde{U} = y(y-2)/2 - \lambda_e$ and $p_{0x}$ is the same pressure gradient as in the SHS channel. The flux is $\tilde{Q} = \int_{z=-P_z}^{P_z} \int_{y=0}^{2} \tilde{U} p_{0x} dy dz = (\tilde{Q} - 2P_z \lambda_e) p_{0x}$, or by integrating over one period $\tilde{Q} = (\tilde{Q} - 2P_z \lambda_e) \Delta p_0 / 2$. Equating the flux of the equivalent channel with the flux of the SHS channel, $\tilde{Q} p_{0x} + Ma \tilde{Q} \Gamma_{0x} = 1$, we find

$$ \lambda_e = \frac{DR_0(\Delta p_R - \Delta p_U)}{P_z \Delta p_R (\Delta p_U DR_0 + \Delta p_R (1 - DR_0))}, \quad (3.39) $$

which can be used to convert results from $DR_0$ to $\lambda_e$.

## 4. Results

In §4.1, we investigate the leading-order drag reduction ($DR_0$), bulk surfactant concentration ($c_0$), interfacial surfactant concentration ($\Gamma_0$) and flow field in the strong-exchange problem, (3.30), (3.3a), (3.3b), by varying the bulk diffusion ($\alpha$), partition coefficient ($\beta$), surfactant strength ($\gamma$) and surface diffusion ($\delta$). In §4.2, we address the moderate-exchange problem, (3.1)–(3.3a–d), by varying $\alpha$, $\beta$, $\gamma$, $\delta$ and the exchange strength ($\nu$). In both strong- and moderate-exchange problems, we identify three primary areas of the parameter space. The Marangoni-dominated ($M$) region, where the interfacial surfactant gradient is sufficient to immobilise the liquid-gas interface (low drag reduction); the advection-dominated region ($A$), where interfacial surfactant has been swept to the downstream stagnation point and the liquid-gas interface is mostly shear free (high drag reduction); and the diffusion-dominated ($D$) region, where the interfacial surfactant gradient has been attenuated by diffusion and the liquid-gas interface is mostly shear free (high drag reduction). Throughout §4, the gas fraction and transverse period width are maintained at $\phi_x = \phi_z = 0.5$ and $P_z = 1$ for simplicity. However, general asymptotic solutions are derived for any $\phi_x$, $\phi_z$ and $P_z$.

### 4.1. Strong exchange

#### 4.1.1. Drag reduction

Figure 4(a) shows how $DR_0$ varies with the bulk diffusion ($\alpha$) and surfactant strength ($\gamma$), for $\beta = 1$. In §4.1, we set $\alpha = \delta$ for simplicity. However, general asymptotic solutions are derived for any $\delta$. The governing equations simplify in different regions of the parameter space, where subsets of the terms in (3.30) are dominant. These distinct physical balances are reflected by limits or transitions of $DR_0$, as well as variations in the concentration profiles. The three primary areas of the parameter space, regions $M$, $D$ and $A$, are separated by black lines in figure 4(a). They are analysed asymptotically in Appendix C. Another more general region ($G$) exists where cross-channel gradients in concentration can be comparable to the cross-channel average (shaded in figure 4a). As region $G$ lies beyond the model predictions, its boundaries delimit the domain of validity of the model (Appendix B.2). To contextualise the changes in $DR_0$, we examine how $c_0$ and $\Gamma_0$ vary across one period in figure 4(b–e). Recall that $D_1$ and $D_2$ represent the subdomains over the plastron and solid ridge, respectively ($D_1$ and $D_2$ are separated by the vertical dotted line at $x = \phi_x = 0.5$), and $c_0 = \Gamma_0$ in $D_1$.
<figure><img src="image_7.png" /><figcaption>Figure 4. The leading-order drag reduction ($DR_0$) and surfactant distribution ($c_0 = Γ_0$) in the strong-exchange problem, for $β = 1$, $φ_x = 0.5$, $ϵ = 0.1$, $φ_z = 0.5$, $Da = 1$ and $P_z = 0.5$, computed using (3.30). (a) Contours of $DR_0$, where $DR_0 = 0$ exhibits a no-slip SHS and $DR_0 = 1$ exhibits a shear-free plastron. The Marangoni (M), advection (A) and diffusion-dominated (D) regions are separated by black lines and $DR_0$ is approximated by (4.1), (4.7) and (4.5) for the M, A and D regions, respectively. The generalised (G) region shaded in brown is discussed in Appendix B.2. The dashed magenta lines describing when $DR_0 ≈ 0.5$ are given by (4.3) and (4.4). (b–e) Plots of $c_0$ for varying surfactant strength ($γ$) and bulk diffusion ($α$), where the asymptotic curves are: (b) $α = δ = 1$, A: -- (4.8), M: ... (4.2) with $γ = 100$, (c) $γ = 0.1$, A: -- (4.8) with $α = δ = 0.01$, D: ... (4.6) with $α = δ = 10$, (d) $α = δ = 0.01$, A: -- (4.8), M: ... (4.2) with $γ = 100$, and (e) $γ = 100$, M: -- (4.2) with $α = δ = 1$, D: ... (4.6) with $α = δ = 1000$. The star identifies the point of the ($α, γ$)-plane where we examine the flow field in § 4.1.2.</figcaption></figure>

To start our journey around the parameter space in figure 4(a), we address region M, in which Marangoni effects are sufficiently strong to render the interface almost immobile. That is, the surfactant is transported along the liquid-gas interface by the flow and accumulates at stagnation points, where the adverse Marangoni stress generated is large
enough to reduce the streamwise velocity at the liquid-gas interface to negligible values along the whole length of the plastron. As shown in Appendix C.1, DR₀ is close to zero in region M and takes the value

$$DR_0 \approx \frac{1}{\gamma} \left( \alpha + \delta + \frac{\epsilon^2 s}{\alpha} + \frac{\phi_x(E+1)}{(E-1)} \right) \quad \text{for } \gamma \gg \max \left( 1, \alpha, \beta, \delta, \frac{\epsilon^2}{\alpha} \right), \quad (4.1)$$

where $E \equiv \exp(2\alpha(1-\phi_x)/(\alpha^2 + \epsilon^2 s_4))$ and $s \equiv s_1 + s_2 + s_3$ (note that $s > 0$ for all $\phi_z$, $P_z$ and $Da$ examined in Appendix B.1). The leading-order drag reduction can be simplified within two sub-regions: $M_D$, in which Marangoni effects start to compete with bulk diffusion, and $M_G$, in which Marangoni effects start to compete with shear dispersion. In the $M_D$ case, (4.1) reduces to $DR_0 \approx \alpha/(\gamma(1-\phi_x)) + \delta/\gamma$ for $\gamma \gg (\alpha, \delta) \gg \max(1, \beta, \epsilon^2/\alpha)$, demonstrating how bulk diffusion (noting that $\alpha/\gamma \propto \hat{\mu}\hat{D}\hat{P}_z/(\hat{A}\hat{G}\hat{H}\hat{L}_d)$) and surface diffusion (noting that $\delta/\gamma \propto \hat{\mu}\hat{D}_I\phi_z\hat{P}_z/(\hat{A}\hat{G}\hat{H}^2)$) weaken the immobilising effects of surfactant. In the $M_G$ case, (4.1) reduces to $DR_0 \approx \epsilon^2(s_4\phi_x + s(1-\phi_x))/(\gamma\alpha(1-\phi_x))$ for $\gamma \gg \epsilon^2/\alpha \gg \max(1, \alpha, \beta, \delta)$, demonstrating how shear dispersion starts to mobilise the interface (noting that $\epsilon^2/(\alpha\gamma) \propto \hat{\mu}\hat{U}^2\hat{H}^5/(\hat{D}\hat{A}\hat{G}\hat{L}_d\hat{P}_z\hat{P}_x^2)$ shows a quadratic dependence on the bulk flow speed).

Representative concentration profiles in region M are shown in figure 4(b,d,e). Here, the leading-order concentration solution in $D_1$ (Appendix C.1) is linear with a shallow gradient,

$$c_0 \approx 1 + \frac{\beta}{\gamma} \left( x - \frac{\phi_x(E+1)}{E-1} \right) \quad \text{for } \gamma \gg \max \left( 1, \alpha, \beta, \delta, \frac{\epsilon^2}{\alpha} \right). \quad (4.2)$$

The gradient of $c_0$ is controlled by both advection and Marangoni effects (noting that $\beta/\gamma \propto 1/Ma = \hat{\mu}\hat{U}/(\hat{A}\hat{G})$, which shows that an advective flow is needed to set up a Marangoni gradient), providing a surface shear stress sufficient to immobilise the liquid-gas interface. From (4.2), the magnitude of the concentration in $D_1$ is weakly regulated by fluxes driven by bulk diffusion or shear dispersion in domain $D_2$, via the factor $\alpha + s_4\epsilon^2/\alpha$ in E that appears due to continuity of concentration across $D_1$ and $D_2$.

The M region transitions into the G region across the GM boundary, where the concentration field starts to develop appreciable cross-channel gradients. In Appendix B.2, we show that a 3D perturbation ($c'_0$) to the cross-channel-averaged surfactant field ($\langle c_0 \rangle$) becomes comparable to $\langle c_0 \rangle$ when Marangoni effects compete with shear dispersion, with $\gamma = O(\epsilon^2/\alpha)$ for $\beta = O(1)$. Similar to streamwise diffusion, shear dispersion increases drag reduction when approaching the GM boundary, as shown in figure 4(a). We conjecture that at 50% drag reduction, shear dispersion becomes large enough to induce an appreciable cross-channel gradient, defining a possible boundary between regions M and G. However, 3D numerical simulations would be required to test this hypothesis and determine precisely the GM boundary, beyond which our asymptotic results (4.1) and (4.2) in region M are no longer valid. Therefore, our asymptotic approximation below to determine the GM boundary based on estimating $DR_0 = 0.5$ should be used qualitatively rather than quantitatively. At the GM boundary, we find (Appendix C.2)

$$DR_0 \approx 0.5 \quad \text{when } \gamma = \frac{\epsilon^2(2s_1(1-\phi_x) + s_4\phi_x(2+\beta))}{2\alpha(1-\phi_x)} \quad \text{for } \alpha \ll 1, \quad (4.3)$$

which gives the leftmost dashed magenta line in figure 4(a). This approximation agrees with the numerical solution of the 1D strong-exchange problem in the limit $\alpha \to 0$. As the immobilising effect of surfactant decreases across the GM boundary, larger streamwise
velocities in the cross-plane imply that $c_0$ decreases from $c_0 \approx 1$ to smaller values along the channel (see figure 4d). This is due to the velocity and surfactant flux being fixed.

The M region transitions into the D region through the DM boundary for $\gamma = O(\alpha)$ or $\gamma = O(\delta)$, which can be defined (see Appendix C.3) by

$$DR_0 \approx 0.5 \quad \text{when } \gamma = \delta + \frac{\alpha(2 + \beta\phi_x)}{2(1 - \phi_x)} \quad \text{for } \min(\alpha, \delta) \gg 1, \quad (4.4)$$

which gives the rightmost dashed magenta line in figure 4(a). The asymptotic approximation (4.4) agrees with the numerical solution of (3.30) (thick black line in figure 4a) for $\min(\alpha, \delta) \gg 1$. Similar to the GM transition, the immobilising effect of surfactant decreases as we move from M to D by increasing the strength of bulk diffusion. By increasing the strength of bulk diffusion, bulk and interfacial surfactant gradients are attenuated because bulk-surface exchange is strong; this allows the streamwise velocity and drag reduction to increase. As the velocity and surfactant flux are fixed, the increase in velocity leads to a decrease in $c_0$, as shown in figure 4(e).

In summary, the factors promoting drag reduction from a state of interfacial immobilisation (in region M) are bulk diffusion (moving into region D), shear dispersion (moving into region G) or reduced surfactant strength (moving into region G, D or A, of which more details are given below). As we move into region D (or G), the surfactant in the bulk of the channel spreads out via diffusion (or shear dispersion), transporting surfactant from areas of high concentration to low concentration. Then, because exchange between the bulk and interfacial surfactant is strong, the concentration fields rapidly equilibrate; the interfacial surfactant gradient is also attenuated by these diffusive processes and the drag reduction $DR_0$ increases.

We next turn our attention to region D, where diffusion is sufficiently strong for surfactant to be distributed almost uniformly along the bulk and interface. This implies that there are almost no Marangoni stresses to increase the drag and the liquid-gas interface is almost shear-free. As shown in Appendix C.4, $DR_0$ is close to unity in region D and is given by

$$DR_0 \approx 1 - \frac{(1 - \phi_x)\gamma}{(1 + \phi_x\beta)\alpha + (1 - \phi_x)\delta} \quad \text{for } \min(\alpha, \delta) \gg \max(1, \gamma), \quad (4.5)$$

for $\beta = O(1)$. Equation (4.5) shows how drag-promoting Marangoni effects are weakened by strong diffusion in the bulk ($\alpha$) or at the interface ($\delta$). Figure 4(c,e) shows the surfactant profiles in region D in the limit of strong diffusion, where $\alpha$ or $\delta$ are large compared with bulk advection and Marangoni effects. The leading-order solution in $D_1$ is uniform along $x$, such that

$$c_0 \approx \frac{\alpha + \delta(1 - \phi_x)}{\alpha(\beta\phi_x + 1) + \delta(1 - \phi_x)} \quad \text{for } \min(\alpha, \delta) \gg \max(1, \gamma), \quad (4.6)$$

demonstrating how diffusion eliminates gradients of $c_0$ throughout $D_1$ and $D_2$ (the weak gradients contributing to $DR_0$ in (4.5) appear at higher order and are detailed in Appendix C.4). Therefore, drag reduction in region D is impeded by decreasing the bulk and surface diffusion (if $\gamma \gg 1$) or increasing the surfactant strength (i.e. moving into region M).

We finally consider region A, in which surfactants can adsorb onto the interface but generate weak Marangoni stresses, allowing the streamwise flow to advect the interfacial surfactant towards the downstream stagnation point. In this advection-dominated region,
we show in Appendix C.5 that DR₀ is once again close to the shear-free value

$$DR_0 \approx 1 - \frac{\gamma}{2\phi_x (\beta + 1)} \quad \text{for } \epsilon^2 \ll \alpha \ll 1, \quad \gamma \ll \min(1, \beta). \quad (4.7)$$

Region A can also be divided into two sub-regions: A<sub>G</sub> which balances advection and shear dispersion, and A<sub>D</sub> which balances advection and diffusion; these balances give the same DR<sub>0</sub> as in (4.7). The corresponding concentration profiles in region A are shown in figure 4(b-d). The leading-order solution in D<sub>1</sub> exhibits a surfactant gradient that increases monotonically towards the downstream end of the plastron

$$c_0 \approx \frac{1}{\beta + 1} + \frac{\beta}{(\beta + 1)} \exp \left( \frac{(1 + \beta)(x - \phi_x)}{\alpha + \delta + \epsilon^2 s_1 / \alpha} \right) \quad \text{for } \epsilon^2 \ll \alpha \ll 1, \quad \gamma \ll \min(1, \beta). \quad (4.8)$$

The downstream boundary layer in (4.8) can be short compared with the plastron length, such that the surface shear stress is negligible almost everywhere on the interface and the bulk flow experiences a largely shear-free boundary.

Region A transitions into region D across the AD boundary, where advection and diffusion balance. The concentration and drag reduction at the AD boundary can be obtained in closed form following the method shown in Appendix C.5. However, in figure 4(a), we see that crossing the AD boundary (shown by the vertical black line at $\alpha = 1$ and $\gamma \ll 1$) does not affect DR<sub>0</sub> to leading order. When $\gamma \ll \min(1, \beta)$ and the surfactant strength is weak, either the streamwise velocity at the surface is large enough to advect most of the interfacial surfactant to the downstream stagnation point ($\alpha \ll 1$) or diffusion is strong enough to attenuate any interfacial surfactant gradient that forms ($\alpha \gg 1$). As the interface is mostly shear-free in both A and D, we do not pursue this limit. Bulk diffusion dominates over bulk advection for $\alpha \gg O(1)$, so we employ $\alpha = 1$ to illustrate the AD boundary in figure 4(a) (via the thick vertical black line on the right of the A region). When bulk diffusion is weak, shear dispersion dominates over bulk advection for $\epsilon^2/\alpha \gg 1$ (Appendix B.2), so we use $\alpha = \epsilon^2$ to illustrate the AG boundary (the thick vertical black line on the left of the A region in figure 4a). Region A transitions into region M across the AM boundary, such that Marangoni effects dominate advection for $\gamma \gg \max(1, \beta)$.

We next investigate the dependence of the drag reduction on the partition coefficient ($\beta$) by evaluating DR<sub>0</sub> for $\beta = 10$ in figure 5(a) and $\beta = 100$ in figure 5(b). As $\beta$ grows, the advective flux of surfactant at the liquid-gas interface increases, sweeping more of the interfacial surfactant towards the downstream stagnation point. This means that a larger portion of the upstream part of the interface is shear-free, however, the surfactant concentration remains unchanging at the downstream end of the interface. Overall, the net (nonlinear) effect with increasing $\beta$ is a reduction in drag. The limit $\beta \gg 1$ also corresponds to the near-insoluble surfactant limit, with nearly all surfactant being transported along the interface in region D<sub>1</sub>. A number of the asymptotic approximations given in (4.1)-(4.8) simplify for $\beta \gg 1$: (4.3) reduces to DR<sub>0</sub> ≈ 0.5 when $\gamma = \epsilon^2 s_4 \beta \phi_x / (2\alpha(1 - \phi_x))$ for $\alpha \ll 1 \ll \beta$, such that shear dispersion in D<sub>2</sub> determines DR<sub>0</sub>; (4.4) reduces to DR<sub>0</sub> ≈ 0.5 when $\gamma = \alpha \beta \phi_x / (2(1 - \phi_x))$ for $\alpha \gg \beta \gg 1$ and $\alpha = O(\delta)$, with bulk diffusion determining DR<sub>0</sub>.

Another feature that emerges from figure 5(a,b) is the flattening of the central contours of DR<sub>0</sub> at the AM boundary as $\beta$ increases for $\beta \gg 1$, due to advection at the interface becoming stronger. To extract this feature asymptotically, in Appendix C.6 we assume $\beta \gg 1 \gg \max(\alpha, \delta, \epsilon^2/\alpha)$ and $\gamma = O(\beta)$ at the AM boundary, showing that
<figure><img src="image_8.png" /><figcaption>Figure 5. The leading-order drag reduction ($DR_0$) and surfactant distribution ($c_0 = \Gamma_0$) in the strong-exchange problem, for $\phi_x = 0.5$, $\epsilon = 0.1$, $\phi_z = 0.5$, $Da = 1$ and $P_z = 0.5$, computed using (3.30). Contours of $DR_0$ for (a) $\beta = 10$ and (b) $\beta = 100$, where $DR_0 = 0$ exhibits a no-slip SHS and $DR_0 = 1$ exhibits a shear-free plastron at the SHS. The Marangoni (M), advection (A) and diffusion-dominated (D) regions are separated by black lines and $DR_0$ is approximated by (4.1), (4.7) and (4.5) for the M, A and D regions, respectively. The dashed magenta lines describing when $DR_0 \approx 0.5$ are given by (4.9) and (4.4) for the AM and DM boundaries, respectively. Plots of $c_0$ for varying surfactant strength ($\gamma$) and $\alpha = \delta = 1$, where (c) $\beta = 10$, A: -- (4.8), M: ... (4.2) with $\gamma = 1000$, and (d) $\beta = 100$, A: -- (4.8), M: ... (4.2) with $\gamma = 10000$.</figcaption></figure>

$DR_0 \approx 1 - \gamma/(2\beta\phi_x)$ for $\gamma/\beta \le 2\phi_x$ and $\epsilon^2 \ll \alpha \ll 1$ (the case where $\gamma/\beta > 2\phi_x$ is also considered in Appendix C.6 and gives $DR_0 \approx -(1/\beta + \gamma \ln(1-2\beta\phi_x/\gamma)/(2\phi_x\beta^2))$ provided $1-2\phi_x\beta/\gamma \gg \exp(-\beta)$). This demonstrates how a large surface advective flux weakens the immobilising effects of surfactant (noting that $\gamma/\beta \propto \hat{A}\hat{G}/\hat{\mu}\hat{U}$). Using this expression for $DR_0$ at the AM boundary, we predict that

$$DR_0 \approx 0.5 \quad \text{when } \gamma = \phi_x \beta \quad \text{for } \beta \gg 1 \gg \max \left( \alpha, \delta, \frac{\epsilon^2}{\alpha} \right), \qquad (4.9)$$
which gives the central dashed magenta lines in figure 5(a,b) (agreement with the 1D numerical solution of the strong-exchange problem improves for increasing $\beta$). At the AM boundary, the surfactant concentration at the interface exhibits a stagnant cap (He et al. 1991), or a piecewise-linear distribution, if $\gamma/\beta \le 2\phi_x$, where the upstream interface is almost shear-free and the downstream interface is effectively no-slip. In Appendix C.6, we show that

$$c_0 \approx \begin{cases} 0 & \text{for } -\phi_x \le x \le x_0, \\ \frac{\beta}{\gamma}(x - \phi_x) + 1 & \text{for } x_0 \le x \le \phi_x, \end{cases} \quad \text{for } \min(\beta, \gamma) \gg 1 \gg \max\left(\alpha, \delta, \frac{\epsilon^2}{\alpha}\right), \quad (4.10)$$

where $x_0 = \phi_x - \gamma/\beta$ (see, e.g., the red curve for $\gamma = 100$ in figure 5d). If $\gamma/\beta > 2\phi_x$, then the interfacial concentration distribution at the AM boundary is linear, $c_0 \approx \beta/\gamma(x - \phi_x) + 1$ at leading order. When the surface advection overcomes Marangoni effects, for $\beta \gg \gamma$, surfactants can accumulate at the downstream stagnation point. The blue curve for $\gamma = 10$ in figure 5(d) demonstrates that $c_0 \approx 0$ throughout most of $D_1$: the large surfactant advective flux at the interface contracts the downstream boundary layer, making the interface almost shear-free. In summary, stronger nonlinearity in the form of stagnant cap profiles in the interfacial concentration can appear with increasing $\beta$ near the transition from region A to M. This results in a sharp decline in $DR_0$ towards zero, which can be seen as the AM transition becomes less smooth, i.e. the vertical distance between $DR_0 = 0.5$ and $DR_0 = 0.05$ contours decreases from figure 5(a) to (b) as the surface advection increases from $\beta = 10$ to $\beta = 100$. However, we note that the location of the transition is captured by the simple expression (4.9) when $\gamma/\beta \le 2\phi_x$.

### 4.1.2. Flow field

Whenever there is partial drag reduction, $0 \le DR_0 < 1$, the surfactant has a non-uniform concentration in $D_1$, reducing the streamwise velocity at the plastron and generating a secondary flow in the bulk. To illustrate, videos of the leading-order streamwise velocity $u_0$, wall-normal velocity $v_1$, transverse velocity $w_1$ and the streamwise gradient of the surfactant distribution $c_{0x}$ are given in supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.264, evaluated at a point in the parameter space where $DR_0 \approx 0.5$, shown by the star in figure 4(a). The corresponding flow field is shown in figure 6. Figure 6(a) shows $u_0$ at the centre of domain $D_1$, where $x = 0$. The streamwise velocity, built from the components shown in figure 2, attains a maximum value at the channel centre, decaying towards either SHS ($y = 0, 2$) to satisfy no-slip at the solid wall ($-P_z \le z \le \phi_z$ and $\phi_z \le z \le P_z$), whilst allowing for slip over the plastron ($-\phi_z \le z \le \phi_z$). The slip velocity at $\mathcal{I}$ reduces as one progresses through $D_1$, in accordance with the rise in $c_{0x}$ for increasing $x$; see figure 6(c) and supplementary movie 1. The associated secondary cross-channel velocity field ($v_1, w_1$) shown in figure 6(b) (recall that the leading-order components are $(v_0, w_0) = (0, 0)$) advects particles from $\mathcal{I}$ towards the centre of the channel ($y = 1$ and $z = \pm P_z$). The absolute maximum of $v_1$ occurs above the centre of the plastron ($z = 0$), whereas the absolute maximum of $w_1$ occurs above the transverse contact lines ($z = \pm \phi_z$). Supplementary movie 1 and figure 6(d) show how $v_1$ and $w_1$ grow in magnitude as one progresses through $D_1$ along the positive x-direction as the interface is immobilised. As mentioned in § 3.2.1, the present long-wave model does not capture rapid adjustments of the flow field near the contact lines at $x = \pm \phi_x$, where domains $D_1$ and $D_2$ meet.
<figure><img src="image_9.png" /><figcaption>Figure 6. Contour maps of the flow field in the strong-exchange problem, for $\alpha = 0.1$, $\beta = 1$, $\gamma = 1.3$, $\delta = 0.1$, $\phi_x = 0.5$, $\phi_z = 0.5$ and $P_z = 0.5$, corresponding to 50% normalised drag reduction given by the star in figure 4(a). (a) Leading-order streamwise velocity $u_0$ and (b) leading-order wall-normal velocity $v_1$ with ($v_1, w_1$) streamlines at the centre of the plastron, $x = 0$. (c) Leading-order streamwise velocity $u_0$ and (d) leading-order transverse velocity $w_1$ with ($u_0, w_1$) streamlines at the interfaces, $y = 0$ or 2. The thick black lines in (a–d) represent the solid regions of the SHS.</figcaption></figure>

## 4.2. Moderate exchange

We now turn our attention to the moderate-exchange problem (3.1)–(3.3a–d), in which the bulk and interfacial concentration fields decouple for $v/\epsilon^2 \sim O(1, \alpha, \beta, \delta)$. At the plastron, surfactant adsorbs onto the interface at the upstream end ($c_0 > \Gamma_0$) and desorbs at the downstream end ($\Gamma_0 > c_0$), as regulated by the bulk–surface exchange parameter $v$ in (3.35). The moderate-exchange model does not take into account shear-dispersion effects, however, we ensure that the bulk diffusion strength is larger than the threshold identified in the strong-exchange problem (§4.1.1), $\alpha \gg \epsilon^2$, where cross-channel concentration gradients first become significant.

### 4.2.1. Drag reduction

Figure 7(a) shows how $DR_0$ varies with surface exchange strength ($v$) and surfactant strength ($\gamma$), for $\alpha = 1$, $\beta = 1$ and $\delta = 1$. Similar to the strong-exchange problem in §4.1.1, the moderate-exchange governing equations (3.1)–(3.2) simplify in different regions of the parameter space, representing distinct physical balances. The analysis in the moderate-exchange problem is simplified by the fact that for $v/\epsilon^2 \gg \max(1, \alpha, \beta, \delta)$ the moderate-exchange problem transitions to the strong-exchange problem, which is discussed in detail in §4.1.1. In the strong-exchange limit, $v/\epsilon^2 \gg \max(1, \alpha, \beta, \delta)$, we reach the AD boundary and region $M$ identified in figure 4(a). Accordingly, the regimes
<figure><img src="image_10.png" /><figcaption>Figure 7. The leading-order drag reduction ($DR_0$), bulk surfactant ($c_0$) and interfacial surfactant distribution ($\Gamma_0$) in the moderate-exchange problem, for $\alpha = 1$, $\beta = 1$, $\delta = 1$, $\phi_x = 0.5$, $\epsilon = 0.1$, $\phi_z = 0.5$ and $P_z = 0.5$, computed using (3.1)–(3.3a–d). (a) Contours of $DR_0$, where $DR_0 = 0$ exhibits a no-slip SHS and $DR_0 = 1$ exhibits a shear-free plastron. The Marangoni ($M$), advection–diffusion ($AD$), Marangoni-exchange ($M_E$) and advection–diffusion-exchange ($AD_E$) regions are separated by black lines and $DR_0$ is approximated by (4.1), (4.7), (4.11a,b) and (4.15) in the $M$, $AD$, $M_E$ and $AD_E$ regions, respectively. The dashed magenta lines describing when $DR_0 \approx 0.5$ are given by (4.4) and (4.13). Plots of $\Gamma_0$ and $c_0$ for varying surfactant strength ($\gamma$) and exchange strength ($\nu$), where the asymptotic curves are: (b) $\nu = 10$, $AD: \dots (4.2)$ with $\gamma = 100$, $M: \dots (4.8)$, (c) $\gamma = 0.03$, $AD_E: \dots (4.16)$, $AD: \dots (4.8)$, (d) $\nu = 0.001$, $AD_E: \dots (4.16)$, $M_E: \dots (4.12a,b)$ with $\gamma = 10$, and (e) $\gamma = 100$, $M_E: \dots (4.12a,b)$, $M: \dots (4.2)$.</figcaption></figure>

$AD$ and $M$ are also identified in figure 7(a) where bulk-surface exchange is strong. New regimes $M_E$ and $AD_E$ appear for weak exchange, $\nu/\epsilon^2 \ll \max(1, \alpha, \beta, \delta)$, and are analysed asymptotically in Appendix D. The four primary areas of the parameter space, namely $M_E$, $AD_E$, $M$ and $AD$, are separated by black lines in figure 7(a). Corresponding surface concentration distributions $\Gamma_0(x)$ (where, unlike the strong-exchange problem, $c_0 \neq \Gamma_0$)
in $D_1$) are illustrated in figure 7(b–e). The primary feature of figure 7(a) is that the drag-reduction transition from large $DR_0$ to small $DR_0$ shifts only modestly, despite the exchange strength $\nu$ varying across many orders of magnitude, for reasons that we explain in the following paragraph.

We begin our exploration of the parameter space in the moderate-exchange problem (figure 7) with region $M_E$, where Marangoni effects are strong and surfactant exchange between the bulk and interface is weak. In Appendix D.1 we show that $DR_0$ remains close to the immobilised value

$$DR_0 \approx \frac{\delta}{\gamma} \quad \text{for } \gamma \gg \max(1, \alpha, \beta, \delta), \quad \frac{\nu}{\epsilon^2} \ll \min(1, \alpha, \beta, \delta). \quad (4.11a,b)$$

Here, the immobilising effects of surfactant are weakened by strong interfacial diffusion ($\delta/\gamma \propto \hat{\mu}D_I\phi_z\hat{P}_z/(\hat{A}\hat{G}\hat{H}^2)$). The corresponding concentration fields $c_0$ and $\Gamma_0$ are shown in figure 7(d,e). The leading-order solution on $\mathcal{I}$ in $M_E$ is given by

$$\Gamma_0 \approx 1 + \frac{\beta x}{\gamma} \quad \text{for } \gamma \gg \max(1, \alpha, \beta, \delta), \quad \frac{\nu}{\epsilon^2} \ll \min(1, \alpha, \beta, \delta). \quad (4.12a,b)$$

The gradient of $\Gamma_0$ in region $M_E$ is controlled by the relative strength of surface advection and Marangoni effects, $\beta/\gamma$, which is sufficiently small to immobilise the liquid-gas interface. Furthermore, there is very little adsorption and desorption at the interface and the bulk and interfacial surfactant concentrations are close to their background values at leading order (i.e. $c_0 \approx 1$ and $\Gamma_0 \approx 1$). In contrast to (4.2) in $M$ and the strong-exchange problem of §4.1.1, where $c_0 = \Gamma_0$ and $\Gamma_0 < 1$ for all $x$, in $M_E$ and the weak-exchange problem, $c_0 \approx 1$ and $\Gamma_0 > 1$ for $x > 0$ (see, e.g., the dashed and blue curves in figure 7e), in order to satisfy the net flux condition (2.26) at leading order. Furthermore, $c_0$ has turning points for $-\phi_x < x < \phi_x$ where the adsorption-desorption fluxes attain a local maximum and then decrease towards the contact lines, as shown by the green and orange curves in the inset in figure 7(d).

Region $M_E$ gives way to region $M$ as the exchange strength increases. It is notable that $DR_0 \approx \delta/\gamma$ (see (4.11a,b)) is regulated mainly by surface diffusion in region $M_E$, whereas both bulk and surface diffusion control $DR_0$ in region $M$, where $DR_0 \approx \alpha/(\gamma(1-\phi_x)) + \delta/\gamma$ (see §4.1.1). This is due to the absence of the strong coupling between $\Gamma_0$ and $c_0$ in region $M_E$, which imposed $\Gamma_0 = c_0$ in region $M$. In figure 7(a), we find that $\nu = 10\epsilon^2$ best captures the centre of this transition (this can also be tested by varying $\epsilon$). As exchange weakens and the bulk and interfacial surfactant move out of equilibrium, variations in $\Gamma_0$ increase, thus making the interface more susceptible to Marangoni effects that increase the drag.

Depending on the strength of diffusion (via the parameters $\alpha$ and $\delta$, which vary orthogonally to the $(\nu, \gamma)$-plane in figure 7a), the lower half of the $(\nu, \gamma)$-plane is composed of either a sub-region of $A$, $D$, or the transition region $AD$ between $A$ and $D$. In figure 7(a), because $\alpha = \delta = 1$, this corresponds to the transition region $AD$. In contrast, in figure 8(a) we plot $DR_0$ for $\alpha = \delta = 0.1$ and $\alpha = \delta = 10$ in figure 8(b). The $M_E$ region transitions into the diffusion-exchange-dominated ($D_E$) region across the $M_E D_E$ boundary for $\min(\alpha, \delta) \gg O(1)$ (figure 8b), whereas the $M_E$ region transitions into the advection-exchange-dominated ($A_E$) region across the $M_E A_E$ boundary for $\max(\alpha, \delta) \ll O(1)$ (figure 8a). At the $M_E D_E$ boundary, the asymptotic analysis in Appendix D.2 shows that

$$DR_0 \approx 0.5 \quad \text{when } \gamma = \delta \quad \text{for } \min(\alpha, \delta) \gg 1, \quad \frac{\nu}{\epsilon^2} \ll 1, \quad (4.13)$$
<figure><img src="image_11.png" /><figcaption>Surfactant-contaminated superhydrophobic channels</figcaption></figure>

Figure 8. The leading-order drag reduction ($DR_0$), bulk surfactant ($c_0$) and interfacial surfactant distribution ($\Gamma_0$) in the moderate-exchange problem, for $\beta = 1$, $\phi_x = 0.5$, $\epsilon = 0.1$, $\phi_z = 0.5$ and $P_z = 0.5$, computed using (3.1)–(3.3a–d). Contours of $DR_0$ for (a) $\alpha = \delta = 0.1$ and (b) $\alpha = \delta = 10$, where $DR_0 = 0$ exhibits a no-slip SHS and $DR_0 = 1$ exhibits a shear-free plastron at the SHS. The Marangoni (M), advection (A), diffusion (D), Marangoni-exchange ($M_E$), advection-exchange ($A_E$), diffusion-exchange ($D_E$) regions are separated by black lines and $DR_0$ is approximated by (4.1), (4.7), (4.5), (4.11a,b), (4.15) and (4.15) in M, A, D, $M_E$, $A_E$ and $D_E$, respectively. The dashed magenta lines describing when $DR_0 \approx 0.5$ are given by (4.4) in (a,b) (right curves), (4.13) in (b) and (4.14) (a). Plots of $\Gamma_0$ and $c_0$ for varying surfactant strength ($\gamma$), where (c) $\alpha = \delta = 0.1$ and $\nu = 10^{-4}$, $A_E: --$ (4.16), $M_E: \cdots$ (4.12a,b) with $\gamma = 10$, and (d) $\alpha = \delta = 10$ and $\nu = 0.01$, $D_E: --$ (4.17), $M_E: \cdots$ (4.12a,b) with $\gamma = 1000$. The star identifies the point of the ($\nu, \gamma$)-plane where we examine the flow field in §4.2.2.

which gives the leftmost dashed line in figure 8(b). This asymptote agrees with the numerical solution presented in figure 8(b) for $\delta \gg 1$ and $\nu \to 0$. The range of validity for the asymptote given in (4.13) is extended to $\delta = O(1)$ to include the $M_E A_D E$ boundary, for which we give partial justification in Appendix D.3. At the $M_E A_E$ boundary, Appendix D.4
demonstrates that

$$DR_0 \approx 0.5 \quad \text{when } \gamma = \frac{\beta\phi_x}{4} + \frac{\delta}{2} \quad \text{for } \max(\alpha, \delta) \ll 1, \quad \frac{\nu}{\epsilon^2} \ll 1, \quad (4.14)$$

which gives the leftmost dashed line in figure 8(a). This asymptote approximates the numerical solution presented in figure 8(a) where $\delta = 0.1$ as $\nu \to 0$. Agreement improves for smaller $\delta$. It is notable that the threshold (4.14) differs from the strong-exchange limit (4.9) by a factor of four.

We next analyse the drag reduction in both the $A_E$ and $D_E$ regions, where advection or diffusion induce a near-uniform distribution of surfactant at the interface. In Appendix D.5 we show that $DR_0$ in regions $A_E$ and $D_E$ is close to the shear-free value

$$DR_0 \approx 1 - \frac{\gamma}{\delta} \quad \text{for } \gamma \ll \min(1, \alpha, \beta, \delta), \quad \frac{\nu}{\epsilon^2} \ll \min(1, \alpha, \beta, \delta). \quad (4.15)$$

Equation (4.15) shows how drag-promoting Marangoni effects are weakened by strong diffusion at the interface, when surfactant exchange between the bulk and interface is weak. Figure 8(c,d) shows the surfactant profiles in $A_E$ and $D_E$, the leading-order solution on $\mathcal{I}$ being given by

$$\Gamma_0 \approx \frac{2\phi_x\beta\exp\left(\frac{\beta(\phi_x+x)}{\delta}\right)}{\delta\left(\exp\left(\frac{2\beta\phi_x}{\delta}\right) - 1\right)} \quad \text{for } \gamma \ll \min(1, \alpha, \beta, \delta), \quad \frac{\nu}{\epsilon^2} \ll \min(1, \alpha, \beta, \delta). \quad (4.16)$$

From (4.16), weak surface diffusion, $\delta \ll O(1)$, decreases (increases) $\Gamma_0$ at the upstream (downstream) end of the interface. A downstream boundary layer forms to satisfy the mass balance condition (2.26), condensing the Marangoni effect to a small region near $x = \phi_x$. The streamwise velocity flows over an almost shear-free boundary at the upstream end of the interface. When surface diffusion is strong, $\delta \gg O(1)$, (4.16) reduces to

$$\Gamma_0 = 1 + \frac{\beta x}{\delta} \quad \text{for } \gamma \ll \min(1, \alpha, \beta, \delta), \quad \frac{\nu}{\epsilon^2} \ll \min(1, \alpha, \beta, \delta). \quad (4.17)$$

Noting that $\beta/\delta \propto \hat{H}^2\hat{U}/(\hat{D}_I\phi_z\hat{P}_z)$, (4.17) demonstrates how the gradient of $\Gamma_0$ is controlled by the ratio of advection to diffusion at the interface, which allows surface diffusion to regulate $DR_0$ in (4.15). In figure 8(a,b), the horizontal boundary of the $M_E$ region (shown with a black solid line) moves downwards when surface diffusion decreases ($\delta = 0.1$ in figure 8a) and upwards when it increases ($\delta = 10$ in figure 8b). Figure 8(c,d) demonstrates how increasing surface diffusion attenuates gradients of surfactant at the interface, increasing the drag reduction from sub-region $A_E$ to $D_E$ in figure 8(a,b).

### 4.2.2. Flow field

Videos of $u_0$, $v_1$, $w_1$ and $\Gamma_{0x}$ are given in supplementary movie 2. Again, we choose an example in figure 9 for which $DR_0 = 0.5$, shown with a star in figure 8(a). Figure 9(a) shows that $u_0$ has a similar structure to the strong-exchange problem at $x = 0$ (see figure 6a). However, when comparing figure 9(c) with figure 6(c), $u_0$ exhibits different behaviour along $\mathcal{I}$. The streamwise velocity decreases slowly at the upstream end of $\mathcal{I}$, $\Gamma_{0x}$ and $\Gamma_{0xx}$ increase slowly with $x$, and there is a uniform distribution of surfactant at the interface which is almost shear-free. The streamwise velocity decreases rapidly at the
<figure><img src="image_12.png" /><figcaption>Figure 9. Contour maps of the flow field in the moderate-exchange problem, for $\alpha = 0.1$, $\beta = 1$, $\gamma = 0.2$, $\delta = 0.1$, $\nu = 10^{-4}$, $\phi_x = 0.5$, $\phi_z = 0.5$ and $P_z = 0.5$, corresponding to 50% normalised drag reduction given by the star in figure 8(a). (a) Leading-order streamwise velocity $u_0$ and (b) leading-order wall-normal velocity $v_1$ with ($v_1, w_1$) streamlines at the centre of the plastron, $x = 0$. (c) Leading-order streamwise velocity $u_0$ and (d) leading-order transverse velocity $w_1$ with ($u_0, w_1$) streamlines at the interfaces, $y = 0$ or 2. The thick black lines in (a–d) represent the solid regions of the SHS.</figcaption></figure>

centre of $\mathcal{I}$, where there is a sudden change of surfactant gradient at the interface. The streamwise velocity then decreases at the downstream end of $\mathcal{I}$, where there is a more linear distribution of surfactant at the interface which is almost no-slip. The cross-channel flow in figure 9(b) closely resembles that observed in figure 6(b). In supplementary movie 2 and figure 9(d), the magnitudes of both $v_1$ and $w_1$ increase until they attain a maximum at the start of the low-slip region, after which, they decrease towards $x = \phi_x$. The streamwise location of the maximum of $v_1$ and $w_1$ is approximately the same as the streamwise location of the maximum of $\Gamma_{0xx}$ on $\mathcal{I}$, the ‘corner’ in the surfactant field.

4.3. *Comparison with numerical simulations*

Finally, we compare our model predictions with the numerical simulations detailed in tables SI and SII and figures S1 and S2 in Temprano-Coleto et al. (2023), which were designed to be representative of microchannel applications. They used finite-element simulations to solve the steady 3D Stokes equations, which were coupled to advection–diffusion equations for bulk and interfacial surfactant. The equation of state and adsorption–desorption kinetics, which link the velocity field to the bulk and interfacial concentrations, were derived from the nonlinear Frumkin isotherm (Chang & Franses 1995). The authors considered a channel with only one SHS at the bottom wall ($y = 0$), so our theory in § 3 is adjusted accordingly for a solid surface at the top wall ($y = 2$).
<figure><img src="image_13.png" /><figcaption>Figure 10. The drag reduction ($DR_0$) and bulk surfactant distribution evaluated at the SHS ($c_0(x, 0, 0)$), using the theory presented here and numerical simulations provided by Temprano-Coleto et al. (2023). (a) Contours of $DR_0$, where $DR_0 = 0$ exhibits a no-slip SHS and $DR_0 = 1$ exhibits a shear-free SHS, for $\beta = 38.6$, $\delta = 1$, $\nu = 0.2$, $\phi_x = 0.99$, $\phi_z = 2/3$ and $P_z = 1$. The Marangoni (M) and diffusion-dominated (D) regions are separated by the black line (along which $DR_0 = 0.5$). The advection-dominated (A) region appears for $\alpha < 0.01$ and $\gamma < 100$. (b) Plot of $c_0$, where $\alpha = 0.4$, $\beta = 38.6$, $\gamma = 1.2 \times 10^4$, $\delta = 1$, $\nu = 0.2$, $\phi_x = 0.99$, $\phi_z = 2/3$ and $P_z = 1$ (corresponding to the stars in a,c), computed using (3.30) (red) and (4.2) (dashed and black). (c) Scatter plot of $DR_0$ using our theory ((3.38) and (3.1)-(3.3a–d)) and the 159 numerical simulations detailed in Temprano-Coleto et al. (2023), the colourbar gives the magnitude of $\epsilon$ for given data point and the dashed line is where $DR_0 = 0.5$. (d) Same data as in (c): scatter plot of $DR_0 \in [0.1, 1]$. The colourbar gives the magnitude of $\alpha$ for given data point and the black points have $\alpha \gg 1.5$. The arrows indicate the region of parameter space for a given data point: orange points are in region D, yellow points are in region A and blue points are in region M.</figcaption></figure>

This amounts to reevaluating $\tilde{U}$, $\tilde{U}$ and $\tilde{U}$, as well as modifying the surface surfactant flux term in (3.10) to account for the contribution of a single interface at $y = 0$. Parameters are detailed in table SI in Temprano-Coleto et al. (2023).

In figure 10(a,b) we compare our model to a representative example from numerical simulations performed in Temprano-Coleto et al. (2023), with parameters $\hat{H} = 6 \times 10^{-5}$, $\hat{P}_z = 6 \times 10^{-5}$ m and $\phi_z = 2/3$, from which we can evaluate $P_z$, $\tilde{Q}$, $\tilde{Q}$, $\tilde{q}$ and $\tilde{q}$ using (3.18a,b), (3.19a-c). Using (3.27) (appropriately adjusted for a single SHS), we can calculate the transport coefficients: $\alpha = 0.4$, $\beta = 38.6$, $\gamma = 1.2 \times 10^4$, $\delta = 1$, $\epsilon = 0.02$, $\nu = 0.2$ and $\phi_x = 0.99$; with $Pe = Pe_l = 20.6$, $Da = 30.1$, $Ma = 1.2 \times 10^4$ and $Bi = 0.2$, as in Temprano-Coleto et al. (2023). Note that Ma is defined differently to Temprano-Coleto et al. (2023); the value provided here uses our definition of Ma. Figure 10(a) shows that the example lies within the Marangoni-dominated region M, where $DR_0$ is close to zero and the liquid-gas interface is immobilised. Using this information, in figure 10(b), we compare the numerically simulated bulk surfactant concentration
$c_0(x, 0, 0)$ from Temprano-Coleto et al. (2023) (blue solid curve) with our prediction for $c_0$ given by the solution to the 1D ODE model solving (3.30) subject to (3.3a,b) (red solid line), and the (indistinguishable) asymptotic solution in $M$ (4.2) (black dashed line). The gradient of $c_0$ at the centre of the plastron computed using our theory ($dc_0/dx(0) \approx 0.19$) and the numerical simulation ($dc/dx(0, 0, 0) \approx 0.16$) of Temprano-Coleto et al. (2023) are similar, noting that our theoretical prediction does not require any fitting parameters. However, the numerical simulation reveals thin boundary layers close to the upstream and downstream contact lines. These boundary layers are not captured by our long-wave theory. Within these inner regions, the 3D Stokes and surfactant transport equations govern the flow and surfactant field near the no-flux stagnation points at the upstream and downstream ends of the interface. This comparison suggests that the inner layers at the upstream and downstream ends are not needed to estimate the leading-order drag reduction in region $M$, but capturing the inner layers may be necessary for more accurate predictions.

In figure 10(c), we compare the drag reduction results from all the 3D numerical simulations of Temprano-Coleto et al. (2023) with the leading-order drag reduction predictions from our 1D model. These simulations span the whole parameter space characteristic of realistic microchannel applications. In order to approximate the drag reduction from the streamwise slip length results provided by Temprano-Coleto et al. (2023), we integrate the streamwise velocity field when $\phi_z = 0$ for the solid-walled and SHS flows. For the flow over a single SHS, we replace the mixed boundary conditions at $y = 0$ with $\lambda_e u_{0y} - u_0 = 0$, to find $u_0$ as a function of $y$ and $\lambda_e$. We then use $DR_0 = (\Delta p_R - \Delta p_0)/(\Delta p_R - \Delta p_U)$ to relate the leading-order drag reduction and streamwise slip length, $DR_0 = 1 - 1/(2\lambda_e + 1)$, in a similar manner to Landel et al. (2020). There is no obvious correlation between the scatter and the size of $\epsilon$ (shown in colour), which varies up to $\epsilon \approx 0.5$ as highlighted by the colourbar in figure 10(c). The scatter could be due to unresolved cross-channel concentration gradients or streamwise boundary layers at the ends of the plastron; both hypotheses require further testing against targeted numerical simulations. Nonetheless, our theory compares well with simulations for $DR_0 > 0.1$ in figure 10(d); the root-mean-squared error between the simulated data from Temprano-Coleto et al. (2023) and our theory is 0.07. The majority of the data points have $DR_0 < 0.5$ and lie in region $M$. Those that have $DR_0 > 0.5$ are classified as belonging to region A when $\alpha < 1$ or D when $\alpha > 1$, as highlighted by the colourbar in figure 10(d). None of the simulations with high drag reduction ($DR > 0.5$) had a bulk diffusion strength small enough to lie on the boundary of region $G$ ($\alpha = O(\epsilon^2)$), where cross-channel concentration gradients become important and shear dispersion decreases the drag. Hence, the effect of shear dispersion we have described cannot be investigated from the simulations of Temprano-Coleto et al. (2023).

## 5. Discussion

The drag-reducing potential of SHSs may be compromised by trace amounts of surfactant (Peaudecerf et al. 2017). In this paper, we have derived an asymptotic theory for 3D laminar flow, in a plane-periodic channel with SHSs made of an array of long but finite-length longitudinal grooves along both walls, which has been contaminated with soluble surfactant. The mass, momentum and surfactant equations are solved in the Stokes flow limit, where the adsorption–desorption kinetics and equation of state are linearised. We have investigated regimes where cross-channel concentration gradients are small, developing a long-wave theory that accounts for a rapidly equilibrating surfactant-driven transverse flow. This results in a 1D model (3.1)–(3.3a–d) for surfactant transport, which incorporates advection, diffusion, Marangoni effects and exchange between the bulk and
the interface. No parameter fitting is required, in contrast to the 3D theory outlined in Temprano-Coleto et al. (2023) and the 2D theory in Landel et al. (2020). Using this theory, we gain access to parts of the parameter space that are unavailable using models that assume uniform shear stress at the interface (Landel et al. 2020), and we make asymptotic predictions for the drag reduction and surfactant concentration distribution that complement expensive numerical simulations of the 3D flow and surfactant equations (Temprano-Coleto et al. 2023).

We have investigated the leading-order drag reduction ($DR_0$) across the parameter space, varying the strength of surface advection ($\beta$), bulk ($\alpha$) and surface ($\delta$) diffusion, Marangoni effects ($\gamma$) and exchange between the bulk and the interface ($\nu$), compared with bulk advection. When exchange of surfactant is strong, the bulk and interfacial concentration are in equilibrium at leading order. We derived and solved a composite equation (3.30) that includes shear dispersion to qualitatively highlight where 3D effects become important (region G in figure 4a). When exchange of surfactant is moderate, the bulk and interfacial concentrations decouple. We have identified three primary regions of the parameter space in both the strong- and moderate-exchange problems (figures 4 and 7). In the Marangoni-dominated ($M$) regime, the interface is immobilised and the drag reduction vanishes to leading order (low-drag-reduction regime, $DR_0 \ll 1$). The interfacial surfactant distribution is linear with a shallow gradient. In the advection- ($A$) and diffusion-dominated ($D$) regimes, the interface is shear-free and the drag reduction is unaffected by the surfactant at leading order (high-drag-reduction regime, $1 - DR_0 \ll 1$). The interfacial surfactant distribution can be non-uniform in $A$, near the $AM$ boundary, where we find both exponential and almost piecewise-linear stagnant cap profiles. To clarify the underlying physics associated with these results, the dependence of the 3D velocity field on surfactant transport at the bulk and interface has been determined in both strong- and moderate-exchange limits (figures 6 and 9; supplementary movies 1 and 2).

Table 1 summarises asymptotic approximations of the leading-order drag reduction $DR_0$ in regions $M$, $A$ and $D$, and the parts of the parameter space that describe them. Starting in region $M$, where there is no drag reduction and the interface is immobilised, we present approximations of the drag reduction when bulk diffusion ($M_D$), shear dispersion ($M_G$) or surface diffusion ($M_E$) allow small surface mobilisation. The drag can be reduced by strengthening diffusion across the $DM$ boundary (when $\min(\gamma/(\alpha\beta), \gamma/\alpha) \sim O(1)$), strengthening shear dispersion across the $GM$ boundary (when $\alpha\gamma/\epsilon^2 \sim O(1)$) or by reducing the surfactant strength relative to advection across the $AM$ boundary (when $\min(\gamma/\beta, \gamma) \sim O(1)$). The quantity $\gamma/\alpha$, which affects $DR_0$ in regions $M$ and $D$, is identified in Temprano-Coleto et al. (2023) as the mobilisation length. The quantity $\gamma/\beta$, which affects $DR_0$ in regions $A$, is identified in Sundin & Bagheri (2022) to determine whether the surfactant concentration is in the stagnant cap regime or not.

A number of dimensionless ratios appear in table 1 that increase drag: $\gamma/\alpha$, $\alpha\gamma/\epsilon^2$, $\gamma/\beta$, $\gamma/(\alpha\beta)$ and $\gamma/\delta$. We give these ratios in terms of dimensional parameters in table 2. All the ratios given in table 2 have the common factor $\hat{A}\hat{C}\hat{L}_d/\hat{\mu}$, where $\hat{A}$ is the surface activity, $\hat{C}$ is the bulk concentration scale, $\hat{L}_d$ is the depletion length and $\hat{\mu}$ is the dynamic viscosity. Here, $\hat{C}\hat{L}_d$ measures the level of surfactant adsorbed on the plastron; $\hat{A}\hat{C}\hat{L}_d$ gives the corresponding surface tension reduction, making $\hat{A}\hat{C}\hat{L}_d/\hat{\mu}$ a velocity scale generated by interfacial Marangoni effects. The factor $\hat{P}_z\hat{D}$, where $\hat{P}_z$ is the transverse pitch and $\hat{D}$ is the bulk diffusivity, decreases the drag across the $DM$ boundary but increases the drag across the $GM$ boundary; this reflects the smoothing effect of diffusion at the $DM$ boundary and shear dispersion at the $GM$ boundary. When exchange is weak, the surface diffusivity
S.D. Tomlinson and others

A number of assumptions have allowed the present model to be derived. When cross-channel diffusion is weak or the long-wave parameter $\epsilon$ is not sufficiently small, the system lies in region G (identified in § 4) where cross-channel concentration gradients become comparable to streamwise variation. In region G, our 1D model breaks down and the full 3D transport equations must be solved to resolve the flow and surfactant fields. Furthermore, a host of higher-order physical effects associated with flows over SHSs may alter surfactant transport, e.g. interface curvature and the gas subphase, that have been neglected in our model but constitute important extensions (for a detailed discussion of these and other effects, see Lee *et al.* 2016). Another application of SHSs is in the thermal management of electronics, where streamwise and spanwise thermocapillary stresses arise due to temperature gradients at the liquid-gas interface (as in considered in Kirk *et al.* 2020). With minor modifications it is possible that the theory outlined herein could also give insight into these diabatic flows.

To summarise, this paper has highlighted the range of physical balances that arise when one considers the effect of soluble surfactant in laminar flows bounded by SHSs. In the appropriate regimes, our results provide a comprehensive analytical framework that can guide the design of surfactant-contaminated SHSs. Our simple closed-form theoretical predictions of the drag reduction can help minimise drag in realistic 3D SHS microchannels and other applications, where surfactant traces may be naturally present.

**Supplementary movies.** Supplementary movies are available at https://doi.org/10.1017/jfm.2023.264.

**Funding.** We acknowledge support from CBET-EPSRC (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 (CCBY) licence to any Author Accepted Manuscript version arising. F.T.-C. acknowledges support from a distinguished postdoctoral fellowship from the Andlinger Center for Energy and the Environment.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.**

● Samuel D. Tomlinson https://orcid.org/0000-0002-7180-9443;

● Frédéric Gibou https://orcid.org/0000-0001-7022-5262;

● Paolo Luzzatto-Fegiz https://orcid.org/0000-0003-3614-552X;

● Fernando Temprano-Coleto https://orcid.org/0000-0002-2179-3148;

● Oliver E. Jensen https://orcid.org/0000-0003-0172-6578;

● Julien R. Landel https://orcid.org/0000-0003-3159-8749.

**Appendix A. Numerical methods**

A.1. Transport equations

We solve the transport equations (3.1), (3.2) and (3.30), subject to the boundary conditions, (3.3a–d), using Chebyshev collocation. An in-depth discussion of the $(N + 1) \times (N + 1)$ differentiation matrix, $D$, that forms the basis of Chebyshev collocation technique can be found in Trefethen (2000). Briefly, it allows us to approximate $f_{\xi}$ using $Df$, where $f = (f(0), \dots, f(N))^{\mathrm{T}}$ is the solution vector on a grid defined by $\xi(i) = \cos(i\pi/N)$ for $i = 0, 1, \dots, N$. We map domains $D_1$, for $x \in [-\phi_x, \phi_x]$, and $D_2$, for $x \in [\phi_x, 2 - \phi_x]$, to a discrete space $\xi(i) \in [-1, 1]$ for $i = 0, 1, \dots, N$. We solve the transport equations using differentiation matrices and then map the numerical solution back to physical space.

The linearised transport equations can be solved analytically to obtain an initial guess (see, e.g., Appendices C.5 and D.5). We concatenate the solution in both domains as $c_0 = (c_{0,1}(0), \dots, c_{0,1}(N), c_{0,2}(0), \dots, c_{0,2}(N))^{\mathrm{T}}$, where $c_{0,i} = (c_{0,i}(0), \dots, c_{0,i}(N))^{\mathrm{T}}$ is
Surfactant-contaminated superhydrophobic channels

the solution in $D_i$ for $i = 1, 2$. The initial guess evolves to the nonlinear state by solving

$$L_1(c_0^{old})c_{0,1}(1:N-1) = 1 \quad \text{where } c_{0,1}(0) = c_{0,2}(N), \qquad (A1a)$$

$$L_2(c_0^{old})c_{0,2}(1:N-1) = 1 \quad \text{where } c_{0,2}(0) = c_{0,1}(N), \qquad (A1b)$$

where $L_i$ are discrete approximations to (3.30), (3.1) and (3.2) in $D_i$ for $i = 1, 2$. We update the solution via $c_0^{new} = c_0^{old} + r(c_0 - c_0^{old})$ where $r$ is a relaxation factor, until $\tau_i = \|L_i(c_{0,i}^{new} - c_{0,i}^{old})\|_{\infty}/r$ falls below a specified tolerance for $i = 1, 2$.

A.2. Leading-order flow field

Problems (3.15a-g) and (3.16a-g) were solved by modifying the framework introduced in Game et al. (2017). We consider the domain $A_s = \{y \in [0, 1]\} \times \{z \in [0, P_z]\}$, constructing the rest of the solution using symmetry arguments. The domain is decomposed into two parts, separated by $z = \phi_z$, where continuity of the variable and its first derivative is enforced. The domains are transformed to facilitate Chebyshev collocation discretisations and the partial differential equations are transformed to discrete space using techniques outlined in Trefethen (2000). The discontinuous boundary conditions at $z = \phi_z$, (3.15b-e)-(3.16b-e), introduce integrable stress singularities into the problem. The leading-order contributions of these singularities are subtracted to produce less singular problems. The unknown strengths of the singularities are determined by imposing a regularity condition, as follows.

Introducing a local polar coordinate system centred at the contact line $(r, \theta) = ((z - \phi_z)^2 + y^2)^{1/2}$, $\tan^{-1}(y, z - \phi_z))$ and assuming $r \ll 1$, (3.15a-c), (3.16a-c) can be used to evaluate the singular part of the solution ($U^s$) via

$$U_{rr}^s + \frac{U_r^s}{r} + \frac{U_{\theta\theta}^s}{r^2} = b, \quad \text{such that} \quad U^s(r, 0) = 0, \quad U_{\theta\theta}^s(r, \pi) = d, \qquad (A2a-c)$$

where $b$ and $d$ are constants that can be chosen to construct (3.15a-g) and (3.16a-g). The solution must remain bounded as $r \to 0$ and (A2a-c) can be solved to give

$$U^s = B_1^s U_1^s + O(r^{3/2}) = B_1^s r^{1/2} \sin(\theta/2) + O(r^{3/2}), \qquad (A3)$$

where $B_1^s$ is an singularity strength that must be evaluated as part of the solution. We substitute $U = B_1^s U_1^s + U^r$ into (3.15a), (3.16a) where $U^r$ is the residual solution with the singularity in (A3) removed, to obtain

$$B_1^s \nabla_{\perp}^2 U_1^s + \nabla_{\perp}^2 U^r = b. \qquad (A4)$$

The system is completed with a condition which requires regularity in first derivatives,

$$U_z^r(0, \phi_z) \quad \text{is constant,} \qquad (A5)$$

where this constant can be set to zero arbitrarily (Game et al. 2017). Equations (A4), (A5) are combined into a matrix problem which determines the singularity strength as part of the solution. The numerical convergence is improved significantly when compared with a single-domain solution without singularity removal.
## Appendix B. Strong exchange and moderate cross-channel diffusion

### B.1. Derivation of shear dispersion terms

We assume here that $Pe^{-1} = Pe_I^{-1} = Bi = O(\epsilon^2)$ for $\epsilon \ll 1$ and rescale $Pe^{-1} = \epsilon^2 \mathcal{P}^{-1}$, $Pe_I^{-1} = \epsilon^2 \mathcal{P}_I^{-1}$ and $Bi = \epsilon^2 \mathcal{B}$, such that cross-channel diffusion is weak and exchange is strong, with $Da = O(1)$. This scaling is chosen so that we can investigate the transition from weak to strong cross-diffusion problems, i.e. for moderate diffusion, where we first deviate from a well-mixed bulk surfactant concentration. We substitute (3.4) into (2.18a-c)-(2.28) and take the $O(1)$ system. In domains $\mathcal{D}_1$ and $\mathcal{D}_2$,

$$ \nabla_{\perp} \cdot \mathbf{u}_{\perp 0} = 0, \quad \nabla_{\perp}^2 \mathbf{u}_0 - \nabla p_0 = 0, \quad \mathbf{u}_{\perp 0} \cdot \nabla_{\perp} c_0 = 0, \qquad (\text{B1a-c}) $$

on the interface $\mathcal{I}$,

$$ \left. \begin{array}{l} \mathbf{n} \cdot \nabla u_0 - Ma\Gamma_{0x} = 0, \quad \mathbf{n} \cdot \nabla w_0 - Ma\Gamma_{0z} = 0, \quad v_0 = 0, \\ \mathbf{n} \cdot \nabla c_0 - Da(c_0 - \Gamma_0) = 0, \quad (w_0\Gamma_0)_z = 0, \end{array} \right\} \qquad (\text{B2a-e}) $$

on the interface contour $\partial\mathcal{I}$,

$$ w_0\Gamma_0 = 0 \quad \text{at } z = \pm \phi_z P_z, \qquad (\text{B3}) $$

between domains $\mathcal{D}_1$ and $\mathcal{D}_2$ we have (3.8a,b), on the ridge $\mathcal{R}$ and solid $\mathcal{S}$ we have (3.11a-d), and at the ends of the transverse period we have (3.12a-d). The total fluid flux is given by (3.9) and total surfactant flux is given by

$$ \int_{z=-P_z}^{P_z} \int_{y=0}^{2} u_0 c_0 \, dy \, dz + \frac{2Da}{\mathcal{B}\mathcal{P}} \int_{z=-P_z}^{P_z} u_0 \Gamma_0 \, dz = 1, \qquad (\text{B4}) $$

and the drag reduction becomes (3.13).

Similar to § 3.2.1, from (B1a-c)-(B4), (3.8a,b)-(3.9), (3.11a-d)-(3.12a-d) and (B4), we have that $\Gamma_0 = \Gamma_0(x)$ where $\Gamma_0 = c_0(x, 0, z) = c_0(x, 2, z)$, $p_0 = p_0(x)$ and $v_0 = w_0 = 0$. The streamwise velocity field is given by $u_0 = \tilde{U}p_{0x} + Ma\tilde{U}\Gamma_{0x}$ in $\mathcal{D}_1$ and $u_0 = \tilde{U}p_{0x}$ in $\mathcal{D}_2$, where $\tilde{U}$, $\tilde{U}$ and $\tilde{U}$ are given by (3.15a-g)-(3.17a-c). Substituting $u_0$ into (3.9), we recover the velocity flux constraints satisfying (3.25a-c), where $\tilde{Q}$, $\tilde{Q}$, $\tilde{q}$ and $\tilde{q}$ are defined in (3.18a,b), (3.19a-c).

We use the $O(\epsilon^2)$ bulk surfactant equation to determine $c_0$. In $\mathcal{D}_1$ and $\mathcal{D}_2$,

$$ \mathcal{P}^{-1} \nabla_{\perp}^2 c_0 - u_0 c_{0x} - \mathbf{u}_{\perp 1} \cdot \nabla_{\perp} c_0 = 0. \qquad (\text{B5}) $$

The third term in (B5) involves the first-order velocity components $v_1$ and $w_1$. However, we shortly assume that cross-channel gradients are small, which means that $v_1$ and $w_1$ are not required to evaluate $c_0$.

Decomposing the bulk surfactant field into a cross-channel average and residual component, we write

$$ c_0(x, y, z) = \langle c_0 \rangle(x) + c'_0(x, y, z), \qquad (\text{B6}) $$

where $\langle c'_0 \rangle = 0$. We can then evaluate the shear-dispersion contributions

$$ \langle \tilde{U}c_0 \rangle = \tilde{Q}\langle c_0 \rangle + \langle \tilde{U}c'_0 \rangle, \quad \langle \tilde{U}c_0 \rangle = \tilde{Q}\langle c_0 \rangle + \langle \tilde{U}c'_0 \rangle, \quad \langle \tilde{U}c_0 \rangle = \tilde{Q}\langle c_0 \rangle + \langle \tilde{U}c'_0 \rangle. \qquad (\text{B7a-c}) $$
Substituting (B7a-c) into (B4) and using (3.14a,b), the surfactant flux constraints become

$$ \left( \frac{2Da}{\mathcal{B}\mathcal{P}} \tilde{q}\Gamma_0 + \tilde{Q}\langle c_0 \rangle + 4P_z \langle \tilde{U}c'_0 \rangle \right) p_{0x} + \left( \frac{2Da}{\mathcal{B}\mathcal{P}} \tilde{q}\Gamma_0 + \tilde{Q}\langle c_0 \rangle + 4P_z \langle \tilde{U}c'_0 \rangle \right) Ma\Gamma_{0x} = 1 \quad \text{in } \mathcal{D}_1, \qquad (\text{B8a}) $$

$$ (\tilde{Q}\langle c_0 \rangle + 4P_z \langle \tilde{U}c'_0 \rangle) p_{0x} = 1 \quad \text{in } \mathcal{D}_2. \qquad (\text{B8b}) $$

In this approach, we are not accounting for variations in $\Gamma_0$ driven by $c'_0$. This will bring a level of refinement that we will address elsewhere.

Substituting (B6) into (B5), (B2d), (3.11d), (3.12d) and assuming $|c'_0| \ll \langle c_0 \rangle$, i.e. the magnitude of the cross-channel concentration gradients is small, we have

$$ \left. \begin{array}{l} \mathcal{P}^{-1} \nabla_{\perp}^2 c'_0 = u_0 \langle c_0 \rangle_x, \quad \text{subject to} \quad \mathcal{D}^+ c'_0(0, z_s) = 0, \quad c'_{0y}(0, z_{ns}) = 0, \\ \mathcal{D}^- c'_0(2, z_s) = 0, \quad c'_{0y}(2, z_{ns}) = 0, \quad c'_{0z}(y, -1) = 0, \quad c'_{0z}(y, 1) = 0 \quad \text{in } \mathcal{D}_1, \end{array} \right\} (\text{B9a-g}) $$

where $\mathcal{D}^{\pm} = \pm \partial_y - Da$ and

$$ \left. \begin{array}{l} \mathcal{P}^{-1} \nabla_{\perp}^2 c'_0 = u_0 \langle c_0 \rangle_x, \quad \text{subject to} \quad c'_{0y}(0, z) = 0, \quad c'_{0y}(2, z) = 0, \\ c'_{0z}(y, -1) = 0, \quad c'_{0z}(y, 1) = 0 \quad \text{in } \mathcal{D}_2. \end{array} \right\} (\text{B10a-e}) $$

Using superposition, we may then write

$$ c'_0 = \mathcal{P} \tilde{C} p_{0x} \langle c_0 \rangle_x + \mathcal{P} Ma \tilde{C} \Gamma_{0x} \langle c_0 \rangle_x \quad \text{in } \mathcal{D}_1, \quad c'_0 = \mathcal{P} \tilde{C} p_{0x} \langle c_0 \rangle_x \quad \text{in } \mathcal{D}_2. \qquad (\text{B11a,b}) $$

Substituting (B11a,b) into (B9a-g), (B10a-e), in order to obtain $c'_0$ we require the solution to the following three boundary-value problems: shear dispersion due to flow driven by a pressure gradient over the plastron,

$$ \left. \begin{array}{l} \nabla_{\perp}^2 \tilde{C} = \tilde{U}, \quad \text{subject to} \quad \mathcal{D}^+ \tilde{C}(0, z_s) = 0, \quad \tilde{C}_y(0, z_{ns}) = 0, \\ \mathcal{D}^- \tilde{C}(2, z_s) = 0, \quad \tilde{C}_y(2, z_{ns}) = 0, \quad \tilde{C}_z(y, -1) = 0, \quad \tilde{C}_z(y, 1) = 0; \end{array} \right\} (\text{B12a-g}) $$

shear dispersion due to flow driven by surface shear,

$$ \left. \begin{array}{l} \nabla_{\perp}^2 \tilde{C} = \tilde{U}, \quad \text{subject to} \quad \mathcal{D}^+ \tilde{C}(0, z_s) = 0, \quad \tilde{C}_y(0, z_{ns}) = 0, \\ \mathcal{D}^- \tilde{C}(2, z_s) = 0, \quad \tilde{C}_y(2, z_{ns}) = 0, \quad \tilde{C}_z(y, -1) = 0, \quad \tilde{C}_z(y, 1) = 0; \end{array} \right\} (\text{B13a-g}) $$

and dispersion due to the flow in $\mathcal{D}_2$,

$$ \tilde{C}_{yy} = \tilde{U}, \quad \text{subject to} \quad \tilde{C}_y(0) = 0, \quad \int_{y=0}^{2} \tilde{C} dy = 0. \qquad (\text{B14a-c}) $$

The solutions to (B12a-g) and (B13a-g) are found numerically, following the procedure described in Appendix A. Equation (B14a-c) can be integrated directly to give $\tilde{C} = 1/5 - y^3/6 + y^4/24$. The shear dispersion contributions become

$$ \left. \begin{array}{l} \langle \tilde{U}c'_0 \rangle = \mathcal{P} \langle \tilde{U}\tilde{C} \rangle p_{0x} \langle c_0 \rangle_x + \mathcal{P} \langle \tilde{U}\tilde{C} \rangle Ma \langle c_0 \rangle_x^2, \\ \langle \tilde{U}c'_0 \rangle = \mathcal{P} \langle \tilde{U}\tilde{C} \rangle p_{0x} \langle c_0 \rangle_x + \mathcal{P} \langle \tilde{U}\tilde{C} \rangle Ma \langle c_0 \rangle_x^2, \quad \langle \tilde{U}c'_0 \rangle = \mathcal{P} \langle \tilde{U}\tilde{C} \rangle p_{0x} \langle c_0 \rangle_x. \end{array} \right\} (\text{B15a-c}) $$

Substituting (B15a-c) into (B8) with $p_{0x}$ from (3.25a) in $\mathcal{D}_1$ and (3.25c) in $\mathcal{D}_2$, yields the steady integrated surfactant transport equations and boundary conditions (3.3a), (3.3b)
<figure><img src="image_14.png" /><figcaption>Figure 11. Plots of (a) $s_1$ defined in (3.29a), (b) $s_2$ defined in (3.29b) and (c) $s_3$ defined in (3.29c), the coefficients multiplying the shear dispersion terms in the moderate cross-channel diffusion and strong-exchange problem, for varying $\phi_z$ and $Da$ with $P_z = 1$.</figcaption></figure>

and (3.28), with $\langle c_0 \rangle$ replaced by $c_0$ and the coefficients $s_1, s_2, s_3$ and $s_4$ defined by (3.29a-d). The coefficients $s_1, s_2$ and $s_3$ are plotted as functions of $\phi_z$ and $Da$ in figure 11; the coefficient $s_4 = 3/35$ for all $\phi_z$ and $Da$, where $P_z = 1$.

## B.2. Validity of the shear dispersion approximation

When bulk diffusion is not sufficiently strong, concentration gradients normal to the SHS become comparable to the streamwise variation, the long-wave theory outlined in § 3 breaks down and 3D numerical simulation of the governing equations is required in order to resolve the coupled flow and surfactant fields. Here, we compare the size of the cross-channel-averaged surfactant concentration $\langle c_0 \rangle$ to the cross-channel variation $c'_0$, introduced in Appendix B.1, assuming that gradients normal to the SHS become important when $c'_0 = O(\langle c_0 \rangle)$. This analysis gives an approximate boundary for the range of validity of our model.

We assume throughout this section that $\beta = O(1)$ (the expressions for $\beta \ll 1$ and $\beta \gg 1$ can be derived using similar arguments). To calculate the boundary between regions $M$ and $G$, expand the cross-channel-averaged surfactant concentration as follows (see Appendix C.1):

$$ \langle c_0 \rangle = 1 + \frac{\beta}{\gamma} \left( x - \frac{\phi_x(E+1)}{(E-1)} \right) + \dots, \quad \text{for } \gamma \gg \max \left( 1, \alpha, \delta, \frac{\epsilon^2}{\alpha} \right). \quad (B16) $$

Substituting (B16) into (B11a,b), using $p_{0x} = 1/\tilde{Q} - Ma\tilde{Q}\Gamma_{0x}/\tilde{Q}$ from (3.25a-c) and noting that $\Gamma_{0x} = \langle c_0 \rangle_x = O(1/\gamma)$ in the strong-exchange regime, we find that $c'_0 = O(\mathcal{P}/\gamma)$ because $Ma = O(\gamma)$. As $\mathcal{P} = O(\epsilon^2/\alpha)$, we have $c'_0 = O(\epsilon^2/(\alpha\gamma))$, so that $c'_0 = O(\langle c_0 \rangle) = O(1)$ when $\gamma = O(\epsilon^2/\alpha)$ in region $M$. Hence, $\gamma = O(\epsilon^2/\alpha)$ defines the boundary of validity of region $M$ near the region $G$, where shear dispersion terms become important (see figure 4).

Similarly, to calculate the boundary between regions $A$ and $G$, expand as follows (see Appendix C.5):

$$ \langle c_0 \rangle = \left. \frac{1}{\beta+1} + \frac{\beta}{(\beta+1)} \exp \left( \frac{(1+\beta)(x-\phi_x)}{\alpha + \delta + \epsilon^2 s_1 / \alpha} \right) + \dots, \right\} \quad (B17) \\ \text{for } \epsilon^2 \ll \alpha \ll 1, \quad \gamma \ll \min(1, \beta). $$

Substituting (B17) into (B11a,b) and using $p_{0x} = 1/\tilde{Q} - Ma\tilde{Q}\Gamma_{0x}/\tilde{Q}$, we find that $c'_0 = O(\mathcal{P})$. As $\mathcal{P} = O(\epsilon^2/\alpha)$, we have that $c'_0 = O(\langle c_0 \rangle)$ when $\alpha = O(\epsilon^2)$ in region $A$.
Appendix C. Asymptotic solutions for strong exchange

C.1. *Strong Marangoni effect: region M*

Assuming that $\beta = O(1)$ and $\gamma \gg \max(1, \alpha, \delta, \epsilon^2/\alpha)$, but retaining the effects of diffusion and shear dispersion, we expand the concentration field from § 3.2 using $c_0 = c_{00} + c_{01}/\gamma + c_{02}/\gamma^2 + \cdots$ in (3.3a), (3.3b), (3.30). At $O(\gamma^2)$, we have

$$
\left.
\begin{array}{l}
c_{00x}^3 = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00} - \alpha c_{00x} - \frac{\epsilon^2 s_4}{\alpha} c_{00x} = 1 \quad \text{in } \mathcal{D}_2, \\
\text{subject to } c_{00}(\phi_x^-) = c_{00}(\phi_x^+), \quad c_{00}(-\phi_x) = c_{00}(2-\phi_x).
\end{array}
\right\}
\quad \text{(C1a-d)}
$$

The system in (C1a-d) requires $c_{00} = 1$ in $D_1 \cup D_2$. At $O(1)$ and $O(\gamma)$, we have

$$
\left.
\begin{array}{l}
c_{01x} = \beta \quad \text{in } \mathcal{D}_1, \quad c_{01} - \alpha c_{01x} - \frac{\epsilon^2 s_4}{\alpha} c_{01x} = 0 \quad \text{in } \mathcal{D}_2, \\
\text{subject to } c_{01}(\phi_x^-) = c_{01}(\phi_x^+), \quad c_{01}(-\phi_x) = c_{01}(2-\phi_x).
\end{array}
\right\}
\quad \text{(C2a-d)}
$$

We can integrate (C2a-d) to show that $c_{01} = \beta(x - \phi_x(E + 1)/(E - 1))$ in $D_1$, where $E = \exp(2\alpha(1 - \phi_x)/(\alpha^2 + \epsilon^2 s_4))$, giving (4.2). At $O(\gamma^{-1})$, we have that

$$
c_{02x} = (\beta + 1)c_{01} - c_{01}c_{01x} - (\alpha + \delta)c_{01x} - \frac{\epsilon^2}{\alpha} \left( s_1 c_{01x} + \frac{s_2}{\beta} c_{01x}^2 + \frac{s_3}{\beta^2} c_{01x}^3 \right) \quad \text{in } \mathcal{D}_1. \tag{C3}
$$

From (C2a-d) and (C3) we can compute the concentration increase over $D_1$,

$$
\Delta c_0 = \frac{2\phi_x\beta}{\gamma} + \frac{2\phi_x\beta}{\alpha\gamma^2(1-E)} \left( \alpha(\phi_x - \alpha - \delta) - \epsilon^2 s + (\alpha(\phi_x + \alpha + \delta) + \epsilon^2 s)E \right), \quad (C4)
$$

where $s \equiv s_1 + s_2 + s_3 > 0$, which is substituted into (3.38) to give (4.1). The corrections to $c_0 \approx c_{00} = 1$ and $DR_0 = 0$ are small provided $\gamma \gg \max(1, \alpha, \delta, \epsilon^2/\alpha)$, defining the boundaries of the asymptotic region $M$.

C.2. *Strong Marangoni effect and strong shear dispersion: the GM boundary*

At the *GM* boundary, assume that $\beta = O(1)$ and $\gamma \gg 1$. Rescale $\epsilon^2/\alpha = \gamma/a$ where $a = O(1)$. Expand the concentration field from § 3.2, using $c_0 = c_{00} + c_{01}/\gamma + \cdots$ in (3.3a), (3.3b) and (3.30). At $O(\gamma^3)$, shear dispersion dominates with

$$
c_{00x}^3 = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00x} = 0 \quad \text{in } \mathcal{D}_2, \qquad \text{(C5a,b)}
$$

subject to (C1c,d). In order to calculate $c_{00}$ we must proceed to the next order. At $O(1)-O(\gamma^2)$, Marangoni effects and advection enter

$$
\left.
\begin{array}{l}
(\beta + 1)c_{00} - c_{00}c_{01x} - \frac{1}{a} \left( s_1 c_{01x} + \frac{s_2}{\beta} c_{01x}^2 + \frac{s_3}{\beta^2} c_{01x}^3 \right) = 1 \quad \text{in } \mathcal{D}_1, \\
\\
c_{00} - \frac{s_4}{a} c_{01x} = 0 \quad \text{in } \mathcal{D}_2,
\end{array}
\right\}
\quad \text{(C6a,b)}
$$

subject to (C1c,d). As $c_{01x} = \beta$ in $M$ (see Appendix C.1), we expect $s_2 c_{01x}^2 / (a\beta) \approx s_2 c_{01x} / a$ and $s_3 c_{01x}^3 / (a\beta^2) \approx s_3 c_{01x} / a$ as we approach the GM boundary, such that we
instead solve the linearised problem

$$ (\beta + 1)c_{00} - c_{00}c_{01x} + \frac{s}{a}c_{01x} \approx 1 \quad \text{in } \mathcal{D}_1, \quad c_{00} - \frac{s_4}{a}c_{01x} = 1 \quad \text{in } \mathcal{D}_2. \quad (C7a,b) $$

Integrating (C7a,b) over the period gives $c_{00}$ as the solution to

$$ (c_{00} - 1)(\phi_x - 1)(s - ac_{00}) = s_4\phi_x(1 - c_{00}(\beta + 1)). \quad (C8) $$

By solving (C8), we can then integrate (C7a) in $\mathcal{D}_1$ to find $\Delta c_0$, which is substituted into (3.38) to give (4.3).

### C.3. Strong diffusion and strong Marangoni effect: the DM boundary

At the DM boundary, assume that $\beta = O(1)$ and $\gamma \gg 1$. Rescale $\alpha = a\gamma$ and $\delta = d\gamma$ where $a \sim d \sim O(1)$. Expand the concentration field using $c_0 = c_{00} + c_{01}/\gamma + \cdots$ in (3.3a), (3.3b), (3.30). At $O(\gamma)$, Marangoni effects, diffusion and shear dispersion dominate with

$$ -c_{00}c_{00x} - (a+d)c_{00x} - \frac{\epsilon^2 s_3}{a\beta^2}c_{00x}^3 = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00x} = 0 \quad \text{in } \mathcal{D}_2, \quad (C9a,b) $$

subject to (C1c,d). Hence, $c_{00x} = 0$ in both $\mathcal{D}_1$ and $\mathcal{D}_2$. In order to calculate $c_{00}$ we must proceed to the next order. At $O(1)$, the advection terms enter

$$ (\beta + 1)c_{00} - c_{00}c_{01x} - (a+d)c_{01x} = 1 \quad \text{in } \mathcal{D}_1, \quad c_{00} - ac_{01x} = 1 \quad \text{in } \mathcal{D}_2, \quad (C10a,b) $$

subject to (C2c,d). Integrating (C10a,b) over $\mathcal{D}_1$ gives $c_{00}$ as the solution to

$$ (c_{00} - 1)(\phi_x - 1)(a + d + c_{00}) = a\phi_x(c_{00}(\beta + 1) - 1). \quad (C11) $$

By solving (C11), we can then integrate (C10a) in $\mathcal{D}_1$ to find $\Delta c_0$, which is substituted into (3.38) to give (4.4).

### C.4. Strong diffusion: region D

Assume that $\beta = O(1)$ and $\min(\alpha, \delta) \gg \max(1, \gamma)$. Let $\delta = d\alpha$, where $d = O(1)$. Expand the concentration field from § 3.2 using $c_0 = c_{00} + c_{01}/\alpha + \cdots$ in (3.3a), (3.3b), (3.30). At $O(\alpha)$, diffusion dominates with

$$ c_{00x} = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00x} = 0 \quad \text{in } \mathcal{D}_2, \quad (C12a,b) $$

subject to (C1c,d). In order to calculate $c_{00}$ we must proceed to the next order. At $O(1)$, advection enters

$$ (\beta + 1)c_{00} - (1 + d)c_{01x} = 1 \quad \text{in } \mathcal{D}_1, \quad c_{00} - c_{01x} = 1 \quad \text{in } \mathcal{D}_2, \quad (C13a,b) $$

subject to (C2c,d). Integrating (C13a,b) over $\mathcal{D}_1 \cup \mathcal{D}_2$ and imposing periodicity, $c_{01x}(-\phi_x) = c_{01}(2 - \phi_x)$, $c_{00}$ is given in (4.6). Using $c_{01x}$, we can obtain the jump in $c_0$ over $\mathcal{D}_1$,

$$ \Delta c_0 \approx \frac{2\phi_x}{\alpha}c_{01x} = \frac{2\phi_x(1 - \phi_x)\beta}{(1 - \phi_x)(\alpha + \delta) + \phi_x\alpha(\beta + 1)}, \quad (C14) $$

which is substituted into (3.38) to give (4.5). The corrections to $c_0 = c_{00}$ and $DR_0 = 1$ are small provided $\min(\alpha, \delta) \gg \max(1, \gamma, \beta)$, defining the boundaries of the asymptotic region D.
### C.5. Weak Marangoni effect: region A

Assume that $\beta = O(1)$ and $\gamma \ll 1$. Expand the concentration field from §3.2 using $c_0 = c_{00} + \gamma c_{01} + \cdots$ in (3.3a), (3.3b), (3.30). At $O(1)$, advection, diffusion and shear dispersion dominate with

$$ (\beta + 1)c_{00} - \left(\alpha + \delta + \frac{\epsilon^2 s_1}{\alpha}\right) c_{00x} = 1 \quad \text{in } \mathcal{D}_1, \quad c_{00} - \left(\alpha + \frac{\epsilon^2 s_4}{\alpha}\right) c_{00x} = 1 \quad \text{in } \mathcal{D}_2, \qquad (\text{C15a,b}) $$

subject to (C1c,d). Assuming $\epsilon^2 \ll \alpha \ll 1$, we integrate (C15a,b) to get (4.8). In either case, we have that

$$ \Delta c_0 \approx \frac{\beta}{\beta + 1} \qquad (\text{C16}) $$

which is substituted into (3.38) to give (4.7). The corrections to $c_0 = c_{00}$ and $DR_0 = 1$ are small provided $\epsilon^2 \ll \alpha \ll 1$ and $\gamma \ll 1$, defining the boundaries of the asymptotic region A. Note that (C15a,b) is linear and has a general solution that can be used to determine the AD boundary. However, as this does not change the drag reduction to leading order (the interface remains shear free for $\gamma \ll \min(1, \beta)$), we do not investigate this limit here.

### C.6. Strong advection and strong Marangoni effect: the AM boundary

Assume that $\beta = O(\gamma)$, $\gamma \gg \max(\alpha, \delta, \epsilon^2/\alpha)$ and $\epsilon^2 \ll \alpha \ll 1$. Rescale $\alpha = a/\gamma$, $\beta = b\gamma$, $\delta = d/\gamma$ and $\epsilon^2/\alpha = e/\gamma$ where $a, b, d$ and $e$ are positive $O(1)$ constants. Expand the concentration field from §3.2, using $c_0 = c_{00} + c_{01}/\gamma + \cdots$ in (3.3a), (3.3b) and (3.30). At $O(\gamma)$, Marangoni effects and advection dominate with

$$ c_{00}(b - c_{00x}) = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00} = 1 \quad \text{in } \mathcal{D}_2, \qquad (\text{C17a,b}) $$

subject to (C1c,d). The system in (C17a,b) gives a linear profile $c_{00} = \beta(x - \phi_x)/\gamma + 1$ for $\gamma/\beta > 2\phi_x$ in $\mathcal{D}_1$, such that $c_{00} \ge 0$ for all $x \in [-\phi_x, \phi_x]$. For $\gamma/\beta > 2\phi_x$, we find that $\Delta c_0 = 2\phi_x\beta/\gamma$ and $DR_0 = 0$ at leading order, hence we proceed to $O(1)$ where we find that $\Delta c_0 = 2\phi_x\beta/\gamma + 2\phi_x/\gamma + \ln(1 - 2\beta\phi_x/\gamma)/\beta$ provided $1 - 2\phi_x\beta/\gamma \gg \exp(-\beta)$. For $\gamma/\beta \le 2\phi_x$ a piecewise-linear solution exists with $c_{00} = 0$ for all $-\phi_x \le x \le \phi_x - \gamma/\beta$ and $c_{00} = \beta(x - \phi_x)/\gamma + 1$ for all $\phi_x - \gamma/\beta \le x \le \phi_x$. This nonlinear profile typically represents the emergence of the stagnant cap profile (He et al. 1991), where surfactant is swept to the downstream end of the plastron where it gives rise to a strong gradient. For $\gamma/\beta \le 2\phi_x$, we find that $\Delta c_0 = 1$, which is substituted into (3.38) to give (4.9).

## Appendix D. Asymptotic solutions for weak exchange

### D.1. Strong Marangoni effect: region $M_E$

Assume that $\beta = O(1)$, $\gamma \gg \max(1, \alpha, \delta)$ and $\nu/\epsilon^2 \ll \min(1, \alpha, \delta)$. Rescale $\nu/\epsilon^2 = n/\gamma$ where $n = O(1)$. Expand the concentration field from §3.3 using $\Gamma_0 = \Gamma_{00} + \Gamma_{01}/\gamma + \Gamma_{02}/\gamma^2 + \cdots$ and $c_0 = c_{00} + c_{01}/\gamma + \cdots$ into (3.1)–(3.3a–d). At $O(\gamma)$, Marangoni
effects, advection and diffusion dominate with

$$
\left.
\begin{array}{l}
(c_{00} - \alpha c_{00x})_x = 0, \quad \Gamma_{00}\Gamma_{00x} = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00} - \alpha c_{00x} = 1 \quad \text{in } \mathcal{D}_2, \\
\text{subject to} \quad c_{00}(\phi_x^-) = c_{00}(\phi_x^+), \quad c_{00}(-\phi_x) = c_{00}(2-\phi_x), \\
\qquad c_{00}(\pm\phi_x) - \alpha c_{00x}(\pm\phi_x) = 1, \quad \Gamma_{00}(\pm\phi_x)\Gamma_{00x}(\pm\phi_x) = 0.
\end{array}
\right\} \text{(D1a-g)}
$$

The system in (D1a-g) requires $\Gamma_{00}$ constant on $\mathcal{I}$. At $O(1)$, surface advection enters

$$
\left.
\begin{array}{l}
(c_{01} - \alpha c_{01x})_x = 0, \quad \beta - \Gamma_{01x} = 0 \quad \text{in } \mathcal{D}_1, \quad c_{01} - \alpha c_{01x} = 0 \quad \text{in } \mathcal{D}_2, \\
\text{subject to} \quad c_{01}(\phi_x^-) = c_{01}(\phi_x^+), \quad c_{01}(-\phi_x) = c_{01}(2-\phi_x), \\
\qquad c_{01}(\pm\phi_x) - \alpha c_{01x}(\pm\phi_x) = 0, \quad \beta - \Gamma_{01x}(\pm\phi_x) = 0.
\end{array}
\right\} \text{(D2a-g)}
$$

We can integrate (D2a-g) to show that $\Gamma_{01} = \beta x + C_1$ on $\mathcal{I}$, using the net flux condition $\int_{x=-\phi_x}^{\phi_x} (\Gamma_{01} - c_{01}) dx = 0$ from (2.26) to find $C_1 = 0$, which gives (4.12a,b). At $O(\gamma^{-1})$, we have

$$
\Gamma_{02x} = \beta \Gamma_{01} - \Gamma_{01} \Gamma_{01x} - \delta \Gamma_{01x} + c_{01} - \alpha c_{01x} \quad \text{on } \mathcal{I}. \qquad (\text{D3})
$$

Integrating (D2a-g) and (D3) on $\mathcal{I}$, altogether we have that $\Delta\Gamma_0 = 2\phi_x\beta/\gamma - 2\phi_x\beta\delta/\gamma^2$, which is substituted into (3.38) to give (4.12a,b). The corrections to $\Gamma_0 = \Gamma_{00} = 1$ and $DR_0 = 0$ are small provided $\gamma \gg \max(1, \alpha, \delta)$ and $\nu/\epsilon^2 \ll \min(1, \alpha, \delta)$, defining the boundaries of the asymptotic region $M_E$.

D.2. *Strong Marangoni effect and strong diffusion: the M<sub>E</sub>D<sub>E</sub> boundary*

At the $M_E D_E$ boundary, assume that $\beta = O(1)$ and $\gamma \gg 1$. Rescale $\alpha = a\gamma$, $\delta = d\gamma$ and $\nu/\epsilon^2 = n/\gamma$ where $a, d$ and $n$ are positive $O(1)$ constants. Expand the concentration field from § 3.3 using $\Gamma_0 = \Gamma_{00} + \Gamma_{01}/\gamma + \dots$ and $c_0 = c_{00} + c_{01}/\gamma \dots$ into (3.1)-(3.3a-d). At $O(\gamma)$, Marangoni effects and diffusion dominate with

$$
\left.
\begin{array}{l}
c_{00xx} = 0, \quad ac_{00x} + \Gamma_{00}\Gamma_{00x} + d\Gamma_{00x} = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00x} = 0 \quad \text{in } \mathcal{D}_2, \\
\text{subject to} \quad c_{00}(\phi_x^-) = c_{00}(\phi_x^+), \quad c_{00}(-\phi_x) = c_{00}(2-\phi_x), \\
\qquad c_{00x}(\pm\phi_x) = 0, \quad \Gamma_{00}(\pm\phi_x)\Gamma_{00x}(\pm\phi_x) + d\Gamma_{00x}(\pm\phi_x) = 0.
\end{array}
\right\} \text{(D4a-g)}
$$

The system in (D4a-g) requires $\Gamma_{00}$ and $c_{00}$ are constant. At $O(1)$, advection appears with

$$
\left.
\begin{array}{l}
c_{01xx} = 0, \quad c_{00} - ac_{01x} + \beta \Gamma_{00} - (\Gamma_{00} + d) \Gamma_{01x} = 1 \quad \text{in } \mathcal{D}_1, \quad c_{01x} = 0 \quad \text{in } \mathcal{D}_2, \\
\text{subject to} \quad c_{01}(\phi_x^-) = c_{01}(\phi_x^+), \quad c_{01}(-\phi_x) = c_{01}(2-\phi_x), \\
\qquad c_{00} - ac_{01x}(\pm\phi_x) = 1, \quad \beta \Gamma_{00} - (\Gamma_{00} + d) \Gamma_{01x}(\pm\phi_x) = 0.
\end{array}
\right\} \text{(D5a-g)}
$$

The system in (D5a-g) requires that $c_{00} = 1$ and $\Gamma_{01x} = (\beta\Gamma_{00})/(d+\Gamma_{00})$. Using information at $O(1/\gamma)$, we can show that $\Gamma_{00} = 1$. We can then evaluate $\Delta\Gamma_0 \approx 2\phi_x\beta/(\gamma(d+1))$ and substitute it into (3.38) to find $DR_0 \approx 0.5$ when $\gamma = \delta$, as in (4.13).

D.3. *Strong Marangoni effect, advection and diffusion: the M<sub>E</sub>AD<sub>E</sub> boundary*

At the $M_E A D_E$ boundary, assume that $\beta = O(1)$ and $\nu/\epsilon^2 \ll 1$. Expand the concentration field from § 3.3 using $\Gamma_0 = \Gamma_{00} + \epsilon^2\Gamma_{01}/\nu + \dots$ and $c_0 = c_{00} + \epsilon^2 c_{01}/\nu \dots$ into
(3.1)–(3.3a–d). At O(1), advection, diffusion and Marangoni effects dominate with

$$
\left.
\begin{array}{l}
(c_{00} - \alpha c_{00x})_x = 0, \quad \beta \Gamma_{00} - \gamma \Gamma_{00} \Gamma_{00x} - \delta \Gamma_{00x} = 0 \quad \text{in } \mathcal{D}_1, \\
c_{00} - \alpha c_{00x} = 1 \quad \text{in } \mathcal{D}_2, \quad \text{subject to} \quad c_{00}(\phi_x^-) = c_{00}(\phi_x^+), \\
c_{00}(-\phi_x) = c_{00}(2-\phi_x), \quad c_{00}(\pm\phi_x) - \alpha c_{00x}(\pm\phi_x) = 1, \\
\beta \Gamma_{00}(\pm\phi_x) - \gamma \Gamma_{00}(\pm\phi_x) \Gamma_{00x}(\pm\phi_x) - \delta \Gamma_{00x}(\pm\phi_x) = 0.
\end{array}
\right\}
\quad (\text{D6a-g})
$$

We can integrate (D6a–g) directly over $\mathcal{I}$ using the no-net-flux condition $\int_{x=-\phi_x}^{\phi_x} \Gamma_{00} dx = 2\phi_x$ as $c_{00} = 1$ to derive

$$2\phi_x\beta - \gamma\Delta\Gamma_{00}(\Gamma_{00}(-\phi_x) + \Gamma_{00}(\phi_x))/2 = \delta\Delta\Gamma_{00}. \quad (D7)$$

From (D7) we can evaluate $\Delta\Gamma_0$ and substitute it into (3.38) to obtain

$$DR_0 = \frac{\delta + \gamma \left( \frac{\Gamma_{00}(-\phi_x) + \Gamma_{00}(\phi_x)}{2} - 1 \right)}{\delta + \gamma \left( \frac{\Gamma_{00}(-\phi_x) + \Gamma_{00}(\phi_x)}{2} \right)}. \quad (D8)$$

As long as $(\Gamma_{00}(-\phi_x) + \Gamma_{00}(\phi_x))/2 \approx 1$ (as in figure 7d), then $DR_0 \approx 1 - \gamma/(\delta + \gamma)$ and $DR_0 = 0.5$ when $\gamma = \delta$. Note that we recover (4.11a,b) for $\gamma \gg \delta$ in $M_E$ and (4.15) for $\gamma \ll \delta$ in $A_E$, $D_E$ and $AD_E$.

### D.4. Strong Marangoni effect and strong advection: the $M_E A_E$ boundary

At the $M_E A_E$ boundary, assume that $\beta = O(1)$ and $\alpha \ll 1$. Rescale, $\delta = d\alpha$ and $\nu/\epsilon^2 = n\alpha$ where $d$ and $n$ are $O(1)$ constants. Expand the concentration field from § 3.3 using $\Gamma_0 = \Gamma_{00} + \alpha\Gamma_{01} + \dots$ and $c_0 = c_{00} + \alpha c_{01} + \dots$ into (3.1)–(3.3a–d). At $O(1)$, advection and Marangoni effects dominate with

$$
\left.
\begin{array}{l}
c_{00x} = 0, \quad \beta\Gamma_{00} - \gamma\Gamma_{00}\Gamma_{00x} = 0 \quad \text{in } \mathcal{D}_1, \quad c_{00} = 1 \quad \text{in } \mathcal{D}_2, \\
\quad \text{subject to} \quad c_{00}(\phi_x^-) = c_{00}(\phi_x^+), \quad c_{00}(-\phi_x) = c_{00}(2-\phi_x), \\
\quad c_{00}(\pm\phi_x) = 1, \quad \beta\Gamma_{00}(\pm\phi_x) - \gamma\Gamma_{00}(\pm\phi_x)\Gamma_{00x}(\pm\phi_x) = 0.
\end{array}
\right\}
\quad (\text{D9a-g})
$$

The system in (D9a–g) requires $\Gamma_{00} = 0$ for $-\phi_x \le x \le x_0$ and $\Gamma_{00} = \beta(x-x_0)/\gamma$ for $x_0 \le x \le \phi_x$, for $\gamma \le \beta\phi_x$. This nonlinear profile typically denotes the presence of the stagnant cap profile (He et al. 1991). The constant $x_0 = \phi_x - 2(\phi_x\gamma/\beta)^{1/2}$ is evaluated by the no-net-flux condition $\int_{x=-\phi_x}^{\phi_x} \Gamma_{00} dx = 2\phi_x$ as $c_{00} = 1$. We then evaluate $\Delta\Gamma_0 = 2(\phi_x\beta/\gamma)^{1/2}$ and substitute it into (3.38) to find $DR_0 \approx 0.5$ when $\gamma = \beta\phi_x/4$, as in (4.14).

To extend the region of validity of the above solution, we retain the surface diffusion term at leading order, such that $c_{00} = 1$ in $\mathcal{D}_1 \cup \mathcal{D}_2$ but now (D9b,g) become

$$\left.
\begin{array}{l}
\beta \Gamma_{00} - \gamma \Gamma_{00} \Gamma_{00x} - \delta \Gamma_{00x} = 0 \quad \text{in } \mathcal{D}_1, \\
\text{subject to} \quad \beta \Gamma_{00}(\pm \phi_x) - \gamma \Gamma_{00}(\pm \phi_x) \Gamma_{00x}(\pm \phi_x) - \delta \Gamma_{00x}(\pm \phi_x) = 0.
\end{array}
\right\}
\quad (\text{D10a,b})$$

We can integrate (D10a,b) term by term, using $\int_{x=-\phi_x}^{\phi_x} \Gamma_{00} dx = 2\phi_x$ and the fact that weak surface diffusion makes $\Gamma_{00}(-\phi_x)$ exponentially small (recall that $\Gamma_{00}(-\phi_x) = 0$ when
there was no surface diffusion), to calculate $\Delta\Gamma_0$, substitute this into (3.38) to obtain

$$DR_0 = 1 + \frac{\delta}{2\phi_x\beta} - \left( \frac{\gamma}{\phi_x\beta} + \left( \frac{\delta}{2\phi_x\beta} \right)^2 \right)^{1/2} . \qquad (D11)$$

We therefore have that $DR_0 = 0.5$ when $\gamma = \beta\phi_x/4 + \delta/2$, as in (4.14); we have now evaluated the leading and first-order correction in the limit where Marangoni effects are strong and diffusion is weak.

### D.5. Weak Marangoni effect: region $A_E$

Assume that $\beta = O(1)$, $\gamma \ll \min(1, \alpha, \delta)$ and $v/\epsilon^2 \ll \min(1, \alpha, \delta)$. Rescale $v/\epsilon^2 = n\gamma$ where $n = O(1)$. Expand the concentration field from § 3.3 using $\Gamma_0 = \Gamma_{00} + \gamma\Gamma_{01} + \dots$ and $c_0 = c_{00} + \gamma c_{01} + \dots$ into (3.1)–(3.3a–d). At $O(1)$,

$$\left.
\begin{array}{l}
(c_{00} - \alpha c_{00x})_x = 0, \quad \beta \Gamma_{00} - \delta \Gamma_{00x} + c_{00} - \alpha c_{00x} = 1 \quad \text{in } \mathcal{D}_1, \\
c_{00} - \alpha c_{00x} = 1 \quad \text{in } \mathcal{D}_2, \quad \text{subject to} \quad c_{00}(\phi_x^-) = c_{00}(\phi_x^+), \\
c_{00}(-\phi_x) = c_{00}(2-\phi_x), \quad c_{00}(\pm\phi_x) - \alpha c_{00x}(\pm\phi_x) = 1, \\
\beta \Gamma_{00}(\pm\phi_x) - \delta \Gamma_{00x}(\pm\phi_x) = 0.
\end{array}
\right\}
\qquad (D12a-g)$$

Integrating (D12a-g) on $\mathcal{I}$ and making use of the no net-flux condition, $\int_{x=-\phi_x}^{\phi_x} \Gamma_{00} dx = 2\phi_x$, we obtain (4.16). Therefore, $\Delta\Gamma_0 \approx 2\beta\phi/\delta$, which is substituted into (3.38) to give (4.15). The corrections to $\Gamma_0 = \Gamma_{00}$ and $DR_0 = 1$ are small provided $\gamma \ll \min(1, \alpha, \delta)$ and $v/\epsilon^2 \ll \min(1, \alpha, \delta)$, defining the boundaries of the asymptotic region $A_E$ and $D_E$.

### REFERENCES

* 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., STEFFES, C. & HARDT, S. 2010 Thermocapillary flow on superhydrophobic surfaces. *Phys. Rev. E* **82**, 037301.

* BIBEN, T. & JOLY, L. 2008 Wetting on nanorough surfaces. *Phys. Rev. Lett.* **100** (18), 186103.

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

* BOND, W.N. & NEWTON, D.A. 1928 Bubbles, drops, and Stokes' law. *Phil. Mag.* **5** (30), 794-800.

* CHANG, C.H. & FRANSES, E.I. 1995 Adsorption dynamics of surfactants at the air/water interface: a critical review of mathematical models, data, and mechanisms. *Colloids Surf.* A **100**, 1-45.

* 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. Engng* **88**, 71-81.

* FACCHINI, M.C., DECESARI, S., MIRCEA, M., FUZZI, S. & LOGLIO, G. 2000 Surface tension of atmospheric wet aerosol and cloud/fog droplets in relation to their organic carbon content and chemical composition. *Atmos. Environ.* **34**, 4853-4857.

* FRUMKIN, A.N. & LEVICH, V.G. 1947 Effect of surface-active substances on movements at the boundaries of liquid phases. *Zhur. Fiz. Khim.* **21**, 1183-1204.

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

* GOLOVIN, K.B., GOSE, J.W., PERLIN, M., CECCIO, S.L. & TUTEJA, A. 2016 Bioinspired surfaces for turbulent drag reduction. *Phil. Trans. R. Soc.* A **374**, 20160189.

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

* 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. *Trans. ASME J. Heat Transfer* **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.

LAUGA, E. & STONE, H.A. 2003 Effective slip in pressure-driven Stokes flow. *J. Fluid Mech.* **489**, 55–77.

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

LEVICH, V. 1962 *Physicochemical Hydrodynamics*. Prentice-Hall.

LEWIS, M.A. 1991 Chronic and sublethal toxicities of surfactants to aquatic animals: a review and risk assessment. *Water Res.* **25**, 101–113.

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. Fluid Mech.* **949**, A18.

MCNAIR, R., JENSEN, O.E. & LANDEL, J.R. 2022 Surfactant spreading in a two-dimensional cavity and emergent contact-line singularities. *J. Fluid Mech.* **930**, A15.

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

OU, J. & ROTHSTEIN, J.P. 2005 Direct velocity measurements of the flow past drag-reducing ultrahydrophobic surfaces. *Phys. Fluids* **17**, 103606.

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. USA* **114**, 7254–7259.

PEREIRA, R., ASHTON, I., SABBAGHZADEH, B., SHUTLER, J.D. & UPSTILL-GODDARD, R.C. 2018 Reduced air–sea CO₂ 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.

ROTHSTEIN, J.P. 2010 Slip on superhydrophobic surfaces. *Annu. Rev. Fluid Mech.* **42**, 89–109.

SBRAGAGLIA, M. & PROSPERETTI, A. 2007 A note on the effective slip properties for microchannel flows with ultrahydrophobic surfaces. *Phys. Fluids* **19**, 043603.

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.

SUNDIN, J., ZALESKI, S. & BAGHERI, S. 2021 Roughness on liquid-infused surfaces induced by capillary waves. *J. Fluid Mech.* **915**, R6.

TAYLOR, G.I. 1953 Dispersion of soluble matter in solvent flowing slowly through a tube. *Proc. R. Soc. A* **219**, 186–203.

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. USA* **120** (3), e2211092120.

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

TREFETHEN, L.N. 2000 *Spectral Methods in MATLAB*. SIAM.
S.D. Tomlinson and others

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.

WEXLER, J.S., JACOBI, I. & STONE, H.A. 2015 Shear-driven failure of liquid-infused surfaces. *Phys. Rev. Lett.* **114**, 168301.

WONG, T., KANG, S.H., TANG, S.K.Y., SMYTHE, E.J., HATTON, B.D., GRINTHAL, A. & AIZENBERG, J. 2011 Bioinspired self-repairing slippery surfaces with pressure-stable omniphobicity. *Nature* **477**, 443–447.