# Unsteady evolution of slip and drag 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

Recognising that surfactants may impede the drag reduction resulting from
superhydrophobic surfaces (SHSs), and that surfactant concentrations can
fluctuate in space and time, we model the unsteady transport of soluble
surfactant in a channel bounded by two SHSs. The flow is laminar, pressure-
driven, and the SHSs are periodic in the streamwise and spanwise directions.
We assume that the channel length is much longer than the streamwise period,
the streamwise period is much longer than the channel height and spanwise
period, and bulk diffusion is sufficiently strong for cross-channel
concentration gradients to be small. By combining long-wave and homogenisation
theories, we derive an unsteady advection–diffusion equation for surfactant
flux transport over the length of the channel, which is coupled to a quasi-
steady advection–diffusion equation for surfactant transport over individual
plastrons. As diffusion over the length of the channel is typically small, the
leading-order surfactant flux is governed by a nonlinear advection equation
that we solve using the method of characteristics. We predict the propagation
speed of a bolus of surfactant and describe its nonlinear evolution via
interaction with the SHS. The propagation speed can fall significantly below
the average streamwise velocity as the surfactant adsorbs and rigidifies the
plastrons. Smaller concentrations of surfactant are therefore advected faster
than larger ones, so that wave-steepening effects can lead to shock formation
in the surfactant-flux distribution. These findings reveal the spatio–temporal
evolution of the slip velocity and enable prediction of the dynamic drag
reduction and effective slip length in microchannel applications.

###### keywords:

Marangoni convection, drag reduction, microfluidics

##  1 Introduction

Surfactants are chemical compounds that are advected and diffuse throughout a
fluid, where they adsorb onto liquid–liquid or liquid–gas interfaces
(Manikantan & Squires, 2020). They have been shown to impair the effective
slip length and drag reduction in superhydrophobic microchannels (Peaudecerf
et al., 2017). Superhydrophobic surfaces (SHSs) use chemically-coated
microscopic structures to suspend a fluid over a series of gas pockets (Lee et
al., 2016). The combination of no-slip structures and shear-free liquid–gas
interfaces generates the drag reduction for microchannel flows. Hence, SHSs
have been considered for applications in biofluidics (Darmanin & Guittard,
2015), heat transfer (Lam et al., 2015) and marine hydrodynamics (Xu et al.,
2020), both in laminar and turbulent flows. Field studies have shown that
surfactant is present in the ocean and that the surfactant concentration can
vary significantly in space and time (Pereira et al., 2018; Frossard et al.,
2019). Traces of surfactant have been measured in rivers, estuaries and fog
(Lewis, 1991; Facchini et al., 2000). They are also present in most industrial
and laboratory environments (Manikantan & Squires, 2020). Surfactants that
have been absorbed onto the liquid–gas interfaces of SHSs are advected
downstream by the flow and accumulate at stagnation points (i.e. liquid–solid
contact lines), generating an adverse Marangoni force that may negate any
drag-reducing effects in laminar (Peaudecerf et al., 2017) or turbulent
(Tomlinson et al., 2023b) flows. In order to better understand how
environmental surfactants may compromise drag reduction by SHSs in
applications, this paper addresses unsteady surfactant transport in a laminar
pressure-driven channel flow bounded between streamwise- and spanwise-periodic
SHSs, investigating how surfactant is advected and diffuses over length scales
and time scales that are large compared to the dimensions of the SHS texture.

Experimental studies first suggested that naturally-occurring surfactants
could affect channel flows bounded by SHSs comprising spanwise ridges (Kim &
Hidrovo, 2012), as well as finite-length streamwise ridges (Bolognesi et al.,
2014). They found that the flow rate and wall shear stress closely resembled a
channel with solid walls, and thus their SHSs offered only a modest drag
reduction (Kim & Hidrovo, 2012; Bolognesi et al., 2014). Schäffel et al.
(2016) showed that experimentally-measured slip lengths on SHSs consisting of
pillars were much smaller than predicted by surfactant-free simulations; this
was true whether surfactant was explicitly added or not, suggesting that
naturally-occurring surfactants played a key role. As noted earlier, a
requirement for surfactant effects to manifest on SHSs is the presence of
stagnation points perpendicular to the flow, at which surfactant can
accumulate to generate a surface tension gradient. Song et al. (2018) showed
that surface tension gradients emerged in their experiments for finite
streamwise ridges, increasing the drag compared to those configurations with
concentric ridges that lack stagnation points perpendicular to the flow.

To investigate the effect of weak surfactant concentrations, Peaudecerf et al.
(2017) introduced simulations inclusive of surfactant dynamics; they showed
that a plastron interface could be immobilised by concentrations below levels
commonly occurring in the environment and in engineered systems. For example,
Hourlier-Fargette et al. (2018) showed that polydimethylsiloxane (PDMS)
surfaces can release uncrosslinked oligomers that behave as surfactants. The
simulations of Peaudecerf et al. (2017) also predicted that surfactant
impairment would decrease as the streamwise plastron interface length
increased; this was confirmed by their experiments (Peaudecerf et al., 2017).
These experiments also showed that, if the driving pressure was suddenly
removed, a reverse flow established at the interface, decaying with time as
$1/t$ at intermediate times. This time scaling was predicted by a similarity
solution driven by surfactant relaxation, assuming advection-dominated flow.
In contrast to these plastron-scale findings, there is presently no theory
that includes the combined effects of solubility, advection and diffusion,
that describes inhomogeneous surfactant transport across multiple plastrons,
or that can model the effects of unsteady surfactant concentration at the
inflow.

Steady scaling theories were constructed for a pressure-driven channel flow
with two-dimensional gratings (approximating long, spanwise-oriented gratings
Landel et al., 2020), as well as for long gratings with finite spanwise
extent, assuming spatially periodic flow (Temprano-Coleto et al., 2023). Both
theories are in agreement with the slip velocity and drag predicted in full
numerical simulations. Temprano-Coleto et al. (2023) further validated their
theory by performing experiments with SHS gratings of various lengths, finding
that surfactant effects decrease with the square of the interface length.
These theories assume that the surfactant concentration is small (as may be
expected when surfactant is not explicitly added) and that the shear stress is
approximately uniform at the liquid-gas interface. They do not consider the
stagnant cap regime, first exhibited in air bubbles rising in surfactant-
contaminated water (Bond & Newton, 1928; Frumkin & Levich, 1947). To examine
effects of non-uniform shear stresses at the liquid–gas interface, Tomlinson
et al. (2023a) considered gratings of finite spanwise extent, and assumed that
bulk diffusion was strong enough to suppress cross-channel concentration
gradients, allowing systematic asymptotic approximations to be developed.
Several dimensionless groups were identified by Temprano-Coleto et al. (2023)
and Tomlinson et al. (2023a) that influence the drag in superhydrophobic
channels. Temprano-Coleto et al. (2023) showed that surfactant impairment in
their simulations and experiments was well predicted by a single dimensionless
group, when the surfactant properties, SHS dimensions and flow velocities are
constrained within physically realizable ranges. Using these physical
constraints, scaling analysis provided the dimensionless group as the ratio
between the streamwise length of the interface and a surfactant-determined
lengthscale, labelled “mobilization length”. Without constraints on surfactant
or flow properties, Tomlinson et al. (2023a) found several other relevant
dimensionless groups by calculating asymptotic solutions for the concentration
field and drag across the whole parameter space; these depend on a velocity
scale generated by interfacial Marangoni effects, the surfactant diffusivity
and the flow rate. Another dimensionless group found by Tomlinson et al.
(2023a) can be used to predict whether the surfactant concentration field is
in the stagnant-cap regime.

Both Temprano-Coleto et al. (2023) and Tomlinson et al. (2023a) assumed that
the velocity and bulk concentration fields are steady and spatially periodic.
That is, they did not allow surfactant to enter the channel with a non-uniform
distribution that varies in space and time (as exhibited in Frossard et al.,
2019) over multiple periods. To address this case, we use multiscale
homogenisation techniques (such as those outlined in Bottaro, 2019) to study
the time- and space-varying effects of surfactant over the whole SHS, without
needing to numerically resolve the small details over each texture period,
which can be computationally very expensive. We show how a time-dependent one-
dimensional asymptotic theory, derived from the three-dimensional Stokes and
surfactant transport equations, can be adapted to describe the unsteady
evolution of slip and drag in a laminar pressure-driven channel flow with
streamwise- and spanwise-periodic grooves, allowing for time-dependent
distributions of surfactant flux at the channel inlet. The problem exhibits
multiple length and time scales, which we exploit to derive and solve a quasi-
steady advection–diffusion problem for surfactant concentration over moderate
length scales (i.e. the streamwise period of the SHS) and an unsteady
advection–diffusion problem for surfactant flux over long length scales (i.e.
the streamwise length of the channel), whilst assuming that bulk diffusion is
strong enough for cross-channel concentration gradients to be small (Tomlinson
et al., 2023a). The surfactant concentration transport equations are nonlinear
and of mixed hyperbolic-parabolic type; the unsteady evolution of the
surfactant flux over the length of the channel is predominantly hyperbolic,
allowing the formation of shocks. The problem possesses a number of distinct
asymptotic regimes, which we exploit to reveal how the shocks forming in the
space- and time-dependent surfactant-flux distribution affect the slip length
and drag reduction. The slip length and drag reduction are key quantities of
interest for practical applications that can be shown to satisfy their own
unsteady partial differential equations over long length scales. We predict
the propagation speed of a disturbance to the surfactant flux and investigate
how excess surfactant can be advected out of the channel to maximise the time-
averaged drag reduction for laminar microchannel applications. Furthermore, by
investigating the unsteady transport of surfactant in laminar channel flows,
the theory developed here is a step towards describing the unsteady transport
of surfactant in unsteady turbulent boundary-layer flows over SHSs, for
applications in marine hydrodynamics, for instance.

The paper is arranged as follows. In §2, the problem is formulated and
homogenised to derive an unsteady advection–diffusion equation for surfactant
flux transport through the channel. At leading order, we derive a purely
advective transport equation for the surfactant flux, valid at the channel
scale. In §3, results are presented for the surfactant flux, drag reduction,
propagation speed, slip velocity and concentration field. We describe the
parameter space and identify regions of high and low drag reduction. We detail
results for two (bell-shaped) canonical distributions of surfactant flux. In
particular, one profile induces a transition between the high and low drag
reduction regions of the parameter space, giving rise to shock formation. We
study these cases using both theoretical and numerical methods, providing
closed-form asymptotic predictions of drag reduction. In §4, we summarise and
discuss our main results. We provide a table with closed-form asymptotic
predictions for the flux propagation speed in the different parts of the
parameter space. These predictions are expressed both using relevant non-
dimensional and dimensional parameters, and are intended as a useful guide for
applications.

##  2 Formulation

###  2.1 Governing equations

Figure 1: A schematic illustrating the multiple length scales present in a
channel flow bounded by streamwise and spanwise-periodic SHSs. Top figure: the
surfactant flux varies in the streamwise direction, $\hat{\chi}$, and time,
$\hat{\tau}$, over the length of the channel, $2\hat{L}_{x}$. Middle figure:
over each period, the concentration of surfactant varies quasi-steadily in the
streamwise direction, $\hat{x}$. Multiple periodic cells are shown, each with
a varying concentration of surfactant. The origin of the Cartesian coordinate
system, $\hat{\boldsymbol{x}}=\boldsymbol{0}$, is placed in the middle of the
interface of the central period, at a distance $\hat{L}_{x}$ along the
channel. Bottom figure: a periodic cell, identifying domains 1 and 2 and, at
the SHS, the ridge, solid wall and liquid–gas interface. For each periodic
cell, the half channel height is $\hat{H}$, the streamwise (spanwise) gas
fraction is $\phi_{z}$ ($\phi_{z}$) and the streamwise (spanwise) period
length is $2\hat{P}_{x}$ ($2\hat{P}_{z}$).

We consider a laminar pressure-driven fluid flow, contaminated with soluble
surfactant, in a channel bounded between two SHSs that are periodic in the
streamwise and spanwise directions, as illustrated in figure 1. We use hats to
indicate dimensional quantities. The streamwise, wall-normal and spanwise
directions are denoted by $\hat{x}$-, $\hat{y}$\- and $\hat{z}$-coordinates,
where $\hat{\boldsymbol{x}}=(\hat{x},\,\hat{y},\,\hat{z})$ is the space vector
and $\hat{t}$ is time. Assuming that the fluid is incompressible and
Newtonian, we define the velocity vector
$\hat{\boldsymbol{u}}=(\hat{u}(\hat{\boldsymbol{x}},\,\hat{t}),\,\hat{v}(\hat{\boldsymbol{x}},\,\hat{t}),\,\hat{w}(\hat{\boldsymbol{x}},\,\hat{t}))$,
pressure field $\hat{p}(\hat{\boldsymbol{x}},\,\hat{t})$, bulk surfactant
concentration field $\hat{c}(\hat{\boldsymbol{x}},\,\hat{t})$ and interfacial
surfactant concentration field $\hat{\Gamma}(\hat{x},\,\hat{z},\,\hat{t})$.
The streamwise length of the channel is $2\hat{L}_{x}$ and the periodic cell
has streamwise (spanwise) 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}$). The channel
height is $2\hat{H}$. The SHSs are made up of $2N+1$ periodic cells in the
streamwise direction. For $n\in\\{-N,\,...,\,N\\}$, the $n$th periodic cell is
split into two subdomains along the streamwise direction, similarly to
Temprano-Coleto et al. (2023),

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

At the SHS, $\hat{y}=0$ and $\hat{y}=2\hat{H}$, we define the $n$th interface,
ridge and solid region, as

$\displaystyle\hat{\mathcal{I}}^{n}$
$\displaystyle=\\{\hat{x}-2n\hat{P}_{x}\in[-\phi_{x}\hat{P}_{x},\,\phi_{x}\hat{P}_{x}]\\}\times\\{\hat{z}\in[-\phi_{z}\hat{P}_{z},\,\phi_{z}\hat{P}_{z}]\\},$
(2a) $\displaystyle\hat{\mathcal{R}}^{n}$
$\displaystyle=\\{\hat{x}-2n\hat{P}_{x}\in[-\phi_{x}\hat{P}_{x},\,\phi_{x}\hat{P}_{x}]\\}\times\\{\hat{z}\in[-\hat{P}_{z},\,-\phi_{z}\hat{P}_{z}]\cup[\phi_{z}\hat{P}_{z},\,\hat{P}_{z}]\\},$
(2b) $\displaystyle\hat{\mathcal{S}}^{n}$
$\displaystyle=\\{\hat{x}-2n\hat{P}_{x}\in[\phi_{x}\hat{P}_{x},\,(2-\phi_{x})\hat{P}_{x}]\\}\times\\{\hat{z}\in[-\hat{P}_{z},\,\hat{P}_{z}]\\}.$
(2c)

The steady equations that govern the fluid and surfactant in each periodic
cell are described in detail by Tomlinson et al. (2023a). Here, we highlight
differences due to the unsteady transport of surfactant, whilst allowing the
concentration of surfactant to vary over the long length scale and slow time
scale associated with the channel. The bulk surfactant is coupled to the
steady incompressible flow through an unsteady advection–diffusion equation.
In $\hat{\mathcal{D}}_{1}^{n}$ and $\hat{\mathcal{D}}_{2}^{n}$,

$\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}-\hat{c}_{\hat{t}}=0,$
(3a–c)

where $\hat{\mu}$ is dynamic viscosity and $\hat{D}$ is the surfactant bulk
diffusivity. The interfacial surfactant is coupled to the flow through an
unsteady advection–diffusion equation and a linear equation of state, such
that
$(\hat{\sigma}_{\hat{x}},\,\hat{\sigma}_{\hat{z}})=(-\hat{A}\hat{\Gamma}_{\hat{x}},\,-\hat{A}\hat{\Gamma}_{\hat{z}})$,
where $\hat{\sigma}$ is the surface tension and $\hat{A}$ is the surface
activity (Manikantan & Squires, 2020). At $\hat{y}=0$ ($\hat{y}=2\hat{H}$) and
along $\hat{\mathcal{I}}^{n}$, the boundary and bulk–interface coupling
conditions for surfactant and flow and the interfacial surfactant transport
equation are

$\hat{\mu}\boldsymbol{n}\cdot\hat{\boldsymbol{\nabla}}\hat{u}-\hat{A}\hat{\Gamma}_{\hat{x}}=0,\quad\hat{v}=0,\quad\hat{\mu}\boldsymbol{n}\cdot\hat{\boldsymbol{\nabla}}\hat{w}-\hat{A}\hat{\Gamma}_{\hat{z}}=0,\quad\hat{D}\boldsymbol{n}\cdot\hat{\boldsymbol{\nabla}}\hat{c}-\hat{K}_{a}\hat{c}+\hat{K}_{d}\hat{\Gamma}=0,\\\
\hat{D}_{I}(\hat{\Gamma}_{\hat{x}\hat{x}}+\hat{\Gamma}_{\hat{z}\hat{z}})+\hat{K}_{a}\hat{c}-\hat{K}_{d}\hat{\Gamma}-(\hat{u}\hat{\Gamma})_{\hat{x}}-(\hat{w}\hat{\Gamma})_{\hat{z}}-\hat{\Gamma}_{\hat{t}}=0,$
(5)

where $\boldsymbol{n}$ is the unit normal to the interface (pointing into the
channel), $\hat{D}_{I}$ is the surfactant interfacial diffusivity,
$\hat{K}_{a}$ is the adsorption rate and $\hat{K}_{d}$ is the desorption rate.
At $\hat{y}=0$ ($\hat{y}=2\hat{H}$) on $\partial\hat{\mathcal{I}}^{n}$, there
is no flux of surfactant

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

At $\hat{y}=0$ ($\hat{y}=2\hat{H}$) along
$\hat{\mathcal{R}}^{n}\cup\hat{\mathcal{S}}^{n}$, the flow and surfactant
boundary conditions are

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

Throughout $\hat{\mathcal{D}}^{n}_{1}$ and $\hat{\mathcal{D}}^{n}_{2}$, we
assume that the flow and concentration fields are periodic in the spanwise
directions,

$\hat{\boldsymbol{q}}(\hat{x},\,\hat{y},\,-\hat{P}_{z})=\hat{\boldsymbol{q}}(\hat{x},\,\hat{y},\,\hat{P}_{z}).$
(8)

Across interfaces between $\hat{\mathcal{D}}_{1}^{n}$,
$\hat{\mathcal{D}}_{2}^{n}$ and $\hat{\mathcal{D}}_{1}^{n+1}$ the flow and
concentration fields are assumed to be continuous between subdomains,

$\displaystyle\hat{\boldsymbol{q}}(((2n+\phi_{x})\hat{P}_{x})^{-},\,\hat{y},\,\hat{z})$
$\displaystyle=\hat{\boldsymbol{q}}(((2n+\phi_{x})\hat{P}_{x})^{+},\,\hat{y},\,\hat{z}),$
(9a)
$\displaystyle\hat{\boldsymbol{q}}(((2(n+1)-\phi_{x})\hat{P}_{x})^{-},\,\hat{y},\,\hat{z})$
$\displaystyle=\hat{\boldsymbol{q}}(((2(n+1)-\phi_{x})\hat{P}_{x})^{+},\,\hat{y},\,\hat{z}),$
(9b)

