A theory for the slip and drag of superhydrophobic surfaces with surfactant

TL;DR

This paper develops a theoretical model to predict how surfactants affect slip and drag on superhydrophobic surfaces in laminar flow, improving accuracy over surfactant-free models and enabling easier inclusion of surfactant effects in simulations.

Contribution

The paper introduces a linearized, scaling-based theory for surfactant effects on slip and drag in superhydrophobic surfaces, validated by numerical simulations.

Findings

Model accurately predicts surfactant influence on slip length.

Surfactant presence can cause overestimation of slip in previous models.

The theory reduces computational complexity by avoiding full simulations.

Abstract

Superhydrophobic surfaces (SHSs) have the potential to reduce drag at solid boundaries. However, multiple independent studies have recently shown that small amounts of surfactant, naturally present in the environment, can induce Marangoni forces that increase drag, at least in the laminar regime. To obtain accurate drag predictions, one must solve the mass, momentum, bulk surfactant and interfacial surfactant conservation equations. This requires expensive simulations, thus preventing surfactant from being widely considered in SHS studies. To address this issue, we propose a theory for steady, pressure-driven, laminar, two-dimensional flow in a periodic SHS channel with soluble surfactant. We linearise the coupling between flow and surfactant, under the assumption of small concentration, finding a scaling prediction for the local slip length. To obtain the drag reduction and interfacial shear, we find a series solution for the velocity field by assuming Stokes flow in the bulk and uniform interfacial shear. We find how the slip and drag depend on the nine dimensionless groups that together characterize the surfactant transport near SHSs, the gas fraction and the normalized interface length. Our model agrees with numerical simulations spanning orders of magnitude in each dimensionless group. The simulations also provide the constants in the scaling theory. Our model significantly improves predictions relative to a surfactant-free one, which can otherwise overestimate slip and underestimate drag by several orders of magnitude. Our slip length model can provide the boundary condition in other simulations, thereby accounting for surfactant effects without having to solve the full problem.

Equations

$$x=\frac{\hat{x}}{\hat{h}},y=\frac{\hat{y}}{\hat{h}},\mathbf{u}=\frac{\hat{\mathbf{u}}}{\hat{U}},\nabla p=\frac{\hat{\nabla }\hat{p}}{\hat{G}},c=\frac{\hat{c}}{{\hat{c}}_0},\Gamma =\frac{\hat{\Gamma }}{{\hat{\Gamma }}_m}.$$

$$\nabla \cdot \mathbf{u}$$

$$Re\nabla \cdot (\mathbf{u}\mathbf{u})$$

$$\nabla \cdot (\mathbf{u}c)$$

$$\frac{\mathrm{d}}{\mathrm{d}x}\left(u_I\Gamma \right)$$

$$S(c_I,\Gamma )=Bi\left(k{c}_I(1-\Gamma )-{\text{e}}^{A\Gamma }\Gamma \right).$$

$$\mathbf{u}(x,y)$$

$$c(x,y)$$

$$\mathrm{and}p(x,y)$$

$$\mathbf{u}$$

$$v$$

$$\frac{\partial c}{\partial y}$$

$$\frac{\mathrm{d}\Gamma }{\mathrm{d}x}$$

$$\frac{\chi k}{Pe}{\frac{\partial c}{\partial y}|}_I=S(c_I,\Gamma )\text{on the interface.}$$

$$\sigma =1+MaCa\left[\ln \left(1-\Gamma \right)-\frac{A{\Gamma }^2}{2}\right],$$

$${\frac{\partial u}{\partial y}|}_I$$

$$\text{that is}{\frac{\partial u}{\partial y}|}_I$$

$$S(c_I,\Gamma )\approx Bi\left(k{c}_I-\Gamma \right),$$

$${\frac{\partial u}{\partial y}|}_I\approx Ma\frac{\mathrm{d}\Gamma }{\mathrm{d}x}\text{on the interface}.$$

$$k^{\ast }\equiv Mak,{\Gamma }^{\ast }\equiv Ma\Gamma .$$

$$\nabla \cdot \mathbf{u}$$

$$Re\nabla \cdot (\mathbf{u}\mathbf{u})$$

$$\nabla \cdot (\mathbf{u}c)$$

$$\frac{\mathrm{d}}{\mathrm{d}x}\left(u_I{\Gamma }^{\ast }\right)$$

$$\mathbf{u}$$

$$v$$

$$\frac{\partial c}{\partial y}$$

$$\frac{\mathrm{d}{\Gamma }^{\ast }}{\mathrm{d}x}$$

$$\frac{\chi k^{\ast }}{Pe}{\frac{\partial c}{\partial y}|}_I$$

$${\frac{\partial u}{\partial y}|}_I$$

$${\frac{\partial c}{\partial y}|}_I=\frac{BiPe}{\chi k^{\ast }}\left(k^{\ast }{c}_I-{\Gamma }^{\ast }\right).$$

$${\int }_{-g/2}^{g/2}{\frac{\partial c}{\partial y}|}_I\mathrm{d}x=\frac{BiPe}{\chi k^{\ast }}{\int }_{-g/2}^{g/2}\left(k^{\ast }{c}_I-{\Gamma }^{\ast }\right)\mathrm{d}x=0,$$

$${\Gamma }^{\ast }\sim k^{\ast }\mp \Delta {\Gamma }^{\ast },c_I\sim 1\mp \Delta c_I,$$

$$\Delta {\Gamma }^{\ast }\sim g{\gamma }_{Ma},$$

$${\gamma }_{Ma}=\frac{1}{g}{\int }_{-g/2}^{g/2}{\frac{\partial u}{\partial y}|}_I\mathrm{d}x$$

$${\frac{\partial c}{\partial y}|}_I\sim \frac{\Delta c_I}{\delta }\sim \frac{BiPe}{\chi k^{\ast }}\left(g{\gamma }_{Ma}-k^{\ast }\Delta c_I\right),$$

$$\frac{\delta }{g}$$

$$\frac{\delta }{g}={\delta }_{0,3}{\left(1+{\delta }_{1,3}{g}^{2}Pe\right)}^{-1/3},$$

$$\Delta c_I\sim {\gamma }_{Ma}\frac{\frac{BiPe}{\chi k^{\ast }}g\delta }{1+\frac{BiPe}{\chi }\delta },$$

Full text

A theory for the slip and drag of superhydrophobic surfaces with surfactant

Julien R. Landela, François J. Peaudecerfb,111Present address: Institute of Environmental Engineering, Department of Civil, Environmental and Geomatic Engineering, ETH Zürich, 8093 Zürich, Switzerland, Fernando Temprano-Coletoc,

Frédéric Gibouc, Raymond E. Goldsteinb, and Paolo Luzzatto-Fegizc,222Email address: fegiz@engineering.ucsb.edu

a Department of Mathematics, University of Manchester,

Oxford Rd, Manchester M13 9PL, UK

b Department of Applied Mathematics and Theoretical Physics,

Centre for Mathematical Sciences, University of Cambridge,

Wilberforce Road, Cambridge CB3 0WA, UK

c Department of Mechanical Engineering, University of California, Santa Barbara,

Santa Barbara CA 93106, USA

March 13, 2024

Abstract

Superhydrophobic surfaces (SHSs) have the potential to reduce drag at solid boundaries. However, multiple independent studies have recently shown that small amounts of surfactant, naturally present in the environment, can induce Marangoni forces that increase drag, at least in the laminar regime. To obtain accurate drag predictions, one must solve the mass, momentum, bulk surfactant and interfacial surfactant conservation equations. This requires expensive simulations, thus preventing surfactant from being widely considered in SHS studies. To address this issue, we propose a theory for steady, pressure-driven, laminar, two-dimensional flow in a periodic SHS channel with soluble surfactant. We linearise the coupling between flow and surfactant, under the assumption of small concentration, finding a scaling prediction for the local slip length. To obtain the drag reduction and interfacial shear, we find a series solution for the velocity field by assuming Stokes flow in the bulk and uniform interfacial shear. We find how the slip and drag depend on the nine dimensionless groups that together characterize the surfactant transport near SHSs, the gas fraction and the normalized interface length. Our model agrees with numerical simulations spanning orders of magnitude in each dimensionless group. The simulations also provide the constants in the scaling theory. Our model significantly improves predictions relative to a surfactant-free one, which can otherwise overestimate slip and underestimate drag by several orders of magnitude. Our slip length model can provide the boundary condition in other simulations, thereby accounting for surfactant effects without having to solve the full problem.

1 Introduction

Superhydrophobic surfaces (SHSs) consist of hydrophobic coatings equipped with micro- or nano-scale textures, such that a layer of air (known as a “plastron”; see e.g. 1) is retained when the surface is submerged in water (see e.g. the reviews of 2, 3, 4, 5, 6). The air layer is held in place by the texture, with the upper edges of the micro- or nano-structures making contact with the water. Since air is approximately 50 times less viscous than water, the plastron has often been approximated as a shear-free surface in analytical models (7, 8, 9, 10, 11, 12, 13, 14, 15), leading to the expectation that SHSs could achieve very large drag reduction. Potential applications include high-Reynolds-number, turbulent flows (e.g. 16, 17, 18, 19, 20, 21), as well as low-Reynolds-number, internal flows, which are the focus of the present paper (e.g. 9, 22, 23, 11, 24, 25, 26, 27, 28). At low Reynolds numbers, the use of SHSs has been proposed to reduce what are otherwise very large pressure differences across microchannels, as is the case in microfluidic devices or in micro-cooling applications (29, 30, 31), as well as to minimize Taylor dispersion in the chemical or biological analysis of species (10).

However, laminar-flow experimental results have been mixed. While early works reported large drag reduction (e.g. 22, 32, 33), several more recent studies found no benefits, even though a plastron was clearly retained on the surface (24, 25, 34). (5) reviewed possible sources of experimental errors that might have affected some of the early measurements.

A key step towards solving this puzzle has come with the realization that surfactants could induce Marangoni stresses that impair drag reduction. More specifically, (24) experimentally examined flow over an SHS consisting of gratings perpendicular to the flow, for which they found no measurable slip at the surface. (25) also found negligible slip for an SHS consisting of gratings aligned with the flow, in contradiction with traditional theoretical and numerical results. (25) hypothesized that surfactant effects could be to blame. Following this hypothesis, surfactants naturally present in water would adsorb onto the air–water interface, as sketched in figure 1($a$). They would then be advected by the flow and accumulate at downstream stagnation points, where the interface terminates in a three-phase contact line. The resulting surfactant gradient would therefore produce a Marangoni stress opposing the fluid motion, thereby decreasing slip and increasing drag (figure 1$b$). Since traditional models of SHSs are surfactant-free, they cannot account for this additional surfactant-induced Marangoni drag.

Motivated by this hypothesis, (26) performed detailed measurements of the interface slip on an SHS comprising posts. They reported slip velocities far smaller than predicted by surfactant-free simulations. The slip pattern also exhibited strong anisotropy, consistently with what may be expected from surfactant-induced Marangoni stresses in their geometry. Deliberately adding surfactant resulted in a further small decrease in slip, although the magnitude of this change was within experimental uncertainty. (27) performed unsteady microchannel experiments over SHS consisting of gratings aligned with the flow. By introducing unsteady forcing, they uncovered complex interfacial responses that could only be explained by surfactant effects. They found that, if the driving pressure difference across the microchannel is suddenly removed, the plastron starts flowing backwards relative to the initial flow due to a surfactant-induced Marangoni force. The reverse flow decays as the inverse power of time, consistently with a similarity solution that assumes advection-dominated surfactant transport at the interface.

Since numerous works (24, 25, 26, 27, 28) observed drastically reduced slip even in nominally clean conditions, (27) performed steady simulations inclusive of surfactants, where they could precisely control surfactant concentrations. They found that surfactant effects can impair drag reduction even at extremely low surfactant concentrations, well below values naturally occurring in the laboratory or the environment. They also found that increasing the streamwise distance between stagnation points on the SHS helped to reduce the surfactant gradient and to increase slip. This explained the large slip achieved in the previous experiments of (32), who used a circular rheometer with annular gratings. Annular gratings are effectively infinitely long, without any stagnation point for surfactants to accumulate, thus avoiding Marangoni stresses. To illustrate this sensitivity of surfactant-induced Marangoni stresses with respect to the interface geometry, (35) devised an experiment whereby a complex maze is solved by a small amount of surfactant, which is introduced at the maze entrance.

More recently, (28) performed detailed experiments on an SHS consisting of a rectangular cavity with small streamwise length. They found that the rectangular gas–liquid interface exhibits recirculation, with reverse flow developing either along the middle or the sides of the plastron, depending on whether the gas–liquid interface is deformed towards the liquid phase (convex) or towards the gas phase (concave), respectively. They performed simulations where a uniform stress was applied to the interface (to approximate a Marangoni stress), showing that the experimentally-observed recirculation pattern could be induced by surfactants.

While the importance of surfactants is an emerging topic in the context of superhydrophobic surfaces, it should be noted that the importance of surfactant effects has been well-established in many other interfacial flows, often after protracted scientific debates that sought explanations for surprising phenomena. Well-known examples can be found for small bubbles rising in water, where the increased drag due to surfactant adsorption has been studied extensively (see e.g. 36, 37, 38, 39, and references therein), as well as in dip-coating problems, where accounting for Marangoni stresses is important to predict the coating thickness (40). In the ocean, the impact of naturally-occurring surfactants is well-established, as they have important effects on wave breaking and gas fluxes (41). Furthermore, steady motions in the bulk (such as internal waves or Langmuir circulations) can cause accumulation of surfactants at the surface. The resulting change in the amplitude of capillary waves affects light scattering, as revealed by satellite photographs (42). In laboratory models of oceanic flows, surfactant accumulation can be disproportionately important, driving stresses that qualitatively change the interior flow (43). Traces of surfactant have also been shown to modify drastically the behavior of the air–water interface of small bubbles probed with atomic force microscopy (see 44, 45). While a free-slip boundary condition would have been expected, force measurements demonstrated a cross-over between free-slip and no slip depending on the approaching speed of the cantilever or its probing frequency. These modified hydrodynamic boundary conditions are well-modelled by theories that include traces of surfactant, at levels undetectable through traditional surface tension measurements (44, 45). These findings further support the notion that surfactant traces can qualitatively alter the hydrodynamics.

Predicting surfactant effects is also important since surface-active molecules are inevitably present in both natural end engineered applications. Indeed, biological or environmental samples have been found to contain large amounts of surface-active compounds, including water from seas, rivers, estuaries and fog (42, 46, 47). For engineered systems, recently (48) used experiments involving insoluble liquid drops in water to demonstrate that uncrosslinked chains of polydimethylsiloxane (PDMS) can act as a surfactant. Since PDMS is one of the most common materials for microchannel fabrication, their results imply that surfactants are commonly present in microfluidic systems.

While this mounting evidence shows the importance of surfactant effects to superhydrophobic surfaces (at least in the laminar regime), there are presently no theoretical models that can predict slip as a function of surfactant properties and flow geometry.

In this paper, we build a scaling theory that describes slip length, and the associated Marangoni shear stress, in surfactant-laden laminar flows over SHS. As noted earlier, these surfactant effects are induced by accumulation of surfactant at stagnation points on the plastron, which are unavoidable in real applications (except in annular gratings in a rotating flow). As a fundamental model of such a flow, we consider a two-dimensional SHS consisting of transverse grooves, such as those considered by (24), (25) and (15). This case also serves as an upper bound for the slip and drag reduction that will be obtained in a three-dimensional flow over finite rectangular gratings. Furthermore, the model developed here constitutes a stepping stone towards a more complex theory for three-dimensional flow over SHSs with surfactant.

