# Laminar drag reduction in surfactant- contaminated superhydrophobic channels

Samuel D. Tomlinson1<EMAIL_ADDRESS>Frédéric Gibou2  Paolo
Luzzatto-Fegiz2  Fernando Temprano-Coleto3  Oliver E. Jensen1  Julien R.
Landel1  1Department of Mathematics, University of Manchester, Oxford Road,
Manchester M13 9PL, UK 2Department of Mechanical Engineering, University of
California, Santa Barbara, CA 93106, USA 3Andlinger Center for Energy and the
Environment, Princeton University, Princeton, NJ 08544, USA

###### Abstract

While 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 impact 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.

###### keywords:

Marangoni convection, drag reduction, microfluidics

##  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 et al., 2016; Park et al., 2021), the reduced
liquid–solid contact area at the SHS reduces viscous drag at the surface when
compared to 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 et al.,
2014; Lee et al., 2016; Landel et al., 2020), high-Reynolds-number turbulent
flows (e.g. Park et al., 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 et al., 2010; Cheng et al., 2015;
Lam et al., 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 impact
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 impact of surfactants on SHS drag reduction was suggested by experimental
studies in the last decade, which reported only a modest drag reduction when
compared to 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 et al. (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 to 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 to 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 impact.
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, whilst 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 et al., 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 as in 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 mobilization 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 mobilization length
depends on the normalised surfactant concentration, Marangoni number,
Damköhler number and Biot number. Since the mobilization length can be much
longer for typical small-scale applications, they concluded that the
mobilization length alone controls the slip in most small-scale applications.
If the grating length is smaller than the mobilization 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 with 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 non-linear 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 non-linear
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 et al., 2015). They can self-repair and are more robust than
SHSs if properly designed. 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}\,$mol/m3 for water–air SHSs and
$\hat{C}_{c}=5\times 10^{-5}\,$mol/m3 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.

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 et al.,
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 impact 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{\boldsymbol{x}}=(\hat{x},\,\hat{y},\,\hat{z})$. The liquid is assumed to
be incompressible and Newtonian with dynamic viscosity $\hat{\mu}$, velocity
$\hat{\boldsymbol{u}}=(\hat{u}(\hat{\boldsymbol{x}}),\,\hat{v}(\hat{\boldsymbol{x}}),\,\hat{w}(\hat{\boldsymbol{x}}))$,
pressure $\hat{p}(\hat{\boldsymbol{x}})$, bulk surfactant distribution
$\hat{c}(\hat{\boldsymbol{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

$\displaystyle\hat{\mathcal{D}}_{1}$
$\displaystyle=\\{\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}]\\},$
(2.1a) $\displaystyle\hat{\mathcal{D}}_{2}$
$\displaystyle=\\{\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}]\\},$
(2.1b)

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

$\displaystyle\hat{\mathcal{I}}$
$\displaystyle=\\{\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}]\\},$
(2.2a) $\displaystyle\hat{\mathcal{R}}$
$\displaystyle=\\{\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}]\\},$
(2.2b) $\displaystyle\hat{\mathcal{S}}$
$\displaystyle=\\{\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}]\\},$
(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.

(a) (b)

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{\boldsymbol{x}}=\boldsymbol{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{\mathcal{R}}$ (outlined in red), solid surfaces, $\hat{\mathcal{S}}$
(yellow), and liquid–gas interfaces, $\hat{\mathcal{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{\mathcal{D}}_{1}$ and
$\hat{\mathcal{D}}_{2}$ (magenta).

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{\mathcal{D}}_{1}$ and
$\hat{\mathcal{D}}_{2}$, we therefore have

$\hat{\boldsymbol{\nabla}}\cdot\hat{\boldsymbol{u}}=0,\quad\hat{\mu}\hat{\nabla}^{2}\hat{\boldsymbol{u}}-\hat{\boldsymbol{\nabla}}\hat{p}=\boldsymbol{0},\quad\hat{D}\hat{\nabla}^{2}\hat{c}-\hat{\boldsymbol{u}}\cdot\hat{\boldsymbol{\nabla}}\hat{c}=0,$
(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{\mathcal{I}}$, we balance the tangential components of the stress,
$\hat{\mathsfbi{T}\,}\cdot\boldsymbol{n}$ where
$\hat{\mathsfbi{T}\,}=-\hat{p}\mathsfbi{I}+\hat{\mu}(\hat{\boldsymbol{\nabla}}\hat{\boldsymbol{u}}+(\hat{\boldsymbol{\nabla}}\hat{\boldsymbol{u}})^{T})$
and $\boldsymbol{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

$\hat{\mu}\boldsymbol{n}\cdot\hat{\boldsymbol{\nabla}}\hat{u}-\hat{A}\hat{\Gamma}_{\hat{x}}=0,\quad\hat{\mu}\boldsymbol{n}\cdot\hat{\boldsymbol{\nabla}}\hat{w}-\hat{A}\hat{\Gamma}_{\hat{z}}=0,\quad\hat{v}=0,\quad\hat{D}\boldsymbol{n}\cdot\hat{\boldsymbol{\nabla}}\hat{c}-\hat{K}_{a}\hat{c}+\hat{K}_{d}\hat{\Gamma}=0,\\\
\quad\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,$
(2.5)

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}.$
(2.6a, 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.$ (2.7a–d)

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

$\hat{\boldsymbol{q}}(-\phi_{x}\hat{P}_{x},\,\hat{y},\,\hat{z})=\hat{\boldsymbol{q}}((2-\phi_{x})\hat{P}_{x},\,\hat{y},\,\hat{z}),\quad\hat{\boldsymbol{q}}(\hat{x},\,\hat{y},\,-\hat{P}_{z})=\hat{\boldsymbol{q}}(\hat{x},\,\hat{y},\,\hat{P}_{z}).$
(2.8a, b)

The bulk flow and concentration are continuous at the boundary between
$\hat{\mathcal{D}}_{1}$ and $\hat{\mathcal{D}}_{2}$

$\hat{\boldsymbol{q}}((\phi_{x}\hat{P}_{x})^{-},\,\hat{y},\,\hat{z})=\hat{\boldsymbol{q}}((\phi_{x}\hat{P}_{x})^{+},\,\hat{y},\,\hat{z}),$
(2.9)

where the superscripts $-$ and $+$ mean that the boundary condition is
evaluated in $\hat{\mathcal{D}}_{1}$ and $\hat{\mathcal{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
(1)–(2.9) across the channel to show how the streamwise bulk surfactant flux
changes as surfactants adsorb and desorb

$\frac{\text{d}}{\text{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}})\,\text{d}\hat{y}\,\text{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})\,\text{d}\hat{z};$
(2.10)

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

$\frac{\text{d}}{\text{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}})\,\text{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})\,\text{d}\hat{z}.$
(2.11)

The boundary-value problem (1)–(2.9) 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}\,\text{d}\hat{y}\,\text{d}\hat{z},$
(2.12)

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}})\,\text{d}\hat{y}\,\text{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}})\,\text{d}\hat{z},$
(2.13)

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})\,\text{d}\hat{z}\,\text{d}\hat{x}=0.$
(2.14)

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

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\,\text{d}\hat{y}\,\text{d}\hat{z}/(4\hat{P}_{z}\hat{H})$
is the cross-channel average. In the limit when $\hat{\mathcal{I}}$ is
immobilized (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}},$
(2.15)

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.15) differs from other definitions of the drag reduction which compare the
flow over a SHS to 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.15) 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\boldsymbol{x}_{\perp}=\frac{\hat{\boldsymbol{x}}_{\perp}}{\epsilon\hat{P}_{x}},\quad
u=\frac{\hat{u}}{\epsilon\hat{U}},\quad\boldsymbol{u}_{\perp}=\frac{\hat{\boldsymbol{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}},$ (2.16a–g)

where $\hat{\boldsymbol{x}}_{\perp}\equiv(\hat{y},\,\hat{z})$ and
$\hat{\boldsymbol{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 et al., 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

$\displaystyle\mathcal{D}_{1}$
$\displaystyle=\\{x\in[-\phi_{x},\,\phi_{x}]\\}\times\\{y\in[0,\,2]\\}\times\\{z\in[-P_{z},\,P_{z}]\\},$
(2.17a) $\displaystyle\mathcal{D}_{2}$
$\displaystyle=\\{x\in[\phi_{x},\,2-\phi_{x}]\\}\times\\{y\in[0,\,2]\\}\times\\{z\in[-P_{z},\,P_{z}]\\},$
(2.17b)

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

$\displaystyle\mathcal{I}$
$\displaystyle=\\{x\in[-\phi_{x},\,\phi_{x}]\\}\times\\{y\in\\{0,\,2\\}\\}\times\\{z\in[-\phi_{z}P_{z},\,\phi_{z}P_{z}]\\},$
(2.18a) $\displaystyle\mathcal{R}$
$\displaystyle=\\{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}]\\},$
(2.18b) $\displaystyle\mathcal{S}$
$\displaystyle=\\{x\in[\phi_{x},\,2-\phi_{x}]\\}\times\\{y\in\\{0,\,2\\}\\}\times\\{z\in[-P_{z},\,P_{z}]\\}.$
(2.18c)

We substitute the non-dimensionalisation (2.2) into the governing equations
(1)–(2.15) to acquire rescaled governing equations in terms of the non-
dimensional variables given in (2.2) and $\epsilon^{2}$. In $\mathcal{D}_{1}$
and $\mathcal{D}_{2}$,

$\epsilon^{2}u_{x}+\boldsymbol{\nabla}_{\perp}\cdot\boldsymbol{u}_{\perp}=0,\quad\epsilon^{2}\boldsymbol{u}_{xx}+\nabla^{2}_{\perp}\boldsymbol{u}-\boldsymbol{\nabla}p=\boldsymbol{0},\\\
\quad\Pen^{-1}(\epsilon^{2}c_{xx}+\nabla^{2}_{\perp}c)-\epsilon^{2}uc_{x}-\boldsymbol{u}_{\perp}\cdot\boldsymbol{\nabla}_{\perp}c=0,$
(2.20)

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

$\boldsymbol{n}\cdot\boldsymbol{\nabla}u-\Ma\Gamma_{x}=0,\quad\boldsymbol{n}\cdot\boldsymbol{\nabla}w-\Ma\Gamma_{z}=0,\quad
v=0,\quad\boldsymbol{n}\cdot\boldsymbol{\nabla}c-Da(c-\Gamma)=0,\\\
\quad\Pen^{-1}_{I}(\epsilon^{2}\Gamma_{xx}+\Gamma_{zz})+\Bi(c-\Gamma)-\epsilon^{2}(u\Gamma)_{x}-(w\Gamma)_{z}=0,$
(2.22)

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,
$\Pen_{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-\Pen^{-1}_{I}\Gamma_{x}=0\quad\text{at}\quad x=\pm\phi_{x},\quad
w\Gamma-\Pen^{-1}_{I}\Gamma_{z}=0\quad\text{at}\quad z=\pm\phi_{z}P_{z}.$
(2.23a, b)

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

$u=0,\quad w=0,\quad v=0,\quad\ c_{y}=0.$ (2.24a–d)

Periodicity and continuity between subdomains require that, for
$\boldsymbol{q}=(\boldsymbol{u},\,p_{x},\,c)$,

$\boldsymbol{q}(-\phi_{x},\,y,\,z)=\boldsymbol{q}(2-\phi_{x},\,y,\,z),\quad\boldsymbol{q}(x,\,y,\,-P_{z})=\boldsymbol{q}(x,\,y,\,P_{z}),$
(2.25a, b)

and

$\boldsymbol{q}(\phi_{x}^{-},\,y,\,z)=\boldsymbol{q}(\phi_{x}^{+},\,y,\,z).$
(2.26)

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\,\text{d}y\,\text{d}z=1,$ (2.27)
$\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\left(uc-\frac{c_{x}}{\Pen}\right)\,\text{d}y\,\text{d}z+\frac{2\Da}{\Bi\Pen}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}\left(u\Gamma-\frac{\Gamma_{x}}{\Pen_{I}}\right)\,\text{d}z=1,$
(2.28)
$\int_{x=-\phi_{x}}^{\phi_{x}}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}(\Gamma-c)\,\text{d}z\,\text{d}x=0,$
(2.29)