where we have defined
$\hat{\boldsymbol{q}}=(\hat{u},\,\hat{v},\,\hat{w},\,\hat{c})$.

We can integrate (1)–(8) across the channel to derive equations relating the
bulk and interfacial flux of fluid and surfactant. The unsteady surfactant
transport equations model how the bulk and interfacial surfactant fluxes
change as surfactants adsorb and desorb at the liquid–gas interface and the
concentration field evolves in time,

$\displaystyle\int_{\hat{\mathscr{A}}_{n}}\hat{c}_{\hat{t}}\,\text{d}\hat{A}+\frac{\text{d}}{\text{d}\hat{x}}\int_{\hat{\mathscr{A}}_{n}}(\hat{u}\hat{c}-\hat{D}\hat{c}_{\hat{x}})\,\text{d}\hat{A}-2\int_{\hat{\mathscr{I}}_{n}}(\hat{K}_{d}\hat{\Gamma}-\hat{K}_{a}\hat{c})\,\text{d}\hat{z}=0\quad\text{in}\quad\mathcal{D}_{1}^{n},$
(10a)
$\displaystyle\int_{\hat{\mathscr{I}}_{n}}\hat{\Gamma}_{\hat{t}}\,\text{d}\hat{z}+\frac{\text{d}}{\text{d}\hat{x}}\int_{\hat{\mathscr{I}}_{n}}(\hat{u}\hat{\Gamma}-\hat{D}_{I}\hat{\Gamma}_{\hat{x}})\,\text{d}\hat{z}+\int_{\hat{\mathscr{I}}_{n}}(\hat{K}_{d}\hat{\Gamma}-\hat{K}_{a}\hat{c})\,\text{d}\hat{z}=0\quad\text{in}\quad\mathcal{D}_{1}^{n},$
(10b)
$\displaystyle\int_{\hat{\mathscr{A}}_{n}}\hat{c}_{\hat{t}}\,\text{d}\hat{A}+\frac{\text{d}}{\text{d}\hat{x}}\int_{\hat{\mathscr{A}}_{n}}(\hat{u}\hat{c}-\hat{D}\hat{c}_{\hat{x}})\,\text{d}\hat{A}=0\quad\text{in}\quad\mathcal{D}_{2}^{n}.$
(10c)

Here,
$\int_{\hat{\mathscr{A}}_{n}}\cdot\,\text{d}\hat{A}\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}$
and
$\int_{\hat{\mathscr{I}}_{n}}\cdot\,\text{d}\hat{z}\equiv\int_{\hat{z}=-\hat{P}_{z}}^{\hat{P}_{z}}\cdot\,\text{d}\hat{z}$
for
$\hat{x}-2n\hat{P}_{x}\in[-\phi_{x}\hat{P}_{x},\,(2-\phi_{x})\hat{P}_{x}]$.
For a steady flow driven in the streamwise direction, the cross-channel
integrated streamwise velocity field, referred to hereafter as the flux of
fluid, $\hat{Q}$, is uniform along the length of the channel,

$\hat{Q}=\int_{\hat{\mathscr{A}}_{n}}\hat{u}\,\text{d}\hat{A}.$ (11)

In contrast, the cross-channel integrated total flux of surfactant, referred
to hereafter as the flux of surfactant, $\hat{K}=\hat{K}(\hat{x},\,\hat{t})$,
can vary along the length of the channel due to unsteady effects, according to

$\int_{\hat{\mathscr{A}}_{n}}\hat{c}_{\hat{t}}\,\text{d}\hat{A}+2\int_{\hat{\mathscr{I}}_{n}}\hat{\Gamma}_{\hat{t}}\,\text{d}\hat{z}+\hat{K}_{\hat{x}}=0\quad\text{in}\quad\hat{\mathcal{D}}_{1}^{n},\quad\int_{\hat{\mathscr{A}}_{n}}\hat{c}_{\hat{t}}\,\text{d}\hat{A}+\hat{K}_{\hat{x}}=0\quad\text{in}\quad\hat{\mathcal{D}}_{2}^{n},$
(12a, b)

where we have reformulated (10) and defined

$\displaystyle\hat{K}$
$\displaystyle=\int_{\hat{\mathscr{A}}_{n}}(\hat{u}\hat{c}-\hat{D}\hat{c}_{\hat{x}})\,\text{d}\hat{A}+2\int_{\hat{\mathscr{I}}_{n}}(\hat{u}\hat{\Gamma}-\hat{D}_{I}\hat{\Gamma}_{\hat{x}})\,\text{d}\hat{z}\quad\text{in}\quad\hat{\mathcal{D}}_{1}^{n},$
(13a) $\displaystyle\hat{K}$
$\displaystyle=\int_{\hat{\mathscr{A}}_{n}}(\hat{u}\hat{c}-\hat{D}\hat{c}_{\hat{x}})\,\text{d}\hat{A}\quad\text{in}\quad\hat{\mathcal{D}}_{2}^{n}.$
(13b)

We also define $\hat{K}_{m}=\max(\hat{K}(\hat{x},\,0))$ to be the maximum
initial surfactant flux along the length of the channel.

Defining the cross-channel-averaged pressure drop per period
$\Delta_{n}\hat{p}(\hat{t})\equiv\langle\hat{p}\rangle((2n-\phi_{x})\hat{P}_{x})-\langle\hat{p}\rangle((2(n+1)-\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, we can define the normalised drag reduction over
the $n$th cell as

${DR}_{n}(\hat{t})=\frac{\Delta_{n}\hat{p}_{I}-\Delta_{n}\hat{p}}{\Delta_{n}\hat{p}_{I}-\Delta_{n}\hat{p}_{U}},$
(14)

where $\Delta_{n}\hat{p}=\Delta_{n}\hat{p}_{I}$ when the liquid–gas interface
is immobilised by surfactant and is no-slip (${DR}_{n}=0$) and
$\Delta_{n}\hat{p}=\Delta_{n}\hat{p}_{U}$ when the liquid–gas interface is
unaffected by surfactant and is shear-free (${DR}_{n}=1$).

###  2.2 Non-dimensionalisation and scalings

In table 1, we summarise the different length, time and velocity scales of
interest in the transport problem described in (1)–(13) and figure 1, assuming
that the channel has an order-one cross-channel aspect ratio
$\hat{H}\sim\hat{P}_{z}$, but small channel-height-to-streamwise-period ratio
$\epsilon=\hat{H}/\hat{P}_{x}\ll 1$ and small streamwise-period-to-channel-
length ratio $\mathcal{E}=\hat{P}_{x}/\hat{L}_{x}\ll 1$. Defining
$\epsilon\hat{U}=\hat{Q}/(\hat{H}^{2})$ as a velocity scale,
$\hat{P}=\hat{\mu}\hat{U}/\hat{H}$ as a pressure scale,
$\hat{C}=\hat{K}_{m}/\hat{Q}$ as a bulk concentration scale and
$\hat{G}=\hat{K}_{a}\hat{C}/\hat{K}_{d}$ as an interfacial concentration
scale, we non-dimensionalise (1)–(14) using multiple time and spatial scales:

$t=\frac{\hat{t}}{\epsilon\hat{P}_{x}/\hat{U}},\quad
T=\frac{\hat{t}}{\hat{P}_{x}/\epsilon\hat{U}},\quad\tau=\frac{\hat{t}}{\hat{P}_{x}/(\epsilon\mathcal{E}\hat{U})},\quad\boldsymbol{x}_{\perp}=\frac{\hat{\boldsymbol{x}}_{\perp}}{\epsilon\hat{P}_{x}},\quad
x=\frac{\hat{x}}{\hat{P}_{x}},\quad\chi=\frac{\hat{x}}{\hat{P}_{x}/\mathcal{E}},\\\
\boldsymbol{u}_{\perp}=\frac{\hat{\boldsymbol{u}}_{\perp}}{\hat{U}},\quad
u=\frac{\hat{u}}{\epsilon\hat{U}},\quad p=\frac{\hat{p}}{\hat{P}},\quad
c=\frac{\hat{c}}{\hat{C}},\quad\Gamma=\frac{\hat{\Gamma}}{\hat{G}},\quad
K=\frac{\hat{K}}{\hat{K}_{m}},$ (15)

where $\hat{\boldsymbol{x}}_{\perp}\equiv(\hat{y},\,\hat{z})$ and
$\hat{\boldsymbol{u}}_{\perp}\equiv(\hat{v},\,\hat{w})$. In this paper, we
focus on the distinguished limit in which $\mathcal{E}=\lambda\epsilon^{2}$ as
$\epsilon\rightarrow 0$, where $\lambda$ is an $O(1)$ constant (this scaling
clarifies the asymptotics and is reasonable from an applications point of
view). This non-dimensionalisation yields a long-wave theory with rapid cross-
channel transport of surfactant over each period, which will be homogenised to
describe the slow transport of surfactant over multiple periods. For
$n\in\\{-N,\,...,\,N\\}$, the $n$th periodic cell becomes, in dimensionless
form (using quantities without hats),

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

where $P_{z}=\hat{P}_{z}/\hat{H}$. At $y=0$ ($y=2$), the regions of the SHS
are given by

$\displaystyle\mathcal{I}^{n}$
$\displaystyle=\\{x-2n\in[-\phi_{x},\,\phi_{x}]\\}\times\\{z\in[-\phi_{z}P_{z},\,\phi_{z}P_{z}]\\},$
(17a) $\displaystyle\mathcal{R}^{n}$
$\displaystyle=\\{x-2n\in[-\phi_{x},\,\phi_{x}]\\}\times\\{z\in[-P_{z},\,-\phi_{z}P_{z}]\cup[\phi_{z}P_{z},\,P_{z}]\\},$
(17b) $\displaystyle\mathcal{S}^{n}$
$\displaystyle=\\{x-2n\in[\phi_{x},\,2-\phi_{x}]\\}\times\\{z\in[-P_{z},\,P_{z}]\\}.$
(17c)

The length of the channel becomes
$2L_{x}=2\hat{L}_{x}/\hat{P}_{x}=2/\mathcal{E}=2/(\lambda\epsilon^{2})$.

Length scale | Time scale | Velocity scale | Interpretation
---|---|---|---
$\epsilon\hat{P}_{x}$ | $\displaystyle\frac{\epsilon\hat{P}_{x}}{\hat{U}}$ | $\hat{U}$ | Spanwise flow across a unit cell
$\hat{P}_{x}$ | $\displaystyle\frac{\hat{P}_{x}}{\epsilon\hat{U}}$ | $\epsilon\hat{U}$ | Streamwise flow along a unit cell
$\displaystyle\frac{\hat{P}_{x}}{\mathcal{E}}$ | $\displaystyle\frac{\hat{P}_{x}}{\epsilon\mathcal{E}\hat{U}}$ | $\epsilon\hat{U}$ | Streamwise flow over multiple cells

Table 1: Summary of the length, time and velocity scales of interest involved
in the problem described in (1)–(13) and figure 1, where $\epsilon\hat{U}$ is
the moderate streamwise velocity scale at the scale of individual cells (owing
to continuity of flux), $\epsilon=\hat{H}/\hat{P}_{x}$ is the ratio of channel
height to streamwise period length and $\mathcal{E}=\hat{P}_{x}/\hat{L}_{x}$
is the ratio of streamwise period length to channel length.

We then assume that the flow and surfactant variables are functions of the
short length scale and rapid time scale ($\boldsymbol{x}_{\perp}$ and $t$
respectively), moderate length scale and intermediate time scale ($x$ and $T$
respectively) and long length scale and slow time scale ($\chi$ and $\tau$
respectively), where these six variables are treated as independent of each
other. In $\mathcal{D}_{1}^{n}$ and $\mathcal{D}_{2}^{n}$, the incompressible
Stokes and surfactant transport equations in (1) become

$\displaystyle\epsilon^{2}(u_{x}+\lambda\epsilon^{2}u_{\chi})+\boldsymbol{\nabla}_{\perp}\cdot\boldsymbol{u}_{\perp}=0,$
(18a)
$\displaystyle\epsilon^{2}(u_{xx}+2\lambda\epsilon^{2}u_{x\chi}+\lambda^{2}\epsilon^{4}u_{\chi\chi})+\nabla^{2}_{\perp}u-p_{x}-\lambda\epsilon^{2}p_{\chi}=0,$
(18b) $\displaystyle\epsilon^{2}(\boldsymbol{u}_{\perp
xx}+2\lambda\epsilon^{2}\boldsymbol{u}_{\perp
x\chi}+\lambda^{2}\epsilon^{4}\boldsymbol{u}_{\perp\chi\chi})+\nabla^{2}_{\perp}\boldsymbol{u}_{\perp}-\boldsymbol{\nabla}_{\perp}p=\boldsymbol{0},$
(18c)
$\displaystyle(\epsilon^{2}(c_{xx}+2\lambda\epsilon^{2}c_{x\chi}+\lambda^{2}\epsilon^{4}c_{\chi\chi})+\nabla^{2}_{\perp}c)/\Pen-\epsilon^{2}u(c_{x}+\lambda\epsilon^{2}c_{\chi})-\boldsymbol{u}_{\perp}\cdot\boldsymbol{\nabla}_{\perp}c\hskip
28.45274pt$ $\displaystyle-
c_{t}-\epsilon^{2}(c_{T}+\lambda\epsilon^{2}c_{\tau})=0,$ (18d)

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}$. At $y=0$ ($y=2$) and
along $\mathcal{I}^{n}$, the boundary conditions for flow and surfactant, the
coupling conditions, and the interfacial surfactant transport equations in
(3a–c) give

$\displaystyle\boldsymbol{n}\cdot\boldsymbol{\nabla}u-\Ma(\Gamma_{x}+\lambda\epsilon^{2}\Gamma_{\chi})=0,$
(19a) $\displaystyle v=0,$ (19b)
$\displaystyle\boldsymbol{n}\cdot\boldsymbol{\nabla}w-\Ma\Gamma_{z}=0,$ (19c)
$\displaystyle\boldsymbol{n}\cdot\boldsymbol{\nabla}c-Da(c-\Gamma)=0,$ (19d)
$\displaystyle(\epsilon^{2}(\Gamma_{xx}+2\lambda\epsilon^{2}\Gamma_{x\chi}+\lambda^{2}\epsilon^{4}\Gamma_{\chi\chi})+\Gamma_{zz})/\Pen_{I}-\epsilon^{2}(u\Gamma)_{x}-\lambda\epsilon^{4}(u\Gamma)_{\chi}\hskip
56.9055pt$
$\displaystyle-(w\Gamma)_{z}-\Gamma_{t}-\epsilon^{2}(\Gamma_{T}+\lambda\epsilon^{2}\Gamma_{\tau})-\Bi(c-\Gamma)=0,$
(19e)

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. At $y=0$ ($y=2$) on
$\partial\mathcal{I}^{n}$, the no-flux interfacial surfactant boundary
conditions in (5) become

$u\Gamma-(\Gamma_{x}+\lambda\epsilon^{2}\Gamma_{\chi})/\Pen_{I}=0\ \
\text{at}\ \ x=\pm\phi_{x},\ \ w\Gamma-\Gamma_{z}/\Pen_{I}=0\ \ \text{at}\ \
z=\pm\phi_{z}P_{z}.$ (20a, b)

At $y=0$ ($y=2$) along $\mathcal{R}^{n}\cup\mathcal{S}^{n}$, the no-flux bulk
flow and bulk surfactant boundary conditions in (6a, b) give

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

Defining $\boldsymbol{q}=(u,\,v,\,w,\,c)$, across $\mathcal{D}_{1}^{n}$ and
$\mathcal{D}_{2}^{n}$, the spanwise (8) and streamwise (9a) continuity
conditions for the flow and surfactant become

$\displaystyle\boldsymbol{q}(x,\,y,\,-P_{z},\,t,\,T,\,\chi,\,\tau)$
$\displaystyle=\boldsymbol{q}(x,\,y,\,P_{z},\,t,\,T,\,\chi,\,\tau),$ (22a)
$\displaystyle\boldsymbol{q}((2n+\phi_{x})^{-},\,y,\,z,\,t,\,T,\,\chi,\,\tau)$
$\displaystyle=\boldsymbol{q}((2n+\phi_{x})^{+},\,y,\,z,\,t,\,T,\,\chi,\,\tau).$
(22b)

The streamwise flow and surfactant continuity condition for $\boldsymbol{q}$
between one cell and the next, i.e. between $\mathcal{D}_{2}^{n}$ and
$\mathcal{D}_{1}^{n+1}$, in (9b) is replaced by a stronger assumption to allow
the use of homogenisation theory (Bottaro, 2019), namely that $\boldsymbol{q}$
is a periodic function of the moderate length scale $x$, such that

$\boldsymbol{q}(2n-\phi_{x},\,y,\,z,\,t,\,T,\,\chi,\,\tau)=\boldsymbol{q}(2(n+1)-\phi_{x},\,y,\,z,\,t,\,T,\,\chi,\,\tau).$
(23)

Slow variations of flow properties from cell to cell will be accommodated via
dependence of the flow and surfactant variables on the long length scale
$\chi$ and slow time scale $\tau$.

The bulk and interfacial surfactant fluxes (10) satisfy

$\displaystyle\int_{\mathscr{A}_{n}}(c_{t}+\epsilon^{2}(c_{T}+\lambda\epsilon^{2}c_{\tau}))\,\text{d}A+\epsilon^{2}\frac{\text{d}}{\text{d}x}\int_{\mathscr{A}_{n}}\left(uc-\frac{c_{x}}{\Pen}-\lambda\epsilon^{2}\frac{c_{\chi}}{\Pen}\right)\,\text{d}A\hskip
39.83368pt$
$\displaystyle+\lambda\epsilon^{4}\frac{\text{d}}{\text{d}\chi}\int_{\mathscr{A}_{n}}\left(uc-\frac{c_{x}}{\Pen}-\lambda\epsilon^{2}\frac{c_{\chi}}{\Pen}\right)\,\text{d}A-\frac{2\Da}{\Pen}\int_{\mathscr{I}_{n}}(\Gamma-c)\,\text{d}z=0\
\ \text{in}\ \ \mathcal{D}_{1},$ (24a)
$\displaystyle\int_{\mathscr{I}_{n}}(\Gamma_{t}+\epsilon^{2}(\Gamma_{T}+\lambda\epsilon^{2}\Gamma_{\tau}))\,\text{d}z+\epsilon^{2}\frac{\text{d}}{\text{d}x}\int_{\mathscr{I}_{n}}\left(u\Gamma-\frac{\Gamma_{x}}{\Pen_{I}}-\lambda\epsilon^{2}\frac{\Gamma_{\chi}}{\Pen_{I}}\right)\,\text{d}z\hskip
39.83368pt$
$\displaystyle+\lambda\epsilon^{4}\frac{\text{d}}{\text{d}\chi}\int_{\mathscr{I}_{n}}\left(u\Gamma-\frac{\Gamma_{x}}{\Pen_{I}}-\lambda\epsilon^{2}\frac{\Gamma_{\chi}}{\Pen_{I}}\right)\,\text{d}z+\Bi\int_{\mathscr{I}_{n}}(\Gamma-c)\,\text{d}z=0\
\ \text{in}\ \ \mathcal{D}_{1},$ (24b)
$\displaystyle\int_{\mathscr{A}_{n}}(c_{t}+\epsilon^{2}(c_{T}+\lambda\epsilon^{2}c_{\tau}))\,\text{d}A+\epsilon^{2}\frac{\text{d}}{\text{d}x}\int_{\mathscr{A}_{n}}\left(uc-\frac{c_{x}}{\Pen}-\lambda\epsilon^{2}\frac{c_{\chi}}{\Pen}\right)\,\text{d}A\hskip
39.83368pt$
$\displaystyle+\lambda\epsilon^{4}\frac{\text{d}}{\text{d}\chi}\int_{\mathscr{A}_{n}}\left(uc-\frac{c_{x}}{\Pen}-\lambda\epsilon^{2}\frac{c_{\chi}}{\Pen}\right)\,\text{d}A=0\
\ \text{in}\ \ \mathcal{D}_{2},$ (24c)

where
$\int_{\mathscr{A}_{n}}\cdot\,\text{d}A\equiv\int_{z=-P_{z}}^{P_{z}}\int_{y=0}^{2H}\cdot\,\text{d}y\,\text{d}z$
and
$\int_{\mathscr{I}_{n}}\cdot\,\text{d}z\equiv\int_{z=-P_{z}}^{P_{z}}\cdot\,\text{d}z$
for $x-2n\in[\phi_{x},\,2-\phi_{x}]$. In $\mathcal{D}_{1}^{n}$ and
$\mathcal{D}_{2}^{n}$, the flux of fluid (11) is given by

$\int_{\mathscr{A}_{n}}u\,\text{d}A=1.$ (25)

The flux of surfactant, $K=K(x,\,t,\,T,\,\chi,\,\tau)$, is related to changes
in the bulk and surface concentration via (11), which becomes

$\displaystyle\int_{\mathscr{A}_{n}}(c_{t}+\epsilon^{2}(c_{T}+\lambda\epsilon^{2}c_{\tau}))\,\text{d}A+\frac{2\Da}{\Bi\Pen}\int_{\mathscr{I}_{n}}(\Gamma_{t}+\epsilon^{2}(\Gamma_{T}+\lambda\epsilon^{2}\Gamma_{\tau}))\,\text{d}z\hskip
54.06006pt$
$\displaystyle+\epsilon^{2}(K_{x}+\lambda\epsilon^{2}K_{\chi})=0\quad\text{in}\quad\mathcal{D}_{1}^{n},$
(26a)
$\displaystyle\int_{\mathscr{A}_{n}}(c_{t}+\epsilon^{2}(c_{T}+\lambda\epsilon^{2}c_{\tau}))\,\text{d}A+\epsilon^{2}(K_{x}+\lambda\epsilon^{2}K_{\chi})=0\quad\text{in}\quad\mathcal{D}_{2}^{n},$
(26b)

where the flux of surfactant (13) is given by

$\displaystyle K$
$\displaystyle=\int_{\mathscr{A}_{n}}\left(uc-\frac{c_{x}}{\Pen}-\lambda\epsilon^{2}\frac{c_{\chi}}{\Pen}\right)\,\text{d}A$
$\displaystyle\hskip
110.96556pt+\frac{2\Da}{\Bi\Pen}\int_{\mathscr{I}_{n}}\left(u\Gamma-\frac{\Gamma_{x}}{\Pen_{I}}-\lambda\epsilon^{2}\frac{\Gamma_{\chi}}{\Pen_{I}}\right)\,\text{d}z\quad\text{in}\quad\mathcal{D}_{1}^{n},$
(27a) $\displaystyle K$
$\displaystyle=\int_{\mathscr{A}_{n}}\left(uc-\frac{c_{x}}{\Pen}-\lambda\epsilon^{2}\frac{c_{\chi}}{\Pen}\right)\,\text{d}A\quad\text{in}\quad\mathcal{D}_{2}^{n},$
(27b)

and $\max(K(x,\,t,\,T,\,\chi,\,\tau))=1$ at $t=T=\tau=0$.

The normalised drag reduction (14) over the $n$th periodic cell becomes

${DR}_{n}(t,\,T,\,\chi,\,\tau)=\frac{\Delta^{n}p_{I}-\Delta^{n}p}{\Delta^{n}p_{I}-\Delta^{n}p_{U}},$
(28)

where $\Delta^{n}{p}\equiv\langle p\rangle(2n-\phi_{x})-\langle
p\rangle(2(n+1)-\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})$.