The problem definition and governing equations are described in §2. In §3, we present the key assumptions which allow us to develop a low-order scaling model for the local slip length at the plastron, as a function of the relevant dimensionless numbers. In §4, a model for the interior flow in a microchannel with a superhydrophobic side is developed, and coupled to the slip-length model to obtain the effective slip length and drag reduction for the overall channel flow. The overall theory is tested against numerical simulations of the full governing equations. The computational setup is described in §5, and results are reported in §6. Each parameter is varied over several orders of magnitude, confirming each aspect of the theory. The performance, key assumptions and potential uses of the theory are discussed in §7, with conclusions presented in §8. To ease adoption and testing of our model, MATLAB codes that automate the theoretical calculations are included as online supplementary material (49).

2 Problem description and governing equations

We study a steady, laminar, two-dimensional liquid flow with a small concentration of surfactant in a channel with a periodic array of flat gas–liquid interfaces on one side, as illustrated in figure 2($a$). This geometry is typical of microchannel experiments, where the smooth side of the channel is made transparent to ensure optical access (see e.g. 23, 25, 26, 27, 28). We use hats to denote dimensional quantities throughout the paper, whilst dimensionless quantities are without hats. The dimensional velocity field is $\hat{\boldsymbol{u}}(\hat{x},\hat{y})=(\hat{u}(\hat{x},\hat{y}),\hat{v}(\hat{x},\hat{y}))$. The surfactant bulk and interfacial concentration fields are $\hat{c}(\hat{x},\hat{y})$ and $\hat{\Gamma}(\hat{x})$, respectively. Owing to the periodicity of the geometry, we can restrict our study to a single periodic cell of total length $\hat{L}$ and height $2\hat{h}$, as shown in figure 2(a). This cell has a centred gas–liquid interface (hereafter designated as “the interface”) of length $\hat{g}$ at $\hat{y}=-\hat{h}$, with solid surfaces on either side of the interface. The solid surfaces have overall combined length $\hat{l}=\hat{L}-\hat{g}$. Opposite to the interface is a solid surface, located at $\hat{y}=\hat{h}$. The flow is driven in the positive $\hat{x}$ direction by a constant streamwise mean pressure drop, per unit length, $\hat{G}=-\Delta\hat{p}/\hat{L}>0$.

We deliberately choose to study the transverse flow over SHS gratings of arbitrary but finite length $\hat{g}$, instead of the longitudinal flow over infinitely long gratings as has been done in many previous theoretical and numerical studies (7, 8, 9, 12, 13, 16, 29, 15, 21). Indeed, as mentioned in §1, the establishment of adverse surfactant-induced Marangoni stresses requires stagnation points, as is necessarily found at the end of real SHS gratings, except in the special case of annular gratings (32). One of the aim of our model is to predict the effect of $\hat{g}$ on the effective slip length, following the observations made by (27) that increasing $\hat{g}$ reduces surfactant-induced Marangoni stresses. As also noted earlier, the two-dimensional flow studied here will yield an upper bound for the slip and drag reduction that can be expected in a three-dimensional flow over finite rectangular gratings.

The governing steady conservation equations for mass, momentum, bulk surfactant, and interfacial surfactant can be found in dimensional form in (27). We non-dimensionalize them using the channel half-height $\hat{h}$ as the length scale, the mean pressure drop per unit length $\hat{G}$ as the scale for pressure gradients, the corresponding velocity $\hat{U}=\hat{G}\hat{h}^{2}/\hat{\mu}$ as the velocity scale (with $\hat{\mu}$ the dynamic viscosity), the background bulk surfactant concentration $\hat{c}_{0}$ as the bulk concentration scale, and the maximum packing concentration of the surfactant at the interface $\hat{\Gamma}_{m}$ (50) as the interfacial concentration scale, such that

[TABLE]

The governing equations are, in dimensionless form,

[TABLE]

where bold symbols are used for vectors, $u_{I}(x)$ designates the velocity at the interface, and $p(x,y)$ is the bulk pressure. The subscript $I$ designates the limit of the bulk quantity considered, as it approaches the interface. In general, this limit is equal to the value taken by the quantity at the interface, except where mentioned explicitly. For instance, we have $u_{I}(x)=\lim_{y\to-1^{+}}u(x,y)=u(x,y=-1)$ for $\left|x\right|<g/2$. The Reynolds number $Re$, and bulk and interfacial Péclet numbers $Pe,Pe_{I}$ are defined below after (17), together with all other dimensionless groups in the problem. A summary is also provided in table 1.

We assume that the source–sink term modelling the flux of surfactants between the bulk and the interface follows kinetics consistent with the Frumkin isotherm, which has been found to model accurately single-component surfactant systems (51, 50),

[TABLE]

with $c_{I}(x)=\lim_{y\to-1^{+}}c(x,y)$ for $\left|x\right|<g/2$. Here $A$ is the Frumkin interaction parameter, which takes negative values $A<0$ for surfactants with attractive intermolecular interactions and positive values $A>0$ in the case of repulsive interactions. This sign convention for $A$ coincides with the one adopted by (50), but the opposite convention can also be found elsewhere in the literature (e.g. in 51). In (6), we note that the bulk concentration near the interface $c_{I}$ is different from the interfacial concentration $\Gamma$. This follows the subsurface layer model, where adsorption and desorption kinetics occur between a bulk subsurface layer and the interface (51). We note that $S>0$ corresponds to an adsorption flux and $S<0$ to a desorption flux, see figure 2(b). By definition, $\boldsymbol{u}$ and $c$ are periodic (with period $L$), while the pressure $p$ has a normalized mean drop per unit of length of $-1$, which is enforced by imposing a net pressure drop of value $L$ across each periodic unit of length $L$, so that

[TABLE]

The boundary conditions include

[TABLE]

Additionally, the continuity of the surfactant fluxes between the bulk and the interface is given by

[TABLE]

The last piece of the model is the balance of forces between the viscous drag from the bulk flow and the surfactant Marangoni force at the interface. The decrease in surface tension $\sigma$ induced by the surfactant is given by an equation of state consistent with the Frumkin isotherm (50), that is

[TABLE]

and the Marangoni shear stress at the interface is given by the gradient of surface tension, yielding the last boundary condition

[TABLE]

The chosen characteristic scales imply the following definitions for the dimensionless groups. The Reynolds number is $Re=\hat{\rho}\hat{h}\hat{U}/\hat{\mu}$, with $\hat{\rho}$ the liquid density. The bulk and interface Péclet numbers are $Pe=\hat{h}\hat{U}/\hat{D}$ and $Pe_{I}=\hat{h}\hat{U}/\hat{D}_{I}$, where $\hat{D}$ and $\hat{D}_{I}$ are the bulk and interface surfactant diffusivities, respectively. The Biot number is $Bi=\hat{\kappa}_{d}\hat{h}/\hat{U}$. The effective bulk concentration is $k=\hat{\kappa}_{a}\hat{c}_{0}/\hat{\kappa}_{d}$, where $\hat{\kappa}_{a}$ and $\hat{\kappa}_{d}$ are the adsorption and desorption coefficients, respectively. Note that, consistently with the canonical definition of Frumkin kinetics, the adsorption and desorption coefficients $\hat{\kappa}_{a}$ and $\hat{\kappa}_{d}$ have different units, so that $k$ is non-dimensional. The surfactant adsorption–desorption kinetics are parameterized by $\chi=\hat{\kappa}_{d}\hat{h}/(\hat{\kappa}_{a}\hat{\Gamma}_{m})$. We note that $\chi k=\hat{c}_{0}\hat{h}/\hat{\Gamma}_{m}$ in (14) is effectively the non-dimensional ratio between the characteristic bulk and interfacial concentration scales. The Marangoni number is $Ma=n_{\sigma}\hat{R}\hat{T}\hat{\Gamma}_{m}/(\hat{\mu}\hat{U})$, where $n_{\sigma}$ is a parameter associated with the Frumkin isotherm (51), $\hat{R}$ is the universal gas constant, and $\hat{T}$ is the absolute temperature. Temperature-driven Marangoni effects are not considered in this study and we assume that temperature is uniform in the domain. Note also that the capillary number $Ca=\hat{\mu}\hat{U}/\hat{\sigma}_{0}$ (where $\hat{\sigma}_{0}$ is the surface tension of a clean interface) has no effect in our model, since it does not appear in the final form of the Marangoni boundary condition (17) and we do not consider any other physical mechanism, such as interface curvature, in which it could play a role (§7.5.3 provides a discussion of this assumption).

The governing equations (2–5) with the periodicity and boundary conditions (7–14 and 17) define a complex nonlinear coupled problem where the unknowns are the two-dimensional velocity field $\boldsymbol{u}$, the pressure $p$, the bulk concentration $c$ and the interfacial concentration $\Gamma$. This transport problem depends on nine non-dimensional numbers, which collectively depend on a combination of flow, liquid and surfactant properties, as well as geometry, namely $Re$, $Pe$, $Pe_{I}$, $Bi$, $k$, $\chi$, $Ma$, $g=\hat{g}/\hat{h}$ and $\phi=g/L=\hat{g}/\hat{L}$. Here $g$ is the normalized interface length, whereas $\phi$ is the gas fraction. According to (17), a surfactant-induced Marangoni shear can develop at the interface when a gradient of interfacial surfactant concentration forms.

The main goal of this study is to determine a low-order model for the interfacial Marangoni shear rate $\left.\partial u/\partial y\right|_{I}\geq 0$ and the interfacial velocity $u_{I}\geq 0$ as a function of the nine non-dimensional numbers above, considering realistic parameter regimes. Such model could be used, for instance, to parameterise a slip-length condition in direct numerical simulations of flow over SHS, without having to solve the full complex coupled problem above.

3 Scaling theory for slip length with surfactant traces

3.1 Introducing the Marangoni concentration $k^{*}\equiv k\,Ma$ for small concentrations

The key assumption we propose is that the normalised interfacial surfactant concentration $\Gamma$ is sufficiently small, such that (6) and (17) can be linearised. The same assumption was made by (52) for the study of air bubbles rising in contaminated water. Hence, we obtain kinetics congruent with the Henry isotherm (51), namely

[TABLE]

In many realistic situations where surfactants are not artificially added, we indeed expect to have low effective bulk concentrations, i.e. $k\ll 1$, which generally lead to small interfacial concentrations $\Gamma$. The interfacial concentration is usually away from saturation, i.e. $\Gamma\ll 1$, because the maximum packing concentration $\hat{\Gamma}_{m}$ used in surfactant models is in fact based on geometrical arguments (53) or on achieving good fit of experimental data based on a specific kinetic model (51). Hence, $\hat{\Gamma}_{m}$ is usually not attained even when the bulk concentration is at the critical micellar concentration. We also have $A\lesssim 1$ for common surfactants (see 50). We discuss further the relevance of our assumption $\Gamma\ll 1$ in the context of applications in §7.

We take advantage of the linearisation of (6) and (17) to propose a parameter reduction in our problem, by introducing the following rescaled effective Marangoni concentration and surface concentration

[TABLE]

Substituting $k=k^{*}/Ma$ and $\Gamma=\Gamma^{*}/Ma$ into (5), (14) and (13), with the Henry kinetics (18) and (19), we obtain a set of equations where $k$ and $Ma$ have been combined to form $k^{*}$, thereby reducing by one the number of dimensionless groups. This can be verified by examining the updated version of the complete set of equations (2)-(13), which becomes

[TABLE]

with boundary conditions

[TABLE]

such that the three quantities $k$, $\Gamma$ and $Ma$ have been replaced by the dimensionless number $k^{*}$ and the variable $\Gamma^{*}$.

3.2 Scaling theory for surfactant dynamics

To make further progress in modelling the shear rate $\left.\partial u/\partial y\right|_{I}$ and velocity $u_{I}$, we perform a scale analysis on the equations in our problem, starting with rearranging (29), which expresses continuity of surfactant fluxes between the bulk and the interface

[TABLE]

For steady flows, adsorption and desorption fluxes between the bulk and the interface are in balance overall, implying

[TABLE]

such that, by the mean value theorem, there is a point on the interface where $\left.{\partial c}/{\partial y}\right|_{I}=0$ and $k^{*}c_{I}=\Gamma^{*}$. With a flow in the positive $x$-direction, interfacial surfactant $\Gamma^{*}$ is advected downstream, such that the beginning of the interface has a lower surfactant concentration, implying that $\Gamma^{*}<k^{*}c_{I}$, and that an adsorption flux exists from the bulk onto the beginning of the interface, such that $\left.\partial c/\partial y\right|_{I}>0$ there, as illustrated in figure 2(b). By the same argument, near the end of the interface, a higher surfactant concentration $\Gamma^{*}>k^{*}c_{I}$ leads to desorption from the interface into the bulk, implying $\left.\partial c/\partial y\right|_{I}<0$. Therefore, somewhere along the interface, we must have $\left.\partial c/\partial y\right|_{I}=0$. We designate by $x_{0}$ this location where the kinetics flux $S=0$, as depicted in figure 2(b).

In addition, at the beginning of the interface, $c_{I}$ is less than the bulk concentration, i.e. $c_{I}<1$ with our nondimensionalization, whereas towards the end of the interface, where surfactants accumulate, $c_{I}>1$. This means that, at a specific location along the interface, the concentration near the interface is equal to the background bulk concentration, that is $c_{I}=1$. Taking $c_{I}\sim 1$ along the interface, we then find that (31) implies that the interfacial concentration scales as $\Gamma^{*}\sim k^{*}$.

Next, assuming that the variations of $c_{I}$ and $\Gamma$ scale in the same way for the adsorption region, $-g/2<x<x_{0}$, and the desorption region, $x_{0}<x<g/2$, we have

[TABLE]

for the adsorption ($-$) and desorption ($+$) regions, respectively (see figure 2b). The quantities $\Delta\Gamma^{*}$ and $\Delta c_{I}$ are the characteristic variations of $\Gamma^{*}$ and $c_{I}$, respectively. We must have $\Delta\Gamma^{*}>k^{*}\Delta c_{I}>0$ to satisfy the direction of the kinetics flux, as described above.

From the relation between Marangoni stress and surfactant gradient (19), we also have

[TABLE]

where

[TABLE]

is the average shear rate induced by Marangoni stresses along the interface, such that $\gamma_{Ma}=0$ corresponds to free-slip at the interface and $\gamma_{Ma}=1$ corresponds to a no-slip interface. Then, a scale analysis of (31) gives

[TABLE]

where $\delta$ is the typical thickness of the diffusive layer of bulk surfactant. To estimate $\delta$, we can use the bulk advection–diffusion equation (4). At high Péclet numbers, $Pe\gg 1$, the diffusive layer of surfactant forms a thin boundary layer. As explained in detail in appendix B, there are two main asymptotic regimes depending on whether there is slip or not at the interface. For large slip and small interfacial shear rate, $\gamma_{Ma}\ll 1$, we can show that the boundary layer thickness scales as (see appendix B)

[TABLE]

where $\delta_{0,1}$, $\delta_{1,1}$, $\delta_{0,2}$, $\delta_{1,2}$ are empirical parameters which need to be determined. We note that the scaling $\delta\sim Pe^{-1/2}$ at large Péclet numbers corresponds to having a uniform velocity in the diffusive boundary layer, consistently with the case $\gamma_{Ma}\ll 1$.

For negligible slip at the interface and $\gamma_{Ma}\sim 1$, we obtain

[TABLE]

for any $g>0$, and with $\delta_{0,3}$ and $\delta_{1,3}$ two empirical parameters which need to be determined. This corresponds to the Lévêque regime (54, 55), giving a power law $\delta\sim Pe^{-1/3}$ at large Péclet numbers owing to a linear shear rate profile in the diffusive boundary layer. The scalings (37)–(39) assume that: (i) the variation of the bulk concentration along the interface is sufficiently smooth; (ii) the boundary layer is not confined vertically, i.e. $\delta\lesssim 1$; and (iii) the diffusive boundary layers between consecutive interfaces are independent. As we will discuss in §6, our scaling prediction remains accurate even for confined diffusive boundary layers $\delta\sim 1$.