where bulk and surface fluxes are related to adsorption–desorption at the
interface via

$\displaystyle\frac{\text{d}}{\text{d}x}\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\left(uc-\frac{c_{x}}{\Pen}\right)\,\text{d}y\,\text{d}z$
$\displaystyle=\frac{2\Da}{\epsilon^{2}\Pen}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}(\Gamma-c)\,\text{d}z,$
(2.30a)
$\displaystyle\frac{\text{d}}{\text{d}x}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}\left(u\Gamma-\frac{\Gamma_{x}}{\Pen_{I}}\right)\,\text{d}z$
$\displaystyle=-\frac{\Bi}{\epsilon^{2}}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}(\Gamma-c)\,\text{d}z.$
(2.30b)

The group $\Da/(\Bi\Pen)=\hat{K}_{a}/(\hat{K}_{d}\hat{H})$ in (2.28) 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/(\Bi\Pen)$ compares $\hat{L}_{d}$ to the channel height
$\hat{H}$. The drag reduction (2.15) becomes

${DR}=\frac{\Delta p_{R}-\Delta p}{\Delta p_{R}-\Delta p_{U}},$ (2.31)

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}\cdot\,\text{d}y\text{d}z/(4P_{z})$.

The non-dimensional governing equations (2.2), boundary conditions
(2.20)–(2.26) and flux constraints (2.27)–(2.29), define a 3D boundary-value
problem for $\boldsymbol{u}$, $p$, $c$ and $\Gamma$. The solution depends on 9
dimensionless groups ($\epsilon$, $\Pen$, $\Ma$, $\Pen_{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 8), 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.31).

##  3 Model

With 9 non-dimensional parameters in the problem, we choose distinguished
limits to reveal different dominant physical balances. We take
$\epsilon\rightarrow 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\rightarrow 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 5
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.30), we assume that
$\Pen=\Pen_{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 will be shown in §3.3, they
reduce to a simpler coupled nonlinear system of second-order ODEs:

$\displaystyle(c_{0}-\alpha
c_{0x})_{x}-\frac{\nu}{\epsilon^{2}}(\Gamma_{0}-c_{0})=0\quad\text{in}\quad$
$\displaystyle\mathcal{D}_{1},$ (3.32a)
$\displaystyle(\beta\Gamma_{0}-\gamma\Gamma_{0}\Gamma_{0x}-\delta\Gamma_{0x})_{x}-\frac{\nu}{\epsilon^{2}}(c_{0}-\Gamma_{0})=0\quad\text{in}\quad$
$\displaystyle\mathcal{D}_{1},$ (3.32b) $\displaystyle(c_{0}-\alpha
c_{0x})_{x}=0\quad\text{in}\quad$ $\displaystyle\mathcal{D}_{2}.$ (3.32c)

The total flux of surfactant (2.28) becomes

$\displaystyle c_{0}-\alpha
c_{0x}+\beta\Gamma_{0}-\gamma\Gamma_{0}\Gamma_{0x}-\delta\Gamma_{0x}$
$\displaystyle=1\quad\text{in}\quad\mathcal{D}_{1},$ (3.33a) $\displaystyle
c_{0}-\alpha c_{0x}$ $\displaystyle=1\quad\text{in}\quad\mathcal{D}_{2}.$
(3.33b)

Continuity of bulk surfactant (2.26) and bulk surfactant flux (2.28) between
unit cells, and no flux of interfacial surfactant (2.22) through contact
lines, are together given by

$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}\quad
x=\pm\phi_{x}.$ (3.35)

Equations (3.32)–(3.33) 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.32)), 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.32)–(3.33), 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.32)–(3.33), 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 $\Pen=\Pen_{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})+...,$
(3.36)

in the governing equations (2.2)–(2.31) and take the $O(1)$ approximation. In
domains $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$,

$\boldsymbol{\nabla}_{\perp}\cdot\boldsymbol{u}_{\perp
0}=0,\quad\nabla^{2}_{\perp}\boldsymbol{u}_{0}-\boldsymbol{\nabla}p_{0}=\boldsymbol{0},\quad\Pen^{-1}\nabla^{2}_{\perp}c_{0}-\boldsymbol{u}_{\perp
0}\cdot\boldsymbol{\nabla}_{\perp}c_{0}=0.$ (3.37a–c)

Along the interface $\mathcal{I}$,

$\boldsymbol{n}\cdot\boldsymbol{\nabla}u_{0}-\Ma\Gamma_{0x}=0,\quad\boldsymbol{n}\cdot\boldsymbol{\nabla}w_{0}-\Ma\Gamma_{0z}=0,\quad
v_{0}=0,\\\
\boldsymbol{n}\cdot\boldsymbol{\nabla}c_{0}-Da(c_{0}-\Gamma_{0})=0,\quad\Pen^{-1}_{I}\Gamma_{0zz}+\Bi(c_{0}-\Gamma_{0})-(w_{0}\Gamma_{0})_{z}=0,$
(3.39)

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