###  2.3 Asymptotic homogenisation

We assume that $\Pen\sim\Pen_{I}\sim\Ma\sim O(1)$ and $\Bi\sim\Da\sim
O(\epsilon^{2})$ in the limit $\epsilon\ll 1$, so that bulk–surface exchange
is comparable to advection, diffusion and Marangoni effects in
$\mathcal{D}_{1}^{n}$ and $\mathcal{D}_{2}^{n}$ for $n\in\\{-N,\,...,\,N\\}$.
As discussed in Tomlinson et al. (2023a), this scaling means that we arrive at
the most general form of the surfactant transport equations with moderate
exchange, whereas, if we had assumed that $\Bi\sim\Da\sim O(1)$, then we would
arrive at a sublimit with strong exchange. In the limit $\epsilon\rightarrow
0$, we rescale $\Bi=\epsilon^{2}\mathscr{B}$ and
$\Da=\epsilon^{2}\mathscr{D}$, where $\mathscr{B}$ and $\mathscr{D}$ are
positive $O(1)$ constants. We then substitute the asymptotic expansion

$\begin{pmatrix}u\\\ v\\\ w\\\ p\\\ c\\\ \Gamma\\\
K\end{pmatrix}=\begin{pmatrix}u_{0}\\\ v_{0}\\\ w_{0}\\\ p_{0}\\\ c_{0}\\\
\Gamma_{0}\\\ K_{0}\end{pmatrix}+\epsilon^{2}\begin{pmatrix}u_{1}\\\ v_{1}\\\
w_{1}\\\ p_{1}\\\ c_{1}\\\ \Gamma_{1}\\\
K_{1}\end{pmatrix}+\epsilon^{4}\begin{pmatrix}u_{2}\\\ v_{2}\\\ w_{2}\\\
p_{2}\\\ c_{2}\\\ \Gamma_{2}\\\ K_{2}\end{pmatrix}+...,$ (29)

into (18)–(27). The leading-order, first-order and second-order problems are
addressed in §2.3.1–2.3.3 respectively.

####  2.3.1 Leading-order problem

In the leading-order problem, we simplify the dependence of the velocity,
pressure, and concentration field on the space and time variables in (2.2). In
$\mathcal{D}_{1}^{n}$ and $\mathcal{D}_{2}^{n}$, streamwise gradients of the
velocity and bulk concentration are small compared to cross-channel gradients.
Hence, cross-channel diffusion balances advection and unsteady effects in the
bulk equation, through the two-dimensional problem

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

At $y=0$ ($y=2$) and along $\mathcal{I}^{n}$, streamwise gradients of the
streamwise velocity and surface concentration are small compared to spanwise
gradients. Hence, spanwise diffusion balances advection and unsteady effects
in the interfacial equation, via

$\boldsymbol{n}\cdot\boldsymbol{\nabla}u_{0}-\Ma\Gamma_{0x}=0,\quad
v_{0}=0,\quad\boldsymbol{n}\cdot\boldsymbol{\nabla}w_{0}-\Ma\Gamma_{0z}=0,\\\
\boldsymbol{n}\cdot\boldsymbol{\nabla}c_{0}=0,\quad\Gamma_{0zz}/\Pen_{I}-(w_{0}\Gamma_{0})_{z}-\Gamma_{0t}=0.$
(32)

At $y=0$ ($y=2$) and on $\partial\mathcal{I}^{n}$,

$u_{0}\Gamma_{0}-\Gamma_{0x}/\Pen_{I}=0\quad\text{at}\quad x=\pm\phi_{x},\quad
w_{0}\Gamma_{0}-\Gamma_{0z}/\Pen_{I}=0\quad\text{at}\quad z=\pm\phi_{z}P_{z}.$
(33a, b)

At $y=0$ ($y=2$) and along $\mathcal{R}^{n}\cup\mathcal{S}^{n}$,

$u_{0}=0,\quad v_{0}=0,\quad w_{0}=0,\quad\ c_{0y}=0.$ (34a–d)

As there are no streamwise gradients of $\boldsymbol{u}_{0}$ in
(2.3.1)–(30a–c), the two-dimensional problem does not capture inner regions
near $x=2n\pm\phi_{x}$, governed by the three-dimensional Stokes equations to
ensure continuity of $\boldsymbol{u}_{0}$ across domains $\mathcal{D}_{1}^{n}$
and $\mathcal{D}_{2}^{n}$.

In $\mathcal{D}_{1}^{n}$ and $\mathcal{D}_{2}^{n}$, the surfactant field
evolves faster in time $t$ than any changes to the flux of surfactant and
bulk–surface exchange at leading-order, so that (24)–(25) give

$\int_{\mathscr{A}_{n}}c_{0t}\,\text{d}A=0,\quad\int_{\mathscr{I}_{n}}\Gamma_{0t}\,\text{d}z=0,\quad\int_{\mathscr{A}_{n}}u_{0}\,\text{d}A=1.$
(35a–c)

According to (27), the flux of surfactant,
$K_{0}=K_{0}(x,\,t,\,T,\,\chi,\,\tau)$, is given by,

$\displaystyle K_{0}$
$\displaystyle=\int_{\mathscr{A}_{n}}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}A+\frac{2\mathscr{D}}{\mathscr{B}\Pen}\int_{\mathscr{I}_{n}}\left(u_{0}\Gamma_{0}-\frac{\Gamma_{0x}}{\Pen_{I}}\right)\,\text{d}z\quad\text{in}\quad\mathcal{D}_{1}^{n},$
(36a) $\displaystyle K_{0}$
$\displaystyle=\int_{\mathscr{A}_{n}}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}A\quad\text{in}\quad\mathcal{D}_{2}^{n}.$
(36b)

The leading-order solution can be expected to decay exponentially fast in time
$t$ to $\Gamma_{0}=\Gamma_{0}(x,\,T,\,\chi,\,\tau)$,
$c_{0}=c_{0}(x,\,T,\,\chi,\,\tau)$, $p_{0}=p_{0}(x,\,T,\,\chi,\,\tau)$,
$v_{0}=w_{0}=0$ and $K_{0}=K_{0}(x,\,T,\,\chi,\,\tau)$ (Mcnair et al., 2022),
satisfying (2.3.1). That is, at moderate ($T$) and slow ($\tau$) time scales,
the concentration field does not vary in the spanwise direction, and there are
no concentration gradients to generate velocities in the cross-plane. Using
linear superposition, we can decompose $u_{0}$ into a contribution from the
streamwise pressure gradient $p_{0x}$ which drives the flow and the streamwise
interfacial concentration gradient $\Gamma_{0x}$ which inhibits it owing to
adverse Marangoni forces, via

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

where the steady velocity profiles $\tilde{U}(y,\,z)$, $\bar{U}(y,\,z)$ and
$\breve{U}(y,\,z)$ are described in Tomlinson et al. (2023a) and satisfy a
Poisson-type problem in the cross-section. Substituting (2.3.1) into (2.3.1c),
we obtain relations between the volume flux, pressure gradient and interfacial
surfactant gradient in $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$,

$\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}^{n},\quad\breve{Q}p_{0x}=1\quad\text{in}\quad\mathcal{D}_{2}^{n},$
(38a–c)

where the fluxes $\tilde{Q}$, $\bar{Q}$, $\breve{Q}$, $\tilde{q}$, $\bar{q}$
and $q$ are described by Tomlinson et al. (2023a).

####  2.3.2 First-order problem

In the first-order problem, we relate the surfactant concentration over
individual plastrons to quantities that vary over the long length scale and
slow time scale associated with the channel. Solvability conditions are
imposed on the first-order problem to constrain $u_{0}$, $c_{0}$ and
$\Gamma_{0}$. These conditions are provided by the conservation arguments that
result in the surfactant transport equations at $O(\epsilon^{2})$. Hence, (24)
gives

$\displaystyle\int_{\mathscr{A}_{n}}c_{1t}\,\text{d}A$
$\displaystyle=-\int_{\mathscr{A}_{n}}c_{0T}\,\text{d}A-\frac{\text{d}}{\text{d}x}\int_{\mathscr{A}_{n}}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}A$
$\displaystyle\hskip
153.6447pt+\frac{2\mathscr{D}}{\Pen}\int_{\mathscr{I}_{n}}(\Gamma_{0}-c_{0})\,\text{d}z\quad\text{in}\quad\mathcal{D}_{1},$
(39a) $\displaystyle\int_{\mathscr{I}_{n}}\Gamma_{1t}\,\text{d}z$
$\displaystyle=-\int_{\mathscr{I}_{n}}\Gamma_{0T}\,\text{d}z-\frac{\text{d}}{\text{d}x}\int_{\mathscr{I}_{n}}\left(u_{0}\Gamma_{0}-\frac{\Gamma_{0x}}{\Pen_{I}}\right)\,\text{d}z\hskip
85.35826pt$ $\displaystyle\hskip
159.3356pt-\mathscr{B}\int_{\mathscr{I}_{n}}(\Gamma_{0}-c_{0})\,\text{d}z\quad\text{in}\quad\mathcal{D}_{1},$
(39b) $\displaystyle\int_{\mathscr{A}_{n}}c_{1t}\,\text{d}A$
$\displaystyle=-\int_{\mathscr{A}_{n}}c_{0T}\,\text{d}A-\frac{\text{d}}{\text{d}x}\int_{\mathscr{A}_{n}}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}A\quad\text{in}\quad\mathcal{D}_{2},$
(39c)

and (27) becomes

$\displaystyle K_{1}$
$\displaystyle=\int_{\mathscr{A}_{n}}\left(u_{0}c_{1}+u_{1}c_{0}-\frac{c_{1x}}{\Pen}-\lambda\frac{c_{0\chi}}{\Pen}\right)\,\text{d}A\hskip
142.26378pt$ $\displaystyle\hskip
76.82234pt+\frac{2\mathscr{D}}{\mathscr{B}\Pen}\int_{\mathscr{I}_{n}}\left(u_{0}\Gamma_{1}+u_{1}\Gamma_{0}-\frac{\Gamma_{1x}}{\Pen_{I}}-\lambda\frac{\Gamma_{0\chi}}{\Pen_{I}}\right)\,\text{d}z\quad\text{in}\quad\mathcal{D}^{n}_{1},$
(40a) $\displaystyle K_{1}$
$\displaystyle=\int_{\mathscr{A}_{n}}\left(u_{0}c_{1}+u_{1}c_{0}-\frac{c_{1x}}{\Pen}-\lambda\frac{c_{0\chi}}{\Pen}\right)\,\text{d}A\quad\text{in}\quad\mathcal{D}^{n}_{2}.$
(40b)

To avoid secular growth of the net mass of surfactant in (39), we require that
the right-hand sides of (39) are zero. Hence, the bulk and interfacial
concentrations evolve over moderate time scales according to

$\displaystyle\int_{\mathscr{A}_{n}}c_{0T}\,\text{d}A+\frac{\text{d}}{\text{d}x}\int_{\mathscr{A}_{n}}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}A-\frac{2\mathscr{D}}{\Pen}\int_{\mathscr{I}_{n}}(\Gamma_{0}-c_{0})\,\text{d}z=0\quad\text{in}\quad\mathcal{D}_{1}^{n},$
(41a)
$\displaystyle\int_{\mathscr{A}_{n}}\Gamma_{0T}\,\text{d}A+\frac{\text{d}}{\text{d}x}\int_{\mathscr{I}_{n}}\left(u_{0}\Gamma_{0}-\frac{\Gamma_{0x}}{\Pen_{I}}\right)\,\text{d}z+\mathscr{B}\int_{\mathscr{I}}(\Gamma_{0}-c_{0})\,\text{d}z=0\quad\text{in}\quad\mathcal{D}_{1}^{n},$
(41b)
$\displaystyle\int_{\mathscr{A}_{n}}c_{0T}\,\text{d}A+\frac{\text{d}}{\text{d}x}\int_{\mathscr{A}_{n}}\left(u_{0}c_{0}-\frac{c_{0x}}{\Pen}\right)\,\text{d}A=0\quad\text{in}\quad\mathcal{D}_{2}^{n}.$
(41c)

Substituting the velocity and flux conditions (2.3.1)–(37a, b) into the
surfactant transport equations (41) gives us ODEs that govern the unsteady
advection, diffusion and exchange of surfactant over one period,

$\displaystyle\theta c_{0T}+c_{0x}-\alpha c_{0xx}-\nu(\Gamma_{0}-c_{0})$
$\displaystyle=0\quad\text{in}\quad\mathcal{D}_{1}^{n},$ (42a)
$\displaystyle\sigma\Gamma_{0T}+\beta\Gamma_{0x}-\gamma(\Gamma_{0}\Gamma_{0x})_{x}-\delta\Gamma_{0xx}-\nu(c_{0}-\Gamma_{0})$
$\displaystyle=0\quad\text{in}\quad\mathcal{D}_{1}^{n},$ (42b)
$\displaystyle\theta c_{0T}+c_{0x}-\alpha c_{0xx}$
$\displaystyle=0\quad\text{in}\quad\mathcal{D}_{2}^{n}.$ (42c)

We describe (42) as the unsteady moderate-exchange equations. The steady-state
problem was solved in Tomlinson et al. (2023a), where the transport
coefficients $\alpha$, $\beta$, $\gamma$, $\delta$ and $\nu$ (specified below)
were defined in terms of physical and geometrical parameters and the fluxes
$\tilde{Q}$, $\bar{Q}$, $\tilde{q}$ and $\bar{q}$. Combining (42) with (36)
gives a set of constraints on the total surfactant flux over one period,

$\displaystyle\theta c_{0T}+\sigma\Gamma_{0T}+K_{0x}$ $\displaystyle=0,\
\text{where}\ K_{0}=c_{0}-\alpha
c_{0x}+\beta\Gamma_{0}-\gamma\Gamma_{0}\Gamma_{0x}-\delta\Gamma_{0x}\
\text{in}\ \mathcal{D}_{1}^{n},$ (43a) $\displaystyle\theta c_{0T}+K_{0x}$
$\displaystyle=0,\ \text{where}\ K_{0}=c_{0}-\alpha c_{0x}\ \text{in}\
\mathcal{D}_{2}^{n}.$ (43b)

We solve (42)–(43) subject to boundary conditions which enforce continuity and
periodicity, of both the surfactant concentration and flux, between subdomains

$\displaystyle c_{0}(2n+\phi_{x}^{-},\,T,\,\chi,\,\tau)$
$\displaystyle=c_{0}(2n+\phi_{x}^{+},\,T,\,\chi,\,\tau),$ (44a) $\displaystyle
c_{0}(2n-\phi_{x},\,T,\,\chi,\,\tau)$
$\displaystyle=c_{0}(2(n+1)-\phi_{x},\,T,\,\chi,\,\tau),$ (44b) $\displaystyle
K_{0}(2n+\phi_{x}^{-},\,T,\,\chi,\,\tau)$
$\displaystyle=K_{0}(2n+\phi_{x}^{+},\,T,\,\chi,\,\tau),$ (44c) $\displaystyle
K_{0}(2n-\phi_{x},\,T,\,\chi,\,\tau)$
$\displaystyle=K_{0}(2(n+1)-\phi_{x},\,T,\,\chi,\,\tau),$ (44d)
$\displaystyle[\beta\Gamma_{0}-\gamma\Gamma_{0}\Gamma_{0x}-\delta\Gamma_{0x}](2n\pm\phi_{x},\,T,\,\chi,\,\tau)$
$\displaystyle=0.$ (44e)