With $\delta$ assumed known in terms of $g$ and $Pe$, we rearrange (36) to solve for $\Delta c_{I}$

[TABLE]

such that, dividing by $\delta$, we obtain a scaling that relates the kinetics flux to the shear

[TABLE]

3.3 Scaling for the interfacial velocity and for the slip length

We now seek a scaling expression for $u_{I}$. We integrate the interfacial advection–diffusion equation (24) from the upstream stagnation point $x=-g/2$ to $x_{0}$. We find

[TABLE]

where we used the no-slip boundary condition (25) at $x=-g/2$ for the left hand side, the no flux boundary condition (28) at $x=-g/2$ for the first term on the right hand side, as well as the continuity of flux condition (29) for the last term. To write the right-hand side in terms of $\gamma_{Ma}$, note that $\left.\Gamma^{*}\right|_{x_{0}}\sim k^{*}$ and $\left.\mathrm{d}\Gamma^{*}/\mathrm{d}x\right|_{x_{0}}\sim\gamma_{Ma}$. For the last term, we use (41) to scale the integral

[TABLE]

Substituting into (42) we obtain a scaling relation between interfacial velocity and shear. Introducing empirical prefactors (to be determined) ahead of each term, we write

[TABLE]

where $a_{1},a_{2}$ are empirical parameters; the choice of writing $2/a_{1}$ for the overall prefactor leads to a more convenient expression for the results later in §6.1.

As also noted in the previous section, scaling expressions for the boundary layer thickness $\delta$ are given by (37), (38) or (39), which depend on $g$, $\phi$ and $\gamma_{Ma}$. Therefore, our scaling is a nonlinear function of the Marangoni shear rate. However, in the comparison of our model with numerical simulations (see §6), we find that the nonlinear dependence of $\gamma_{Ma}$ with $\delta$ is actually weak. Consequently, we adopt only (39) in our model. This proves to be a good approximation and allows us to regard $\delta$ as independent from $\gamma_{Ma}$.

Furthermore, we note that a first-order linear expansion of the concentrations $c$ and $\Gamma$ near $x_{0}$ predicts $a_{2}=1/8$ since $x_{0}=0$ (that is, $x_{0}$ is at the mid-gap location) due to the balance of desorption and adsorption fluxes along the interface.

A characteristic scale for the slip length near $x_{0}$, which corresponds to the mid-gap of the interface under our assumptions, is therefore simply

[TABLE]

This scaling prediction shows that the local slip length $\lambda_{x_{0}}$ depends strongly on the Marangoni concentration $k^{*}=k\,Ma$ and the normalised gap length $g$. It is intuitive that increasing the gap length tends to increase the slip length, since it would reduce the concentration gradient at the interface and thus the opposing Marangoni stress. In contrast, increasing the effective bulk surfactant concentration $k$ or the Marangoni number tends to reduce the slip length, as expected. We also find increasing the bulk or interfacial Péclet numbers, $Pe$ or $Pe_{I}$, reduces $\lambda_{x_{0}}$. Increasing the Biot or $\chi$ numbers has a positive effect on the slip length.

However, we note that (45) is only a local measure of the characteristic slip length near the middle of the interface, where $S(x_{0})=0$. In order to have an effective or global slip length over the entire SHS which takes into account all interfaces and solid ridges, we also need to model the channel flow over the SHS. In the next section, we analyse the remaining governing equations for the flow, i.e. the continuity and Navier–Stokes equations (21) and (22), to study how the flow is affected by a SHS with a surfactant-induced Marangoni stress over the interfaces.

4 Complete model for effective slip in channel flows with one-sided periodic transverse ridges

4.1 Stokes flow model for SHS channels with surfactant contamination

According to equation (19), interfacial surfactant concentration gradients can generate a Marangoni shear rate at the interface $\left.\partial u/\partial y\right|_{I}\approx\textrm{d}\Gamma^{*}/\textrm{d}x\geq 0$. In this section, we derive an expression for how interfacial stresses with arbitrary profile can affect the flow over a periodic SHS. The geometry follows the same schematic presented in figure 2. Such a periodic SHS arrangement was studied in detail by (9) and (13) for a shear-free interface, i.e. with $\left.\partial u/\partial y\right|_{I}=0$ along the interface, at low Reynolds number. Here we generalize their approach to also study the case where $\left.\partial u/\partial y\right|_{I}\geq 0$. We also assume $Re\ll 1$, such that (22) simplifies to the Stokes flow equation

[TABLE]

Taking the curl of (46) and using the continuity equation (21), we find that the pressure field $p$ and the vorticity field, $\boldsymbol{\omega}=\nabla\times\boldsymbol{u}$ are both solutions of Laplace’s equation. Using the superposition principle to solve Laplace’s equation for the vorticity, we decompose it as the sum of the two-dimensional Poiseuille flow component, which is a pressure driven flow in a channel with full solid walls on both sides (denoted by a subscript $p$), and a deviating component (denoted by a subscript $d$), such that

[TABLE]

where $\omega_{p}=y$. As the flow is incompressible, we can also use the streamfunction $\boldsymbol{\Psi}$, defined such that $\boldsymbol{u}=\nabla\times\boldsymbol{\Psi}$, and which is the solution of the biharmonic equation $\nabla^{4}\boldsymbol{\Psi}=\boldsymbol{0}$. Note that $\boldsymbol{\Psi}=(0,0,\Psi)$ for two-dimensional flows. The solution for the deviating component of the vorticity is obtained using separation of variables considering the periodicity of the flow with wavelength $L$. Integrating twice, we then obtain the deviation streamfunction (9, 13). Noting that the mean pressure gradient imposed by the deviating field is zero, neglecting the constant of integration, and using the no-flow boundary condition $v=0$ in (25) and (26), the deviating component of the streamfunction is

[TABLE]

where $k_{n}=2\pi n/L$, and $B$, $E$, $e_{n}$ and $d_{n}$ are unknowns to be determined using the other boundary conditions. The streamfunction for the Poiseuille component is

[TABLE]

Up to this point, (4.1) and (49) are general solutions for any arrangement and geometry of SHS in a two-dimensional Stokes flow channel: i.e. they are not limited to one-sided SHS, symmetric patterns, or a particular shear rate profile at the interface.

With our geometry, using the no-slip boundary condition on the solid wall side at $y=1$ for all $x$, we find $B=-E$ and

[TABLE]

Hence, the deviating streamfunction simplifies to

[TABLE]

To determine the unknowns $E$ and $d_{n}$ for $n\geq 1$, we can use the no-slip boundary condition on the SHS side. At $y=-1$ for $g/2<\left|x\right|<L/2$, we have the condition

[TABLE]

where

[TABLE]

We then apply the last boundary condition on the interface, where we assume that there is an arbitrary shear rate profile $\left.\partial u/\partial y\right|_{I}(x)\geq 0$. Hence, we obtain the general condition

[TABLE]

for $\left|x\right|<g/2$, $y=-1$, and with

[TABLE]

To make further progress and obtain a relationship between the interfacial shear rate and the interfacial velocity, we now assume that the interfacial shear rate is uniform along the interface: $\left.\partial u/\partial y\right|_{I}=\gamma_{Ma}$, where $0\leq\gamma_{Ma}\leq 1$ corresponds to the interface-averaged surfactant-induced Marangoni shear rate, as introduced previously in (34). This assumption is consistent with having a uniform concentration gradient, following the linearised coupling condition (19). In the context of air bubbles rising in surfactant-contaminated water, this assumption is also consistent with the ‘uniformly retarded’ regime described for instance by (39). In the context of SHS, a similar assumption was made by (56) to model viscous effects from a gas phase trapped inside the cavities of the SHS. This allowed them to decouple the flow above the interface from the flow in the cavity of the SHS. We discuss further the relevance of this assumption in applications in §7. Hence, (54) becomes

[TABLE]

If ${\gamma}_{Ma}=0$ in the equation above, the interface is stress free and the surfactant concentration gradient at the interface vanishes. The surface is completely immobilized if ${\gamma}_{Ma}=1$, and the flow follows a channel Poiseuille flow.

Following (9), we can compute an approximation of the solution by truncating the series in equations (52) and (56) at $n=N-1$, multiplying (52) and (56) by $\cos(2\pi mr)$ for $m\in[0,N-1]$ (with $r=x/L$) and integrating them for $r\in(\phi/2,1/2)$ and $r\in(0,\phi/2)$, respectively, where $\phi$ is the gas fraction. Summing together the results for each $m$ in one single equation, we finally obtain a linear system of $N$ equations for the $N$ unknown coefficients $E$ and $d_{n}$ for $n\in[1,N-1]$, which we can solve numerically. The linear system in matrix form is, for $m\in[0,N-1]$ and $n\in[0,N-1]$,

[TABLE]

with $U_{0}=E$ and $U_{n}$ = $d_{n}$. The square matrix $\boldsymbol{A}_{m,n}$ has coefficients

[TABLE]

and the vector $\boldsymbol{B}_{m}$ has coefficients

[TABLE]

Care must be taken at large $n$, where the system is not well conditioned, as pointed out by (13). We provide, as supplementary material, MATLAB routines solving the linear system (57) (49).

4.2 Interfacial slip velocity

Once all the coefficients $E$ and $d_{n}$ are computed, the non-dimensional slip velocity at the interface $u_{I}$ can be determined to machine precision, depending on the size $N$ of the matrix $\boldsymbol{A}$, such that

[TABLE]

Through its coefficients $E$ and $d_{n}$, (65) is a function of the uniform Marangoni interfacial shear rate ${\gamma}_{Ma}$ and of the two non-dimensional geometrical parameters $g$ and $\phi$. Hence, we have

[TABLE]

The function $\mathcal{G}$ is only known implicitly through the solution of the linear system (57). In practice, it would be useful to obtain an explicit analytical solution, or at least a scaling expression for $\mathcal{G}$ which can give an approximate solution to the coupled surfactant–flow transport problem in combination with (44). In the linear system (57), we can factorize all the coefficients of $\boldsymbol{B}_{m}$ by $(1-{\gamma}_{Ma})$. This means that $E$ and $d_{n}$ are proportional to $(1-\gamma_{Ma})$ for all $n\geq 1$. Thus, the velocity at the interface is such that

[TABLE]

where, again, $\mathcal{F}$ is an implicit function. Now, $\mathcal{F}$ is decoupled from the surfactant transport problem since it does not depend on $\gamma_{Ma}$. It can thus be computed to arbitrary numerical precision for each couple of geometrical non-dimensional parameters $(g,\phi)$ and for all $x$ by solving the linear system (57) in the surfactant-free case, i.e. setting $\gamma_{Ma}=0$ in (63) and (64).

Based on this observation, $\mathcal{F}(g,\phi,x)=u_{I}(x)/(2(1-\gamma_{Ma}))$ is a normalised interfacial velocity. In figure 3(a), we plot on a log–log scale this normalised interfacial velocity at the middle of the gap, $x=0$:

[TABLE]

as a function of $g$ and for different $\phi$ (shown with different colors, see legend). For $\phi\ll 1$ and $g\lesssim 1$, the normalised interfacial velocity follows a linear asymptotic trend

[TABLE]

plotted with a black dotted line in figure 3(a). We can see that for $\phi=0.99$ the interfacial velocity still follows a linear scaling, although with a higher slope than in the asymptotic limit (69), as shown by the black dot-dashed line in figure 3(a), which was computed using (120). At large gap length, $g\gg 1$, and for low gas fraction, $\phi\ll 1$, the interfacial velocity collapses on the asymptotic plateau

[TABLE]

plotted with a black dashed line in figure 3(a). More details about the behaviour of the interfacial velocity $u_{Ic}$ with $\phi$ and $g$ and the two asymptotic limits (69) and (70) can be found in appendix C. The transition observed at $g\sim 1$ from a linear trend towards a plateau is due to the importance of the opposite wall at $y=1$ through viscous effects.

We note that the behaviour of $u_{I,c}/(2(1-\gamma_{Ma}))$ is similar across all $g$ and for any $\phi$. This function goes from a linear behaviour for $g\lesssim 1$ to a plateau for $g\gtrsim 1$, and with simple asymptotics in the case $\phi\ll 1$. Most of the data in figure 3(a) follows these limiting regimes, suggesting that asymptotic results are sufficiently accurate in many applications.

This common behavior of the interfacial velocity might also suggest that the velocity field follows a closed analytical form. However, we have not been able to demonstrate this theoretically from the biharmonic equation. As far as we are aware, the case of Stokes flow in a transverse channel with mixed boundary conditions changing twice (on one or both channel sides), which is reminiscent of the longitudinal-channel work of (7), has not been shown to have a closed analytical form in the literature. It would be valuable to re-examine the present problem with conformal mapping tools similar to those used by (15, 57).

Figure 3(b) plots curves of $u_{I,c}/(2(1-\gamma_{Ma}))$ versus gas fraction $\phi$, with $g$ as a parameter. As the gas fraction $\phi$ increases towards $1$, the normalised interfacial velocity increases rapidly at any fixed $g$. In the limit $\phi\to 1$ we have

[TABLE]

which can be predicted from the velocity field with uniform boundary conditions at the top and bottom sides, that is $u(y=1)=0$ and $\textrm{d}u_{u}/\textrm{d}y(y=-1)={\gamma_{Ma}}$, respectively. The solution to the Stokes problem (46) with these uniform boundary conditions is independent of $x$:

[TABLE]

and can be used to yield the limit of $u_{Ic}\to u_{u}(y=-1)$ for $\phi\to 1$. At $g\lesssim 1$, it is also clear from figure 3(b) that $u_{Ic}/(2(1-\gamma_{Ma}))\to 1$ only for gas fraction very close to 1, i.e. in the limit $\phi\to 1$, as already observed in figure 3(a). This result confirms the range of validity of the first scaling (37) for the diffusive boundary layer thickness $\delta$. Then, we can show (see appendix C) that in the limit of large gap length, $g\gg 1$, the normalised interfacial velocity follows the asymptotic hyperbolic trend

[TABLE]

plotted with a black dashed line in figure 3(b). The asymptotic result (73) is valid for any $\phi$. This result is consistent with (70) and (71).

4.3 Predictions of the interfacial shear rate, effective slip length and drag reduction

We now have two independent expressions relating the interfacial velocity $u_{I}$ and the Marangoni shear $\gamma_{Ma}$. The scaling (44) was found based on near-interface surfactant dynamics, whereas (67) was derived from a Stokes flow solution. Eliminating the interface velocity, we deduce a scaling expression for the average Marangoni shear rate,

[TABLE]

where $a_{1}$ and $a_{2}$ are the empirical parameters that were introduced in §3.2. This predictive scaling depends only on the properties of the flow, fluid and surfactant through the non-dimensional numbers $k^{*}=k\,Ma$, $Pe_{I}$, $Bi$, $Pe$ and $\chi$, and on the two geometrical parameters $g$ and $\phi$. As noted earlier, it assumes a sufficiently small concentration of surfactant and a small Reynolds number in the flow, and the diffusive boundary layer thickness $\delta$ depends only weakly on $\gamma_{Ma}$ following (37), (38) or (39). The parameters $a_{1}$, $a_{2}$, as well as $\delta_{0,i}$ and $\delta_{1,i}$ (with $i=1$, 2 or 3 for the scaling predictions (37), (38) or (39), respectively) for $\delta$, are determined empirically by fitting to our numerical simulations in §6.

We can also compute a global effective slip length $\lambda_{e}$ as defined by (9), which corresponds to the value $\lambda_{e}$ such that an equivalent channel flow under the same pressure gradient, but with a uniform Navier slip boundary condition $u(y=-1)=\lambda_{e}\left.\partial{u}/\partial{y}\right|_{y=-1}$ replaces the mixed conditions of the SHS at the bottom boundary. We can show that the contribution of the effective slip length $\lambda_{e}$ is such that the total volume flux in the channel is the sum of the background Poiseuille volume flux, $Q_{p}=2/3$, and the volume flux of the deviating flow,