$w_{0}\Gamma_{0}-\Pen^{-1}_{I}\Gamma_{0z}=0\quad\text{at}\quad
z=\pm\phi_{z}P_{z}.$ (3.40)

There are inner regions near $x=\pm\phi_{x}$, within which
$u_{0}\Gamma_{0}-\Pen^{-1}_{I}\Gamma_{0x}=0$ is imposed on
$\partial\mathcal{I}$, and the boundary conditions
$\boldsymbol{u}_{0}(\phi_{x}^{-},\,y,\,z)=\boldsymbol{u}_{0}(\phi_{x}^{+},\,y,\,z)$
and
$\boldsymbol{u}_{0}(-\phi_{x},\,y,\,z)=\boldsymbol{u}_{0}(2-\phi_{x},\,y,\,z)$
are imposed between $\mathcal{D}_{1}$ and $\mathcal{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 $\mathcal{D}_{1}$ and
$\mathcal{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),$ (3.41a, b)

and continuity of volume and surfactant flux between $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$, such that

$\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}u_{0}\,\text{d}y\,\text{d}z=1,$ (3.42)
$\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}y\,\text{d}z+\frac{2\Da}{\Bi\Pen}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z},}\left(u_{0}\Gamma_{0}-\frac{\Gamma_{0x}}{\Pen_{I}}\right)\,\text{d}z=1,$
(3.43)

for all $-\phi_{x}\leq x\leq 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.$ (3.44a–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}\quad
z=\pm P_{z}.$ (3.45a–d)

The drag reduction becomes

${DR}_{0}=\frac{\Delta p_{R}-\Delta p_{0}}{\Delta p_{R}-\Delta p_{U}}.$ (3.46)

(a) (b)

Figure 2: (a) Contour plot of $\tilde{U}$ defined by (3.47a, b), the leading-
order contribution to the streamwise flow due to the pressure gradient,
$p_{0x}$. (b) Contour plot of $\bar{U}$ defined by (3.49), the leading-order
contribution to the streamwise flow due to the surfactant gradient,
$\Gamma_{0x}$. Since $p_{0x}<0$, $\tilde{U}$ contributes positively to the
leading-order streamwise velocity component $u_{0}$, whilst $\bar{U}$
contributes negatively, since $\Gamma_{0x}>0$, following (2). 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$.

From (3.36)–(3.44a–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.36b), and $\Gamma_{0x}$, from
(3.37a–ca). Using linear superposition, we can write

$u_{0}=\tilde{U}p_{0x}+\Ma\bar{U}\Gamma_{0x}\quad\text{in}\
\mathcal{D}_{1}\quad\text{and}\quad u_{0}=\breve{U}p_{0x}\quad\text{in}\
\mathcal{D}_{2},$ (3.47a, b)

where $\tilde{U}(y,\,z)$ and $\bar{U}(y,\,z)$ are velocity contributions
related to the bulk pressure and surfactant gradient, respectively, in
$\mathcal{D}_{1}$, and $\breve{U}(y)$ is the velocity contribution related to
the bulk pressure gradient in $\mathcal{D}_{2}$. Hence, we must solve the
following boundary-value problems: streamwise flow driven by a pressure
gradient over $\mathcal{D}_{1}$,

$\nabla^{2}_{\perp}\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,\\\
\tilde{U}(2,\,z_{ns})=0,\quad\tilde{U}_{z}(y,\,-P_{z})=0,\quad\tilde{U}_{z}(y,\,P_{z})=0;$
(3.49)

streamwise flow driven by a surfactant-induced Marangoni shear stress over
$\mathcal{D}_{1}$,

$\nabla^{2}_{\perp}\bar{U}=0,\quad\text{subject
to}\quad\bar{U}_{y}(0,\,z_{s})=1,\quad\bar{U}(0,\,z_{ns})=0,\quad\bar{U}_{y}(2,\,z_{s})=-1,\\\
\bar{U}(2,\,z_{ns})=0,\quad\bar{U}_{z}(y,\,-P_{z})=0,\quad\bar{U}_{z}(y,\,P_{z})=0;$
(3.51)

and streamwise flow driven by a pressure gradient over $\mathcal{D}_{2}$,

$\breve{U}_{yy}=1,\quad\text{subject to}\quad\breve{U}(0)=0,\
\quad\breve{U}(2)=0;$ (3.52a–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.47a, b, 3.49) 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.51) can be integrated to give $\breve{U}=y(y-2)/2$ and
$\breve{Q}=\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\breve{U}\,\text{d}y\,\text{d}z=-4P_{z}/3$.
We define the volume and surface fluxes

$\tilde{Q}=\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\tilde{U}\,\text{d}y\,\text{d}z,\quad\bar{Q}=\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\bar{U}\,\text{d}y\,\text{d}z,$
(3.53a, b)
$q=\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}u_{0}\,\text{d}z,\quad\tilde{q}=\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}\tilde{U}\,\text{d}z,\quad\bar{q}=\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}\bar{U}\,\text{d}z.$
(3.54a–c)

Bulk ($\tilde{Q}$, $\bar{Q}$) and surface fluxes ($\tilde{q}$, $\bar{q}$) are
plotted as functions of $\phi_{z}$ and $P_{z}$ in figure 3.

(a) (b) (c) (d)

Figure 3: Plot of (a) $\tilde{Q}$ defined in (2a) and (b) $\tilde{q}$ defined
in (3.53a, bb), the contributions to the bulk and surface flux due to the
pressure gradient $p_{0x}<0$, for varying $\phi_{z}$ and $P_{z}$. Plot of (c)
$\bar{Q}$ defined in (2b) and (d) $\bar{q}$ defined in (3.53a, bc), the
contributions to the bulk and surface flux due to the surfactant gradient
$\Gamma_{0x}>0$, for varying $\phi_{z}$ and $P_{z}$. Since $p_{0x}<0$,
$\tilde{Q}$ and $\tilde{q}$ contribute positively to the leading-order
velocity flux, whilst $\bar{Q}$ and $\bar{q}$ contribute negatively, since
$\Gamma_{0x}>0$, following (3a).

Next, we substitute the expansions (3.36) into the governing equations
(2.2)–(2.24a–d) 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 concentration field $\Gamma_{1}$ driven by streamwise gradients of the
leading-order quantities ($\boldsymbol{u}_{0}=(u_{0}(x,y,z),\,0,\,0)$,
$c_{0}(x)$ and $\Gamma_{0}(x)$), and to close the leading-order problem. In
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$,

$\boldsymbol{\nabla}_{\perp}\cdot\boldsymbol{u}_{\perp
1}=-u_{0x},\quad\nabla^{2}_{\perp}\boldsymbol{u}_{\perp
1}-\boldsymbol{\nabla}_{\perp}p_{1}=\boldsymbol{0},\quad\Pen^{-1}\nabla^{2}_{\perp}c_{1}=-\Pen^{-1}c_{0xx}+u_{0}c_{0x}.$
(3.55a–c)

Along $\mathcal{I}$,

$\boldsymbol{n}\cdot\boldsymbol{\nabla}w_{1}-\Ma\Gamma_{1z}=0,\quad
v_{1}=0,\quad\boldsymbol{n}\cdot\boldsymbol{\nabla}c_{1}-Da(c_{1}-\Gamma_{1})=0,\\\
\Pen^{-1}_{I}\Gamma_{1zz}+\Bi(c_{1}-\Gamma_{1})-\Gamma_{0}w_{1z}=-\Pen_{I}^{-1}\Gamma_{0xx}+(u_{0}\Gamma_{0})_{x},$
(3.57)

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

$w_{1}\Gamma_{0}-\Pen^{-1}_{I}\Gamma_{1z}=0\quad\text{at}\quad
z=\pm\phi_{z}P_{z}.$ (3.58)

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

$v_{1}=0,\quad w_{1}=0,\quad c_{1y}=0.$ (3.59a–c)

Transverse periodicity can be rewritten as the symmetry conditions

$v_{1z}=0,\quad w_{1}=0,\ \quad c_{1z}=0\quad\text{at}\quad z=\pm P_{z}.$
(3.60a–c)