In (42)–(44), we have introduced the following transport coefficients:

$\displaystyle\alpha$ $\displaystyle=\frac{4P_{z}}{\Pen}\quad\quad\quad\quad$
$\displaystyle\text{(bulk diffusion)};$ (45a) $\displaystyle\beta$
$\displaystyle=\frac{2\mathscr{D}\tilde{q}}{\mathscr{B}\Pen\tilde{Q}}\quad\quad\quad\quad$
$\displaystyle\text{(partition coefficient)};$ (45b) $\displaystyle\gamma$
$\displaystyle=\frac{2\Ma\mathscr{D}(\tilde{q}\bar{Q}/\tilde{Q}-\bar{q})}{\mathscr{B}\Pen}\quad\quad\quad\quad$
$\displaystyle\text{(surfactant strength)};$ (45c) $\displaystyle\delta$
$\displaystyle=\frac{4\phi_{z}P_{z}\mathscr{D}}{\mathscr{B}\Pen\Pen_{I}}\quad\quad\quad\quad$
$\displaystyle\text{(surface diffusion)};$ (45d) $\displaystyle\nu$
$\displaystyle=\frac{4\phi_{z}P_{z}\mathscr{D}}{\Pen}\quad\quad\quad\quad$
$\displaystyle\text{(exchange strength)};$ (45e) $\displaystyle\theta$
$\displaystyle=4P_{z}\quad\quad\quad\quad$ $\displaystyle\text{(bulk
capacitance)};$ (45f) $\displaystyle\sigma$
$\displaystyle=\frac{4\phi_{z}P_{z}\mathscr{D}}{\mathscr{B}\Pen}\quad\quad\quad\quad$
$\displaystyle\text{(surface capacitance)}.$ (45g)

The bulk (surface) diffusion coefficient $\alpha>0$ ($\delta>0$) compares the
strength of bulk (interfacial) streamwise diffusion to advection. The
partition coefficient $\beta>0$ characterises the distribution of the
surfactant flux, where for $\beta\gg 1$ ($\beta\ll 1$) the interfacial (bulk)
surfactant flux dominates. The surfactant strength $\gamma>0$ characterises
the impact of Marangoni stresses on the interfacial surfactant flux. The
exchange strength $\nu>0$ compares the rate of adsorption to advection. The
remaining parameters, $\theta$ and $\sigma$, are associated with time-
dependent variations and were not reported in Tomlinson et al. (2023a). The
bulk capacitance coefficient $\theta>0$ characterises the transverse aspect
ratio of the channel and specifies the bulk response to time-dependent changes
in the surfactant flux. The surface capacitance $\sigma>0$ is the rescaled (by
$4\phi_{z}P_{z}$) surfactant depletion depth
$L_{d}=\mathscr{D}/(\mathscr{B}\Pen)$; $\sigma$ captures the manner in which
solubility regulates the surface response to gradients in the surfactant flux.
The dependence of the transport coefficients on dimensional parameters will be
discussed later in §4.

We solve (42)–(44) subject to the initial condition
$c_{0}(x,\,0,\,\chi,\,0)=1$ to illustrate the dependence of the bulk
concentration field on the intermediate time scale $T$ in figure 2(a);
convergence to a steady state for different values of $\theta=\sigma$ is
illustrated using $c_{0}(\phi_{x},\,T,\,\chi,\,\tau)$ in figure 2(b). The
initially uniform concentration falls to an equilibrium state, periodic over
the unit cell, in which the positive gradient ($c_{0x}>0$) in
$\mathcal{D}_{1}$ generates an interfacial stress opposing the mean flow. The
time taken to reach a steady state increases with $\theta$ and $\sigma$.
However, as our primary objective is to investigate surfactant transport over
the full length of the channel, we assume that the leading-order solution is
close to equilibrium and the concentration field no longer depends on the
intermediate time $T$, i.e. $\Gamma_{0}=\Gamma_{0}(x,\,\chi,\,\tau)$ and
$c_{0}=c_{0}(x,\,\chi,\,\tau)$. As the concentration field does not depend on
$T$, from (43), we have that $K_{0x}=0$ in $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$, and therefore the problem in each period simplifies to
finding $c_{0}=c_{0}(x;\,K_{0})$ and $\Gamma_{0}=\Gamma_{0}(x;\,K_{0})$ for a
given surfactant flux $K_{0}=K_{0}(\chi,\,\tau)$. Hence, we solve the steady
moderate-exchange equations from Tomlinson et al. (2023a), given by

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

subject to the steady surfactant flux conditions

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

and boundary conditions given in (44).

(a)(b)

Figure 2:  The unsteady bulk concentration field ($c_{0}$) in a unit cell
($\mathcal{D}_{1}\cup\mathcal{D}_{2}$), computed using (42)–(44) with
$c_{0}(x,\,0,\,\chi,\,0)=1$. When bulk–surface exchange is strong
($c_{0}\approx\Gamma_{0}$), the drag reduction depends on $\Delta
c_{0}=c_{0}(\phi_{x},\,T,\,\chi,\,\tau)-c_{0}(-\phi_{x},\,T,\,\chi,\,\tau)$
and a large gradient of $c_{0}$ in $\mathcal{D}_{1}$ can partially immobilise
the liquid–gas interface. Plot of (a) $c_{0}(x,\,T,\,\chi,\,\tau)$ at
different times $T$ for $\theta=\sigma=1$; and (b) leading-order bulk
concentration at the end of the plastron $c_{0}(\phi_{x},\,T,\,\chi,\,\tau)$
for varying bulk capacitance $\theta=\sigma$, computed using (46) where
$\alpha=\delta=10$ and $\beta=\gamma=1$.

The solution to the surfactant concentration transport equations (44, 46, 46)
exhibits multiple asymptotic regimes, which are discussed in detail in
Tomlinson et al. (2023a). Briefly, we distinguish a strong-exchange problem
($\nu\gg\max(1,\,\alpha,\,\delta)$), where the $c_{0}$ and $\Gamma_{0}$ fields
are in equilibrium ($c_{0}\approx\Gamma_{0}$), from a moderate-exchange
problem ($\nu=O(1,\,\alpha,\,\delta)$), where $c_{0}$ and $\Gamma_{0}$ are
distinct. In the strong-exchange problem, we identify three primary areas of
parameter space and two significant boundaries between them; these are
summarised in figure 3(a). In the Marangoni-dominated (M) region (analysed in
Appendix A.1), the interfacial surfactant gradient immobilises the liquid–gas
interface (leading to low drag reduction); in the advection-dominated (A)
region (Appendix A.2), the interfacial surfactant is swept to the downstream
stagnation point of each plastron and the liquid–gas interface is mostly
shear-free (high drag reduction); and in the diffusion-dominated (D) region
(Appendix A.3), the surfactant gradient is attenuated by diffusion and the
liquid–gas interface is mostly shear-free (high drag reduction). Across the
advection–Marangoni (AM) (Appendix A.4) and the diffusion–Marangoni (DM)
boundaries (Appendix A.5), these effects compete to partially immobilise the
liquid–gas interface (moderate drag reduction). Each of these regions has an
analogue when exchange is weak ($\nu\ll\min(1,\,\alpha,\,\delta)$; see figure
3b). These sub-regions have the same leading-order physics as regions M, D and
A and are referred to as Marangoni–exchange (M${}_{\text{E}}$),
diffusion–exchange (D${}_{\text{E}}$) and advection–exchange
(A${}_{\text{E}}$) sub-regions.

(a)(b)

Figure 3:  Contours of the leading-order drag reduction (${DR}_{0}$) in the
(a) $(\alpha,\,\gamma)$-plane for $\nu=100$ (strong-exchange problem) and (b)
in the $(\nu,\,\gamma)$-plane for $\alpha=10$ (illustrating the transition
from weak to strong exchange as $\nu$ increases), computed using (46) where
$\alpha=\delta$, $\beta=1$, $K_{0}=1$, $\phi_{x}=0.5$ and $P_{z}=1$ in both
panels. When ${DR}_{0}=0$ the liquid-gas interface is no-slip and when
${DR}_{0}=1$ the liquid-gas interface is shear-free. The Marangoni (M),
advection (A) and diffusion (D) regions, the advection–Marangoni (AM) and
diffusion–Marangoni (DM) boundaries (taken to lie in between the ${DR}=0.1$
and ${DR}_{0}=0.9$ contours), the Marangoni–exchange (M${}_{\text{E}}$) and
diffusion–exchange (D${}_{\text{E}}$) sub-regions, and the
diffusion–Marangoni–exchange boundary (DM${}_{\text{E}}$), are separated by
black lines. The advection–diffusion and moderate-exchange boundaries are
illustrated using black lines at $\alpha=0.3$ and $\nu=100$ in (a) and (b),
respectively. The white lines illustrate the relationship between $K_{0}$ and
$\gamma$ that is discussed at the end of §2.3.2.

The link between surfactant flux $K_{0}$ and surfactant concentration is
evident from (46), by noting that $K_{0}$ can be scaled to unity under the
mapping $c_{0}\rightarrow K_{0}c^{*}_{0}$, $\Gamma_{0}\rightarrow
K_{0}\Gamma^{*}_{0}$ and $\gamma\rightarrow\gamma^{*}/K_{0}$. Equivalently, by
solving the surfactant transport equations (46, 46) with $K_{0}=1$, we can
capture variations in the surfactant flux parametrically through variations in
${\gamma}$. For instance, increasing $K_{0}$ from $1/2$ to $1$ for fixed
$\alpha=1$ and $\gamma=100$ is equivalent (for the rescaled concentrations
$c^{*}_{0}$ and $\Gamma^{*}_{0}$) to setting $K_{0}=1$ and increasing $\gamma$
from 50 to 100 (illustrated by the right white line in figure 3a), thus moving
away from the DM boundary and further into the M region. Similarly, increasing
$K_{0}$ from $0.01$ to $1$ for fixed $\alpha=0.1$ and $\gamma=10$ has the more
dramatic effect of moving from the A region (high drag reduction) to the M
region (low drag reduction), by varying $\gamma$ from $0.1$ to $10$ with
$K_{0}=1$ (illustrated by the left white line in figure 3a). While we could
reduce the number of parameters in the plastron-scale problem by using the
rescaled concentrations $c^{*}_{0}$ and $\Gamma^{*}_{0}$ and by subsuming the
parameter $K_{0}$ into the rescaled surfactant strength parameter
$\gamma^{*}=\gamma K_{0}$, we choose to retain $K_{0}$ explicitly and use the
concentrations $c_{0}$ and $\Gamma_{0}$, because it provides a crucial link
between the plastron-scale and large channel-scale problems. We will return to
these examples in §3.

####  2.3.3 Second-order problem

In the second-order problem, we allow quantities such as the slip length and
drag reduction to vary over the long length scale of the channel and slow time
scale. Solvability conditions are imposed on the second-order problem to
constrain $u_{1}$, $c_{1}$ and $\Gamma_{1}$ appearing in (39)–(40). These
conditions are provided by the conservation arguments that result in the
surfactant transport equations at $O(\epsilon^{4})$. In $\mathcal{D}_{1}^{n}$,
(24a, b) give bulk and interfacial equations, which can be combined into

$\int_{\mathscr{A}_{n}}(c_{2t}+c_{1T})\,\text{d}A+\frac{2\mathscr{D}}{\mathscr{B}\Pen}\int_{\mathscr{I}_{n}}(\Gamma_{2t}+\Gamma_{1T})\,\text{d}z+K_{1x}=-\lambda\bigg{(}\int_{\mathscr{A}_{n}}c_{0\tau}\,\text{d}A\\\
+\frac{2\mathscr{D}}{\mathscr{B}\Pen}\int_{\mathscr{I}_{n}}\Gamma_{0\tau}\,\text{d}z+K_{0\chi}-\lambda\epsilon^{2}\bigg{(}\int_{\mathscr{A}_{n}}\frac{c_{0\chi\chi}}{\Pen}\,\text{d}A+\frac{2\mathscr{D}}{\mathscr{B}\Pen}\int_{\mathscr{I}_{n}}\frac{\Gamma_{0\chi\chi}}{\Pen_{I}}\,\text{d}z\bigg{)}\bigg{)},$
(48)

and in $\mathcal{D}_{2}$, (24c) gives

$\int_{\mathscr{A}_{n}}(c_{2t}+c_{1T})\,\text{d}A+K_{1x}=-\lambda\bigg{(}\int_{\mathscr{A}_{n}}c_{0\tau}\,\text{d}A+K_{0\chi}-\lambda\epsilon^{2}\int_{\mathscr{A}_{n}}\frac{c_{0\chi\chi}}{\Pen}\,\text{d}A\bigg{)},$
(49)

using the definition of the leading- and first-order surfactant fluxes $K_{0}$
and $K_{1}$ given in (36) and (40), respectively. In (48)–(49), we have
retained $O(\epsilon^{2})$ diffusion terms with respect to $\chi$ in order to
regularise any shocks that may arise because of the nonlinear dependence of
$c_{0}$ and $\Gamma_{0}$ on $K_{0}$ (discussed further in §3 below).

To avoid secular growth of the net mass of surfactant, we require that the
combined right-hand sides (i.e. source/sink terms) of the surfactant transport
equations, (48)–(49), integrate to zero along one period. We know from §2.3.2
that $c_{0}=c_{0}(x;\,K_{0})$ and $\Gamma_{0}=\Gamma_{0}(x;\,K_{0})$ where
$K_{0}=K_{0}(\chi,\,\tau)$. As $c_{0}$ is assumed to be periodic in $x$, we
can use the cell with $n=0$ as representative of all others. Hence,
integrating the surfactant transport equations (48)–(49) over one period,
using the velocity fields and fluxes from (2.3.1)–(37a, b) and using the
definition of the transport coefficients in (45), we obtain

$\frac{\text{d}C_{0}}{\text{d}\tau}+\frac{\text{d}A_{0}}{\text{d}\chi}-\frac{\text{d}M_{0}}{\text{d}\chi}-\frac{\text{d}D_{0}}{\text{d}\chi}-\lambda\epsilon^{2}\frac{\text{d}^{2}D_{1}}{\text{d}\chi^{2}}=0,$
(50)

where

$\displaystyle C_{0}(K_{0})$
$\displaystyle=\theta\int_{x=-\phi_{x}}^{2-\phi_{x}}c_{0}\,\text{d}x+\sigma\int_{x=-\phi_{x}}^{\phi_{x}}\Gamma_{0}\,\text{d}x,\quad$
$\displaystyle\text{(total weighted concentration)},$ (51a) $\displaystyle
A_{0}(K_{0})$
$\displaystyle=\int_{x=-\phi_{x}}^{2-\phi_{x}}c_{0}\,\text{d}x+\beta\int_{x=\phi_{x}}^{\phi_{x}}\Gamma_{0}\,\text{d}x\quad$
$\displaystyle\text{(advective flux)},$ (51b) $\displaystyle M_{0}(K_{0})$
$\displaystyle=\gamma\bigg{[}\frac{\Gamma_{0}^{2}}{2}\bigg{]}_{x=-\phi_{x}}^{\phi_{x}}\quad$
$\displaystyle\text{(Marangoni flux)},$ (51c) $\displaystyle D_{0}(K_{0})$
$\displaystyle=\alpha\bigg{[}c_{0}\bigg{]}_{x=-\phi_{x}}^{2-\phi_{x}}+\delta\bigg{[}\Gamma_{0}\bigg{]}_{x=-\phi_{x}}^{\phi_{x}}\quad$
$\displaystyle\text{(primary diffusive flux)},$ (51d) $\displaystyle
D_{1}(K_{0})$
$\displaystyle=\alpha\int_{x=-\phi_{x}}^{2-\phi_{x}}c_{0}\,\text{d}x+\delta\int_{x=-\phi_{x}}^{\phi_{x}}\Gamma_{0}\,\text{d}x,\quad$
$\displaystyle\text{(secondary diffusive flux)}.$ (51e)

We can express (50) as a nonlinear advection–diffusion equation for the
leading-order surfactant flux:

$\frac{\partial C_{0}}{\partial K_{0}}\frac{\partial
K_{0}}{\partial\tau}+\left(\frac{\partial A_{0}}{\partial
K_{0}}-\frac{\partial M_{0}}{\partial K_{0}}-\frac{\partial D_{0}}{\partial
K_{0}}\right)\frac{\partial
K_{0}}{\partial\chi}-\lambda\epsilon^{2}\frac{\partial}{\partial\chi}\left(\frac{\partial
D_{1}}{\partial K_{0}}\frac{\partial K_{0}}{\partial\chi}\right)=0.$ (52)

Equation (52) describes the spatio–temporal evolution of a disturbance to the
flux of surfactant. It predicts how such disturbances are advected and spread
by the flow over the long length scale ($\chi$) and slow time scale ($\tau$)
that are characteristic of the channel flow. This equation is motivated by
environmental surfactant concentrations that can vary significantly in space
and time across the large length scales involved in applications (Frossard et
al., 2019). As the initial surfactant flux evolves with respect to $\chi$ and
$\tau$, the surfactant flux transport (52) is coupled to the surfactant
concentration transport (46, 46) over a given streamwise period. We illustrate
the relationship between the surfactant flux and the resulting surfactant bulk
concentration in figure 4, which shows how they both vary over $500$
plastrons, in a case which is representative of the numerical simulations we
perform in this study. Where $K_{0}$ is elevated, stronger adsorption can be
expected to lead to interfacial rigidification, reducing the proportion of net
flux $K_{0}$ carried by the interface and increasing the local drag. The
impact of these changes on the evolution of the flux field is described by
(52). While $K_{0}$ varies smoothly with $\chi$, the underlying concentration
field has a wavy multiscale structure (inset). We now aim to solve this
coupled transport problem to compute the time- and space-varying leading-order
drag reduction and slip associated for specific initial and boundary
conditions of relevance to applications.

Figure 4: A plot illustrating the relationship between the total flux of
surfactant ($K_{0}=(1+\exp(-9\chi^{2}))/2$) over the length of the SHS
$(\chi\in[-1,\,1])$ and the corresponding concentration field ($c_{0}$ and
$\Gamma_{0}$) over $2N=500$ periods $(x\in[-100,\,100])$, computed using (65a)
where $\alpha=\delta=0.1$, $\beta=1$, $\gamma=10$, $\phi_{x}=0.5$ and
bulk–surface exchange is strong ($c_{0}=\Gamma_{0}$).

We can solve (52) using analytical methods and numerical techniques.
Analytically, we will neglect $O(\epsilon^{2})$ terms. This yields a
hyperbolic problem for which the solution is found using the method of
characteristics, generating simple formulae that can be used by
experimentalists and practitioners. Numerically, we will retain secondary
diffusion terms to regularise any shocks that may arise in hyperbolic problem
and to validate the analytical results.