[TABLE]

where the maximum value for the deviating flux is $Q_{d}\to 2$, as $\lambda_{e}\to\infty$. The effective slip length as a function of the deviating flux is

[TABLE]

From equation (4.1), the deviating streamfunction is

[TABLE]

and substituting into (76) yields

[TABLE]

Following from the linearity of the governing equations, and of the boundary conditions, $E$ also scales linearly with $(1-\gamma_{Ma})$. Accordingly, we can find the explicit dependence of $\lambda_{e}$ with the Marangoni shear rate $\gamma_{Ma}$,

[TABLE]

where $E_{0}$ is the first coefficient of the vector $\boldsymbol{U}_{n}$ (see (57)) in the surfactant-free case, i.e. $E_{0}=E$ for $\gamma_{Ma}=0$, and $\gamma_{Ma}$ is expressed by (74). As expected, $0\leq\lambda_{e}<\infty$, since $0\leq\gamma_{Ma}\leq 1$ and $0\leq E_{0}\leq 1$.

The corresponding drag reduction due to the presence of the SHS in our pressure-driven channel flows, inclusive of surfactant, can be computed as

[TABLE]

where $C_{f}=\left<\hat{\tau}_{s}\right>/(\hat{\rho}(\hat{Q}/2)^{2})$ is the laminar friction coefficient for a pressure-driven flow through a SHS channel with surfactants and $C_{f,p}=\left<\hat{\tau}_{s}\right>_{p}/(\hat{\rho}(\hat{Q}_{p}/2)^{2})$ is the laminar friction coefficient for the equivalent Poiseuille channel flow driven with the same pressure gradient and for the same channel height. The quantities $\left<\tau_{s}\right>$ and $\left<\tau_{s}\right>_{p}$ are the surface stresses averaged along both top and bottom surfaces for an SHS channel flow and a Poiseuille channel flow with the same geometry, respectively. Since, by construction, the pressure gradient is the same for the flow in the SHS channel and the Poiseuille channel flow, we have $\left<\hat{\tau}_{s}\right>=\left<\hat{\tau}_{s}\right>_{p}$. Then, using (75) and (76) into (80) we find

[TABLE]

The maximum possible drag reduction is $DR\to 15/16$ as $\lambda_{e}\to\infty$. We can compute $\lambda_{e}$ in (81) using (79) and (74). We also provide, as supplementary materials (49), MATLAB routines computing $\lambda_{e}$, $DR$ and $\gamma_{Ma}$ for any specified flow-related, surfactant or geometrical parameters.

5 Surfactant-laden numerical simulations

To test the validity of our theoretical model and its predictions for the surfactant-induced Marangoni shear $\gamma_{Ma}$ in (74) and for the effective slip length $\lambda_{e}$ in (79), we performed 137 surfactant-laden numerical simulations of the full governing equations (2)–(13).

We varied the nine dimensionless numbers independently over several orders of magnitude to comprehensively explore the parameter space. As introduced in §2, these dimensionless groups are the Reynolds number $Re=\hat{\rho}\hat{h}\hat{U}/\hat{\mu}$, the bulk and interface Péclet numbers $Pe=\hat{h}\hat{U}/\hat{D}$ and $Pe_{I}=\hat{h}\hat{U}/\hat{D}_{I}$, the Biot number $Bi=\hat{\kappa}_{d}\hat{h}/\hat{U}$, the non-dimensional bulk concentration $k=\hat{\kappa}_{a}\hat{c}_{0}/\hat{\kappa}_{d}$, the surfactant adsorption–desorption kinetics number $\chi=\hat{\kappa}_{d}\hat{h}/(\hat{\kappa}_{a}\hat{\Gamma}_{m})$, the Marangoni number $Ma=n_{\sigma}\hat{R}\hat{T}\hat{\Gamma}_{m}/(\hat{\mu}\hat{U})$, the gas fraction $\phi=\hat{g}/\hat{L}$ and the non-dimensional interfacial length $g=\hat{g}/\hat{h}$. The Frumkin interaction parameter, used in equation (6), is kept constant at $A=-1$ for all our simulations. Since this parameter has a weak influence on the surfactant-induced Marangoni shear rate, we chose a value for $A$ corresponding to moderate attractive interactions between the adsorbed surfactant molecules. This value is close to the measured value for the common surfactant sodium dodecyl sulfate in de-ionised water: $A=-2.4$ (51, 50). The aim is also to obtain values for the empirical parameters $a_{1}$, $a_{2}$ in (74) and $\delta_{0,i}$ and $\delta_{1,i}$ (with $i=1$, 2 or 3) in the uniform shear regime.

The model described by the dimensional form of equations (2) to (13) was implemented in COMSOL Multiphysics 5.2® in two-dimensional finite-element numerical simulations. The SHS channel geometry shown in figure 2(a) was used for the simulation domain, where the range of values for the gap length $\hat{g}$, the ridge length $\hat{l}$, the channel half-height $\hat{h}$ and the streamwise mean pressure drop per unit length $\hat{G}$ are presented in Supplementary Table S1.

When designing the mesh of the domain, we were particularly careful to ensure we could capture strong possible variations of some variables near the stagnation points at the beginning and end of the interface ($x=\pm g/2$), and in the vicinity of the interface. For each simulation, the maximum size of the mesh elements at the stagnation points, on the interface, and in the bulk, is detailed in Supplementary Table S1. Across all the simulations, the maximum density of elements close to the two stagnation points of the interface is 200 per micron, while the lowest density of elements at the middle of the interface is 20 per micron.

To implement the model in COMSOL, we combine the Laminar Flow module with a Dilute Species Transport module for the transport equations in the bulk (2–4). The equation for the transport of surfactant on the interface (5) is implemented through a General Form Boundary PDE, with a source term corresponding to the Frumkin kinetics flux $S$ (6). This flux also serves to implement the condition for the continuity of the diffusive flux and the kinetics flux (14) at the interface for the Dilute Species Transport module. The non-uniform distribution of surfactants at the interface yield Marangoni forces, which modify the Laminar Flow module, as stated in (17), through a weak contribution at the interface coupled to a free-slip boundary condition, resulting in the required partial slip at the interface.

The flow is forced by a mean pressure drop per unit length, which is implemented through a Periodic Flow Condition between inlet and outlet following (9), also enforcing velocity field periodicity between inlet and outlet. A gauge for the pressure is imposed through a pressure point constraint at a corner of the domain. The initial guess for the velocity, for the stationary solver, is set to the reference Poiseuille profile $u_{p}=(1-y^{2})/2$ in the entire chamber, corresponding to the stream-function (49). Periodic boundary conditions between inlet and outlet following (8) are also imposed in the Dilute Species Transport module for the bulk surfactant concentration $c$.

To ensure the accuracy and stability of the numerical simulations, we discretize the fluid flow with quadratic elements for the velocity field and linear elements for the pressure field (Taylor-Hood elements), as well as quadratic elements for the concentration fields in the bulk and on the interface. We use the MUMPS solver of COMSOL to solve for the steady state of the system, with a relative tolerance of $10^{-5}$. All our 137 COMSOL numerical simulations were fully converged, satisfying this strict relative tolerance.

The surfactant properties correspond to the well-characterized surfactant sodium dodecyl sulfate (SDS), which are well described by Frumkin kinetics (50). The physical parameters were chosen in order to explore a large range of the key non-dimensional numbers. Variations by four to six orders of magnitude were explored, as summarized in table 1 in this section, as well as in figure 8 in appendix A. In five simulations, we explored the limit of high Reynolds number with $Re\geq 1,000$, for which the flow should physically be at or above the transition to a turbulent regime. However, we imposed the flow to remain laminar in these simulations, since we are not interested in the effect of inertial instabilities or turbulence in this study. We will return to this point in §6, when discussing results at large Reynolds numbers under laminar conditions. All other relevant physical and kinetics parameters of the 137 performed simulations are presented in Supplementary Table S1.

6 Results and model performance

6.1 Effective slip length

In figure 4, we compare our scaling predictions for the effective slip length $\lambda_{e}^{\mathrm{theory}}$ with the numerical results $\lambda_{e}^{\mathrm{data}}$. We compute $\lambda_{e}^{\mathrm{theory}}$ using (79), where $\gamma_{Ma}$ follows (74) and the coefficients $E_{0}$ are computed by solving the linear problem (57) in the surfactant-free case for each couple of geometrical parameters $(g,\phi)$. The empirical parameters $a_{1}$, $a_{2}$ in (74) and $\delta_{0,i}$, $\delta_{1,i}$, with $i=1$, 2 or 3 for $\delta$ (see equations (37)–(39)) can be determined using a least-squares fitting approach and the Trust Region Reflective algorithm, as implemented in the package optimize.least_squares of Scipy (58).

First, we determine $\delta_{1,3}$ in $\delta$ by fitting a measure of the characteristic diffusive boundary layer thickness in our numerical simulations, calculated using (41), with the scaling model given in (39). We have only used the Lévêque scaling (39) for $\delta$ in (74). In our numerical simulations, the diffusive boundary layer mostly follows the Lévêque regime, which assumes a background linear shear flow, since the slip velocity $u_{I}$ is small. Moreover, as also noted earlier, the scaling model (74) for $\gamma_{Ma}$ depends weakly on $\delta$. Hence, the choice of scaling for $\delta$, which can vary between (37), (38) or (39) depending on the geometry and the slip, does not appear to be critical. The fit gives

[TABLE]

from the minimization of the sum of the squares of the relative distance of theory from data, i.e. $(\delta^{\mathrm{theory}}-\delta^{\mathrm{data}})^{2}/(\delta^{\mathrm{data}})^{2}$. This prior independent determination of $\delta_{1,3}$ reduces the number of fitting parameters to three in (74): $a_{1}$, $a_{2}$ and $\delta_{0,3}$. This ensures a more accurate and robust fit for $\lambda_{e}^{\mathrm{theory}}$, less sensitive on the actual fitting technique used.

Then, using $\delta_{1,3}=0.0528$ in (74), we fit the effective slip length $\lambda_{e}^{\mathrm{theory}}$ given by (79) to $\lambda_{e}^{\mathrm{data}}$ computed via the deviating flux $Q_{d}$ using (76). Incidentally, computing $\lambda_{e}^{\mathrm{data}}$ using (76) gives an accurate and robust estimation of the effective slip length in our numerical simulations, as it relies solely on the integral quantity $Q_{d}=Q-Q_{p}$ (see 75). Minimising the sum of the squares of the absolute distance between $\lambda_{e}^{\mathrm{theory}}$ and $\lambda_{e}^{\mathrm{data}}$, we obtain

[TABLE]

As we can see in figure 4(a,b), the scaling predictions for $\lambda_{e}^{\mathrm{theory}}$, using the values for $a_{1}$, $a_{2}$, $\delta_{0,3}$ and $\delta_{1,3}$ stated in the previous paragraph, show an excellent agreement with $\lambda_{e}^{\mathrm{data}}$ over a very large range: $${10}^{-12}$\lesssim\lambda_{e}\lesssim 1$.

The nine data points at non-negligible Reynolds numbers, $1\leq Re\leq${10}^{5}$$ (identified with blue circles in figure 4(a,b)), also exhibit good agreement despite violating the low Reynolds number assumption made in our flow model (see §4.1). As explained previously in §5, although the full steady nonlinear Navier–Stokes equation (3) was used in the simulations, the flow remained in the laminar regime for all Reynolds numbers tested.

At large non-dimensional background concentrations, $1\leq k\leq 100$ (identified with green triangles in figure 4(b)), the scaling predictions underestimate slightly the slip length. This is due to the fact that the model assumes a low concentration of surfactant. However, the model still provides a practically useful prediction of the boundary condition at the interface, which can be effectively considered as no-slip for all our simulations with $k\geq 1$. We also find that the maximum boundary layer thickness is $\delta=1.20$, which suggests that our scaling prediction is accurate even if the diffusive boundary layer is vertically confined.

We indicate in figure 4 (as well as in figures 5 and 6) data where the interface properties are strongly nonuniform, which are labeled by vermilion squares and orange diamonds. Qualitatively similar interface non-uniformities have been studied extensively in the context of air bubbles rising in surfactant-contaminated water (e.g. 36, 37, 38, 59), where they correspond to the ‘stagnant cap regime’. In this regime, an upstream part of the interface has a negligible surfactant gradient and can be considered as shear-free ($\gamma_{Ma}\to 0$), whilst the rest of the interface downstream has a large Marangoni shear ($\gamma_{Ma}\to 1$), leading to a no-slip condition over a portion of the bubble known as the ‘stagnant cap’ (hereafter designated as SC). In the SC regime, advection of surfactant along the interface dominates relative to surfactant transport between the interface and the bulk. This makes possible highly non-uniform interfacial concentrations. Since transport between the interface and the bulk is mediated by both the diffusive boundary layer flux and the surfactant kinetics, the SC regime requires that advection along the interface must be large compared to either diffusive or kinetic fluxes (or both).

We briefly summarize here the bubble-flow analysis of (39), and translate it to SHS flow. For a bubble, the SC regime is found when the characteristic interfacial Péclet number is large, and either the adsorption–desorption kinetics flux $S$ is small, or the diffusive flux through the boundary layer is small compared with the interfacial advective flux. Denoting with a superscript ‘bubble’ the results of (39), they showed that this implies $Pe_{I}^{\text{bubble}}\gg 1$, and $\mathcal{K}_{I}=Bi^{\text{bubble}}(1+k)\ll 1$ or $\mathcal{D}_{I}=\chi^{\text{bubble}}(1+k)/(Pe^{\text{bubble}})^{1/2}\ll 1$. For a bubble, the characteristic length and velocity scales are the bubble radius and interfacial velocity in the surfactant-free case. In order to translate these canonical bubble results to SHSs, note that the bubble radius is analogous to the grating length $\hat{g}$. For the SHS, the characteristic velocity scale for these non-dimensional numbers is the mid-gap interfacial velocity in the surfactant-free case, namely $\hat{u}_{Ic}(\gamma_{Ma}=0)$, which differs from the bulk characteristic velocity, such that $\hat{u}_{Ic}(\gamma_{Ma}=0)=2\mathcal{F}_{0}\hat{U}$, according to (67). This contrasts slightly with contaminated air bubbles in water, where the characteristic interfacial velocity in the surfactant-free case scales as the far-field bulk velocity, owing to the absence of rigid no-slip walls. As shown in figure 3 and explained in detail in appendix C, we have $\mathcal{F}_{0}\sim 1$ for $g\gtrsim 1$ (as for bubbles) and $\mathcal{F}_{0}\sim g$ for $g\lesssim 1$.

Therefore, using our dimensionless group definitions of §2, and using a ‘$g$’ subscript to characterize dimensionless groups where we use the lengthscale $\hat{g}$, rather than $\hat{h}$, we have $Pe_{I}^{\text{bubble}}\mapsto\mathcal{F}_{0}gPe_{I}=\mathcal{F}_{0}Pe_{I,g}$ and $\mathcal{K}_{I}=Bi^{\text{bubble}}(1+k)\mapsto\mathcal{K}_{I,g}=Bi_{g}(1+k)/\mathcal{F}_{0}$, as well as $\mathcal{D}_{I}=\chi^{\text{bubble}}(1+k)/(Pe^{\text{bubble}})^{1/2}\mapsto\mathcal{D}_{I,g}=\chi_{g}(1+k)/(\mathcal{F}_{0}Pe_{g})^{1/2}$.

The ranges spanned by the quantities $\mathcal{F}_{0}Pe_{I,g}$, $\mathcal{K}_{I,g}$ and $\mathcal{D}_{I,g}$ are reported in table 1.