Solvability conditions imposed on (3)–(3.59a–c) constrain $u_{0}$, $c_{0}$ and
$\Gamma_{0}$ at $O(1)$. The forcing on the right-hand sides of (3, 3.55a–c)
must be orthogonal to each vector in the null space of the linear operator
that is adjoint to the left-hand sides of (3, 3.55a–c). Conveniently, these
conditions are provided by the conservation arguments that result in the
leading-order velocity and surfactant flux conditions (3.42, 3.43).
Substituting the leading-order streamwise velocity (2) into the leading-order
velocity flux condition (3.42), we obtain

$\tilde{Q}p_{0x}+\Ma\bar{Q}\Gamma_{0x}=1,\quad
q=\tilde{q}p_{0x}+\Ma\bar{q}\Gamma_{0x}\quad\text{in}\quad\mathcal{D}_{1},\quad\breve{Q}p_{0x}=1\quad\text{in}\quad\mathcal{D}_{2}.$
(3.61a–c)

Substituting (3) into the leading-order surfactant flux condition (3.43) and
remembering that $c_{0}(x)=\Gamma_{0}(x)$, we obtain the strong cross-channel
diffusion problem

$\displaystyle\left(\beta+1\right)c_{0}-\gamma
c_{0}c_{0x}-\left(\alpha+\delta\right)c_{0x}$
$\displaystyle=1\quad\text{in}\quad\mathcal{D}_{1},$ (3.62a) $\displaystyle
c_{0}-\alpha c_{0x}$ $\displaystyle=1\quad\text{in}\quad\mathcal{D}_{2},$
(3.62b)

subject to continuity of bulk surfactant (3.33a, b). In (3.62), we have
introduced

$\displaystyle\alpha$
$\displaystyle=\frac{4P_{z}}{\Pen}=\frac{4\hat{D}\hat{P}_{z}}{\hat{H}^{2}\hat{U}}\quad\text{(bulk
diffusion)},$ (3.63a) $\displaystyle\beta$
$\displaystyle=\frac{2\Da\tilde{q}}{\Bi\Pen\tilde{Q}}=\frac{2\hat{K}_{a}\tilde{q}}{\hat{H}\hat{K}_{d}\tilde{Q}}\quad\text{(partition
coefficient)},$ (3.63b) $\displaystyle\gamma$
$\displaystyle=\frac{2\Ma\Da}{\Bi\Pen}\left(\frac{\tilde{q}\bar{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}\bar{Q}}{\tilde{Q}}-\bar{q}\right)\quad\text{(surfactant
strength)},$ (3.63c) $\displaystyle\delta$
$\displaystyle=\frac{4\phi_{z}P_{z}\Da}{\Bi\Pen\Pen_{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)}.$ (3.63d)

The bulk diffusion coefficient $\alpha>0$ is a rescaled (by
$4\hat{P}_{z}/\hat{H}$) inverse Péclet number
$1/\Pen=\hat{D}/(\hat{H}\hat{U})$, nominally assumed $O(1)$. The coefficient
$\alpha$ characterizes 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.62a)
determines the portion of the surfactant flux transported by the pressure-
driven flow along the interface in comparison to 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}\bar{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$ characterizes the nonlinear impact 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/\Pen_{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.36)–(3.44a–d) has been reduced to the
5-parameter problem (3.62) subject to (3.33a, b), a special case (where
$\nu/\epsilon^{2}\gg O(1,\,\alpha,\,\beta,\,\delta)$) of the 6-parameter
problem (3.32)–(3.33), requiring the solution of a first-order ODE to
determine $c_{0}$ and hence ${DR}_{0}$ (explained in §3.4 below). Then, the
leading-order streamwise velocity (2) and first-order problem (3)–(3.59a–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}^{\prime}(x,\,y,\,z)$, we assume that
$\|c_{0}^{\prime}\|\equiv\sqrt{\langle c_{0}^{\prime 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/\Pen=1/\Pen_{I}=\Bi=O(\epsilon^{2})$
and $\Da=\Ma=O(1)$ for $\epsilon\ll 1$, substituting the asymptotic expansion
(3.36) into the governing equations (2.2)–(2.30) (as in §3.2.1). We then use
the first-order system to eliminate $c_{0}^{\prime}$ 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.62) accounting for shear dispersion,

$\displaystyle\left(\beta+1\right)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)$
$\displaystyle=1\quad\text{in}\quad\mathcal{D}_{1},$ (3.64a) $\displaystyle
c_{0}-\frac{\epsilon^{2}}{\alpha}\left(s_{4}c_{0x}\right)$
$\displaystyle=1\quad\text{in}\quad\mathcal{D}_{2},$ (3.64b)

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

$s_{1}=-\frac{16P_{z}^{2}\langle\tilde{U}\tilde{C}\rangle}{\tilde{Q}^{2}},\quad
s_{2}=-\frac{16P_{z}^{2}\tilde{q}(\tilde{Q}(\langle\bar{U}\tilde{C}\rangle+\langle\tilde{U}\bar{C}\rangle)-2\bar{Q}\langle\tilde{U}\tilde{C}\rangle)}{\tilde{Q}^{2}(\tilde{q}\bar{Q}-\bar{q}\tilde{Q})},\\\
s_{3}=-\frac{16P_{z}^{2}\tilde{q}^{2}(\bar{Q}^{2}\langle\tilde{U}\tilde{C}\rangle-\bar{Q}\tilde{Q}(\langle\bar{U}\tilde{C}\rangle+\langle\tilde{U}\bar{C}\rangle)+\tilde{Q}^{2}\langle\bar{U}\bar{C}\rangle)}{\tilde{Q}^{2}(\tilde{q}\bar{Q}-\bar{q}\tilde{Q})^{2}},\quad
s_{4}=-\frac{16P_{z}^{2}\langle\breve{U}\breve{C}\rangle}{\breve{Q}^{2}},$
(3.66)

where $\tilde{C}$, $\bar{C}$ and $\breve{C}$ are solutions of the boundary-
value problems (B.112a, b)–(B.116), respectively (see Appendix B.1), which
depend on geometrical factors and $\Da$. Here, ($\tilde{U},\,\bar{U}$) are
given by (3.47a, b, 3.49), ($\tilde{Q},\,\bar{Q})$ by (2), and
($\tilde{q},\,\bar{q}$) by (3.53a, b). 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$, since
$\breve{U}=y(y-2)/2$, $\breve{C}=1/5-y^{3}/6+y^{4}/24$ and
$\breve{Q}=-4P_{z}/3$ in $\mathcal{D}_{2}$. Comparing the strong (3.62) and
moderate cross-channel diffusion problems (3.64), streamwise diffusion terms
proportional to $\alpha$ and $\delta$ in (3.62) have been replaced by shear
dispersion terms proportional to $\epsilon^{2}/\alpha$ in (3.64), some of
which are nonlinear due to Marangoni effects.

####  3.2.3 Composite equation

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

$\displaystyle(\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)$
$\displaystyle=1\ \text{in}\ \mathcal{D}_{1},$ (3.67a) $\displaystyle
c_{0}-\alpha c_{0x}-\frac{\epsilon^{2}}{\alpha}\left(s_{4}c_{0x}\right)$
$\displaystyle=1\ \text{in}\ \mathcal{D}_{2},$ (3.67b)

subject to continuity of bulk surfactant (3.33a, b). The dependence of
$\alpha$, $\beta$, $\gamma$ and $\delta$ on $\Ma$, $\Pen$, $\Pen_{I}$, $\Bi$,
$\Da$, $\phi_{z}$ and $P_{z}$ is given in (3.63). The dependence of $s_{1}$,
$s_{2}$, $s_{3}$ and $s_{4}$ on $\Da$, $\phi_{z}$ and $P_{z}$ is given in
(3.64) and discussed in Appendix B.1. Together, the solution to (3.67) depends
on nine dimensionless parameters. The numerical method used to solve (3.33a,
3.33b, 3.67) is outlined in Appendix A.1. Incorporation of shear dispersion
effects in (3.67), through perturbation of the cross-channel-averaged
concentration ($c_{0}=\langle c_{0}\rangle(x)+c_{0}^{\prime}(x,\,y,\,z)$),
enables the strong-exchange model to explore regimes where cross-channel
concentration gradients first become significant
($\|c_{0}^{\prime}\|=O(\langle c_{0}\rangle)$ for some $x$), enabling us to
estimate approximately the boundary of validity of the 1D model (3.32)–(3.33)
(Appendix B.2). When bulk diffusion is strong enough for shear dispersion to
become negligible, (3.67) reduces to (3.62), or (3.33) 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}\mathscr{B}$
and $\Da=\epsilon^{2}\mathscr{D}$ with $\mathscr{B}=\mathscr{D}=O(1)$, whilst
retaining $\Pen=\Pen_{I}=\Ma=O(1)$. Substituting the expansion (3.36) into the
governing equations (2.2)–(2.30), then the leading-order strong cross-channel
diffusion equations (3.36), (3.37a–ca–c), (3.40)–(3.42) and (3.43)–(3.44a–d)
are recovered. Along the interface $\mathcal{I}$