We solve (52) subject to an initial condition, $K_{0}(\chi,\,0)$, satisfying
the constraint $\max(K_{0}(\chi,\,0))=1$. However, bar this constraint, we can
choose any initial condition for the distribution of surfactant flux in the
channel. Here, we choose a Gaussian distribution, rather than a step or a ramp
function, as it constitutes a classical example for which behaviours in purely
advective transport systems are observed, such as: wave steepening, wave
expansion or shock formation (Strauss, 2007). We solve (52) subject to the
initial and boundary conditions

$K_{0}(\chi,\,0)=K_{b}+(1-K_{b})\exp(-(10\chi+15/2)^{2})),\quad K_{0}(\pm
L_{x},\,\tau)=K_{b},$ (53a, b)

where $L_{x}$ is sufficiently large for (4b) to be approximately valid and
$K_{b}$ is the background surfactant flux. Taking $K_{b}=1/2$, the
distribution defined in (4) is characteristic of a channel that is
contaminated with a bolus of surfactant that locally doubles the background
surfactant flux. In this case, the drag-reduction values remain in region M
(see figure 3), as we will show in §3. We also take $K_{b}=0.01$, representing
an almost clean channel that is contaminated with a bolus of surfactant. We
investigate $K_{b}=0.01$, as an alternative to $K_{b}=1/2$, as in this case
values of the drag reduction can transition between region A (or D) and region
M (figure 3). As we show in §3, the surfactant-flux distribution can exhibit
shocks with such initial conditions.

(a) (b)

(c) (d)

Figure 5: Coefficients in the nonlinear advection–diffusion equation (52) for
the leading order surfactant flux $K_{0}$: (a) $\partial C_{0}/\partial K_{0}$
defined in (51a) for different $\theta=\sigma$, $\alpha=\delta=0.01$,
$\beta=10$ and $\gamma=5$; (b) $\partial A_{0}/\partial K_{0}$ defined in
(51b) for different $\beta$, $\alpha=\delta=0.01$, $\theta=\sigma=1$ and
$\gamma=5$; (c) $\partial M_{0}/\partial K_{0}$ defined in (51c) for different
$\gamma$, $\alpha=\delta=0.01$, $\theta=\sigma=1$ and $\beta=10$; (d)
$\partial D_{0}/\partial K_{0}$ defined in (51d) for different
$\alpha=\delta$, $\beta=10$, $\gamma=5$ and $\theta=\sigma=1$, where we vary
the surfactant flux ($K_{0}$) and $\phi_{x}=0.5$. Numerical solutions are
plotted using solid coloured lines and asymptotic solutions are plotted using
dashed black lines using (68), (74a, b) and (82a–c) in regions M, A and D
respectively.

The dependence of $C_{0}$, $A_{0}$, $M_{0}$ and $D_{0}$ in (51) on $K_{0}$,
for different values of $\alpha$, $\beta$, $\gamma$, $\delta$, $\theta$ and
$\sigma$, is illustrated in figure 5. We choose parameter values for $\alpha$,
$\beta$, $\gamma$ and $\delta$ such that drag reduction values are generally
at the AM boundary in the parameter space (see figure 3). By increasing
$K_{0}$, drag reduction values transition from regions A to M when
$K_{0}=O(2\phi_{x}\beta/\gamma)$ (see Appendix A.4). In figure 5(a, b), the
relationship between $C_{0}$ and $A_{0}$ with $K_{0}$ is linear in regions M
(see (66a–ca, b)) and A (see (74a, ba, b)), and nonlinear in between. In
figure 5(c, d), the relationship between $M_{0}$ and $D_{0}$ with $K_{0}$ can
be nonlinear in regions M and A, however, this nonlinearity does not affect
the leading-order surfactant-flux distribution and the drag reduction in
regions M and A, because $\gamma\gg\max(1,\,\alpha,\,\delta)$ and
$\max(\alpha,\,\delta,\,\gamma)\ll 1$, respectively (see Appendices A.1 and
A.2). In figure 5(d), when $\alpha=100$, we transition to region D because
$\min(\alpha,\,\delta)\gg\max(1,\,\gamma)$ (see Appendix A.3). In figure 5, we
see that all the coefficients $C_{0}$, $A_{0}$, $M_{0}$ and $D_{0}$ have a
nonlinear dependence on $K_{0}$ at the AM boundary, which will affect the
leading-order surfactant-flux distribution and drag reduction. We will discuss
this further in §3.

###  2.4 Solving the surfactant flux evolution equation

First, we neglect the $O(\epsilon^{2})$ terms in (52) and seek closed-form
expressions for the surfactant-flux distribution and drag reduction using the
method of characteristics. $\text{d}K_{0}/\text{d}\tau=0$ on the
characteristic curves of (52), which are solutions of the ODE (Strauss, 2007)

$\frac{\text{d}\chi}{\text{d}\tau}=a(K_{0}(\chi,\,\tau))=\frac{A^{\prime}_{0}-M^{\prime}_{0}-D^{\prime}_{0}}{C^{\prime}_{0}},$
(54)

where primes denote derivatives of the functions defined in (51) with respect
to $K_{0}$. The propagation speed, $a$, characterises how fast changes in the
surfactant-flux distribution will be transported in space and time along the
length of the channel. The characteristic curves $\chi=\chi(\tau)$ are
straight lines and $K_{0}$ is constant along each characteristic. We can solve
(52) subject to (4) provided that the characteristics do not intersect. The
characteristic through $(\chi,\,\tau)$ and $(\xi,\,0)$ has gradient
$(\chi-\xi)/\tau=\text{d}\chi/\text{d}\tau=a(K_{0}(\chi,\,\tau))=a(K_{0}(\xi,\,0))$,
and therefore,

$\chi=\xi+a(K_{0}(\xi,\,0))\tau\quad\text{for}\quad\xi\in\mathbb{R},$ (55)

which gives $\xi$ implicitly as a function of $\chi$ and $\tau$, i.e.
$\xi=\xi(\chi,\,\tau)$. Hence, the solution to (52) subject to (4) is given by
$K_{0}(\chi,\,\tau)=K_{0}(\xi,\,0)$. For $\xi_{1},\,\xi_{2}\in\mathbb{R}$, the
characteristic curves $\chi=\xi_{1}+a(K_{0}(\xi_{1},\,0))\tau$ and
$\chi=\xi_{2}+a(K_{0}(\xi_{2},\,0))\tau$ intersect when
$\tau=(\xi_{2}-\xi_{1})/(a(K_{0}(\xi_{1},\,0))-a(K_{0}(\xi_{2},\,0)))$. In the
limit $\xi_{2}\rightarrow\xi_{1}$, the shock time converges to
$\tau=-1/a_{\xi}(K_{0}(\xi_{1},\,0))$. Hence, the time $\tau_{b}$ when a shock
first forms is given by

$\tau_{b}=\min_{\xi\in\mathbb{R}}\left\\{\frac{-1}{a_{\xi}(K_{0}(\xi,\,0))}\right\\},$
(56)

and (55) gives the streamwise location $\chi_{b}$ where a shock first forms.
The integral form of (50) can be used to derive the Rankine–Hugoniot
condition, $u_{s}=[[A_{0}-M_{0}-D_{0}]]/[[C_{0}]]$, where $u_{s}$ is the shock
speed, with the jump bracket defined as
$[[q]]=q(\chi_{s}^{+},\,\tau)-q(\chi_{s}^{-},\,\tau)$ and $\chi_{s}$ is the
location of the shock (Strauss, 2007). The jump condition can then be
integrated to determine the location of the shock for times $\tau>\tau_{b}$,
for admissible shocks that satisfy the entropy condition,

$\chi_{s}(\tau)=u_{s}\tau+B,$ (57)

where the integration constant $B$ is determined using $\chi_{s}=\chi_{b}$ at
$\tau=\tau_{b}$ from (55)–(56).

Second, we retain $O(\epsilon^{2})$ diffusive terms in (52) and solve it
numerically subject to the initial and boundary conditions (4). We use the
method of lines and a backwards-in-time and centered-in-space scheme. As
discussed earlier, retaining small diffusive terms avoids numerical
difficulties associated with shock formation and regularises the shocks
through a small amount of diffusion. This procedure is outlined in detail in
Appendix B.

###  2.5 Quantities of interest for applications

As discussed in §1, the main quantities of interest in SHS applications,
ranging from microchannel to marine transport, are the effective slip length
and drag reduction. Following Tomlinson et al. (2023a), the leading-order drag
reduction (over a plastron) depends on the total flux of surfactant via

${DR}_{0}=1-\frac{\gamma\Delta\Gamma_{0}}{2\phi_{x}\beta},\quad\text{where}\quad\Delta\Gamma_{0}\equiv\Gamma_{0}(\phi_{x};\,K_{0})-\Gamma_{0}(-\phi_{x};\,K_{0}).$
(58)

The drag reduction inherits a dependence on $\tau$ and $\chi$ from $K_{0}$; we
therefore define the space- and time-averaged drag reduction as

$\langle{DR}_{0}\rangle_{\chi}(\tau)=\frac{1}{2L_{x}}\int_{\chi=-L_{x}}^{L_{x}}{DR}_{0}\,\text{d}\chi,\quad\overline{DR}_{0}(\chi)=\frac{1}{\mathcal{T}}\int_{\tau=0}^{\mathcal{T}}{DR}_{0}\,\text{d}\tau,$
(59a, b)

where $[-L_{x},\,L_{x}]$ is the length of the channel and $[0,\,\mathcal{T}]$
is the time interval over which the drag reduction is measured. Given a
surfactant-flux distribution $K_{0}$, we can calculate the corresponding drag-
reduction distribution ${DR}_{0}$ by solving the surfactant transport
equations (44, 46, 46) for each $K_{0}$ to get $\Gamma_{0}(x;\,K_{0})$, and
then use (58) to calculate ${DR}_{0}$. To evaluate the effective slip length
$\lambda_{e}$ over a plastron (Tomlinson et al., 2023a), we integrate the
leading-order streamwise momentum equation (2.3.1b) for an equivalent channel
with the mixed boundary conditions (30a–ca, 33a, ba) replaced by
$\lambda_{e}u_{0y}-u_{0}=0$. We obtain $u_{0}=\check{U}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 volume 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 volume
flux of the equivalent channel with the volume flux of the SHS channel,
$\tilde{Q}p_{0x}+\Ma\bar{Q}\Gamma_{0x}=1$, we find

$\lambda_{e}=\frac{{DR}_{0}(\Delta p_{I}-\Delta p_{U})}{P_{z}\Delta
p_{I}(\Delta p_{U}{DR}_{0}+\Delta p_{I}(1-{DR}_{0}))},$ (60)

which can be used to convert results from ${DR}_{0}$ to $\lambda_{e}$.
Therefore, with $\lambda_{e}$ and ${DR}_{0}$ being directly related to
$K_{0}$, the governing equations, (52)–(4), can also be rewritten as initial
boundary value problems for either the effective slip length or drag
reduction, that vary over the long length scale and slow time scale of the
channel.

##  3 Results

In §3.1–3.4, we investigate how the surfactant flux ($K_{0}$), drag reduction
(${DR}_{0}$), propagation speed ($a$), streamwise velocity ($u_{0}$) and
surfactant concentration ($c_{0}$ and $\Gamma_{0}$) vary with the bulk
diffusion ($\alpha$), partition coefficient ($\beta$), surfactant strength
($\gamma$), surface diffusion ($\delta$), exchange strength ($\nu$), bulk
capacitance ($\theta$), surface capacitance ($\sigma$) and streamwise gas
fraction ($\phi_{x}$). We discuss the quantities that vary over the length of
the channel ($K_{0}$, ${DR}_{0}$ and $a$) and their effect on the flow and
surfactant transport ($u_{0}$, $c_{0}$ and $\Gamma_{0}$) over each period.
Using asymptotic (detailed in Appendix A) and numerical (Appendix B)
techniques, we evaluate the solution to (52, 4) in the main regions and
boundaries illustrated in figure 3 and discussed in §2.3.3: the Marangoni-
dominated (M) region (§3.1); the advection-dominated (A) region (§3.2); the
diffusion-dominated (D) region (§3.2); the advection–Marangoni (AM) boundary
(§3.3); and the diffusion–Marangoni (DM) boundary (§3.4). Shocks in the
surfactant-flux and drag-reduction distribution can arise in regions AM and
DM. Throughout §3, we fix the length of the channel $L_{x}=1$, channel-height-
to-streamwise-period ratio $\epsilon=0.1$, spanwise gas fraction
$\phi_{z}=0.5$ and spanwise period width $P_{z}=1$. We construct asymptotic
solutions for any $L_{x}$, $\phi_{z}$ and $P_{z}$ in Appendix A where
$\epsilon\ll 1$.

###  3.1 Marangoni–dominated region

####  3.1.1 Flow and surfactant flux transport at the channel scale

Figure 6(a) shows how a bolus of surfactant flux, using initial and boundary
conditions (4), is advected along the length of the channel at increasing
times. Taking $K_{b}=1/2$ in (4), $K_{0}$ remains sufficiently large for the
flow to be everywhere in the Marangoni-dominated (M) region. As the surfactant
flux and surfactant concentration increase, the leading-order drag reduction
in figure 6(b) decreases, implying that the liquid–gas interface is more
immobilised. We plot in figure 6(a, b) asymptotic solutions for the leading-
order surfactant flux and drag reduction, (68), derived in the limit where
both Marangoni effects and bulk–surface exchange are strong compared to
advection and diffusion (see region M in figure 3a). As the flux of surfactant
is conserved in (52), the space-averaged drag reduction
($\langle{{DR}_{0}}\rangle_{\chi}$) defined in (58a) is constant provided the
bolus of surfactant flux remains in the channel for the given time interval.
However, if the bolus of surfactant flux is advected out of the channel, then
the space- and time-averaged drag reduction
$\langle\overline{{DR}_{0}}\rangle_{\chi}$ defined in (58) is maximised by
maximising the propagation speed ($a_{\text{M}}$) in region M. In figure 6(b),
$\langle\overline{{DR}_{0}}\rangle_{\chi}=\langle{DR}_{0}\rangle_{\chi}=0.06$
for $\chi\in[-1,\,1]$ and $\tau\in[0,\,1]$; the drag reduction is small as the
liquid–gas interface is mostly immobilised with the background surfactant
flux.

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

Figure 6:  Quantities of interest as a bolus of surfactant flux, (4), passes
through the channel in the Marangoni-dominated region (M), where
$\alpha=\beta=\delta=1$, $\gamma=10$, $\nu=100$, $\phi_{x}=0.5$, $L_{x}=1$ and
$K_{b}=1/2$. (a) Surfactant flux $K_{0}$ and (b) drag reduction ${DR}_{0}$
versus distance along the channel $\chi$ at different times $\tau=0$, $0.2$,
$0.4$, $0.6$, $0.8$ and $1$, computed using (58, 4, 68) for $\theta=\sigma=1$.
The dashed horizontal line in panel (b) is the space-averaged drag reduction
$\langle{DR}_{0}\rangle_{\chi}$, calculated using (58). (c) Bulk $c_{0}$
(dotted lines) and interfacial $\Gamma_{0}$ (solid lines) concentration fields
across a unit cell $\mathcal{D}_{1}\cup\mathcal{D}_{2}$ and (d) streamwise
slip velocity $u_{0}(x,\,0,\,0)$ across a plastron (domain $\mathcal{D}_{1}$),
located at $\chi=-0.5$ (dotted vertical line in panels a and b) when
$\tau=0.2$, $0.3$, $0.4$, $0.5$, $0.6$ and $0.7$, and calculated using (2.3.1,
46).

In region M, the bolus of surfactant flux, (4), is advected downstream by the
flow with a constant propagation speed, $a_{\text{M}}\approx
1/(\theta+\sigma\phi_{x})$ (derived in Appendix A.1), and therefore, the
distribution of surfactant flux and drag reduction in figure 6(a, b) do not
change shape as the bolus moves through the channel. For $\sigma\ll\theta$, we
have that $a_{\text{M}}\approx 1/\theta$, which corresponds to a dimensional
speed $\hat{U}_{m}=\hat{Q}/(4\hat{P}_{z}\hat{H})$. This is the cross-channel-
averaged bulk propagation speed and indicates bulk-dominated surfactant
transport. For a fixed $\theta$, adsorption at the liquid–gas interface causes
the propagation speed of the bolus of surfactant flux to fall significantly
compared to the cross-channel-averaged bulk propagation speed, reducing to
$a_{\text{M}}\approx 1/(\sigma\phi_{x})$ for $\sigma\gg\theta$. Dimensionally,
this corresponds to a reduction in the bulk propagation speed by a factor
$\phi_{x}\phi_{z}L_{d}$ and indicates surface-dominated transport, where
$L_{d}=\hat{K}_{a}/(\hat{H}\hat{K}_{d})$ is the normalised surfactant
depletion length and $\phi_{x}\phi_{z}$ is the area gas fraction of the SHS.
Hence, the propagation of disturbances to the surfactant concentration field
is significantly slower for more insoluble surfactants (large $L_{d}$) and
when the area of adsorption, i.e. the liquid–gas interface
$0<\phi_{x}\phi_{z}<1$, is maximised. For either bulk or surface-dominated
transport, the propagation speed in region M decreases with the bulk
($\theta=4P_{z}$) and surface capacitance ($\sigma=4\phi_{z}P_{z}L_{d}$). For
a fixed volume flux of fluid, an increase in the cross-channel area (spanwise
surface area) will reduce the streamwise bulk (surface) velocity, implying
that the bolus of surfactant flux will be advected slower throughout both
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$.

####  3.1.2 Flow and surfactant transport at the period scale

The magnitude of the background surfactant flux, $K_{b}=1/2$, means that the
liquid–gas interface is almost immobile along the entire SHS and the leading-
order bulk and interfacial concentrations ($\Gamma_{0}$ and $c_{0}$) in each
period are linear with a shallow gradient, as shown in equation (65) and
figure 6(c). In (65), we see that $c_{0}$ depends linearly on $K_{0}$ at
leading-order, however, $\Delta c_{0}$ does not, as the liquid–gas interface
is already immobilised when $K_{0}=O(1)$. As the bolus of surfactant flux
passes over an individual plastron and $K_{0}$ varies from $1/2$ up to $1$ and
back down to $1/2$, the concentration rises (from times $\tau=0.2$ to
$\tau=0.4$) and then falls (from times $\tau=0.4$ to $\tau=0.7$). We observe
adsorption and desorption inside boundary layers around $x=\pm\phi_{x}=\pm
0.5$ where the bulk and interfacial concentrations deviate from each other,
generating local surfactant gradients that reduce the streamwise slip velocity
($u_{0}$) close to the stagnation points ($x=\pm\phi_{x}=\pm 0.5$) in figure
6(d). The streamwise slip velocity inherits a dependence on $K_{0}$ through
$c_{0}$, as more surfactant increases the amount of immobilisation at the
liquid–gas interface. Thus, the streamwise velocity at the interface falls and
then rises as the bolus passes over an individual plastron (see curves from
$\tau=0.2$ to $\tau=0.7$ in the graph in figure 6d). As mentioned in §2.3.1,
the present long-wave theory does not capture inner regions close to the
stagnation points where the streamwise velocity satisfies $u_{0}=0$,
explaining the non-zero value exhibited by $u_{0}(x,0,0)$ in figure 6(d) near
$x=\pm\phi_{x}=\pm 0.5$. As described in (32a), we have instead imposed no
flux of surfactant in the streamwise direction at these stagnation points.