The distinction between the partial SC regime, where the SC fills only part of the interface, and the full SC regime, where the SC fills all the interface, is revealed by an inspection of the shear rate profiles along the interface (not shown here). In the partial SC regime, the shear rate increases abruptly from negligible values to $\gamma_{Ma}\sim 1$ at a particular location along the interface. In the partial SC regime, the non-dimensional numbers in our simulations range approximately: 2.5\text{\times}{10}^{3}$\leq\mathcal{F}_{0}Pe_{I,g}\leq$2.5\text{\times}{10}^{5}, 9.9\text{\times}{10}^{-4}$\leq\mathcal{K}_{I,g}\leq 0.4$ and $0.04\leq\mathcal{D}_{I,g}\leq 0.4$ (see also figure [8](#A1.F8), appendix [A](#A1), for the variations of these numbers across all our numerical simulations and for the different regimes, as well as Supplementary Table S1 for the value of each parameter for each simulation). In the full SC regime, the non-dimensional numbers range approximately: $52\leq\mathcal{F}_{0}Pe_{I,g}\leq$2.5\text{\times}{10}^{4}, 2\text{\times}{10}^{-2}$\leq\mathcal{K}_{I,g}\leq 50$ and 4.0\text{\times}{10}^{-5}$\leq\mathcal{D}_{I,g}\leq 1.3$. The interfacial Péclet number is mostly higher in the partial SC regime than in the full SC regime, which is intuitively expected. We can see in figure 4(a) that the four data points in the partial SC regime (plotted with vermilion squares) are the only data points where the scaling predictions significantly underestimates the effective slip length with \lambda_{e}^{\mathrm{theory}}\leq{5.2\text{\times}{10}^{-5}}, whereas $\lambda_{e}^{\mathrm{data}}\geq$3.5\text{\times}{10}^{-4}$$. This discrepancy is due to the strong non-uniformity of the shear rate profile in the SC regime, not taken into account by our scaling model which is based on the assumption that the shear rate is approximately uniform along the interface (see (34)). The predictions $\lambda_{e}^{\mathrm{theory}}$ in the full SC regime, plotted with orange diamonds, are in reasonable agreement with the data $\lambda_{e}^{\mathrm{data}}$. We can see that $\lambda_{e}^{\mathrm{theory}}$ underestimates slightly the data, although by less than one order of magnitude for all our results in the full SC regime, with ${0.25}\leq\lambda_{e}^{\mathrm{theory}}/\lambda_{e}^{\mathrm{data}}\leq{2.4}$.

The data plotted with black pluses in figure 4, i.e. not in the SC regime, are in a state analogous to the ‘uniformly retarded regime’ described by (39) in their study of air bubbles rising in contaminated water, where they make the case that this regime exists for $\mathcal{K}_{I,g}\sim 1$ and $\mathcal{D}_{I,g}\sim 1$. However, in our simulations we find that the interfacial shear rate is in the ‘uniform’ regime, and thus satisfies our modelling assumption, over a range of $\mathcal{K}_{I,g}$ and $\mathcal{D}_{I,g}$ that spans several orders of magnitude, implying that the vast majority of the simulations satisfy our modelling assumptions. More specifically, we find that simulations in the ‘uniformly retarded regime’ have parameters that satisfy approximately 2.8\text{\times}{10}^{-3}$\leq\mathcal{D}_{I,g}\leq$4.4\text{\times}{10}^{3} and 1.9\text{\times}{10}^{-2}$\leq\mathcal{K}_{I,g}\leq$3.2\text{\times}{10}^{3}. This is most likely due to the fact that some of our simulations are in an intermediate or transition regime between the SC regime and the uniformly retarded regime, and for which $\lambda_{e}$ still follows our scaling prediction, though perhaps with slightly more scatter, as shown by some of the black pluses in figure 4.

In figure 4(c), we show the relative error between the scaling predictions $\lambda_{e}^{\mathrm{theory}}$ and the numerical results $\lambda_{e}^{\mathrm{data}}$ for the effective slip length. The error remains relatively small across all values of the average interfacial shear rate $\gamma_{Ma}$. It is less than approximately 33% for $\gamma_{Ma}\leq 0.7$, except for the four simulations in the partial SC regime plotted with vermilion squares. The relative error is less than 1.7 for $0.7\leq\gamma_{Ma}\leq 1$.

For comparison, we also show with red crosses in figure 4(c) the prediction from a surfactant-free model, which is obtained using (79) with $\gamma_{Ma}=0$. Our model provides consistently better predictions than the one that neglects surfactant effect. In particular, the error made by neglecting surfactant effects becomes very large when the interfacial shear rate increases towards the Poiseuille value $\gamma_{p}=1$. At low shear rate, $\gamma_{Ma}\leq 0.3$ we can see that the two models have comparable (small) relative errors.

Overall, we find that our scaling model for $\lambda_{e}$ provides excellent quantitative predictions across a large range of non-dimensional numbers, beyond the strict range of validity based on our modelling assumptions. Although our model predictions can underestimate the slip length in some cases (at large concentrations, and in the stagnant cap regime), our model remains practically useful as both theory and simulation yield negligible slip in those instances.

6.2 Drag reduction

We compare the drag reduction predicted by our theory ($DR^{\mathrm{theory}}$) with the numerical results from our simulations ($DR^{\mathrm{data}}$), as shown in figure 5. The value of $DR^{\mathrm{theory}}$ is obtained from (81), where the corresponding values of $\lambda_{e}^{\mathrm{theory}}$ are shown earlier in figure 4. Similarly, $DR^{\mathrm{data}}$, is calculated using $\lambda_{e}^{\mathrm{data}}$, whose values are also shown in figure 4. Using a log–log scale, we only plot data for $DR\geq${10}^{-4}$$, which correspond to the more meaningful range for practical applications. The predictions from our scaling model are in very good agreement with the numerical results. Data at even lower drag reductions (not shown here) still exhibit a very good agreement with our theoretical prediction.

In figure 5, we also plot, using red crosses, the drag reduction computed using a surfactant-free model. This is obtained by substituting the values for the surfactant-free $\lambda_{e}^{\mathrm{theory}}$ (plotted with red crosses in figure 4c) into (81). As may be expected, the surfactant-free theory almost always incorrectly predicts a larger drag reduction, with values often more than an order of magnitude larger that the actual ones. This clearly shows that the drag reduction potential of SHSs can be significantly overestimated in conditions where surfactants are important. This is consistent with the findings of (27), who showed that, for SHSs with rectangular longitudinal gratings, surfactant effects become important at very low concentrations, similar to background levels found in the environment. As may be expected, the few surfactant-free predictions in figure 5 that show better agreement with the numerical simulations correspond to lower values of $\gamma_{Ma}$, when the surfactant-free predictions converge towards our model predictions (see figure 4c).

6.3 Interfacial shear rate

We compare in figure 6 the numerical results for the average interfacial shear rate $\gamma_{Ma}^{\mathrm{data}}$ with the theoretical predictions, $\gamma_{Ma}^{\mathrm{theory}}$ computed using (74) using the four empirical parameters optimized for $\lambda_{e}$ in §6.1: $a_{1}={2.30}$ and $a_{2}={0.319}$ for $\gamma_{Ma}$, and $\delta_{0,3}={1.68}$ and $\delta_{1,3}={0.0528}$ for $\delta$ based on (39). The numerical results for $\gamma_{Ma}^{\mathrm{data}}$ have been computed by taking the spatial average of the interfacial shear rate in the interior of the interface $-g/2\leq x\leq g/2$.

In figure 6(a), we show $(1-\gamma_{Ma})$ in a log–log plot to focus on the no-slip limit $\gamma_{Ma}\to 1$. Over the limited range shown on this graph, we find good agreement between our scaling predictions and the data for all our numerical simulations where the interfacial shear rate is found approximately uniform along the interface (see the uniform regime, plotted with black plusses). Similar to $\lambda_{e}$ shown in figure 4, we can see that the four data points in the partial SC regime (vermilion squares) with $(1-\gamma_{Ma}^{\mathrm{data}})\geq 0.4$ are the only ones where the predictions underestimate the data. As discussed earlier, this is due to the strong non-uniformity of the shear rate profile in the SC regime, in contradiction with the uniform assumption made in our model (see (34)). Nevertheless, the theory remains practically useful, as both model and simulation yield shear that is essentially indistinguishable from that of a no-slip boundary.

In figure 6(b) (which is the inset of figure 6a), we plot $(1-\gamma_{Ma}^{\mathrm{data}})$ over the full range of values tested. As the average shear rate tends to the maximum Poiseuille value, $\gamma_{Ma}^{\mathrm{data}}\to\gamma_{p}=1$ or equivalently $(1-\gamma_{Ma}^{\mathrm{data}})\to 0$, $(1-\gamma_{Ma}^{\mathrm{theory}})$ underestimates the data. The difference becomes significant for $(1-\gamma_{Ma}^{\mathrm{theory}})\lesssim${10}^{-3}$$. This is due to the singularity at the two stagnation points and the difficulty associated with resolving it numerically. The shear rate exhibits extreme variations very close to the stagnation points, whilst the shear rate remains flat in the interior of the interface with values very close to the Poiseuille shear rate. We note however that the effect of the singularity appears only in the limit $\gamma_{Ma}^{\mathrm{data}}\to\gamma_{p}=1$, at values practically equivalent to a no-slip boundary condition at the interface.

This can also be seen in figure 6(c), where we plot $\gamma_{Ma}$ directly, for $\gamma_{Ma}\geq 0.7$. The scaling predictions consistently predict a no-slip boundary condition $\gamma_{Ma}^{\mathrm{theory}}\to 1$, as $\gamma_{Ma}^{\mathrm{data}}\to 1$. This shows that the actual error between $\gamma_{Ma}^{\mathrm{theory}}$ and $\gamma_{Ma}^{\mathrm{data}}$ is actually very small in this limit, where we find the simulations at large Reynolds numbers $Re\geq 1$ (blue circles), large non-dimensional concentrations $k\geq 1$ (green triangles) and in the full SC regime (orange diamonds). Predictions at intermediate values (shown by plusses in figure 6c), for $0.7\leq\gamma_{Ma}^{\mathrm{theory}}\leq 1$, show a good agreement with $\gamma_{Ma}^{\mathrm{data}}$ although with a slight overestimation.

Therefore, our scaling model also provides reasonable predictions across the whole range of interfacial shear rate values, even though the model has been fitted for $\lambda_{e}$ and not for $\gamma_{Ma}$. An agreement is found from intermediate to large values, provided the interface is not in a partial SC regime. Our scaling model remains accurate across a broad range of non-dimensional numbers (see table 1 and figure 8, appendix A) and in the full SC regime.

7 Discussion

7.1 Verifying the validity of our main assumptions

The first key assumption in our scaling model is that the non-dimensional interfacial surfactant concentration $\Gamma$ is sufficiently small so that the adsorption–desorption kinetics flux $S$ in (6) and the coupling condition (17) between the viscous stress and the surfactant-induced Marangoni stress can be linearised (see §3.1). To test the validity of the assumption $\Gamma\ll 1$, at least a posteriori, we can note that it implies $\Gamma\sim k\ll 1$, which results from applying (18) at $S=0$ along the interface. As mentioned before, we expect that $k$ should remain low in many applications where surfactants are not artificially added. (27) estimated typical ranges of $k$, depending on whether one considers ‘weak’ or ‘strong’ surfactant. (27) calculated that, for ‘weak’ types of surfactants, the non-dimensional concentration range is {10}^{-9}$\lesssim k\lesssim${10}^{-2}, which supports our hypothesis. Note that the upper bound of this range is given at the critical micellar concentration for the bulk concentration $\hat{c}_{0}$, implying that the worst-case scenario corresponds to water that is saturated with surfactant. Only for ‘strong’ types of surfactant they indicated that the $k\ll 1$ assumption could potentially be invalid, since {10}^{-6}$\lesssim k\lesssim${10}^{3}. Strong surfactants are likely to be found only in applications where they have been artificially added. Nevertheless, the model presented here performed well even at large $k$, as seen for example in figure 4.

The second key assumption made in our scaling model is that the surfactant-induced Marangoni shear rate along the interface is approximately uniform. This is related to having a uniform concentration gradient, following the linearised coupling condition (19). From the broad range of parameters tested, see table 1 and figure 8, appendix A, we find that this assumption is invalid only in the partial stagnant cap (SC) regime, where the concentration gradient presents an abrupt increase at some point along the interface, separating the no-shear and no-slip regions. As we saw in figure 4, in the partial SC regime (see vermilion squares) our model underestimates the slip length. However, it is noteworthy that our scaling model provides reasonably accurate predictions for the full SC regime, where the no-slip region spans the whole interface. Furthermore, our scaling model remains practically useful in both the partial and full SC regimes, since it correctly predicts an essentially negligible effective slip length.

If wishing to strictly determine whether our model applies, we must therefore distinguish the parameter ranges between the full and partial SC regimes. As explained in §6.1, the SC regime exists when the Péclet number at the interface, $\mathcal{F}_{0}Pe_{I,g}$, is large and either $\mathcal{D}_{I,g}$ or $\mathcal{K}_{I,g}$ are small. From our simulations, we cannot find any clear distinction between the partial and full SC regimes based only on $\mathcal{D}_{I,g}$ or $\mathcal{K}_{I,g}$. However, we noted already that the partial SC regime was generally found at larger Péclet numbers, for $\mathcal{F}_{0}Pe_{I,g}\gtrsim${10}^{3}, whilst the full SC regime was found for $1\ll\mathcal{F}_{0}Pe_{I,g}\lesssim${10}^{4}. This is physically intuitive as increasing the external flow velocity would eventually overcome the Marangoni stress at the interface. This would lead to a compression of the finite amount of surfactant adsorbed onto the interface towards the downstream end, thereby freeing the upstream part of the interface from any shear.

Since $\mathcal{F}_{0}Pe_{I,g}\propto\mathcal{F}_{0}\hat{U}$, $\mathcal{K}_{I,g}\propto 1/(\mathcal{F}_{0}\hat{U})$ and $\mathcal{D}_{I,g}\propto 1/(\mathcal{F}_{0}\hat{U})^{1/2}$, we expect to find the partial SC regime in applications where the characteristic velocity near the interface $\mathcal{F}_{0}\hat{U}$ is large. We emphasize again that the characteristic velocity in these dimensionless numbers is the local characteristic velocity near the interface, $\mathcal{F}_{0}\hat{U}$, where the bulk velocity $\hat{U}$ is modulated by the geometrical function $0\leq\mathcal{F}_{0}\leq 1$, which scales as $\mathcal{F}_{0}\sim g$ for $g\lesssim 1$, otherwise $\mathcal{F}_{0}\sim 1$. Hence, our model is valid for applications at sufficiently low $\hat{U}$ or if $g$ is sufficiently small such that the SHS is away from the partial SC regime. Microfluidic applications, such as lab-on-a-chip systems or micro-cooling, where $\hat{U}$ is small would be typical applications for our model. For instance, we can consider a typical micro-fluidic channel with $\hat{h}=50$ microns, a flow of water with characteristic speed ranging $0.1$ to $10\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$, and SHS gratings of length $\hat{g}=$1\text{,}\mathrm{m}\mathrm{m} with gas fraction $\phi\approx 0.95$. If surfactants similar to sodium dodecyl sulfate are present at a concentration of approximately ${10}^{-3}\text{\,}\mathrm{m}\mathrm{M}$ (equivalent to traces naturally present in the water), then we obtain: $800\leq\mathcal{F}_{0}Pe_{I,g}\leq$8\text{\times}{10}^{4}, $12\leq\mathcal{K}_{I,g}\leq$120, $0.7\leq\mathcal{D}_{I,g}\leq$7, and $k=${10}^{-3}$$. This shows that for this geometry with this range of flow speeds, the SHS would be in the uniform regime, far from the stagnant cap regime, such that our model would predict accurate estimates of the impact of surfactant on the slip length, drag reduction and average Marangoni shear rate.

7.2 Comparison to experimental studies of surfactant effects

Another application of our model is to analyze experimental studies reporting degradation of the performance of SHSs where surfactant contamination could be the cause. For example, two recent studies by (27) and (28) identified surfactant as the cause for the reduced or negligible slip measured near SHSs in laminar channel flows. As we discuss in detail in appendix D, the main difficulty in applying our model to predict the reduction in slip in their experiments is the absence of information regarding the surfactant properties and concentration. This is due to the fact that the surfactants were not introduced artificially, but were present as unwanted and unknown contaminants in their experiments.

Nevertheless, we can use our model to analyze a posteriori the impact of surfactant in the studies of (27) and (28). Assuming different possible surfactant types, we find that our theoretical model can predict physically sensible concentrations $\hat{c}_{0}$, which would lead to the reduced slip measured in their experiments. As detailed in appendix D, our model predicts that for instance a ‘strong’ surfactant (see 27) would only require minute traces, far below typical environmental concentrations, to reduce the slip velocity $u_{I}$ as measured by (27) and (28). A weak surfactant, e.g. SDS, would require large $\hat{c}_{0}$ close to the critical micellar concentration, whilst an intermediate surfactant (see appendix D) would require small $\hat{c}_{0}$ at or below typical environmental conditions. Hence, our model predictions are consistent with the experimental results of (27) and (28) attributing their reduced performance to surfactant contaminant traces. Our model can also provide a rational a posteriori explanation for other experimental and field measurements that have reported poor SHS drag reduction performance, in contradiction with surfactant-free theoretical or numerical predictions. Therefore, our model could help design future SHSs to mitigate or avoid surfactant effect a priori, for instance by identifying the optimal geometry and flow conditions for a given surfactant contaminant.

7.3 Analytical limits for slip and drag

It can be instructive and useful to examine practically relevant limits where our results simplify. We discuss effects of key dimensionless groups, and derive expressions in the limits of insoluble surfactant, and of very long gratings. In the latter case, it is possible to immediately predict the drag reduction without the need for solving the full Stokes flow problem. In any other case, we recommend using the MATLAB codes provided as supplementary materials (49).

To model insoluble surfactant, consider the interfacial advection–diffusion equation (24), setting the kinetics term on the right-hand side to zero. Integrating from the upstream stagnation point $x=-g/2$ to $x_{0}$, and dividing through by $\left.\Gamma^{*}\right|_{x_{0}}$, we obtain

[TABLE]

Dividing by $\gamma_{Ma}$ we find the plastron slip length, in the insoluble limit

[TABLE]

Where we assume that $\left.\Gamma\right|_{x_{0}}\simeq\Gamma_{s}$, where $\Gamma_{s}$ is the (uniform) interfacial concentration found in static conditions, and $Ma_{\text{ins}}$ is a Marangoni number for insoluble surfactant, namely

[TABLE]

Therefore, in the insoluble case, the plastron slip length is simply inversely proportional to $Ma_{\text{ins}}$, such that $Ma_{\text{ins}}\rightarrow\infty$ yields zero slip ($\lambda_{x_{0}}\rightarrow 0$), whereas $Ma_{\text{ins}}\rightarrow 0$ allows the plastron to be free-slip ($\lambda_{x_{0}}\rightarrow\infty$), analogously to the soluble case.

In general, computing the effective slip length (or equivalently the drag) requires going thorugh the Stokes flow calculation described in §4. However, in the $g\gg 1$ limit, it is also possible to evaluate analytically the effective slip length $\lambda_{e}$, and therefore the drag reduction. We start from (79), which expresses $\lambda_{e}$ as a function of $\gamma_{Ma}$ and $E_{0}$. Note that, based on (105),

[TABLE]

yielding $\lambda_{e}$ in terms of $(\gamma_{Ma},\phi)$

[TABLE]

To calculate $\gamma_{Ma}$ in the insoluble limit, use (84) to eliminate $u_{I}$ in (67), and obtain

[TABLE]

For $\mathcal{F}_{0}=u_{Ic}/[2(1-\gamma_{Ma})]$, use the large-$g$ approximation (73), that is $\mathcal{F}_{0}\simeq 1/(4-3\phi)$. Substituting into (89) and then into (90), the effective slip length for insoluble surfactant over a long grating is found explicitly as

[TABLE]

For long gratings, analytical expressions for $\lambda_{x_{0}}$ and $\lambda_{e}$ are also possible in the case of soluble surfactant. If $g\gg 1$, we expect the diffusive boundary layer to be limited by the channel height, and therefore $\delta$ will approach a constant. From our simulations, we find $\delta\simeq 1.20$ in this limit. For the plastron slip length, if $g$ is sufficiently large, we expect the second term in (45) to be dominant, yielding

[TABLE]

To find $\lambda_{e}$, we calculate $\gamma_{Ma}$ using (74), where we again set $\delta\simeq 1.20$ and $\mathcal{F}_{0}\simeq 1/(4-3\phi)$. Without further approximation we obtain

[TABLE]

Recalling that $a_{1}=2.30$ and $a_{2}=0.319$, equations (92) and (88) together provide explicitly the effective slip length as a function of surfactant properties and geometry, without the need to solve the full Stokes flow problem. The drag reduction is then found from (81), as before.

7.4 Tentative deductions for turbulent regimes

Applications of our model in turbulent regimes might be possible if a sufficiently thick viscous sublayer exists. If the surfactant transport occurs within the viscous sublayer, where the flow is laminar, the viscous sublayer height would be the appropriate length scale instead of $\hat{h}$, and the flow velocity at the edge of the viscous sublayer would be the relevant velocity scale $\hat{U}$. In that case, the local characteristic velocity $\mathcal{F}_{0}\hat{U}$ at the interface may be sufficiently small to avoid the partial SC regime. While the resulting predictions based on our model would be at best qualitative, it is of great practical interest to explore this tentative application to turbulent flows. Here we restrict ourselves to examining the plastron slip length $\hat{\lambda}_{x_{0}}$, defined in (45) and which does not depend on whether the flow is internal (e.g. channel flow) or external (e.g. a boundary layer).

In a turbulent boundary layer, with dimensional wall shear stress $\hat{\tau}_{w}$, the canonical scales are the shear velocity $\hat{u}_{\nu}=\sqrt{\hat{\tau}_{w}/\hat{\rho}}$ and the viscous length scale $\hat{\delta}_{\nu}=\hat{\nu}/\hat{u}_{\nu}$ (60). The height of the viscous sublayer is of order $10\,\hat{\delta}_{\nu}$. At this distance from a smooth wall, the flow velocity is of order $10\,\hat{u}_{\nu}$ (60). We replace $\hat{h}$ and $\hat{U}$ in our analysis with these turbulent scales and set a representative wall shear stress $\hat{\tau}_{w}=50\,$N/m2.

In practical applications, detection of a specific surfactant type is challenging. However, dimensional surface tension $\hat{\sigma}$ has been measured for both clean ‘synthetic’ seawater (labelled ‘$\hat{\sigma}_{0}$’ below), as well as for seawater samples collected through cruises (61, and references therein). For the purpose of estimating the order of magnitude of $k$ in this example, we use the Langmuir isotherm. While this is less accurate than the Frumkin isotherm, it does not require $k\ll 1$, yet it provides a relation between surfactant and surface tension that can be analytically inverted. In a liquid at equilibrium (51)

[TABLE]

(62) find that seawater that is away from major surfactant sources (such as seasonal blooms of phytoplankton, oil seeps or wastewater treatment facilities) has $\hat{\sigma}_{0}-\hat{\sigma}\sim 10^{-4}$ N/m (see also 63). Setting $n_{\sigma}\approx 2$, $\hat{R}=8.314$ kg m${}^{2}/$(s2 K mol), $\hat{T}\approx 300$ K, $\hat{\Gamma}_{m}\approx 3.9\times 10^{-6}$ mol$/$m2 and rearranging (93) for $k$, we obtain, for low-surfactant oceanic conditions

[TABLE]

Incidentally, this example yields $k\ll 1$, consistently with our set of assumptions. Note that substantially higher $k$ values can occur in oceans and lakes. In order to set up a well-defined calculation, we consider SDS with concentrations $\hat{c}_{0}=(0.01,0.1,1)$ mM, corresponding to $k=1.79\times(10^{-3},10^{-2},10^{-1})$, which bracket the value of $k$ found in (94). We change the length of the grating $\hat{g}$ from $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ to $2.5\text{\,}\mathrm{c}\mathrm{m}$, the latter being the grating length in (17). We use (45) to calculate $\hat{\lambda}_{x_{0}}$, as shown in figure 7.

To interpret figure 7, we note that one needs the effective slip length to be comparable to the thickness of the viscous sublayer in order to achieve meaningful drag reduction in turbulent flow (3). For this substantial effective slip to be possible, one needs the plastron to have an even larger slip length, since of course the solid walls will have no-slip. (For context, recall that, in canonical surfactant-free theories and simulations, the plastron is assumed to have infinite slip length.) Since the viscous sublayer thickness is $10\,\hat{\delta_{\nu}}$, we propose that a plastron slip length of around $100\,\hat{\delta_{\nu}}$ is a tentative relevant threshold for useful drag reduction. This value is marked by a dashed line in figure 7.

Note that, at small grating lengths $\hat{g}$, the first term in the right-hand side of (45) dominates. This is independent of $\hat{g}$. At larger $\hat{g}$, the second term in (45) dominates, eventually following a scaling of $\hat{g}^{5/3}$, as shown in figure 7($b$). The slip length also increases with the inverse of $\hat{c}_{0}$. These results suggest that useful drag reduction may be possible provided the surfactant concentration is not too strong and the plastron is sufficiently long in the streamwise direction. Our conclusions are consistent with the experimental results of (17), who found strong drag reduction for gratings in laboratory experiments, indicating that traces of surfactant were not sufficient to negate drag reduction. However, our theory also indicates that a large drag increase may occur for a ship equipped with SHS, when it navigates through surfactant-rich waters, which are common in the coastal ocean, rivers, and lakes. Finally, we emphasize, once again, that these are tentative deductions, and that our model will require additional work to provide quantitative drag predictions in turbulent flow.

7.5 Relative importance of effects neglected in the present model

7.5.1 Surface rheology

It is also worth discussing some physical effects not considered in our model. For example, surface rheology could play a role in the boundary condition (17) if viscous surface stresses at the interface were comparable to viscous stresses in the bulk. The relevant dimensionless groups accounting for this balance are the Boussinesq numbers $Bo_{\mu}=\hat{\mu}_{s}/\hat{\mu}\hat{g}$ and $Bo_{\kappa}=\hat{\kappa}_{s}/\hat{\mu}\hat{g}$, with $\hat{\mu}_{s}$ and $\hat{\kappa}_{s}$ the surface shear and the surface dilatational viscosities of the surfactant-laden interface, respectively. The precise measurement of $\hat{\mu}_{s}$ and $\hat{\kappa}_{s}$ is itself a challenging problem with many open questions (64).

A recent experimental study by (65), who employed a technique of unprecedented precision, concludes that soluble surfactants can be regarded as surface shear inviscid, with values of $\hat{\mu}_{s}$ below their experimental sensitivity of ${10}^{-8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$. In our problem, we can expect a negligible effect from surface shear viscous stresses, even at the smallest practical SHS length $\hat{g}$. Indeed, assuming a worst-case scenario with $\hat{\mu}_{s}=${10}^{-8}\text{,}\mathrm{kg}\text{,}{\mathrm{s}}^{-1} in water ($\hat{\mu}\approx${10}^{-3}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{s}}^{-1}) we find $Bo_{\mu}\ll{1}$ for $\hat{g}\gg{\hat{\mu}_{s}/\hat{\mu}}=${10}^{-5}\text{,}\mathrm{m}$$, which is the case in practical applications.

Surface dilatational viscosities are even more challenging to measure, since dilatational rheology and Marangoni stresses are necessarily coupled, and therefore hard to distinguish, at an interface subject to compression or expansion ((66, 67)). However, a natural (although unverified) assumption for soluble surfactant is to assume $\hat{\kappa}_{s}\sim\hat{\mu}_{s}$ (68), leading to $Bo_{\kappa}\sim{Bo}_{\mu}$. Thus, surface dilatational viscous stresses can also be considered negligible for SHS geometries with practical gap lengths $\hat{g}\gg${10}^{-5}\text{,}\mathrm{m}$$. Note also that it is common to find an effective surface dilatational viscosity in the literature which can be much larger than $\hat{\mu}_{s}$. However, unlike the true intrinsic viscosity $\hat{\kappa}_{s}$, the effective surface dilatational viscosity actually accounts for dissipation from other non-rheological effects such as adsorption–desorption fluxes, which are already accounted for explicitly in our study.

7.5.2 Viscous stresses in the gas phase

In this model, we have neglected viscous stresses from a gas phase inside the grating compared with the other stresses, namely the driving viscous stress from the liquid phase and the surfactant-induced Marangoni stress at the interface. To assess the validity of this assumption, we compare an order-of-magnitude estimate of the characteristic gas viscous stress with order-of-magnitude estimates of the other two stresses.

Let us consider the condition of continuity of stress at a surfactant-free interface of a SHS, where a viscous gas phase fills two-dimensional rectangular gratings of depth $\hat{H}_{g}$ (see e.g. 56, 69). The gas viscous stress at the interface, normalised by the characteristic driving stress from the liquid phase, is at most of the order of $\epsilon u_{Ic}/H_{g}$, where $u_{Ic}=\hat{u}_{Ic}/\hat{U}=2\mathcal{F}_{0}(g,\phi)$ is the maximum shear-free interfacial velocity computed using (68) at $x=0$ and with $\gamma_{Ma}=0$, $\epsilon=\hat{\mu}_{g}/\hat{\mu}$ is the dynamic viscosity ratio between the gas and liquid phases, and $H_{g}=\hat{H}_{g}/\hat{h}$ is the normalised depth of the grating. From all our simulations presented in figures 4–6, we estimate that the gas viscous stress is negligible compared with the driving stress from the liquid phase, $\epsilon u_{Ic}/H_{g}\ll 1$, for all $\hat{H}_{g}\gtrsim${10}^{-5}\text{,}\mathrm{m}, except at high viscosity ratio $\epsilon\gtrsim 1$. To calculate $\epsilon$ we have assumed that the gas in the grating is air, $\hat{\mu}_{g}=$1.81\text{\times}{10}^{-5}\text{\,}\mathrm{k}\mathrm{g}\mathrm{m}^{-1}\mathrm{s}^{-1}, whilst the liquid viscosity varies over a broad range such that 1.81\text{\times}{10}^{-4}$\leq\epsilon\leq$1.81\text{\times}{10}^{5}.

Compared with the Marangoni stress measured in our simulations, which can be estimated as $k^{*}/g$ when normalised with the driving viscous stress from the liquid phase (see (19) and (33)), the gas viscous stress is also negligible in all our simulations, $(\epsilon u_{Ic}/H_{g})/(k^{*}/g)\ll 1$ for all $\hat{H}_{g}\gtrsim${10}^{-5}\text{,}\mathrm{m} and all $\epsilon$. If we assume $\hat{H}_{g}\sim${10}^{-6}\text{\,}\mathrm{m}, we find that the gas viscous stress is of the same order of magnitude as the Marangoni stress in a small number of simulations only.

We have also studied the effect of air viscosity in the experiments of (27) and (28), who measured the velocity profile near the air–water interface of SHSs made of longitudinal rectangular gratings in laminar channel flows. We find that the normalised air viscous stress at the interface is approximately $\epsilon u_{Ic}/H_{g}\approx 0.05$ in the experiments of (27) and of the order of $0.001$ to 0.01 in the experiments of (28). This ratio falls by at least an order of magnitude when using their measured (reduced) slip velocity $u_{I}$, instead of our theoretical shear-free prediction (68). Since air viscous stresses are typically several orders of magnitude smaller than the characteristic driving force from the water phase, air viscous effects alone cannot explain the negligible or reduced slip velocity measured in their experiments. As we have shown above, the presence of surfactants, as modelled in this study, provides a consistent explanation for the reduced or negligible slip they measured.

The experimental and numerical study of (26) also provides compelling evidence that viscous effects from a gas phase are generally negligible or second order effects. They report experimental local and effective slip lengths in microfluidic channel flows over a SHS made of pillars. Since their geometry differs from the rectangular gratings considered in our model, the effect of Marangoni stresses due to the presence of surfactant is more difficult to estimate. Any surfactants in their experiments are transported over a complex two-dimensional interface with multiple local stagnation points, rather than a one-dimensional interface with two clear stagnation points. Nevertheless, their comparison between experimentally measured slip lengths and the slip lengths obtained from numerical simulations is revealing. Using their notation, the local experimental slip lengths, named $b_{local,exp}$ (see their figure 2b), is approximately 7% to 93% lower (depending on the location at the interface) than the slip lengths obtained numerically, which already account for viscous effect from the gas phase (named $b_{cp}$, see their figure 4b). Moreover, they find that if viscous effects from the gas phase are neglected in the numerical simulations, the effective (global) slip length increases only slightly, from $b_{eff,th}=$4.0\text{,}\mathrm{\SIUnitSymbolMicro m} to $4.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$, compared with $b_{eff}=$1.7\text{\,}\mathrm{\SIUnitSymbolMicro m} as measured experimentally. They attribute the 58% reduction in the experimental effective slip length to ‘interface contamination’, i.e. surfactant, explaining that viscous effects from the gas phase cannot explain the discrepancy with their numerical simulations.

Based on our own simulations and the studies of (26), (27), and (28), we find that viscous effects from a gas phase inside SHS gratings can generally be neglected for most practical applications, as intuitively expected and commonly assumed in the SHS literature. Indeed, in many applications or experimental studies on SHS, the liquid and gas phases are often water and air, respectively, such that $\epsilon\approx 0.02$ is very small. Moreover, the grating depth $\hat{H}_{g}$ can often be made sufficiently large so as to minimize viscous effects from the gas phase. The criterion found based on our simulations is $\hat{H}_{g}\gtrsim${10}^{-5}\text{,}\mathrm{m}$$, which is technically feasible in many applications and often necessary in experiments to prevent collapse of the plastron during the filling of the chamber. We note that in general this criterion depends on the geometry, such as the ratios $H_{g}$ and $H_{g}/g$, and whether the flow is confined in a channel or unbounded in a semi-infinite domain. For further detail about viscous effects from the gas phase in surfactant-free SHS flow, we refer the reader to the theoretical and numerical studies of (56), (70) and (69). These studies also show that gas viscous effects are mainly important at large $\epsilon$ or small $H_{g}$, consistently with our findings. In these particular regimes, both viscous effects from the gas phase and surfactant Marangoni stresses would need to be modelled in order to assess their respective contribution on the drag reduction of the SHS.

7.5.3 Interface deformation

Another physical mechanism not considered in the present study is the effect of interface deformation. Many studies have investigated the effect of lateral or longitudinal curvature of the air–water interface of SHSs (see e.g. 69, 31, 71, and references therein). They have found positive or negative impact depending on the curvature sign (whether it points towards the liquid phase or the gas phase), geometry (transverse or longitudinal SHSs), whether the flow is bounded or unbounded (31), or the Reynolds number (71).

The deformation of the interface could be due to the gravity force, viscous forces or a pressure difference across the interface. These forces must be compared with the surface tension, which resists deformations associated with an increase in surface area of the interface, i.e. flattening the plastron in this specific problem. In general, gravity can be neglected since the smallest length scale in microfluidic applications is much smaller than the capillary length, which in this case has a typical value of $\hat{l}_{c}=\nobreak\sqrt{\hat{\sigma}_{0}/(\Delta\hat{\rho}\,\hat{a}_{g})}\approx$2.7\text{\times}{10}^{-3}\text{,}\mathrm{m}. For this estimate we have chosen the representative values of $\hat{\sigma}_{0}\approx\nobreak$7.2\text{\times}{10}^{-2}\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1} for the surface tension, $\hat{a}_{g}=$9.81\text{,}\mathrm{m}\mathrm{/}\mathrm{s}^{2} for the gravitational acceleration, and $\Delta\hat{\rho}\approx$1000\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3} for the air-water density difference. Similarly, viscous forces are neglected in most applications due to small capillary numbers near the interface $Ca_{I}=u_{I}Ca=\hat{\mu}\hat{u}_{I}/\hat{\sigma}$. We typically find $Ca_{I}\lesssim${10}^{-3}\text{,}$$ for a surfactant-free air–water interface across our range of parameters. We note the effect of surfactant would on the one hand tend to reduce $\sigma$, thus enhancing interfacial deformation. On the other, as we have shown in this study, surfactant would reduce $u_{I}$, thus limiting interfacial deformation. The reduction of $u_{I}$ can occur at concentrations much smaller than concentrations necessary to change the surface tension noticeably (51). Hence, we intuitively expect that viscous forces would have negligible effect, even when combined with surfactant, in regimes where surfactants affect $u_{I}$. The capillary length and the capillary number depend on the properties of the fluids on either side of the interface. Although air–water systems, as assumed above, are the most common across real applications, laboratory and field experiments, these characteristic numbers would need to be examined carefully in more specialised applications (e.g. liquid metals for micro-cooling, 30).