$c_{0y}=0,\quad\Pen^{-1}_{I}\Gamma_{0zz}-(w_{0}\Gamma_{0})_{z}=0.$ (3.68a, b)

The total flux of surfactant (2.28) is given by

$\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}y\,\text{d}z+\frac{2\mathscr{D}}{\mathscr{B}\Pen}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}\left(u_{0}\Gamma_{0}-\frac{\Gamma_{0x}}{\Pen_{I}}\right)\,\text{d}z=1.$
(3.69)

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

$\displaystyle\frac{\text{d}}{\text{d}x}\left(\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}y\,\text{d}z\right)$
$\displaystyle=\frac{2\mathscr{D}}{\Pen}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}(\Gamma_{0}-c_{0})\,\text{d}z,$
(3.70a)
$\displaystyle\frac{\text{d}}{\text{d}x}\left(\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}\left(u_{0}\Gamma_{0}-\frac{\Gamma_{0x}}{\Pen_{I}}\right)\,\text{d}z\right)$
$\displaystyle=-\mathscr{B}\int_{z=-\phi_{z}P_{z},\,y=0}^{\phi_{z}P_{z}}(\Gamma_{0}-c_{0})\,\text{d}z.$
(3.70b)

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\bar{U}\Gamma_{0x}$ in $\mathcal{D}_{1}$ and
$u_{0}=\breve{U}p_{0x}$ in $\mathcal{D}_{2}$, where $\tilde{U}$, $\bar{U}$ and
$\breve{U}$ are given by (3.47a, b)–(3.51).

At $O(\epsilon^{2})$, the first-order strong cross-channel diffusion equations
(3), (3.55a–ca, b) and (3.58)–(3.59a–c) are recovered. On $\mathcal{I}$,

$c_{1y}=\mathscr{D}(c_{0}-\Gamma_{0}),\quad\Pen^{-1}_{I}\Gamma_{1zz}-\Gamma_{0}w_{1z}=-\Pen^{-1}_{I}\Gamma_{0xx}+(u_{0}\Gamma_{0})_{x}-\mathscr{B}(c_{0}-\Gamma_{0}).$
(3.71a, 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 (2) into the velocity flux condition
(3.42), we recover $\tilde{Q}p_{0x}+\Ma\bar{Q}\Gamma_{0x}=1$ in
$\mathcal{D}_{1}$ and $\breve{Q}p_{0x}=1$ in $\mathcal{D}_{2}$, where
$\tilde{Q}$ is given by (2a), $\bar{Q}$ is given by (2b) and
$\breve{Q}=-4P_{z}/3$. Substituting the streamwise velocity (2) into (3.70) we
derive the steady leading-order surfactant transport equations (3.32).
Substituting (2) into the surfactant transport equations (3.69) (or adding
(3.32a) and (3.32b) and integrating) we recover integral constraints on the
surfactant transport equations (3.32), which are given by surfactant flux
constraints (3.33). The coefficients $\alpha$, $\beta$, $\gamma$ and $\delta$
are defined in (3.63) and

$\nu=\frac{4\phi_{z}P_{z}\Da}{\Pen}=\frac{4\hat{K}_{a}\phi_{z}\hat{P}_{z}}{\hat{H}\hat{U}}\quad\text{(surface
exchange)}.$ (3.72)

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.33a, b), continuity of bulk surfactant
flux between $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ (3.33c) and no flux of
interfacial surfactant (3.33d). For $\nu/\epsilon^{2}\gg
O(1,\,\alpha,\,\beta,\,\delta)$ then $c_{0}=\Gamma_{0}$ at leading-order,
recovering the transport equations (3.62) in the strong cross-channel
diffusion and strong exchange limit. The numerical method used to solve the 1D
surfactant transport equations (3.32)–(3.33) 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 $\mathcal{D}_{1}$ and
$\breve{Q}p_{0x}=1$ in $\mathcal{D}_{2}$, such that integrating $-p_{0x}$
across the period gives

$\Delta
p_{U}=-\int_{x=-\phi_{x}}^{\phi_{x}}p_{0x}\,\text{d}x-\int_{x=\phi_{x}}^{2-\phi_{x}}p_{0x}\,\text{d}x=-\frac{2\phi_{x}}{\tilde{Q}}-\frac{2(1-\phi_{x})}{\breve{Q}}.$
(3.73)

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

$\Delta
p_{0}=-\int_{x=-\phi_{x}}^{\phi_{x}}p_{0x}\,\text{d}x-\int_{x=\phi_{x}}^{2-\phi_{x}}p_{0x}\,\text{d}x=-\frac{2\phi_{x}}{\tilde{Q}}-\frac{2(1-\phi_{x})}{\breve{Q}}+\frac{\Ma\bar{Q}\Delta\Gamma_{0}}{\tilde{Q}},$
(3.74)

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.63b, 3.63c) and substituting $\Delta
p_{U}$, $\Delta p_{R}$ and $\Delta p_{0}$ into the leading-order drag
reduction (3.46), we have

${DR}_{0}=1-\frac{\gamma\Delta\Gamma_{0}}{2\phi_{x}\beta}.$ (3.75)

In (3.75), the first term gives the drag reduction when the interface is
shear-free and the second term measures the impact of surfactant. In order to
derive the expression in (3.75), we used the fact that
$\breve{Q}\bar{Q}\tilde{q}/((\tilde{q}\bar{Q}-\bar{q}\tilde{Q})(\breve{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.36b) for an
equivalent channel with the mixed boundary conditions (3.37a–c, 3.43) replaced
by $\lambda_{e}u_{0y}-u_{0}=0$ on $\mathcal{I}$, $\mathcal{R}$ and
$\mathcal{S}$. We obtain $u_{0}=\check{U}\Delta p_{0x}$ where
$\check{U}=y(y-2)/2-\lambda_{e}$ and $p_{0x}$ is the same pressure gradient as
in the SHS channel. The flux is
$\check{Q}=\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2}\check{U}p_{0x}\text{d}y\text{d}z=(\breve{Q}-2P_{z}\lambda_{e})p_{0x}$,
or by integrating over one period
$\check{Q}=(\breve{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\bar{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}))},$ (3.76)

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.67, 3.33a,
3.33b), 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.32)–(3.33), 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.67) 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 $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$
represent the subdomains over the plastron and solid ridge, respectively,
($\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ are separated by the vertical dotted
line at $x=\phi_{x}=0.5$), and $c_{0}=\Gamma_{0}$ in $\mathcal{D}_{1}$.

(a)

(b)

(c)

(d)

(e)