###  3.2 Advection and diffusion–dominated regions

We also briefly discuss asymptotic results (derived in Appendix A.2 and A.3)
for the advection-dominated (A) and diffusion–dominated (D) regions (see
regions A and D in figure 3a). When bulk–surface exchange is strong, the bolus
of surfactant flux in (4) propagates in a similar manner to figure 6(a), but
is advected with speeds
$a_{\text{A}}\approx(\beta+1)/(\phi_{x}\left(\sigma+\theta\right)+\theta(1-\phi_{x})(\beta+1))$
and
$a_{\text{D}}\approx(\alpha(1+\beta)+\delta(1-\phi_{x}))/((\theta+\sigma\phi_{x})(\alpha+\delta(1-\phi_{x})))$
in the A and D regions, respectively, and with a new parameter dependency on
the partition coefficient $\beta=2L_{d}\tilde{q}/\tilde{Q}$. In regions A and
D, the propagation speed increases with $\beta$. For $\beta\gg 1$ ($\beta\ll
1$), the flux of surfactant along the liquid–gas interface is greater
(smaller) than the flux of surfactant in the bulk, which is non-
dimensionalised to unity in (46). Hence, as $L_{d}$ grows, the localised
concentration of surfactant will be advected faster along the liquid–gas
interface, and therefore, the surfactant will be advected faster throughout
$\mathcal{D}_{1}$. When bulk–surface exchange is weak, e.g. regions
M${}_{\text{E}}$ and D${}_{\text{E}}$ depicted in figure 3(b), the propagation
speed is the same in all regions M, A and D (see Appendices A.1, A.2 and A.3),
such that $a_{\text{M}}=a_{\text{D}}=a_{\text{A}}\approx
1/(\theta+\sigma\phi_{x})$.

###  3.3 Advection–Marangoni boundary

####  3.3.1 Flow and surfactant flux transport at the channel scale

In figure 7(a), we use asymptotic (92) and numerical solutions (outlined in
Appendix B) to examine the spatio–temporal evolution of ${DR}_{0}$ along the
length of the channel at the advection–Marangoni (AM) boundary, where
Marangoni effects, interfacial advection and bulk–surface exchange are strong
compared to diffusion (see the AM boundary in figure 3a). Taking $K_{b}=0.01$,
we focus on the case where $\gamma K_{0}\leq 2\phi_{x}\beta$, which places the
surfactant profile at the plastron scale in the stagnant-cap regime (Tomlinson
et al., 2023a). The flow advects the bolus of surfactant flux (4) through the
channel with a propagation speed $a_{\text{AM}}\approx
2\beta/(\gamma(\sigma+\theta)K_{0}+2\beta\theta(1-\phi_{x}))$ that depends on
the local flux $K_{0}$. As the bolus propagates in $(\chi,\,\tau)$-space, the
wave steepens at the upstream edge (rearside of the wave) and ultimately a
shock forms at some location and time along the channel, which we discuss
further in §3.3.2. At the AM boundary, the liquid–gas interface is partially
immobilised and the space-averaged drag reduction in figure 7(a) is given by
$\langle\overline{{DR}_{0}}\rangle_{\chi}=\langle{DR}_{0}\rangle_{\chi}=0.95$
for $\chi\in[-1,\,1]$ and $\tau\in[0,\,0.29]$, which is close to the shear-
free value for reasons that we now discuss. In Appendix A.4, we derive that
${DR}_{0}=0.5$ when $\gamma=\beta\phi_{x}/K_{0}$. Hence, for $\beta=10$,
$\gamma=5$ and $\phi_{x}=0.5$ in figure 7(a), we transition from region A (as
${DR}_{0}>0.5$ for $K_{0}<1$) to region M only when $K_{0}=1$, and therefore,
the flow is mostly dominated by advection.

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

Figure 7:  Quantities of interest as a bolus of surfactant flux, (4), forms a
shock while passing through the channel at the advection–Marangoni boundary
(AM), where $\alpha=\delta=0.1$, $\beta=10$, $\epsilon=0.1$, $\gamma=5$,
$\nu=1000$, $\phi_{x}=0.5$, $L_{x}=1$ and $K_{b}=0.01$. (a) Drag reduction
${DR}_{0}$ versus distance along the channel $\chi$ at different times
$\tau=0$, $0.04$, $0.08$, $0.12$, $0.17$, $0.21$, $0.25$ and $0.29$, computed
using (58, 4), (92) (solid lines, no streamwise diffusion) and (46) (dashed
lines, weak streamwise diffusion), for $\theta=\sigma=1$. The dashed
horizontal line in panel (a) is the space-averaged drag reduction
$\langle{DR}_{0}\rangle_{\chi}$, calculated using (58). (b) Propagation speed
$a_{\text{AM}}$ for varying surfactant flux $K_{0}$ and capacitance
$\theta=\sigma$, computed using (92). (c) Characteristic curves
$\chi=\chi(\tau;\,\xi)$ intersecting and forming a shock at
$(\chi_{b},\,\tau_{b})$ that propagates along $\chi=\chi_{s}(\tau)$, computed
using (4, 55, 56, 57). (d) Bulk $c_{0}$ (dotted lines) and interfacial
$\Gamma_{0}$ (solid lines) concentration fields across a unit cell
$\mathcal{D}_{1}\cup\mathcal{D}_{2}$ and (e) streamwise slip velocity
$u_{0}(x,\,0,\,0)$ across a plastron (domain $\mathcal{D}_{1}$), located at
$\chi=-0.7$ (dotted line in panel a) when $\tau=0$, $0.02$, $0.05$, $0.07$,
$0.08$ and $0.1$, and calculated using (2.3.1, 46).

Figure 7(b) shows how the propagation speed depends on the spatio–temporal
evolution of the bolus of surfactant at the AM boundary, using (92). The wave
steepening in the drag-reduction distribution observed in figure 7(a) occurs
because there is wave steepening in the surfactant-flux distribution, and a
small surfactant flux is advected faster than a large surfactant flux.
Physically speaking, greater surfactant concentrations at a given location in
the channel imply that the liquid–gas interface is more immobilised than it is
for small surfactant concentrations. This implies that the streamwise slip
velocity and thus the propagation speed decreases with increasing $K_{0}$. The
dependence of $a_{\text{AM}}$ on $\beta$, $\theta$ and $\sigma$ is similar as
for $a_{\text{M}}$, which is discussed in §3.1. When $K_{0}\ll 1$, we have
that $a_{\text{AM}}\approx 1/(\theta(1-\phi_{x}))$, the cross-channel-averaged
bulk propagation speed enhanced by a factor proportional to the streamwise
groove length $1-\phi_{x}$.

####  3.3.2 Shock formation and regularisation via streamwise diffusion

The time taken for the shock to form in the surfactant-flux distribution,
$\tau_{b}$, increases as the partition coefficient ($\beta$) decreases and the
bulk and surface capacitance parameters ($\theta$ and $\sigma$) increase. The
location in space where the shock in the surfactant-flux distribution forms,
$\chi_{b}$, does not vary with respect to $\theta$ and $\sigma$; however,
$\chi_{b}$ increases with decreasing $\beta$, for reasons that we now explain.
Recall that, as $\beta$ increases, the propagation speed increases because the
interfacial surfactant flux increases and the bulk surfactant flux is fixed,
and as $\theta$ and $\sigma$ decrease, the propagation speed increases because
the cross-sectional area reduces and the volume flux is fixed. A larger
propagation speed means that the difference between the total flux of
surfactant when the concentration is small and large is greater (see figure
7b), and therefore, the larger difference causes the wave to steepen faster
and the shock forms earlier. Taking $K_{b}=0.01$ in the initial distribution
of surfactant flux prescribed in (4), we can use (55)–(56) to calculate that a
shock will form when $(\chi_{b},\,\tau_{b})\approx(0.34,\,0.13)$; as
illustrated by the red dot in figure 7(c). We can evaluate the shock speed and
location $\chi=\chi_{s}(\tau)$ for times $\tau>\tau_{b}$ using the
Rankine–Hugoniot condition in (57), as shown by the green curve in figure
7(c). The solid curves for $\tau=0.17$, $0.21$, $0.25$ and $0.29$ in figure
7(a) are composed of the solution found using the method of characteristics
combined with the above results for $\chi=\chi_{s}(\tau)$. For $\tau>0.25$,
the minimum in the drag reduction starts to increase owing to nonlinear
interaction with the SHS.

Alternatively, we can solve the advection–diffusion equation in (50)
numerically using the method in Appendix B. Diffusion over the length of the
channel regularises the flow in the vicinity of the shock. This methodology
allows us to compute the distribution of $DR_{0}$ for $\tau\geq\tau_{b}$, as
shown by the dotted curves for $\tau\geq 0.17$ in figure 7(a). The solutions,
inclusive of a small amount of diffusion (dotted lines), remain fairly close
to the solutions without diffusion (solid lines). However, diffusion
marginally reduces the maximum (minimum) amplitude of the surfactant-flux
(drag-reduction) distribution and widens the surfactant-flux (drag-reduction)
distribution as time progresses.

####  3.3.3 Flow and surfactant transport at the period scale

As the bolus of surfactant flux passes over a given plastron at the AM
boundary, the proportion of the liquid–gas interface that is shear-free at the
upstream end and no-slip at the downstream end varies; as seen from the
concentration field in figure 7(d) and the streamwise slip velocity in figure
7(e). From (87), the length of the shear-free region increases as $K_{0}$
decreases, $\beta$ increases and $\gamma$ decreases, as less surfactant is
held at higher concentrations at the downstream stagnation point,
$x=\phi_{x}=0.5$. The maximum concentration at the downstream end of the
interface increases with $K_{0}$ and $\Delta c_{0}=K_{0}$. Therefore,
${DR}_{0}$ attains a minimum when $K_{0}$ is at its maximum, where the
amplitude (difference between the maximum and minimum value of $c_{0}$) and
length (distance between the end of the plastron and the start of no-slip
region on the plastron) of the stagnant cap are greatest. At the leading edge
of the bolus of surfactant flux (where the expansion wave is), large stagnant
caps slowly transition into small ones over a long range in $\chi$, and at the
trailing edge of the bolus of surfactant flux (on the rear side, where the
shock is), small stagnant caps rapidly transition into large ones over a
shorter range in $\chi$.

###  3.4 Diffusion–Marangoni boundary

We derive asymptotic solutions at the diffusion–Marangoni (DM) boundary in
Appendix A.5, where both diffusion and Marangoni effects dominate over
interfacial advection (see the DM boundary in figure 3a). The asymptotic
solution for $a_{\text{DM}}$ has a complex dependence on $\alpha$, $\beta$,
$\delta$, $\theta$, $\sigma$ and $K_{0}$ that is different to $a_{\text{AM}}$.
Nevertheless, we find that $a_{\text{DM}}$ exhibits similar trends as
$a_{\text{AM}}$ with respect to $\alpha$, $\beta$, $\delta$, $\theta$,
$\sigma$ and $K_{0}$, and furthermore, because of its nonlinear dependence on
$K_{0}$, wave-steepening effects can lead to shock formation in the
surfactant-flux and drag-reduction distribution. When bulk–surface exchange is
weak, e.g. at the DM${}_{\text{E}}$ boundary depicted in figure 3(b), we find
that $a_{\text{AM}}=a_{\text{DM}}\approx 1/(\theta+\sigma\phi_{x})$.

##  4 Discussion

Superhydrophobic surfaces (SHSs) can be contaminated by trace amounts of
surfactant, which may reduce their drag-reducing potential for applications in
microchannels and marine hydrodynamics (Peaudecerf et al., 2017; Tomlinson et
al., 2023b). Field studies have shown that environmental surfactant
concentrations can vary significantly in space and time (Frossard et al.,
2019). Some of these variations can occur over length scales comparable to the
size of the SHS, causing local variations in slip and drag reduction. To
explore the impact of these local spatio–temporal variations in surfactant
concentration, we have derived an asymptotic theory to model the unsteady
transport of soluble surfactant in a laminar pressure-driven channel flow
bounded between two SHSs. The SHSs are textured with grooves that are periodic
in the streamwise and spanwise directions. Exploiting the multiple length
scales in the problem, we have derived and solved a quasi-steady nonlinear
advection–diffusion equation for surfactant concentration transport over
moderate length scales, which is coupled to an unsteady nonlinear
advection–diffusion equation for surfactant flux transport over large length
scales. The governing partial differential equations for surfactant flux
transport can be rewritten in terms of the slip length or drag reduction, key
quantities of interest for practical applications. When there is a disturbance
to the surfactant flux that varies over a large number of periods (similar to
those measured in the ocean by Frossard et al., 2019) but over a length
smaller or similar to the size of the SHS, our model allows us to make
predictions about the propagation speed, the shape of the disturbance, and the
evolving distribution of slip length and drag reduction. Furthermore, in
certain regions of parameter space, higher surfactant concentrations can lead
to surface immobilisation and therefore slower surfactant flux transport,
leading to wave-steepening and shock formation in the distribution of
surfactant flux, slip length and drag reduction.

We have investigated the transport of the surfactant flux and the
corresponding reduction in drag along the channel length in distinct
asymptotic regimes (figure 3), defined by the relative strength of the bulk
diffusion ($\alpha$), partition coefficient ($\beta$), surfactant strength
($\gamma$), surface diffusion ($\delta$), exchange strength ($\nu$) and
background surfactant flux ($K_{b}$). Extending the results of Tomlinson et
al. (2023a), we also identify the bulk capacitance ($\theta$) and surface
capacitance ($\sigma$) as key parameters, which quantify the bulk and surface
response to time-dependent changes in the surfactant flux, respectively. The
speed of propagation of disturbances to the surfactant flux across different
asymptotic regimes is summarised in tables 3 and 3. If a bolus of surfactant
flux is injected into the channel, the propagation speed is constant in the
Marangoni- (M), advection- (A) and diffusion-dominated (D) regimes, and the
shape of the bolus of surfactant flux does not change appreciably along the
length of the channel (figure 6). In region M, the interfacial concentration
profile is linear, the liquid–gas interface is immobilised and there is
negligible drag reduction (${DR}_{0}\ll 1$); in region A, the interfacial
concentration is constant and then increases in a downstream boundary layer,
the liquid–gas interface is shear-free and there is substantial drag reduction
($1-{DR}_{0}\ll 1$); and in region D, interfacial concentration profile is
uniform, the liquid–gas interface is shear-free and there is substantial drag
reduction ($1-{DR}_{0}\ll 1$). However, at the advection–Marangoni (AM) and
diffusion–Marangoni (DM) boundaries, the values of $DR_{0}$ can span the whole
range from 0 to 1 along the surface, depending on the local surfactant flux.
Here, a bolus of surfactant flux steepens at its upstream end (rearside of the
distribution) as smaller concentrations of surfactant are advected faster than
larger concentrations, because the liquid–gas interface is more mobile at
lower concentrations, resulting in the formation of a shock (figure 7). We
anticipate that the results observed here for a Gaussian initial distribution
in the surfactant flux would be applicable to other distribution profiles due
to the dominant advective nature of the transport at the channel scale.
Increasing profiles in the surfactant flux would behave similar to the rear
part of the Gaussian distribution (the increasing part) that we studied,
exhibiting wave steepening and potential shocks. Decreasing profiles would
behave similar to the front of the Gaussian distribution, exhibiting wave
expansion. These effects arise only in the strong-exchange limit (large
$\nu$). In contrast, when exchange is weak, the propagation speed given in
table 3 is the same for all $\alpha$, $\beta$, $\gamma$, $\delta$, $\theta$,
$\sigma$ and $\phi_{x}$.

Strong exchange $\displaystyle\nu\gg\max(1,\,\alpha,\,\beta,\,\delta)$ |
---|---
Region | Parameter space | Propagation speed ($a$) | Regime
M |  $\displaystyle\gamma\gg\max\left(1,\,\alpha,\,\beta,\,\delta\right)$, $K_{b}>2\phi_{x}\beta/\gamma$ | $\displaystyle\frac{1}{\theta+\sigma\phi_{x}}$ | ${DR}_{0}\ll 1$
A | $\max(\alpha,\,\delta)\ll 1,\quad\gamma\ll\min(1,\,\beta)$ | $\displaystyle\frac{(\beta+1)}{\phi_{x}(\sigma+\theta)+\theta(1-\phi_{x})(\beta+1)}$ | $1-{DR}_{0}\ll 1$
D | $\min(\alpha,\,\delta)\gg\max(1,\,\gamma)$ | $\displaystyle\frac{\alpha(1+\beta\phi_{x})+\delta(1-\phi_{x})}{(\theta+\sigma\phi_{x})(\alpha+\delta(1-\phi_{x}))}$ | $1-{DR}_{0}\ll 1$
AM |  $\min(\beta,\,\gamma)\gg 1\gg\max(\alpha,\,\delta)$, $K_{b}\leq 2\phi_{x}\beta/\gamma$ | $\displaystyle\frac{2\beta}{\gamma(\sigma+\theta)K_{0}+2\beta\theta(1-\phi_{x})}$ | $0<{DR}_{0}<1$
Weak exchange $\displaystyle\nu\ll\min(1,\,\alpha,\,\beta,\,\delta)$ |
Region | Parameter space | Propagation speed ($a$) | Regime
$\text{M}_{\text{E}}$ | $\gamma\gg\max(1,\,\alpha,\,\beta,\,\delta)$ | $\displaystyle\displaystyle\frac{1}{\theta+\sigma\phi_{x}}$ | ${DR}_{0}\ll 1$
$\text{A}_{\text{E}}$, $\text{D}_{\text{E}}$ | $\gamma\ll\min(1,\,\alpha,\,\beta,\,\delta)$ | $\displaystyle\displaystyle\frac{1}{\theta+\sigma\phi_{x}}$ | $1-{DR}_{0}\ll 1$

Table 2: Summary of the propagation speed and drag reduction regime in the
main regions analysed in the strong-exchange problem: the Marangoni-dominated
region (M), the advection-dominated region (A), the diffusion-dominated region
(D) and the advection-Marangoni (AM) boundary; and their analogues in the
weak-exchange problem: the $\text{M}_{\text{E}}$, $\text{A}_{\text{E}}$ and
$\text{D}_{\text{E}}$ regions. The propagation speed is expressed in terms of
the transport coefficients $\alpha$, $\beta$, $\gamma$, $\delta$, $\nu$,
$\theta$ and $\sigma$ given in (45).