The effect of pressure difference is one of the most common cause of interfacial deformation (e.g. 71). Interfacial curvature typically depends on the ratio of the pressure difference and the surface tension, $\Delta p/\sigma$, following the Young–Laplace law. Thus, surfactant could enhance curvature by reducing $\sigma$, thereby affecting the performance of the SHS. Similar to what we noted for viscous effect, we expect the negative impact of surfactant on $u_{I}$ via Marangoni effects to be generally more important than via interface deformation. If the pressure difference is large enough, there can exist some regimes where both interface deformation and Marangoni stresses are important. The combined effects on SHS performance of negative Marangoni effects and positive or negative interfacial deformation effects would be an important topic for future research.

7.5.4 Three-dimensional effects

Although the geometry used in our model is two-dimensional, we expect the model to give a reasonable estimate of the impact of surfactants for flows above three-dimensional rectangular longitudinal SHS gratings, similar to those used by (27) and many other studies. For three-dimensional gratings with small aspect ratio $w/g=1/15$, (28) observed three-dimensional flows with recirculations along the side boundaries or via the interior, depending on whether the interface was convex or concave. Overall, they found significant reduction of the slip velocity at the interface due to surfactant contamination, which shows that these three-dimensional recirculation flows are secondary effects compared to the mean two-dimensional effects due to the surfactant-induced Marangoni forces. For cases without this recirculation pattern, we expect surfactants to be advected along the grating, forming a longitudinal surfactant gradient which is approximately uniform in the spanwise direction (i.e. across the grating width). Owing to spanwise viscous friction, we can also note that our model would give a lower bound prediction on the surfactant-induced Marangoni shear, or conversely, an upper bound for the effective slip length and maximum drag reduction.