Figure 4:  The leading-order drag reduction (${DR}_{0}$) and surfactant
distribution ($c_{0}=\Gamma_{0}$) in the strong-exchange problem, for
$\beta=1$, $\phi_{x}=0.5$, $\epsilon=0.1$, $\phi_{z}=0.5$, $\Da=1$ and
$P_{z}=0.5$, computed using (3.67). (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.77), (4.83)
and (4.81) 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}\approx 0.5$ are given by (4.79, 4.80). (b–e) Plots
of $c_{0}$ for varying surfactant strength ($\gamma$) and bulk diffusion
($\alpha$), where the asymptotic curves are: (b) $\alpha=\delta=1$, A: $--$
(4.84), M: $\cdot\cdot\cdot$ (4.78) with $\gamma=100$, (c) $\gamma=0.1$, A:
$--$ (4.84) with $\alpha=\delta=0.01$, D: $\cdot\cdot\cdot$ (4.82) with
$\alpha=\delta=10$, (d) $\alpha=\delta=0.01$, A: $--$ (4.84), M:
$\cdot\cdot\cdot$ (4.78) with $\gamma=100$, and (e) $\gamma=100$, M: $--$
(4.78) with $\alpha=\delta=1$, D: $\cdot\cdot\cdot$ (4.82) with
$\alpha=\delta=1000$. The star identifies the point of the
$(\alpha,\,\gamma)$-plane where we examine the flow field in §4.1.2.

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}_{0}$ 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)}{\left(E-1\right)}\right)\quad\text{for}\quad\gamma\gg\max\left(1,\,\alpha,\,\beta,\,\delta,\,\frac{\epsilon^{2}}{\alpha}\right),$
(4.77)

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: $\text{M}_{\text{D}}$, in which Marangoni
effects start to compete with bulk diffusion, and $\text{M}_{\text{G}}$, in
which Marangoni effects start to compete with shear dispersion. In the
$\text{M}_{\text{D}}$ case, (4.77) 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 immobilizing effects of surfactant. In the $\text{M}_{\text{G}}$
case, (4.77) 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 $\mathcal{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}\quad\gamma\gg\max\left(1,\,\alpha,\,\beta,\,\delta,\,\frac{\epsilon^{2}}{\alpha}\right).$
(4.78)

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.78), the magnitude of the concentration in
$\mathcal{D}_{1}$ is weakly regulated by fluxes driven by bulk diffusion or
shear dispersion in domain $\mathcal{D}_{2}$, via the factor
$\alpha+s_{4}\epsilon^{2}/\alpha$ in $E$ that appears due to continuity of
concentration across $\mathcal{D}_{1}$ and $\mathcal{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}^{\prime}$) 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
region 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.77) and (4.78) 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}\quad\gamma=\frac{\epsilon^{2}(2s_{1}(1-\phi_{x})+s_{4}\phi_{x}(2+\beta))}{2\alpha\left(1-\phi_{x}\right)}\quad\text{for}\quad\alpha\ll
1,$ (4.79)

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\rightarrow 0$. As the immobilizing 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}\quad\gamma=\delta+\frac{\alpha(2+\beta\phi_{x})}{2(1-\phi_{x})}\quad\text{for}\quad\min(\alpha,\,\delta)\gg
1,$ (4.80)

which gives the rightmost dashed magenta line in figure 4(a). The asymptotic
approximation (4.80) agrees with the numerical solution of (3.67) (thick black
line in figure 4a) for $\min(\alpha,\,\delta)\gg 1$. Similar to the GM
transition, the immobilizing 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. Since 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}\quad\min(\alpha,\,\delta)\gg\max(1,\,\gamma),$
(4.81)

for $\beta=O(1)$. Equation (4.81) 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 to
bulk advection and Marangoni effects. The leading-order solution in
$\mathcal{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}\quad\min(\alpha,\,\delta)\gg\max(1,\,\gamma),$
(4.82)

demonstrating how diffusion eliminates gradients of $c_{0}$ throughout
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ (the weak gradients contributing to
${DR}_{0}$ in (4.81) 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).

(a) (b)

(c) (d)

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.67). 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.77), (4.83) and (4.81) for the M, A and D
regions, respectively. The dashed magenta lines describing when
${DR}_{0}\approx 0.5$ are given by (4.85) and (4.80) 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.84), M:
$\cdot\cdot\cdot$ (4.78) with $\gamma=1000$, and (d) $\beta=100$, A: $--$
(4.84), M: $\cdot\cdot\cdot$ (4.78) with $\gamma=10000$.

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}_{0}$ is
once again close to the shear-free value

${DR}_{0}\approx
1-\frac{\gamma}{2\phi_{x}\left(\beta+1\right)}\quad\text{for}\quad\epsilon^{2}\ll\alpha\ll
1,\quad\gamma\ll\min(1,\,\beta).\\\ $ (4.83)

Region A can also be divided into two sub-regions: $\text{A}_{\text{G}}$ which
balances advection and shear dispersion, and $\text{A}_{\text{D}}$ which
balances advection and diffusion; these balances give the same ${DR}_{0}$ as
in (4.83). The corresponding concentration profiles in region A are shown in
figure 4(b, c, d). The leading-order solution in $\mathcal{D}_{1}$ 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)\
\text{for}\ \epsilon^{2}\ll\alpha\ll 1,\ \gamma\ll\min(1,\,\beta).$ (4.84)

The downstream boundary layer in (4.84) can be short compared to 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}_{0}$ 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)$.

(a) (b)

(c) (d)

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$% normalized 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.

We next investigate the dependence of the drag reduction on the partition
coefficient ($\beta$) by evaluating ${DR}_{0}$ 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,
whilst 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 $\mathcal{D}_{1}$. A number of the asymptotic
approximations given in (4.77)–(4.84) simplify for $\beta\gg 1$: (4.79)
reduces to ${DR}_{0}\approx 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 $\mathcal{D}_{2}$ determines
${DR}_{0}$; (4.80) reduces to ${DR}_{0}\approx 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}_{0}$.

Another feature that emerges from figure 5(a, b) is the flattening of the
central contours of ${DR}_{0}$ 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 ${DR}_{0}\approx 1-\gamma/(2\beta\phi_{x})$ for
$\gamma/\beta\leq 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 immobilizing 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}\quad\gamma=\phi_{x}\beta\quad\text{for}\quad\beta\gg
1\gg\max\left(\alpha,\,\delta,\,\frac{\epsilon^{2}}{\alpha}\right),$ (4.85)

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\leq 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}\displaystyle 0\ \text{for}\ -\phi_{x}\leq x\leq
x_{0},\\\ \displaystyle\frac{\beta}{\gamma}(x-\phi_{x})+1\ \text{for}\
x_{0}\leq x\leq\phi_{x},\end{cases}\text{for}\ \min(\beta,\,\gamma)\gg
1\gg\max\left(\alpha,\,\delta,\,\frac{\epsilon^{2}}{\alpha}\right),$ (4.86)

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 5d demonstrates that $c_{0}\approx 0$ throughout most of
$\mathcal{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.85) when $\gamma/\beta\leq
2\phi_{x}$.

####  4.1.2 Flow field