Quantity | $\displaystyle\alpha$ | $\displaystyle\beta$ | $\displaystyle\gamma$ | $\displaystyle\delta$ | $\displaystyle\nu$ | $K_{b}$ | $\displaystyle\theta$ | $\displaystyle\sigma$
---|---|---|---|---|---|---|---|---
Prop. to | $\displaystyle\frac{\hat{D}\hat{P}_{z}\hat{H}}{\hat{P}_{x}\hat{Q}}$ | $\displaystyle\frac{\hat{K}_{a}}{\hat{H}\hat{K}_{d}}$ | $\displaystyle\frac{\hat{K}_{a}^{2}\hat{A}\hat{H}^{2}\max(\hat{K}(\hat{x},\,0))}{\hat{P}_{x}\hat{K}_{d}^{2}\hat{\mu}\hat{Q}^{2}}$ | $\displaystyle\frac{\hat{D}_{I}\hat{K}_{a}\hat{P}_{z}}{\hat{P}_{x}\hat{K}_{d}\hat{Q}}$ | $\displaystyle\frac{\hat{K}_{a}\hat{P}_{z}\hat{H}^{2}}{\hat{P}_{x}\hat{Q}}$ | $\displaystyle\frac{\min(\hat{K}(\hat{x},\,0))}{\max(\hat{K}(\hat{x},\,0))}$ | $\displaystyle\frac{\hat{P}_{z}}{\hat{H}}$ | $\displaystyle\frac{\hat{K}_{a}\hat{P}_{z}}{\hat{H}^{2}\hat{K}_{d}}$

Table 3: Summary of the dimensionless parameters appearing in table 3 that
affect the leading-order propagation speed, and their dependence on the
dimensional quantities characterising the flow, surfactant properties and
geometry (outlined in §2).

As a practical illustration, we now evaluate the propagation speeds ($a$)
presented in tables 3 and 3 using parameters characteristic of microchannel
applications. In the analysis that follows, the transport coefficients in (45)
have been appropriately adjusted for the geometry employed in Temprano-Coleto
et al. (2023), which is bounded by one SHS and solid wall rather than the two
SHSs considered in §2. In regions M, $\text{M}_{\text{E}}$,
$\text{A}_{\text{E}}$ and $\text{D}_{\text{E}}$, $\theta\approx 4$,
$\sigma\approx 20.8$, and therefore, the dimensionless propagation speed is
predicted to be $a\approx 0.04$. Using $\epsilon\hat{U}=2.4\times
10^{-4}$\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$ as the velocity scale, a
surfactant flux perturbation is advected out of the channel at approximately
$9.6\times 10^{-6}$\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$ when Marangoni
effects dominate, significantly slower than the fluid itself. For regions A
and D, $\alpha\approx 0.4$, $\beta\approx 3.7$, $\delta\approx 1$, and
therefore, $a\approx 0.19$. Thus, when advection or diffusion dominates, we
find that the surfactant flux disturbance is advected out of the channel at
approximately $4.5\times 10^{-5}$\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$,
approximately five times faster than the propagation speed in region M. The
difference in propagation speed is because shear-free liquid–gas interfaces
(regions A and D) that lack surfactant gradients give rise to greater
streamwise velocities than immobilised liquid–gas interfaces (region M) with
substantial surfactant gradients. We can then vary these parameters to
maximise the propagation speed in regions M, A and D, and therefore, if we
suppose a bolus of surfactant enters and leaves the channel in the measurement
time interval, we can evaluate and then minimise the space and time-averaged
drag reduction for microchannel applications.

Our theory rests on several assumptions, which we summarise below in order to
also suggest possible extensions to this study. First, the asymptotic
expansion requires $\epsilon=\hat{H}/\hat{P}_{x}\ll 1$ and
$\mathcal{E}=P_{x}/L_{x}\ll 1$. This seems reasonable based on the
microchannel configurations considered herein, e.g. $\hat{H}\approx 1\times
10^{-4}$\text{\,}\mathrm{m}$$, $\hat{P}_{x}\approx 1\times
10^{-2}$\text{\,}\mathrm{m}$$ and $\hat{L}_{x}\approx 1\times
10^{-1}$\text{\,}\mathrm{m}$$ in Temprano-Coleto et al. (2023). However, this
may need to be revised in other applications, e.g. marine hydrodynamics, where
the boundary layer grows over the surface of the vessel and $\hat{L}_{x}\gg
1\times 10^{-1}$\text{\,}\mathrm{m}$$. Second, we have only considered the
case where diffusion is strong enough to eliminate cross-channel concentration
gradients. The reader is referred to Tomlinson et al. (2023a) for a discussion
of the parameter regimes where cross-channel concentration gradients first
become important. Third, several potential physical complications can arise
when considering surfactant-contaminated superhydrophobic channels, such as
liquid–gas interface deformation, the interaction of the interior flow with
the external gas subphase (as discussed in Lee et al., 2016) or turbulence in
the outer flow field (as discussed in Tomlinson et al., 2023b), which have
been neglected in this study.

To summarise, we have shown how a disturbance to the surfactant concentration
field can undergo wave-steepening as it propagates under a laminar channel
flow bounded by SHSs. This nonlinear evolution is shared by the distributions
of the effective slip and drag reduction in microchannel applications,
emphasising the importance of treating these as dynamic quantities in time-
evolving flows.

## 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 Asymptotic solutions

###  A.1 Strong Marangoni effect: region M

Assume that $\beta=O(1)$, $\gamma\gg\max(1,\,\alpha,\,\delta)$ and
$\nu\gg\max(1,\,\alpha,\,\delta)$, so that exchange is strong,
$c_{0}=\Gamma_{0}$, and expand using $c_{0}=c_{00}+c_{01}/\gamma+...$. At
$O(\gamma)$, Marangoni effects are comparable to bulk advection and diffusion,
and (46) reduces to

$c_{00}c_{00x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00}-c_{00x}=K_{0}\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject
to}\quad c_{00}(\phi_{x}^{-})=c_{00}(\phi_{x}^{+}),\quad
c_{00}(-\phi_{x})=c_{00}(2-\phi_{x}),$ (61)

which gives $c_{00}=K_{0}$ in $\mathcal{D}_{1}\cup\mathcal{D}_{2}$. At $O(1)$,
Marangoni effects are comparable to advection and bulk diffusion, and (46)
gives

$\beta-c_{01x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{01}-c_{01x}=0\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject to}\quad
c_{01}(\phi_{x}^{-})=c_{01}(\phi_{x}^{+}),\quad
c_{01}(-\phi_{x})=c_{01}(2-\phi_{x}),$ (63)

such that $c_{01}=\beta(x-\phi_{x}(E+1)/(E-1))$ in $\mathcal{D}_{1}$ where
$E\equiv\exp(2(1-\phi_{x})/\alpha)$. Similarly, at $O(1/\gamma)$, (46) reduces
to

$K_{0}c_{02x}=(\beta+1)c_{01}-c_{01}c_{01x}-(\alpha+\delta)c_{01x}\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{02}-c_{02x}=0\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject to}\quad
c_{02}(\phi_{x}^{-})=c_{02}(\phi_{x}^{+}),\quad
c_{02}(-\phi_{x})=c_{02}(2-\phi_{x}),$ (65)

which gives $c_{02}=\beta x^{2}/(2K_{0})+\beta
x(\phi_{x}-\alpha-\delta-2\phi_{x}E/(E-1)))/K_{0}+M_{1}$ where $M_{1}$ is an
integration constant. Hence, we have that

$c_{0}=K_{0}+\beta(x-\phi_{x}(E+1)/(E-1))/\gamma+...,\quad\Delta
c_{0}=2\beta\phi_{x}/\gamma+....$ (66a–c)

Substituting (65) into (51), at leading order we have that

$C_{0}\approx 2(\theta+\sigma\phi_{x})K_{0},\quad A_{0}\approx
2(1+\beta\phi_{x})K_{0},\quad M_{0}\approx 2\beta\phi_{x}K_{0},\quad
D_{0}\approx 0,$ (67a–d)

where $K_{0}=K_{0}(\chi,\,\tau)$ and $D_{0}=O(1/\gamma)$. Then the advection
equation, (52) with $\lambda=0$, has the solution

$K_{0}=K_{0}(\chi-a_{\text{M}}\tau,\,0)\quad\text{where}\quad
a_{\text{M}}\approx\frac{1}{\theta+\sigma\phi_{x}}.$ (68)

Note that when $\sigma\rightarrow 0$, $a_{\text{M}}\rightarrow 1/\theta$.

Next, assume that $\nu\ll\min(1,\,\alpha,\,\delta)$ so that exchange is weak,
and expand using $c_{0}=c_{00}+c_{01}/\gamma+...$ and
$\Gamma_{0}=\Gamma_{00}+\Gamma_{01}/\gamma+...$. At $O(\gamma)$, Marangoni
effects are comparable to bulk advection and bulk diffusion, and (46) reduces
to

$(c_{00}-\alpha
c_{00x})_{x}=0,\quad\Gamma_{00}\Gamma_{00x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00}-\alpha c_{00x}=K_{0}\quad\text{in}\quad\mathcal{D}_{2},\\\
\text{subject to}\quad c_{00}(\phi_{x}^{-})=c_{00}(\phi_{x}^{+}),\quad
c_{00}(-\phi_{x})=c_{00}(2-\phi_{x}),\\\ c_{00}(\pm\phi_{x})-\alpha
c_{00x}(\pm\phi_{x})=K_{0},\quad\Gamma_{00}(\pm\phi_{x})\Gamma_{00x}(\pm\phi_{x})=0,$
(69)

which gives $\Gamma_{00}=c_{00}=K_{0}$ as
$\int_{x=-\phi_{x}}^{\phi_{x}}(\Gamma_{00}-c_{00})\,\text{d}x$ = 0. At $O(1)$,
Marangoni effects are comparable to advection and bulk diffusion, and (46)
gives

$(c_{01}-\alpha
c_{01x})_{x}=0,\quad\beta-\Gamma_{01x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{01}-\alpha c_{01x}=0\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject
to}\quad c_{01}(\phi_{x}^{-})=c_{01}(\phi_{x}^{+}),\quad
c_{01}(-\phi_{x})=c_{01}(2-\phi_{x}),\\\ c_{01}(\pm\phi_{x})-\alpha
c_{01x}(\pm\phi_{x})=0,\quad\beta-\Gamma_{01x}(\pm\phi_{x})=0,$ (71)

such that $\Gamma_{01}=\beta x$ and $c_{01}=0$ as
$\int_{x=-\phi_{x}}^{\phi_{x}}(\Gamma_{01}-c_{01})\,\text{d}x=0$. Hence,

$c_{0}=K_{0}+...,\quad\Gamma_{0}=K_{0}+\beta
x/\gamma...,\quad\Delta\Gamma_{0}=2\beta\phi_{x}/\gamma+....$ (72a, b)

Substituting (71) into (51) we recover (66a–c) and the propagation speed in
(68).

###  A.2 Strong advection: region A

In the advection–dominated (A) region, assume that $\beta=O(1)$, $\gamma\ll 1$
and $\nu\gg\max(1,\,\alpha,\,\delta)$, such that $c_{0}=\Gamma_{0}$ and expand
using $c_{0}=c_{00}+\gamma c_{01}+...$. At $O(1)$, advection is comparable to
diffusion, and (46) reduces to

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

subject to (A.1c, d), which gives
$c_{00}=K_{0}/(\beta+1)+K_{0}\beta\exp((\beta+1)(x-\phi_{x})/(\alpha+\delta))/(\beta+1)$
in $\mathcal{D}_{1}$ and $c_{00}=K_{0}$ in $\mathcal{D}_{2}$, for
$\max(\alpha,\,\delta)\ll 1$. Hence, we have that

$c_{0}=\frac{K_{0}}{\beta+1}\left(1+\beta\exp\left(\frac{(\beta+1)(x-\phi_{x})}{\alpha+\delta}\right)\right)+...,\quad\Delta
c_{0}=\frac{K_{0}\beta}{\beta+1}+....$ (74a, b)

From (73a, b), $c_{0}$ is constant over the upstream end of the liquid–gas
interface and increases exponentially in a boundary layer close to the
downstream stagnation point. The surfactant gradient, size of the boundary
layer and drag reduction increase with decreasing $K_{0}$, as the channel and
liquid–gas interface becomes less contaminated with surfactant, or with
increasing $\alpha$ and $\delta$, as diffusion eliminates the concentration
gradient. Substituting (73a, b) into (51), at leading order we have that

$C_{0}\approx\frac{2\phi_{x}(\sigma+\theta)K_{0}}{\beta+1}+2\theta(1-\phi_{x})K_{0},\quad
A_{0}\approx 2K_{0},\quad M_{0}\approx 0,\quad D_{0}\approx 0,$ (75a–d)

where $K_{0}=K_{0}(\chi,\,\tau)$ and $M_{0}=D_{0}=O(\gamma)$. Then the
advection equation, (52) with $\lambda=0$, has the solution

$K_{0}=K_{0}(x-a_{\text{A}}\tau,\,0)\quad\text{where}\quad
a_{\text{A}}\approx\frac{(\beta+1)}{\phi_{x}(\sigma+\theta)+\theta(1-\phi_{x})(\beta+1)}.$
(76)

Note that when $\beta\rightarrow 0$ and $\sigma\rightarrow 0$,
$a_{\text{A}}\rightarrow 1/\theta$.

Next, assume that $\nu\ll\min(1,\,\alpha,\,\delta)$ and expand using
$c_{0}=c_{00}+\gamma c_{01}+...$ and
$\Gamma_{0}=\Gamma_{00}+\gamma\Gamma_{01}+...$. At $O(1)$, diffusion is
comparable to advection, and (46) reduces to

$(c_{00}-\alpha
c_{00x})_{x}=0,\quad\beta\Gamma_{00}-\delta\Gamma_{00x}+c_{00}-\alpha
c_{00x}=K_{0}\quad\text{in}\quad\mathcal{D}_{1},\\\ c_{00}-\alpha
c_{00x}=K_{0}\quad\text{in}\quad\mathcal{D}_{2},\quad\text{subject to}\quad
c_{00}(\phi_{x}^{-})=c_{00}(\phi_{x}^{+}),\quad
c_{00}(-\phi_{x})=c_{00}(2-\phi_{x}),\\\ c_{00}(\pm\phi_{x})-\alpha
c_{00x}(\pm\phi_{x})=K_{0},\quad\beta\Gamma_{00}(\pm\phi_{x})-\delta\Gamma_{00x}(\pm\phi_{x})=0.$
(77)

which gives $c_{00}=K_{0}$ and $\Gamma_{00}=2\phi_{x}\beta
K_{0}\exp((\beta(\phi_{x}+x))/\delta)/(\delta(\exp((2\beta\phi_{x})/\delta)-1))$
as $\int_{x=-\phi_{x}}^{\phi_{x}}(\Gamma_{00}-c_{00})\,\text{d}x=0$. Hence,

$c_{0}=K_{0}+...,\quad\Gamma_{0}=\frac{2\phi_{x}\beta
K_{0}\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\Delta\Gamma_{0}=\frac{2\phi_{x}\beta
K_{0}}{\delta}+....$ (78a, b)

Substituting (77) into (51), we have

$C_{0}\approx 2(\theta+\sigma\phi_{x})K_{0},\quad A_{0}\approx
2(1+\beta\phi_{x})K_{0},\quad M_{0}\approx 0,\quad D_{0}\approx 2\phi_{x}\beta
K_{0},$ (79a–d)

and we recover (68).

###  A.3 Strong diffusion: region D

Assume that $\beta=O(1)$, $\min(\alpha,\,\delta)\gg\max(1,\,\gamma)$ and
$\nu\gg\max(1,\,\alpha,\,\delta)$ such that $c_{0}=\Gamma_{0}$. Let
$\delta=d\alpha$ where $d=O(1)$ and expand using
$c_{0}=c_{00}+c_{01}/\alpha+...$. At $O(\alpha)$, diffusion dominates and (46)
reduces to

$c_{00x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00x}=0\quad\text{in}\quad\mathcal{D}_{2},$ (80a, b)

subject to (A.1c, d). At $O(1)$, diffusion is comparable to advection, and
(46) gives

$(\beta+1)c_{00}-(1+d)c_{01x}=K_{0}\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00}-c_{01x}=K_{0}\quad\text{in}\quad\mathcal{D}_{2},$ (81a, b)

subject to (61c, d). Integrating (80a, b) over $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$ gives
$c_{00}=(K_{0}(\alpha+\delta(1-\phi_{x})))/(\alpha(\beta\phi_{x}+1)+\delta(1-\phi_{x})).$
Hence, we have that

$c_{0}=\frac{(\alpha+\delta(1-\phi_{x}))K_{0}}{(\alpha(\beta\phi_{x}+1)+\delta(1-\phi_{x}))}+...,\quad\Delta
c_{0}=\frac{2\beta\phi_{x}(1-\phi_{x})K_{0}}{\alpha(\beta\phi_{x}+1)+\delta(1-\phi_{x})}+....$
(82a–c)

From (81a, b), the surfactant concentration and drag increase linearly as the
flux of surfactant increases over the SHS, and the interface becomes
completely shear-free as $K_{0}\rightarrow 0$. Substituting (81a, b) into
(51), at leading order we have that

$C_{0}\approx\frac{2(\theta+\sigma\phi_{x})(\alpha+\delta(1-\phi_{x}))K_{0}}{\alpha(\beta\phi_{x}+1)+\delta(1-\phi_{x})},\quad
A_{0}\approx\frac{2(1+\beta\phi_{x})(\alpha+\delta(1-\phi_{x}))K_{0}}{\alpha(\beta\phi_{x}+1)+\delta(1-\phi_{x})},\\\
M_{0}\approx 0,\quad
D_{0}\approx\frac{2\beta\phi_{x}\delta(1-\phi_{x})K_{0}}{\alpha(\beta\phi_{x}+1)+\delta(1-\phi_{x})},$
(84)

where $K_{0}=K_{0}(\chi,\,\tau)$ and $M_{0}=O(1/\alpha)$. Then the advection
equation, (52) with $\lambda=0$, has the solution

$K_{0}=K_{0}(x-a_{\text{D}}\tau,\,0)\quad\text{where}\quad
a_{\text{D}}\approx\frac{\alpha(1+\beta\phi_{x})+\delta(1-\phi_{x})}{(\theta+\sigma\phi_{x})(\alpha+\delta(1-\phi_{x}))}.$
(85)

Note that when $\beta\rightarrow 0$, $\delta\rightarrow 0$ and
$\sigma\rightarrow 0$, $a_{\text{D}}\rightarrow 1/\theta$.

When $\nu\ll\min(1,\,\alpha,\,\delta)$, the expansion and solution are the
same as in region A, such that we recover (68).

###  A.4 Strong advection and strong Marangoni effect: the AM boundary

At the advection–Marangoni boundary, assume that $\beta=O(\gamma)$,
$\gamma\gg\max(\alpha,\,\delta)$, $\max(\alpha,\,\delta)\ll 1$ and
$\nu\gg\max(1,\,\alpha,\,\delta)$, such that $c_{0}=\Gamma_{0}$. Rescale
$\alpha=a/\gamma$, $\beta=b\gamma$ and $\delta=d/\gamma$, where $a$, $b$ and
$d$ are positive $O(1)$ constants. Expand using
$c_{0}=c_{00}+c_{01}/\gamma+...$. At $O(\gamma)$, Marangoni effects are
comparable to advection, and (46) reduces to

$c_{00}(b-c_{00x})=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00}=K_{0}\quad\text{in}\quad\mathcal{D}_{2},$ (86a, b)

subject to (A.1c, d). For $K_{0}/b\leq 2\phi_{x}$, (A.4) gives