8 Conclusions

In this study, we present a reduced-order scaling model to account for the impact of soluble surfactants in channel flows with superhydrophobic surfaces. The drag reduction potential of superhydrophobic surfaces can be severely reduced if surfactants adsorbed onto the plastron induce Marangoni forces opposed to the flow. These Marangoni forces develop when a gradient of surfactant establishes along the interface.

To simplify the governing equations of this problem, we first linearised the kinetics source terms for the surfactant flux between the bulk and the interface, as well as the coupling condition balancing the viscous force and the surfactant-induced Marangoni force. This linearisation holds for small surfactant concentration $\Gamma\ll 1$, which is a reasonable assumption for most applications where surfactants are not artificially added. Then, integrating the transport equations in the bulk and at the interface, we find a linear relationship between the interfacial slip velocity at mid-gap and the interface-averaged surfactant-induced Marangoni shear, given by (45). This relationship depends explicitly on the non-dimensional numbers $k^{*}=kMa$, which combines both the non-dimensional bulk background surfactant concentration $k$ and the Marangoni number $Ma$, as well as $Pe$, $Pe_{I}$, $g$, $Bi$ and $\chi$.

To obtain a global effective slip length and predict how surfactant transport can affect the flow rate and the drag reduction potential of the SHS, we solve the continuity and momentum conservation equations for low Reynolds number flow. Using a technique based on the work of (9) for surfactant-free SHS flow, we solve Stokes’ equation with mixed boundary conditions and a prescribed shear profile at the interface. In the case of a uniform interfacial shear $\gamma_{Ma}$, the interfacial velocity relates linearly to $1-\gamma_{Ma}$, where the coefficient of proportionality depends on the geometric non-dimensional parameters of the SHS, namely the grating length $g$ and the gas fraction $\phi$. We close the problem and eliminate the interface velocity by using our earlier result, based on the surfactant problem, that also related interface velocity to shear. Hence, we find that the average Marangoni shear $\gamma_{Ma}$ depends on seven non-dimensional parameters: $k^{*}$, $Pe$, $Pe_{I}$, $Bi$, $\chi$, $g$ and $\phi$, following (74). The dependence on the geometry is implicit through the function $\mathcal{F}_{0}(g,\phi)$, which can be solved from the linear problem (57) assuming a surfactant-free Stokes’ flow in the same geometry. We find that the effective slip length is $\lambda_{e}=2(1-\gamma_{Ma})E_{0}/(1-(1-\gamma_{Ma})E_{0})$, see (79), where $E_{0}=Q_{d,0}/2$ with $Q_{d,0}$ the added volume flow rate in an SHS channel flow without any surfactant. The corresponding added flow rate $Q_{d}$ and drag reduction $DR$ due to the SHS, in the general case of a surfactant-contaminated flow, can be determined from the effective slip length following (76) and (81), respectively. These equations show how the slip length, the added flow rate and the drag reduction are affected by the surfactant-induced Marangoni shear rate at the interface.

In order to test the regime of validity and the accuracy of our model, we performed 137 finite-element numerical simulations of the full governing equations in steady, pressure-driven, laminar channel flows, inclusive of soluble surfactants following (2)–(13). We varied the governing non-dimensional groups across a broad range of values to explore the vast parameter space of this problem (see figure 8, appendix A, table 1 and the Supplementary Table S1). The model predictions for $\lambda_{e}$, $DR$ and $\gamma_{Ma}$ follow well the numerical results across almost all the parameter space explored. The model coefficients are determined through a least-squares fit for $\lambda_{e}$, yielding $a_{1}={2.30}$, $a_{2}={0.319}$, $\delta_{0,3}={1.68}$ and $\delta_{1,3}={0.0528}$. The flows that are least well captured by our model corresponds to the ‘partial stagnant cap regime’, which is also found in air bubbles rising in surfactant-contaminated water. This regime occurs at very large $\mathcal{F}_{0}Pe_{I,g}$, and low $\mathcal{D}_{I,g}$ or low $\mathcal{K}_{I,g}$. The partial SC regime exhibits a sharp increase in the shear rate at the transition between a shear-free upstream part and a no-slip downstream part of the interface, which differs from our assumption of a uniform Marangoni shear along the interface. Nevertheless, at least for the simulations performed here, our model predictions are sufficiently accurate for practical purposes. It will be important to test the accuracy of our model also in more complex flows.

Canonical SHS models, which completely neglect surfactant effects, can yield a large error in the prediction of the slip length and of the drag reduction, as shown in figures 4 and 5. In particular, the error is very large, by several orders of magnitude, at large Marangoni stresses. Hence, models neglecting surfactant can significantly overestimate the drag reduction potential of the SHS. This is particularly important in applications where small background environmental surfactant traces are sufficient to induce strong Marangoni forces, as previously found by (27).

Overall, the model we present provides a useful quantitative estimate of the effect of surfactants on the drag reduction potential of SHSs, across a vast part of the parameter space except in the partial stagnant cap regime. Our scaling predictions can be used directly in numerical simulations of flow over SHS in realistic conditions where surfactants cannot be neglected. The effective slip length $\lambda_{e}$ can be used in a Navier-slip boundary condition on the SHS side, without having to solve the full coupled nonlinear surfactant transport problem. This will reduce considerably the computational burden associated with realistic simulations of SHS flows. We also note that our model can be easily adapted for a two-sided SHS channel, via changes in the boundary conditions in the Stokes’ flow problem (see §4.1). This change in boundary conditions will modify the geometric function $\mathcal{F}_{0}$.

Future work will investigate how the model can be modified for more complex three-dimensional flows over SHSs, such as pillars or disordered SHSs. Apart from annular flows (32, 28) or very long air–water interfaces (27), accumulation of surfactant at stagnation points in these three-dimensional problems can also lead to surfactant-induced Marangoni stresses. Predicting the magnitude of these forces and the overall effect on the effective slip length or the drag reduction is a complex problem. Many applications operate at larger Reynolds numbers, where the effect of turbulence on the surfactant Marangoni stresses may be important. At intermediate Reynolds numbers, where the viscous sub-layer forming at the SHS is sufficiently thick compared with the surfactant diffusive boundary layer thickness, our scaling model may still be applicable, though the empirical parameters may differ from those found here. At very large Reynolds numbers, turbulence is likely to enhance the diffusion of surfactant in the bulk and at the interface, which could change the concentration gradients and result in intermittent localised Marangoni forces at the interface. These problems have a direct impact on the performance of SHSs in many applications, and constitute important topics for future studies.

9 Acknowledgments

We gratefully acknowledge financial support from the Raymond and Beverly Sackler Foundation, the Engineering and Physical Sciences Research Council, the European Research Council Grant 247333, Mines ParisTech, the Schlumberger Chair Fund, the California NanoSystems Institute through a Challenge Grant, ARO MURI W911NF-17-1-0306 and ONR MURI N00014-17-1-2676.

Appendix A Key dimensionless numbers across all numerical simulations

To help provide a visual overview of the simulations performed, figure 8 plots the value of each dimensionless group on the vertical axis, with the horizontal axis indicating different simulations. Ranges for each parameters are also reported earlier in table 1. Detailed values are included in table S1 of the Supplementary Materials.

Appendix B Diffusive boundary layer thickness

To determine an estimate of the boundary layer thickness $\delta$ for the surfactant concentration, we builds on the result in §3.2 and perform a scale analysis of the bulk advection–diffusion equation (23), which is expanded below as

[TABLE]

In the surfactant adsorption (resp. desorption) boundary layer forming above the interface, we denote the characteristic variation of the bulk concentration as $\Delta c$. In the streamwise direction, we expect the change $\Delta c$ to take place between between $x=-g/2$ and $x=x_{0}$ (resp. $x_{0}$ and $g/2$), as sketched in figure 2. As explained in §3.2, $x_{0}$ is defined as the interface location where the kinetics flux $S$ vanishes. Under the assumption of low interfacial concentration (see (33) and text above), we previously found that the adsorption and desorption diffusive boundary layers are approximately anti-symmetric and of characteristic streamwise length scale $\sim g$, as depicted in figure 2(b).

If we focus on the adsorption region of the interface, $c$ at the interface is denoted as $c_{I}$, which varies by a scale $\Delta c_{I}$ between $x=-g/2$ and $x=x_{0}$, where $S=0$ implies $c_{I}\sim 1$. In addition, the characteristic cross-stream variation, across the boundary layer, is from $c_{I}$ to $1$, implying that this variation in $c$ also scales as $\Delta c_{I}$. Therefore, in both the $x-$ and $y-$ directions, $\Delta c=\Delta c_{I}$ over characteristic distances $g$ and $\delta$, respectively.