Whenever there is partial drag reduction, $0\leq{DR}_{0}<1$, the surfactant
has a non-uniform concentration in $\mathcal{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, 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
$\mathcal{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}\leq z\leq\phi_{z}$ and $\phi_{z}\leq z\leq P_{z}$), whilst allowing
for slip over the plastron ($-\phi_{z}\leq z\leq\phi_{z}$). The slip velocity
at $\mathcal{I}$ reduces as one progresses through $\mathcal{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 $\mathcal{D}_{1}$
along the positive $x$-direction as the interface is immobilized. 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
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ meet.

###  4.2 Moderate exchange

We now turn our attention to the moderate-exchange problem (3.32)–(3.33), in
which the bulk and interfacial concentration fields decouple for
$\nu/\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 $\nu$ in (3.72). 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

(a)

(b)

(c)

(d)

(e)

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.32)–(3.33). (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 ($\text{M}_{\text{E}}$) and
advection-diffusion-exchange ($\text{AD}_{\text{E}}$) regions are separated by
black lines and ${DR}_{0}$ is approximated by (4.77), (4.83), (4.87) and
(4.91) in the M, AD, $\text{M}_{\text{E}}$ and $\text{AD}_{\text{E}}$ regions,
respectively. The dashed magenta lines describing when ${DR}_{0}\approx 0.5$
are given by (4.80, 4.89). 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: $--$ (4.78) with $\gamma=100$, M:
$\cdot\cdot\cdot$ (4.84), (c) $\gamma=0.03$, $\text{AD}_{\text{E}}$: $--$
(4.92), AD: $\cdot\cdot\cdot$ (4.84), (d) $\nu=0.001$, $\text{AD}_{\text{E}}$:
$--$ (4.92), $\text{M}_{\text{E}}$: $\cdot\cdot\cdot$ (4.88) with $\gamma=10$,
and (e) $\gamma=100$, $\text{M}_{\text{E}}$: $--$ (4.88), M: $\cdot\cdot\cdot$
(4.78).

Figure 7(a) shows how ${DR}_{0}$ varies with surface exchange strength ($\nu$)
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.32)–(3.33) 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
$\nu/\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,
$\nu/\epsilon^{2}\gg\max(1,\,\alpha,\,\beta,\,\delta)$, we reach the AD
boundary and sub-region $\text{M}_{\text{D}}$ identified in figure 4(a).
Accordingly, the regimes AD and $\text{M}_{\text{D}}$ are also identified in
figure 7(a) where bulk–surface exchange is strong. New regimes
$\text{M}_{\text{E}}$ and $\text{AD}_{\text{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 $\text{M}_{\text{E}}$, $\text{AD}_{\text{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 $\mathcal{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 below.

We begin our exploration of the parameter space in the moderate-exchange
problem (figure 7) with region $\text{M}_{\text{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}\quad\gamma\gg\max(1,\,\alpha,\,\beta,\,\delta),\quad\frac{\nu}{\epsilon^{2}}\ll\min(1,\,\alpha,\,\beta,\,\delta).$
(4.87)

Here, the immobilizing effects of surfactant are weakened by strong
interfacial diffusion
($\delta/\gamma\propto\hat{\mu}\hat{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
$\text{M}_{\text{E}}$ is given by

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

The gradient of $\Gamma_{0}$ in region $\text{M}_{\text{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.78) in M and the strong-exchange
problem of §4.1.1, where $c_{0}=\Gamma_{0}$ and $\Gamma_{0}<1$ for all $x$, in
$\text{M}_{\text{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.29) 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 $\text{M}_{\text{E}}$ gives way to region $\text{M}_{\text{D}}$ as the
exchange strength increases. It is notable that ${DR}_{0}\approx\delta/\gamma$
(see (4.87)) is regulated mainly by surface diffusion in region
$\text{M}_{\text{E}}$, whilst both bulk and surface diffusion control $DR_{0}$
in region $\text{M}_{\text{D}}$, 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 $\text{M}_{\text{E}}$, which imposed $\Gamma_{0}=c_{0}$ in region
$\text{M}_{\text{D}}$. 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),
since $\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 $\text{M}_{\text{E}}$ region
transitions into the diffusion-exchange-dominated ($\text{D}_{\text{E}}$)
region across the $\text{M}_{\text{E}}\text{D}_{\text{E}}$ boundary for
$\min(\alpha,\,\delta)\gg O(1)$ (figure 8(b)), whereas the
$\text{M}_{\text{E}}$ region transitions into the advection-exchange-dominated
($\text{A}_{\text{E}}$) region across the
$\text{M}_{\text{E}}\text{A}_{\text{E}}$ boundary for
$\max(\alpha,\,\delta)\ll O(1)$ (figure 8(a)). At the
$\text{M}_{\text{E}}\text{D}_{\text{E}}$ boundary, the asymptotic analysis in
Appendix D.2 shows that

${DR}_{0}\approx
0.5\quad\text{when}\quad\gamma=\delta\quad\text{for}\quad\min(\alpha,\,\delta)\gg
1,\quad\frac{\nu}{\epsilon^{2}}\ll 1,$ (4.89)

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\rightarrow 0$. The range of validity for the asymptote given in (4.89) is
extended to $\delta=O(1)$ to include the
$\text{M}_{\text{E}}\text{AD}_{\text{E}}$ boundary, for which we give partial
justification in Appendix D.3. At the $\text{M}_{\text{E}}\text{A}_{\text{E}}$
boundary, Appendix D.4 demonstrates that

${DR}_{0}\approx
0.5\quad\text{when}\quad\gamma=\frac{\beta\phi_{x}}{4}+\frac{\delta}{2}\quad\text{for}\quad\max(\alpha,\,\delta)\ll
1,\quad\frac{\nu}{\epsilon^{2}}\ll 1,$ (4.90)

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\rightarrow 0$. Agreement improves for smaller $\delta$.
It is notable that the threshold (4.90) differs from the strong-exchange limit
(4.85) by a factor of 4.

(a) (b)

(c) (d)

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.32)–(3.33). 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 ($\text{M}_{\text{E}}$), advection-exchange
($\text{A}_{\text{E}}$), diffusion-exchange ($\text{D}_{\text{E}}$) regions
are separated by black lines and ${DR}_{0}$ is approximated by (4.77), (4.83),
(4.81), (4.87), (4.91) and (4.91) in M, A, D, $\text{M}_{\text{E}}$,
$\text{A}_{\text{E}}$ and $\text{D}_{\text{E}}$, respectively. The dashed
magenta lines describing when ${DR}_{0}\approx 0.5$ are given by (4.80) in
(a,b) (right curves), (4.89) in (b) and (4.90) (a). Plots of $\Gamma_{0}$ and
$c_{0}$ for varying surfactant strength ($\gamma$), where (c)
$\alpha=\delta=0.1$ and $\nu=10^{-4}$, $\text{A}_{\text{E}}$: $--$ (4.92),
$\text{M}_{\text{E}}$: $\cdot\cdot\cdot$ (4.88) with $\gamma=10$, and (d)
$\alpha=\delta=10$ and $\nu=0.01$, $\text{D}_{\text{E}}$: $--$ (4.93),
$\text{M}_{\text{E}}$: $\cdot\cdot\cdot$ (4.88) with $\gamma=1000$. The star
identifies the point of the $(\nu,\,\gamma)$-plane where we examine the flow
field in §4.2.2.

We next analyse the drag reduction in both the $\text{A}_{\text{E}}$ and
$\text{D}_{\text{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 $\text{A}_{\text{E}}$ and $\text{D}_{\text{E}}$ is
close to the shear-free value

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

Equation (4.91) 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
$\text{A}_{\text{E}}$ and $\text{D}_{\text{E}}$, the leading-order solution on
$\mathcal{I}$ being given by

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

From (4.92), 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.29), condensing
the Marangoni effect to a small region near $x=\phi_{x}$. The streamwise
velocity flows over a almost shear-free boundary at the upstream end of the
interface. When surface diffusion is strong, $\delta\gg O(1)$, (4.92) reduces
to

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

Noting that
$\beta/\delta\propto\hat{H}^{2}\hat{U}/(\hat{D}_{I}\phi_{z}\hat{P}_{z})$,
(4.93) 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.91). In figure 8(a, b), the horizontal
boundary of the $\text{M}_{\text{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
$\text{A}_{\text{E}}$ to $\text{D}_{\text{E}}$ in figure 8(a, b).

####  4.2.2 Flow field

(a) (b)

(c) (d)

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%
normalized 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.

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

(a) (b)

(c) (d)

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 panel a and c),
computed using (3.67) (red) and (4.78) (dashed and black). (c) Scatter plot of
${DR}_{0}$ using our theory ((3.75) and (3.32)–(3.33)) and the 159 numerical
simulations detailed in Temprano-Coleto et al. (2023), the colorbar 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.

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$). This amounts to reevaluating $\tilde{U}$, $\bar{U}$ and $\breve{U}$,
as well as modifying the surface surfactant flux term in (3.43) 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}\,$m, $\hat{P}_{z}=6\times 10^{-5}\,$m and
$\phi_{z}=2/3$, from which we can evaluate $P_{z}$, $\tilde{Q}$, $\bar{Q}$,
$\tilde{q}$ and $\bar{q}$ using (2, 3.53a, b). Using (3.63) (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 $\Pen=\Pen_{I}=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.67) subject to (3.33a, b) (red solid line), and
the (indistinguishable) asymptotic solution in M (4.78) (black dashed line).
The gradient of $c_{0}$ at the centre of the plastron computed using our
theory ($\text{d}c_{0}/\text{d}x(0)\approx 0.19$) and the numerical simulation
($\text{d}c/\text{d}x(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 colorbar 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 colorbar 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 superhydrophobic surfaces (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.32)–(3.33) 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).

Strong exchange $\displaystyle\frac{\nu}{\epsilon^{2}}\gg\max(1,\,\alpha,\,\beta,\,\delta)$ |
---|---
Region | Parameter space | ${DR}_{0}$ | Regime
M | $\displaystyle\gamma\gg\max\left(1,\,\alpha,\,\beta,\,\delta,\,\frac{\epsilon^{2}}{\alpha}\right)$ | $\displaystyle\frac{1}{\gamma}\left(\alpha+\delta+\frac{\epsilon^{2}s}{\alpha}+\frac{\phi_{x}(E+1)}{\left(E-1\right)}\right)$ | ${DR}_{0}\ll 1$
$\text{M}_{\text{D}}$ | $\displaystyle\gamma\gg(\alpha,\,\delta)\gg\max\left(1,\,\beta,\,\frac{\epsilon^{2}}{\alpha}\right)$ | $\displaystyle\frac{\alpha}{\gamma(1-\phi_{x})}+\frac{\delta}{\gamma}$ | ${DR}_{0}\ll 1$
$\text{M}_{\text{G}}$ | $\displaystyle\gamma\gg\frac{\epsilon^{2}}{\alpha}\gg\max(1,\,\alpha,\,\beta,\,\delta)$ | $\displaystyle\frac{\epsilon^{2}(s_{4}\phi_{x}+s(1-\phi_{x}))}{\gamma\alpha(1-\phi_{x})}$ | ${DR}_{0}\ll 1$
A | $\epsilon^{2}\ll\alpha\ll 1,\quad\gamma\ll\min(1,\,\beta)$ | $\displaystyle 1-\frac{\gamma}{2\phi_{x}\left(\beta+1\right)}$ | $1-{DR}_{0}\ll 1$
D | $\min(\alpha,\,\delta)\gg\max(1,\,\gamma)$ | $\displaystyle 1-\frac{(1-\phi_{x})\gamma}{(1+\phi_{x}\beta)\alpha+(1-\phi_{x})\delta}$ | $1-{DR}_{0}\ll 1$
Weak exchange $\displaystyle\frac{\nu}{\epsilon^{2}}\ll\min(1,\,\alpha,\,\beta,\,\delta)$ |
Region | Parameter space | ${DR}_{0}$ | Regime
$\text{M}_{\text{E}}$ | $\gamma\gg\max(1,\,\alpha,\,\beta,\,\delta)$ | $\displaystyle\frac{\delta}{\gamma}$ | ${DR}_{0}\ll 1$
$\text{A}_{\text{E}}$, $\text{D}_{\text{E}}$ | $\gamma\ll\min(1,\,\alpha,\,\beta,\,\delta)$ | $\displaystyle 1-\frac{\gamma}{\delta}$ | $1-{DR}_{0}\ll 1$

Table 1: Summary of the asymptotic predictions of the leading-order drag reduction $DR_{0}$ in the main asymptotic regions analysed in the strong-exchange problem with: the Marangoni-dominated region (M) with sub-regions $\text{M}_{\text{D}}$ and $\text{M}_{\text{G}}$, the advection-dominated region (A) and the diffusion-dominated region (D); and their analogues in the weak-exchange problem: the $\text{M}_{\text{E}}$, $\text{A}_{\text{E}}$ and $\text{D}_{\text{E}}$ regions. The drag reduction $DR_{0}$ is expressed in terms of the transport coefficients $\alpha$, $\beta$, $\gamma$, $\delta$ and $\nu$ given in (3.63) and (3.72) and constants ($s_{1}>0$, $s_{2}$, $s_{3}$ and $s_{4}>0$) given in (3.64) where $s\equiv s_{1}+s_{2}+s_{3}$ and $E\equiv\exp(2\alpha(1-\phi_{x})/(\alpha^{2}+\epsilon^{2}s_{4}))$. Quantity | $\displaystyle\frac{\gamma}{\alpha}$ | $\displaystyle\frac{\alpha\gamma}{\epsilon^{2}}$ | $\displaystyle\frac{\gamma}{\beta}$ | $\displaystyle\frac{\gamma}{\alpha\beta}$ | $\displaystyle\frac{\gamma}{\delta}$
---|---|---|---|---|---
Proportional to | $\displaystyle\frac{\hat{A}\hat{C}\hat{L}_{d}^{2}\hat{H}}{\hat{\mu}\hat{D}\hat{P}_{z}}$ | $\displaystyle\frac{\hat{A}\hat{C}\hat{L}_{d}^{2}\hat{D}\hat{H}\hat{P}_{z}}{\hat{\mu}\hat{Q}^{2}}$ | $\displaystyle\frac{\hat{A}\hat{C}\hat{L}_{d}\hat{H}^{3}}{\hat{\mu}\hat{Q}\hat{P}_{x}}$ | $\displaystyle\frac{\hat{A}\hat{C}\hat{L}_{d}\hat{H}^{2}}{\hat{\mu}\hat{D}\hat{P}_{z}}$ | $\displaystyle\frac{\hat{A}\hat{C}\hat{L}_{d}\hat{H}^{2}}{\hat{\mu}\hat{D}_{I}\phi_{z}\hat{P}_{z}}$
Table 2: Summary of the dimensionless ratios appearing in table 1 that affect
the leading-order drag reduction, and their dependence on the dimensional
quantities characterising the flow and surfactant properties and the geometry
(outlined in §2).

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 to 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.67) 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 ($\text{M}_{\text{D}}$), shear dispersion
($\text{M}_{\text{G}}$) or surface diffusion ($\text{M}_{\text{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 $\hat{D}_{I}$ instead decreases the drag across the
$\text{D}_{\text{E}}\text{M}_{\text{E}}$ boundary (figure 8). The velocity
flux $\hat{Q}$ decreases the drag across the AM and GM boundaries,
quadratically in the latter case (table 2). The approximations of the leading-
order drag reduction in table 2 can also be divided into those with a linear
or quadratic dependence on $\hat{L}_{d}$. Recall from §2 that
$\hat{L}_{d}/\hat{H}$ is the normalised surfactant depletion length in §2,
such that for $\hat{L}_{d}/\hat{H}\ll 1$ surfactant is essentially insoluble,
and for $\hat{L}_{d}/\hat{H}\gg 1$ nearly all surfactant is adsorbed to the
interface. Accordingly, bulk diffusion and shear dispersion reduce drag
through parameters that are quadratic in $\hat{L}_{d}$, requiring solubility
for them to be effective.

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 highlights 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 minimize drag in realistic 3D SHS microchannels and other applications,
where surfactant traces may be naturally present.

## Acknowledgements

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.

##  Appendix A Numerical methods

###  A.1 Transport equations

We solve the transport equations, (3.32, 3.33, 3.67), subject to the boundary
conditions, (3.33), using Chebyshev collocation. An in-depth discussion of the
$(N+1)\times(N+1)$ differentiation matrix, $\boldsymbol{D}$, that forms the
basis of Chebyshev collocation technique can be found in Trefethen (2000).
Briefly, it allows us to approximate $f_{\xi}$ using
$\boldsymbol{D}\boldsymbol{f}$, where $\boldsymbol{f}=(f(0),\,...,f(N))^{T}$
is the solution vector on a grid defined by $\xi(i)=\cos(i\pi/N)$ for
$i=0,\,1,\,...,\,N$. We map domains $\mathcal{D}_{1}$, for
$x\in[-\phi_{x},\,\phi_{x}]$, and $\mathcal{D}_{2}$, for
$x\in[\phi_{x},\,2-\phi_{x}]$, to a discrete space $\xi(i)\in[-1,\,1]$ for
$i=0,\,1,\,...,\,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., Appendix C.5 and D.5). We concatenate the solution
in both domains as
$\boldsymbol{c}_{0}=(c_{0,\,1}(0),\,...,\,c_{0,\,1}(N),\,c_{0,\,2}(0),\,...,\,c_{0,\,2}(N))^{T}$,
where $\boldsymbol{c}_{0,\,i}=(c_{0,\,i}(0),\,...,\,c_{0,\,i}(N))^{T}$ is the
solution in $\mathcal{D}_{i}$ for $i=1,\,2$. The initial guess evolves to the
nonlinear state by solving

$\displaystyle\mathsfbi{L}_{1}(\boldsymbol{c}_{0}^{\text{old}})\boldsymbol{c}_{0,\,1}(1:N-1)=\boldsymbol{1}\quad\text{where}\quad\boldsymbol{c}_{0,\,1}(0)=\boldsymbol{c}_{0,\,2}(N),$
(A.94a)
$\displaystyle\mathsfbi{L}_{2}(\boldsymbol{c}_{0}^{\text{old}})\boldsymbol{c}_{0,\,2}(1:N-1)=\boldsymbol{1}\quad\text{where}\quad\boldsymbol{c}_{0,\,2}(0)=\boldsymbol{c}_{0,\,1}(N),$
(A.94b)

where $\mathsfbi{L}_{i}$ are discrete approximations to (3.67, 3.32, 3.33) in
$\mathcal{D}_{i}$ for $i=1,\,2$. We update the solution via
$\boldsymbol{c}_{0}^{\text{new}}=\boldsymbol{c}_{0}^{\text{old}}+r(\boldsymbol{c}_{0}-\boldsymbol{c}_{0}^{\text{old}})$
where $r$ is a relaxation factor, until
$\tau_{i}=\|\mathsfbi{L}_{i}(\boldsymbol{c}_{0,\,i}^{\text{new}}-\boldsymbol{c}_{0,\,i}^{\text{old}})\|_{\infty}/r$
falls below a specified tolerance for $i=1,\,2$.

###  A.2 Leading-order flow field

Problems (3.47a, b, 3.49) were solved by modifying the framework introduced in
Game et al. (2017). We consider the domain
$\mathcal{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 PDEs are transformed
to discrete space using techniques outlined in Trefethen (2000). The
discontinuous boundary conditions at $z=\phi_{z}$, (3.47a, bb–e)–(3.49b–e),