$\displaystyle c_{00}$ $\displaystyle=0\quad\text{for all}\quad-\phi_{x}\leq
x\leq\phi_{x}-K_{0}/b,$ (87a) $\displaystyle c_{00}$
$\displaystyle=b(x-\phi_{x})+K_{0}\quad\text{for all}\quad\phi_{x}-K_{0}/b\leq
x\leq\phi_{x}.$ (87b)

At $O(1)$, Marangoni effects are comparable to advection, and (46) gives

$bc_{01}+c_{00}-c_{00}c_{01x}-c_{01}c_{00x}=K_{0}\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{01}=0\quad\text{in}\quad\mathcal{D}_{2},$ (88a, b)

subject to (61c, d). For $K_{0}/b\leq 2\phi_{x}$, (87) gives

$\displaystyle c_{01}$ $\displaystyle=K_{0}/b\quad\text{for
all}\quad-\phi_{x}\leq x\leq\phi_{x}-K_{0}/b,$ (89a) $\displaystyle c_{01}$
$\displaystyle=x-\phi_{x}-K_{0}\log(b(x-\phi_{x})/K_{0}+1)/b\quad\text{for
all}\quad\phi_{x}-K_{0}/b\leq x\leq\phi_{x}.$ (89b)

The solution for $K_{0}/b>2\phi_{x}$ is outlined in Tomlinson et al. (2023a)
and is not included here as it gives a similar propagation speed to region M.
Substituting (87)–(89) into (51), we have

$C_{0}\approx\frac{(\sigma+\theta)K_{0}^{2}}{2b}+\theta(2-2\phi_{x})K_{0},\quad
A_{0}\approx\frac{\gamma K_{0}^{2}}{2}+2K_{0},\\\ M_{0}\approx\frac{\gamma
K_{0}^{2}}{2},\quad D_{0}\approx 0,$ (91)

where $K_{0}=K_{0}(\chi,\,\tau)$. Then the advection equation, (52) with
$\lambda=0$, has the solution

$K_{0}=K_{0}(\chi-a_{\text{AM}}\tau,\,0)\quad\text{where}\quad
a_{\text{AM}}\approx\frac{2b}{K_{0}(\sigma+\theta)+2b\theta(1-\phi_{x})}.$
(92)

As $b\rightarrow K_{0}/(2\phi_{x})$, the bulk concentration $c_{00}\rightarrow
K_{0}$ and we recover (66a–c). As $b\rightarrow\infty$, the bulk concentration
$c_{00}\rightarrow 0$ everywhere except at $x\approx\phi$ where $c_{00}=K_{0}$
and we recover (74a, b). Furthermore, as $K_{0}\rightarrow 2\phi_{x}b$, the
propagation speed $a_{\text{AM}}\rightarrow a_{\text{M}}$, and as
$K_{0}\rightarrow 0$, the propagation speed $a_{\text{AM}}\rightarrow
a_{\text{A}}$ for $\beta\gg 1$.

Next, assume that $\nu/\epsilon^{2}\ll\min(1,\,\alpha,\,\delta)$ and expand
using $c_{0}=c_{00}+\gamma c_{01}+...$ and
$\Gamma_{0}=\Gamma_{00}+\gamma\Gamma_{01}+...$. At $O(\gamma)$, Marangoni
effects are comparable to advection, and (46) reduces to

$c_{00x}=0,\quad
b\Gamma_{00}-\Gamma_{00}\Gamma_{00x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00}=K_{0}\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject to}\quad
c_{00}(\phi_{x}^{-})=c_{00}(\phi_{x}^{+}),\quad
c_{00}(-\phi_{x})=c_{00}(2-\phi_{x}),\\\ c_{00}(\pm\phi_{x})=K_{0},\quad
b\Gamma_{00}(\pm\phi_{x})-\Gamma_{00}(\pm\phi_{x})\Gamma_{00x}(\pm\phi_{x})=0.$
(93)

For $K_{0}\gamma\leq\beta\phi_{x}$, (A.4) gives

$\displaystyle\Gamma_{00}$ $\displaystyle=0\quad\text{for
all}\quad-\phi_{x}\leq x\leq\phi_{x}-2(\phi_{x}K_{0}/b)^{1/2},$ (94a)
$\displaystyle\Gamma_{00}$
$\displaystyle=b(x-\phi_{x})+2(b\phi_{x}K_{0})^{1/2}\quad\text{for
all}\quad\phi_{x}-2(\phi_{x}K_{0}/b)^{1/2}\leq x\leq\phi_{x},$ (94b)

as $\int_{x=-\phi_{x}}^{\phi_{x}}(\Gamma_{00}-c_{00})\,\text{d}x=0$ and
$c_{00}=K_{0}$. At $O(1)$, Marangoni effects are comparable to advection, and
(46) gives

$c_{01x}=0,\quad
b\Gamma_{01}-\Gamma_{00}\Gamma_{01x}-\Gamma_{01}\Gamma_{00x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{01}=0\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject to}\quad
c_{01}(\phi_{x}^{-})=c_{01}(\phi_{x}^{+}),\quad
c_{01}(-\phi_{x})=c_{01}(2-\phi_{x}),\quad c_{01}(\pm\phi_{x})=0,\\\
b\Gamma_{01}(\pm\phi_{x})-\Gamma_{00}(\pm\phi_{x})\Gamma_{01x}(\pm\phi_{x})-\Gamma_{01}(\pm\phi_{x})\Gamma_{00x}(\pm\phi_{x})=0.$
(96)

For $K_{0}\gamma\leq\beta\phi_{x}$, (A.4) gives $\Gamma_{01}=0$ for all
$-\phi_{x}\leq x\leq\phi_{x}$ as
$\int_{x=-\phi_{x}}^{\phi_{x}}(\Gamma_{01}-c_{01})\,\text{d}x=0$ and
$c_{01}=0$. Substituting (94) into (51), we have

$C_{0}\approx 2(\theta+\sigma\phi_{x})K_{0},\quad A_{0}\approx
2(1+\beta\phi_{x})K_{0},\quad M_{0}\approx 2\beta\phi_{x}K_{0},\quad
D_{0}\approx 0,$ (97a–d)

and we recover (68).

###  A.5 Strong diffusion and strong Marangoni effect: the DM boundary

Assume that $\beta=O(1)$, $\gamma\gg 1$ and $\nu\gg\max(1,\,\alpha,\,\delta)$,
so that $c_{0}=\Gamma_{0}$. Rescale $\alpha=\mathcal{A}\gamma$ and
$\delta=d\gamma$, where $\mathcal{A}$ and $d$ are positive $O(1)$ constants.
Expand using $c_{0}=c_{00}+c_{01}/\gamma+...$. At $O(\gamma)$, Marangoni
effects are comparable to diffusion, and (46) reduces to

$c_{00}c_{00x}+(\mathcal{A}+d)c_{00x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00x}=0\quad\text{in}\quad\mathcal{D}_{2},$ (98a, b)

subject to (A.1c, d), which means that $c_{00}$ is constant. At $O(1)$,
Marangoni effects are comparable to advection and diffusion, and (46) gives

$(\beta+1)c_{00}-c_{00}c_{01x}-(\mathcal{A}+d)c_{01x}=K_{0}\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00}-\mathcal{A}c_{01x}=K_{0}\quad\text{in}\quad\mathcal{D}_{2},$ (99a, b)

subject to (61c, d). Integrating (A.5) over $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$ gives $c_{00}=c_{00}(K_{0})$ as the solution to the
quadratic equation

$(c_{00}-K_{0})(\phi_{x}-1)(\mathcal{A}+d+c_{00})=\mathcal{A}\phi_{x}(c_{00}(\beta+1)-K_{0}).$
(100)

We can then substitute $c_{00}$ back into (98a, b), which gives
$c_{01}=(x((\beta+1)c_{00}-K_{0}))/(\mathcal{A}+d+c_{00})+D_{1}$ in
$\mathcal{D}_{1}$, where $D_{1}$ is an integration constant. Hence,

$c_{0}=c_{00}(K_{0})+...,\quad\Delta
c_{0}=\frac{2\phi_{x}((\beta+1)c_{00}-K_{0}))}{\mathcal{A}+d+c_{00}}+....$
(101a, b)

At the DM boundary, the shear stress at the liquid-gas interface is not
sufficient to completely immobilise the interface and there is partial drag
reduction, as seen in (100). The concentration field is approximately constant
with a shallow gradient in $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$.
Substituting (100) into (51) gives

$C_{0}\approx 2(\theta+\sigma\phi_{x})c_{00},\quad\quad A_{0}\approx
2(1+\beta\phi_{x})c_{00},\\\
M_{0}\approx\frac{4\phi_{x}c_{00}((\beta+1)c_{00}-K_{0})}{\mathcal{A}+d+c_{00}},\quad
D_{0}\approx\frac{2\phi_{x}d((\beta+1)c_{00}-K_{0}))}{\mathcal{A}+d+c_{00}},$
(103)

where $c_{00}=c_{00}(K_{0})$. Then the advection equation, (52) with
$\lambda=0$, has the solution

$K_{0}=K_{0}(\chi-a_{\text{DM}}\tau,\,0)\quad\text{where}\quad
a_{\text{DM}}(K_{0})\approx\frac{A^{\prime}_{0}-M^{\prime}_{0}-D^{\prime}_{0}}{C^{\prime}_{0}},$
(104)

where primes denote derivatives of the functions defined in (101a, b). As
$\mathcal{A}\rightarrow 0$ and $d\rightarrow 0$, the bulk concentration
$c_{00}\rightarrow K_{0}$ and we recover (66a–c). As
$\mathcal{A}\rightarrow\infty$ and $d\rightarrow\infty$, the bulk
concentration
$c_{00}\rightarrow(\alpha+\delta(1-\phi_{x}))K_{0}/(\alpha(\beta\phi_{x}+1)+\delta(1-\phi_{x}))$
and we recover (82a–c). Furthermore, as $K_{0}\rightarrow\infty$, we need
$c_{00}=K_{0}$ to satisfy (100) and the propagation speed
$a_{\text{DM}}\rightarrow a_{\text{M}}$, and as $K_{0}\rightarrow 0$, we
linearise (100) (neglecting terms $O(c_{00}^{2},\,c_{00}K_{0})$) and the
propagation speed $a_{\text{DM}}\rightarrow a_{\text{D}}$.

Next, assume that $\nu\ll\min(1,\,\alpha,\,\delta)$ and expand using
$c_{0}=c_{00}+c_{01}/\gamma+...$ and
$\Gamma_{0}=\Gamma_{00}+\Gamma_{01}/\gamma+...$. At $O(\gamma)$, Marangoni
effects are comparable to diffusion, and (46) reduces to

$c_{00xx}=0,\quad\mathcal{A}c_{00x}+\Gamma_{00}\Gamma_{00x}+d\Gamma_{00x}=0\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{00x}=0\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject to}\quad
c_{00}(\phi_{x}^{-})=c_{00}(\phi_{x}^{+}),\quad
c_{00}(-\phi_{x})=c_{00}(2-\phi_{x}),\\\
c_{00x}(\pm\phi_{x})=0,\quad\Gamma_{00}(\pm\phi_{x})\Gamma_{00x}(\pm\phi_{x})+d\Gamma_{00x}(\pm\phi_{x})=0,$
(105)

which implies that $\Gamma_{00}$ and $c_{00}$ are constant. At $O(1)$,
Marangoni effects are comparable to surface advection and diffusion, and (46)
gives

$c_{01xx}=0,\quad
c_{00}-\mathcal{A}c_{01x}+\beta\Gamma_{00}-(\Gamma_{00}+d)\Gamma_{01x}=K_{0}\quad\text{in}\quad\mathcal{D}_{1},\quad
c_{01x}=0\quad\text{in}\quad\mathcal{D}_{2},\\\ \text{subject to}\quad
c_{01}(\phi_{x}^{-})=c_{01}(\phi_{x}^{+}),\quad
c_{01}(-\phi_{x})=c_{01}(2-\phi_{x}),\\\
c_{00}-\mathcal{A}c_{01x}(\pm\phi_{x})=K_{0},\quad\beta\Gamma_{00}-(\Gamma_{00}+d)\Gamma_{01x}(\pm\phi_{x})=0.$
(107)

such that $c_{00}=\Gamma_{00}=K_{0}$, $\Gamma_{01}=(\beta
xK_{0})/(d+K_{0})+\mathcal{C}$ and $c_{01}=\mathcal{C}$, as
$\int_{x=-\phi_{x}}^{\phi_{x}}(\Gamma_{00}-c_{00})\,\text{d}x=0$ and
$\int_{x=-\phi_{x}}^{\phi_{x}}(\Gamma_{01}-c_{01})\,\text{d}x=0$, where
$\mathcal{C}$ is an integration constant. Hence,

$c_{0}=K_{0}+...,\quad\Gamma_{0}=K_{0}+...,\quad\Delta\Gamma_{0}=\frac{2\phi_{x}\beta}{\gamma(d+K_{0})}+....$
(108a, b)

Substituting (107) into (51) we have

$C_{0}\approx 2(\theta+\sigma\phi_{x})K_{0},\ \ A_{0}\approx
2(1+\beta\phi_{x})K_{0},\ \ M_{0}\approx\frac{2\phi_{x}\beta
K_{0}^{2}}{d+K_{0}},\ \ D_{0}\approx\frac{2\phi_{x}\beta dK_{0}}{d+K_{0}},$
(109a–d)

and we recover (68).

##  Appendix B Numerical solution to the advection–diffusion equation

In §3, we solve (52) whilst retaining the $O(\epsilon^{2})$ secondary-
diffusion operator, partly for numerical convenience and partly to provide a
rational regularisation of shock-like structures that may arise. The unsteady
advection–diffusion equation is solved numerically using the method of lines
and a backwards-in-time and centered-in-space (BTCS) scheme. At each timestep,
we iterate $C_{0}$, $A_{0}$, $M_{0}$, $D_{0}$ and $D_{1}$, using $c_{0}$,
$\Gamma_{0}$ and $K_{0}$ at the current timestep (which are solved using the
Chebyshev collection technique described in Appendix A of Tomlinson et al.,
2023a). Discretising space such that $\chi_{i}=i\Delta\chi$ for $i=0$, $1$, …,
$N_{\chi}=2L_{x}\Delta\chi$, where $2L_{x}$ is length of the channel, we write
(52) at each $\chi_{i}$ for $i=1$, $2$, …, $N-1$:

$\left(\frac{\partial C_{0}}{\partial
K_{0}}\right)_{i}\left(\frac{\text{d}K_{0}}{\text{d}t}\right)_{i}=\left(\frac{\partial
A_{0}}{\partial K_{0}}-\frac{\partial M_{0}}{\partial K_{0}}-\frac{\partial
D_{0}}{\partial K_{0}}-\lambda\epsilon^{2}\frac{\partial^{2}D_{1}}{\partial
K_{0}^{2}}\frac{\partial
K_{0}}{\partial\chi}\right)_{i}\left(\frac{K_{0,\,i+1}-K_{0,\,i-1}}{2\Delta\chi}\right)\\\
-\lambda\epsilon^{2}\left(\frac{\partial D_{1}}{\partial
K_{0}}\right)_{i}\left(\frac{K_{0,\,i+1}-2K_{0,\,i}+K_{0,\,i-1}}{\Delta\chi^{2}}\right)+O(\Delta\chi^{3}),$
(110)

with periodic boundary conditions applied at interior ($i=1$, $N-1$), boundary
interior ($i=0$, $N$) and ghost nodes ($i=-1$, $N+1$):

$K_{0,\,0}=K_{0,\,N},\quad\frac{K_{0,\,1}-K_{0,\,-1}}{2\Delta\chi}=\frac{K_{0,\,N+1}-K_{0,\,N-1}}{2\Delta\chi}+O(\Delta\chi^{3}).$
(111)

Assembling (110)–(111) into a matrix problem, the advection–diffusion equation
in (52) reduces to solving the system of ODEs

$\frac{\text{d}\boldsymbol{K}_{0}}{\text{d}t}=\mathsfbi{A}(\boldsymbol{K}_{0})\boldsymbol{K}_{0}.$
(112)

We solve the problem in (112) using an implicit Euler scheme. Hence, defining
$\tau^{n}=n\Delta\tau$ for $n=1$, $2$, …, $N$, we have that

$(\mathsfbi{I}-\Delta\tau\mathsfbi{A}(\boldsymbol{K}_{0}^{n+1}))\boldsymbol{K}_{0}^{n+1}=\boldsymbol{K}_{0}^{n}.$
(113)

Note that (113) is nonlinear, therefore we initialise (113) with the solution
at the previous step $\boldsymbol{K}_{0}^{n}$ such that
$\mathsfbi{A}=\mathsfbi{A}(\boldsymbol{K}_{0}^{n})$, we then solve (113) for
$\boldsymbol{K}_{0}^{n+1}$ and substitute the new solution into
$\mathsfbi{A}=\mathsfbi{A}(\boldsymbol{K}_{0}^{n+1})$, until
$\boldsymbol{K}_{0}^{n+1}$ varies less than some specified tolerance.

## References

  * Bolognesi et al. (2014) Bolognesi, G, Cottin-Bizonne, C & Pirat, C 2014 Evidence of slippage breakdown for a superhydrophobic microchannel. Phys. Fluids 26, 082004.
  * Bond & Newton (1928) Bond, W N & Newton, D A 1928 Bubbles, drops, and Stokes’ law. Phil. Mag. 5 (30), 794–800.
  * Bottaro (2019) Bottaro, A 2019 Flow over natural or engineered surfaces: an adjoint homogenization perspective. J. Fluid Mech. 877, P1.
  * Darmanin & Guittard (2015) Darmanin, T & Guittard, F 2015 Superhydrophobic and superoleophobic properties in nature. Mater. Today 18 (5), 273–285.
  * Facchini et al. (2000) Facchini, M C, Decesari, S, Mircea, M, Fuzzi, S & Loglio, G 2000 Surface tension of atmospheric wet aerosol and cloud/fog droplets in relation to their organic carbon content and chemical composition. Atmos. Environ. 34, 4853–4857.
  * Frossard et al. (2019) Frossard, A A, Gérard, V, Duplessis, P, Kinsey, J D, Lu, X, Zhu, Y, Bisgrove, J, Maben, J R, Long, M S, Chang, R Y-W & others 2019 Properties of seawater surfactants associated with primary marine aerosol particles produced by bursting bubbles at a model air–sea interface. Environ. Sci. Technol. 53 (16), 9407–9417.
  * Frumkin & Levich (1947) Frumkin, A N & Levich, V G 1947 Effect of surface-active substances on movements at the boundaries of liquid phases. Zhur. Fiz. Khim. 21, 1183–1204.
  * Hourlier-Fargette et al. (2018) Hourlier-Fargette, A, Dervaux, J, Antkowiak, A & Neukirch, S 2018 Extraction of silicone uncrosslinked chains at air–water–polydimethylsiloxane triple lines. Langmuir 34, 12244–12250.
  * Kim & Hidrovo (2012) Kim, T J & Hidrovo, C 2012 Pressure and partial wetting effects on superhydrophobic friction reduction in microchannel flow. Phys. Fluids 24, 112003.
  * Lam et al. (2015) Lam, L S, Hodes, M & Enright, R 2015 Analysis of galinstan-based microgap cooling enhancement using structured surfaces. J. Heat Transf. 137, 091003.