We denote the characteristic streamwise and cross-stream velocities in the diffusive boundary layer as $U_{\delta}$ and $V_{\delta}$, respectively. Hence, a scale analysis of equation (95) gives

[TABLE]

where we can divide throughout by $\Delta c_{I}$. Thus, $\delta$ is a function of the Péclet number $Pe$, the interface length $g$, as well as $U_{\delta}$ and $V_{\delta}$. These velocity scales are expected to depend on the geometrical parameters $g$ and $\phi$, as well as on the interfacial velocity $u_{I}$ and the characteristic shear rate profile $\gamma_{Ma}$. We now seek an explicit dependence of $U_{\delta}$ and $V_{\delta}$ on these parameters.

We assume that the diffusive boundary layer thickness is not affected by the channel height, such that $\delta<1$. We also assume that a diffusive boundary layer above a particular interface is independent from the other interfaces, and thus independent of the gas fraction $\phi$. We retain the dependence on the interface length $g$. We can distinguish two main limits influencing $U_{\delta}$ and $V_{\delta}$, depending on the boundary condition at the interface. This can either consist of a finite slip and negligible shear ($u_{I}>0$ and $\gamma_{Ma}\ll 1$), or of no-slip and finite shear ($u_{I}=0$ and $\gamma_{Ma}\sim 1$). Hence, in general, $U_{\delta}\sim u_{I}+\gamma_{Ma}\delta$. The two cases are analyzed further below.

• •

First, for $u_{I}>0$ and $\gamma_{Ma}\ll 1$, according to the analysis in appendix C, we have $U_{\delta}\sim u_{I}\sim g$ for $g\lesssim 1$, and $U_{\delta}\sim u_{I}\sim 1$ for $g\gtrsim 1$. We determine the scale for $V_{\delta}$ using the continuity equation (21), which gives $V_{\delta}\sim U_{\delta}\delta/g$. Replacing these velocity scales into (96), we find

[TABLE]

and

[TABLE]

where $\delta_{0,1}$, $\delta_{1,1}$, $\delta_{0,2}$, $\delta_{1,2}$ are empirical parameters.

• •

Second, for $u_{I}$ negligible and $\gamma_{Ma}\sim 1$, $U_{\delta}$ depends only on the ratio of the thickness of the diffusive boundary layer and the channel height: $U_{\delta}\sim\delta\gamma_{Ma}=\delta$ for $\delta<1$. This regime is also known as the Lévêque regime (54, 55). Note that $V_{\delta}\sim 0$ in this case. Replacing these velocity scales into (96), we find the asymptotic behaviour

[TABLE]

for any $g>0$, and with $\delta_{0,3}$ and $\delta_{1,3}$ two empirical parameters.

As noted before, the results (97)–(99) are valid provided $\delta<1$, which is satisfied for large enough Péclet numbers or small enough gap length. For an intermediate regime with partial slip and partial shear, i.e. $U_{\delta}\sim u_{I}+\gamma_{Ma}\delta$, we expect that the boundary layer thickness has an exponent between $-1/2$ and $-1/3$. The transition between the slip dominated regime, with scaling (97) or (98), and the shear dominated regime, with scaling (99), should be smooth at low Reynolds numbers.

Appendix C Asymptotic limits for the slip velocity

The computation of the slip velocity yields distinctive simplified behaviours in the limits of large and small gap length $g$, as evidenced in figure 3. In this section, we analytically derive asymptotic limits for the slip velocity profile $u_{I}(x)$, and confirm their agreement with the numerically computed values at mid-gap from figure 3.

We start by considering the so-called dual series comprised of equations (52) and (56)

[TABLE]

with $\alpha_{n}$ and $\beta_{n}$ defined in equations (53) and (55), respectively. From this set of expressions, it is possible to obtain a closed form of the asymptotic behavior of the slip velocity by considering only the leading order of $\alpha_{n}$ and $\beta_{n}$ in the relevant limits. This is done in a similar fashion to (9) and (13), who derived expressions for the effective slip length from the asymptotic behavior of the first coefficient $E$. However, the slip velocity depends on the whole set of coefficients, and in this case it is not enough to derive an expression for the first coefficient only. Indeed, recall the form of $u_{I}(x)$ from (65)

[TABLE]

The derivations of $u_{I}(x)$ and its value at mid-gap $u_{Ic}=u_{I}(x=0)$ in the limits of large and small $g$ are presented in the next two subsections.

C.1 Limit of large gap length

Consider the limit $g\rightarrow\infty$, with the gas fraction $\phi$ fixed. Since $L=g/\phi$, note that this case necessarily implies $L\rightarrow\infty$ as well. Note that, due to our choice of the channel height $\hat{h}$ as the length scale for the nondimensionalization, this limit corresponds to a “narrow” channel with the top wall close to the plastron. Consequently, we have

[TABLE]

and in this limit $\alpha_{n}$ and $\beta_{n}$ can be expanded as

[TABLE]

Taking into account that $k_{n}\sim L^{-1}$, the expressions (100) and (102) lead to the following expansions of the unknown Fourier coefficients

[TABLE]

Substituting (102) and (103) into (100) we arrive at the leading-order dual series for $E^{(0)}$ and $d_{n}^{\,(0)}$. After introducing the changes of variable $\widehat{d}^{\,\,(0)}_{n}=(2\pi{n})^{3}\,d^{\,(0)}_{n}$ and $z=2\pi{x}/L$, this dual series yields

[TABLE]

and its coefficients can be obtained exactly. Indeed, after integrating (104a) from $\phi\pi$ to $\pi$ and (104b) from [math] to $\phi\pi$, one can then sum the two expressions and obtain

[TABLE]

The rest of the coefficients $\widehat{d}^{\,\,(0)}_{n}$ can be retrieved multiplying (104) by harmonics of the form $\cos(mz)$ with $m\in\mathbb{N}$ and $m\geq 1$. Then, applying the same procedure of integration and summation and invoking orthogonality between the functions, we arrive at

[TABLE]

Using the obtained set of coefficients, one can now evaluate the slip velocity. Substituting (102a) and (103) in (101), we have

[TABLE]

Applying (105) and (106), one subsequently obtains

[TABLE]

First, note that for $\phi=0$ the above expression (107) yields $u_{I}(x)=0$ at leading order, as one would expect. We then observe that the expression in brackets in (107) is the Fourier cosine series of a square wave with value 1 for $|z-2j\pi|<\phi\pi$ and 0 for $\phi\pi<|z-2j\pi|\leq\pi$, where $j\in\mathbb{N}$. Consequently, by virtue of the uniqueness of a Fourier series one has, after undoing the change of variables

[TABLE]

The fact that the slip velocity tends to a constant value as $g\rightarrow\infty$ is expected, due to the confinement effect of the top wall. Indeed, the disparity of horizontal and vertical length scales ($g\gg{1}$) leads to a lubrication regime in which the slip velocity asymptotically tends to a constant along the plastron. In such a regime, the velocity field can be approximated as unidirectional in the central “core region” following the thin-gap approximation (see for instance the examples in 72). From (108), the value of the slip velocity at mid-gap $u_{Ic}=u_{I}(x=0)$ would then yield at leading order

[TABLE]

where $u_{Ic}$ has been normalized with $2(1-\gamma_{Ma})$ following Section 4. Notice that this normalization implicitly assumes $0\leq\gamma_{Ma}<1$, however in the case $\gamma_{Ma}=1$ it is straightforward from (108) that $u_{Ic}=0$ at leading order.

The expression (109) is plotted in figure 3($b$), confirming the trend of the values $u_{Ic}$ computed numerically. Moreover, note that within this asymptotic regime $g\rightarrow\infty$, the expression (109) leads to the two following limits

[TABLE]

which are corroborated as well by the asymptotic behavior in figure 3(a).

C.2 Limit of small gap length

Consider now the limit $g\rightarrow{0}$, with the gas fraction $\phi$ fixed. Then, like in the previous case, $g\rightarrow{0}$ necessarily implies $L\rightarrow{0}$ as well. This case corresponds to a “tall” channel with distant top walls. We have

[TABLE]

and therefore $\alpha_{n}$ and $\beta_{n}$ can be expanded as

[TABLE]

Given the functional form of the leading order terms in (111), we introduce the change of variable $\widehat{d}_{n}=e^{k_{n}}\,d_{n}$ and seek the expansions

[TABLE]

After re-scaling the spatial variable $z=2\pi{x}/L$, we insert (111) and (112) into (100) and group the $O(1)$ terms to arrive at the leading-order dual series for $E^{(0)}$ and $\widehat{d}_{n}^{\,\,(0)}$

[TABLE]

which leads to $E^{(0)}=0$ and $\widehat{d}_{n}^{\,\,(0)}=0$. The terms of order $O(L)$ can then be grouped into the following dual series for $E^{(1)}$ and $\widehat{d}_{n}^{\,\,(1)}$

[TABLE]

The coefficients in (114) can be obtained exactly following the procedure of (73). However, in this case it is not necessary to explicitly obtain $E^{(1)}$ and $\widehat{d}_{n}^{\,\,(1)}$ in order to obtain the slip velocity, since the left-hand side of equation (114a) can be determined exactly for $|z|<\phi\pi$ (see page 161 of 73)

[TABLE]

Here $h(t)$ can be retrieved from equation (5.4.60) in page 162 of (73), which in our case simplifies to

[TABLE]

where it is worth noting that the closed form of the integral

[TABLE]

has been used in the derivation above, obtained from (12).

The desired slip velocity for $0\leq z<\phi\pi$ can now be retrieved at leading order introducing the expansions (111) and (112) into (101) and using equation (115)

[TABLE]

Making use of (C.2), integrating and undoing the change of variable we obtain the velocity profile

[TABLE]

The above formula (119), after setting $\gamma_{Ma}=0$ and a change in the variables normalization, is at leading order exactly half of the slip velocity obtained by (12) for a configuration with longitudinal no-shear infinite gaps in a semi-infinite domain. This result is consistent with the analysis of (74), who conclude that the slip velocity profile in such a configuration should be larger than that of the equivalent transverse case by exactly a factor of two.

From (119), we can finally obtain the normalized slip velocity at mid-gap at leading order

[TABLE]

where we have substituted $L=g/\phi$.

From (120) we can corroborate the validity of the linear scaling $u_{I}\sim{g}$ for $g\lesssim\nobreak{1}$. Indeed, the asymptote (120) is plotted for $\phi=0.99$ in figure 3($a$), showing good agreement with the numerically computed slip velocity.

Consider now the limit $\phi\rightarrow{0}$ within the regime of small gap length $g\rightarrow{0}$ investigated in this subsection. Then (120) yields, to leading order in $g$,

[TABLE]

This is congruent with the linear asymptote for small $g$ followed by the values calculated numerically, which is also shown in figure 3($a$).

The coefficient $E$, which appears in the expressions for the effective slip length (79) and drag reduction (81), can also be obtained in this limit from (114) following (73). We make use of (C.2) to arrive at

[TABLE]

Appendix D Application of our model to experimental studies in the literature showing reduced slip

We study the experimental results of (27) and (28) to analyse with our theoretical model how surfactant affected their SHS performance. The slip velocities extrapolated from the measurements of (27) on the interface ($z=0$) at mid-gap ($y=0$) are: $u_{I}\approx$4\text{\times}{10}^{-3}\text{,}$\pm$4\text{\times}{10}^{-3}\text{,} for 2 mm long lanes (see their figure 3D), which is practically negligible; and $u_{I}\approx$5\text{\times}{10}^{-2}\text{\,}$\pm$9\text{\times}{10}^{-3}\text{\,} for 30 mm long lanes (figure 3E), which is significantly reduced compared with the theoretical (surfactant-free) prediction. (Note that we have non-dimensionalised these velocities using the characteristic velocity $U$, following the convention used in the present study.) Similarly, (28) report: $u_{I}\approx$8\text{\times}{10}^{-2}\text{,} for 5 mm long lanes (see their figure 3*b*), which is significantly reduced compared with the theoretical (surfactant-free) prediction; and $u_{I}\approx$8\text{\times}{10}^{-3}\text{\,} for 15 mm long lanes (figure 5), which is practically negligible.

The main difficulty in applying our theoretical model, for instance to predict the reduced slip velocities measured experimentally by (27) and (28), is that the surfactant properties and their concentrations are completely unknown in their experiments. Instead, we use our model to predict the concentration of surfactant, for three different possible surfactant types, which could lead to the measured $u_{I}$ reported in (27) and (28). The three surfactants we choose are: a ‘strong’ poorly soluble surfactant with properties described in (27), a ‘weak’ highly soluble surfactant, namely Sodium Dodecyl Sulfate (SDS), and an ‘intermediate’ type with similar weak properties as SDS but rendered almost insoluble in water by reducing its desorption coefficient to $\hat{\kappa}_{d}=$1\text{,}\mathrm{s}^{-1} (instead of $\hat{\kappa}_{d}=$500\text{\,}\mathrm{s}^{-1} for SDS in water). Surfactants have a large number of parameters ($\hat{\kappa}_{d}$, $\hat{\kappa}_{a}$, $\hat{D}$, $\hat{D}_{I}$, $\hat{\Gamma}_{m}$, $A$, $n_{\sigma}$), in addition to their bulk background concentration $\hat{c}_{0}$, which are almost all used in our theoretical model (see (44) and (74), which we use to compute $u_{I}$. Thus, by choosing only three different types of surfactants from the vast parameter space, the analysis in this section is primarily qualitative. The aim is to show that our theoretical model provides physically meaningful explanations regarding the impact of surfactant in experimental studies showing reduced SHS performance such as those of (27) and (28).

Assuming a strong surfactant, our model predicts that a bulk surfactant concentration $\hat{c}_{0}\sim${10}^{-13}\text{,}\mathrm{m}\mathrm{M} can reduce $u_{I}$ in the same extent and under the same conditions as reported by ([27](#bib.bib27)) for both short and long lanes; and $\hat{c}_{0}\sim${10}^{-15}\text{\,} to ${10}^{-14}\text{\,}\mathrm{m}\mathrm{M}$ for the experiments reported by (28). Assuming the weak SDS surfactant, our model predicts $\hat{c}_{0}\sim 1$ to $10\text{\,}\mathrm{m}\mathrm{M}$ (i.e. near the critical micellar concentration) to obtain the results of (27) and $\hat{c}_{0}\sim 0.1$ to $3\text{\,}\mathrm{m}\mathrm{M}$ for the experimental results of (28). Assuming an intermediate surfactant, $\hat{c}_{0}\sim${10}^{-5}\text{,} to ${10}^{-4}\text{\,}\mathrm{m}\mathrm{M}$ would lead to the results of ([27](#bib.bib27)), and $\hat{c}_{0}\sim$4\text{\times}{10}^{-7}\text{\,} to ${10}^{-5}\text{\,}\mathrm{m}\mathrm{M}$ for the results of (28). These theoretical predictions show that: (i) a very strong surfactant would require only minute traces, unavoidable in normal environmental conditions, to strongly affect slip; (ii) whilst at the other extreme, a weak surfactant such as SDS would require a concentration of the order of the critical micellar concentration to lead to a no-slip or reduced slip condition. Then, an intermediate surfactant would require small concentration at or below typical environmental background concentration to lead to no-slip or reduced slip condition.

Therefore, our theoretical model provides physically sensible predictions with regard to surfactant types and concentrations that may have contaminated the experiments of (27) and (28). This is consistent with their conclusions. We note that our theoretical model assumes a two-dimensional channel geometry with one-dimensional interfaces, whereas the experiments are three-dimensional with two-dimensional (flat) or three-dimensional (curved) interfaces bounded laterally by no-slip walls. Hence, we expect our model to over-predict the interfacial slip velocity (for a given surfactant type and concentration) or over-predict the background surfactant concentration (for a given surfactant type and interfacial slip velocity). This means that even lower surfactant concentrations could have affected the experimental results of (27) and (28).