# A single parameter can predict surfactant impairment of superhydrophobic drag reduction

Fernando Temprano-Coleto, Scott Smith, François Peaudecerf, Julien Landel, Frédéric Gibou, Paolo Luzzatto-Fegiz

► To cite this version:

Fernando Temprano-Coleto, Scott Smith, François Peaudecerf, Julien Landel, Frédéric Gibou, et al. A single parameter can predict surfactant impairment of superhydrophobic drag reduction. Proceedings of the National Academy of Sciences of the United States of America, 2023, 120 (3), pp.e2211092120. 10.1073/pnas.2211092120. hal-03943393

HAL Id: hal-03943393

https://hal.science/hal-03943393

Submitted on 17 Jan 2023

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

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

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

# A single parameter can predict surfactant impairment of superhydrophobic drag reduction

Fernando Temprano-Coletoª ID, Scott M. Smithª, François J. Peaudecerfᵇ ID, Julien R. Landelᶜ ID, Frédéric Gibouª ID, and Paolo Luzzatto-Fegizª,¹ ID

Edited by Joanna Aizenberg, Harvard University, Cambridge, MA; received June 29, 2022; accepted November 21, 2022

Recent experimental and computational investigations have shown that trace amounts of surfactants, unavoidable in practice, can critically impair the drag reduction of superhydrophobic surfaces (SHSs), by inducing Marangoni stresses at the air-liquid interface. However, predictive models for realistic SHS geometries do not yet exist, which has limited the understanding and mitigation of these adverse surfactant effects. To address this issue, we derive a model for laminar, three-dimensional flow over SHS gratings as a function of geometry and soluble surfactant properties, which together encompass 10 dimensionless groups. We establish that the grating length *g* is the key geometric parameter and predict that the ratio between actual and surfactant-free slip increases with *g*². Guided by our model, we perform synergistic numerical simulations and microfluidic experiments, finding good agreement with the theory as we vary surfactant type and SHS geometry. Our model also enables the estimation, based on velocity measurements, of a priori unknown properties of surfactants inherently present in microfluidic systems. For SHSs, we show that surfactant effects can be predicted by a single parameter, representing the ratio between the grating length and the interface length scale beyond which the flow mobilizes the air-water interface. This mobilization length is more sensitive to the surfactant chemistry than to its concentration, such that even trace-level contaminants may significantly increase drag if they are highly surface active. These findings advance the fundamental understanding of realistic interfacial flows and provide practical strategies to maximize superhydrophobic drag reduction.

superhydrophobic surface | drag reduction | surfactant | Marangoni stress | plastron

Superhydrophobic surfaces (SHSs) have the potential to yield enormous technological benefits in fields ranging from microfluidics to maritime transportation, primarily due to their ability to reduce drag (1). Through a combination of hydrophobic chemistry and microscopic surface patterning, these substrates are able to retain a superficial layer of air, thereby producing an apparent slip when in contact with a liquid flow (2). Early theoretical work (3–5) modeled the air pockets trapped within these textures as flat boundaries with no shear, predicting large drag reductions in the laminar regime. Although early experiments found promising levels of drag reduction (6–9), subsequent studies measured a reduced or even nonexistent slip (10–12), pointing at the interfacial stresses induced by surface-active contaminants as one possible cause of this discrepancy. Recently, independent experimental studies have reported time-dependent and spatially complex interfacial dynamics that unequivocally demonstrate the importance of surfactant-induced stresses on SHSs (13, 14). Theoretical and computational works have confirmed the extent to which trace amounts of these surface-active contaminants can reduce slip (15–18). This slip reduction is also consistent with broader findings for small-scale multiphase flows, where environmental levels of surfactants, often extremely difficult to avoid or control, play a central role (19); prominent examples are given by small bubbles rising in water (e.g., refs. 20–23 and references therein) or bubbles probed by atomic force microscopy (24, 25). These flows have been understood through models that include surfactants, sometimes at trace levels that are undetectable by traditional surface tension measurements.

For SHS textures, the concentration gradients that induce Marangoni stresses appear in the streamwise direction, owing to stagnation points at the downstream ends of the interfaces, where advected surfactants accumulate (Fig. 1A and B). Modeling this physical mechanism for realistic SHS geometries is challenging. In addition to the four coupled partial differential equations governing the physics and the ten associated dimensionless numbers (detailed below), there is a major difficulty stemming from the alternating slip/no-slip boundary conditions at the edges where the fluid interface meets the solid substrate. The resulting spatially complex flows constitute a challenge for analytical

## Significance

Trace surfactants, unavoidable in applications, can impair the drag reduction achieved by superhydrophobic surfaces (SHS) as Marangoni stresses immobilize the air-water interface. It is not known how SHS impairment depends on surfactant type and concentration, flow velocity, and SHS geometry; as a result, mitigation strategies are still needed. We introduce a model of this phenomenon and perform simulations and experiments. We find that the interface can be mobilized if it is longer than a critical length scale, which is determined by the surfactant properties, essentially independently of flow velocity. SHS impairment is thereby predicted from a single parameter, namely the ratio of interface length and mobilization scale, providing fundamental insight and practical guidance to achieve superhydrophobic drag reduction.

Author contributions: F.T.-C., F.J.P., J.R.L., F.G., and P.L.-F. designed research; F.T.-C., S.M.S., F.J.P., J.R.L., and P.L.-F. performed research; F.T.-C., S.M.S., F.J.P., J.R.L., and P.L.-F. analyzed data; and F.T.-C., F.J.P., J.R.L., F.G., and P.L.-F. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

Copyright © 2023 the Author(s). Published by PNAS. This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

¹To whom correspondence may be addressed. Email: pfgiz@ucsb.edu.

This article contains supporting information online at http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2211092120/-/DCSupplemental.

Published January 12, 2023.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** (A) Diagram showing the plane channel flow arrangement studied, with an array of slender periodic streamwise gratings on the bottom wall. (B) Unit cell of the SHS (also shown with dashed lines in panel A), periodic in $\hat{x}$ and $\hat{z}$, illustrating the downstream accumulation of surfactant. (C) Streamwise cross-section at mid-grating ($\hat{z} = 0$, also shown with dashed lines in panel B), showing a field of bulk surfactant from a representative simulation, with adsorption/desorption regions at the upstream/downstream ends, respectively. (D) Interface concentration (blue) and shear stress (red), for the same $\hat{z} = 0$ cross-section. Throughout the article, hats denote dimensional quantities.</figcaption></figure>

progress. For this reason, models with surfactants considered only two-dimensional (2D) flows over transverse SHS gratings as this is the simplest geometry that captures detrimental surfactant effects (15, 17). A more realistic configuration is that of streamwise gratings, which are widely used (6, 11, 13, 26), owing to their potential for very high drag reduction in surfactant-free conditions. Gratings have been modeled as infinitely long in surfactant-free theories (3, 5, 27–29); however, modeling surfactant effects requires considering finite streamwise gratings with stagnation points, leading to a three-dimensional (3D)

flow. Theories of realistic gratings inclusive of surfactant are still needed.

Here, we introduce a theory for 3D flow over streamwise SHS gratings with surfactants, by coupling a hydrodynamic solution (for slender, finite gratings with arbitrary shear at the interface) with a scaling analysis of surfactant dynamics (for soluble, dilute surfactants). We use our model to design experiments and simulations where the slip velocity varies across three orders of magnitude, relative to surfactant-free conditions, and thereby achieve a direct comparison between theory and experiments for
realistic SHSs. Our theory can also use velocity measurements to estimate physicochemical parameters of unknown, trace-level surfactants, which are inevitable both in natural and artificial settings. Although the general problem comprises ten dimensionless groups, we show that impairment by trace surfactant is approximately controlled by a single parameter, which depends on surfactant type and concentration, and is independent of flow velocity. Since surface-active molecules are naturally released by polymers widely used in microfabrication (30–33), we expect these results to be valuable over a broad range of fundamental and applied microfluidic research.

## SHS Model for 3D Flow with Surfactants

We consider steady, laminar flow driven by a mean pressure gradient $\hat{G}$ across a channel of half-height $\hat{h}$, where hats denote dimensional quantities. The bottom of the channel is lined with a periodic pattern of slender, rectangular gratings (Fig. 1A). Each gas-liquid interface (the “plastron”) is assumed flat. Due to the periodicity of the array in the streamwise and spanwise directions, we focus on a unit cell consisting of one grating and its surrounding ridges, as depicted in Fig. 1A. The streamwise, wall-normal, and spanwise directions are $\hat{x}$, $\hat{y}$, and $\hat{z}$, respectively, with the coordinate origin at the center of the unit cell (Fig. 1B).

We leverage the disparity of scales between the length $\hat{L}$ and the half-height $\hat{h}$ (Fig. 1B) and define a small parameter $\epsilon = \hat{h}/\hat{L} \ll 1$, which is dimensionless and thus written without hats. Differently from the classic Hele-Shaw flow approximation (28), here we do not assume that the spanwise length scale (the pitch $\hat{P}$) is much larger than $\hat{h}$, since, in microfluidic applications, $\hat{h}$ and $\hat{P}$ are of the order of tens of micrometers, whereas $\hat{L}$ ranges in the millimeter or centimeter scale (6, 7, 11, 26). Consequently, we define the nondimensional coordinates $x = \hat{x}/\hat{L}$, $y = \hat{y}/\hat{h} = \hat{y}/(\epsilon\hat{L})$, and $z = \hat{z}/(\epsilon\hat{L})$. Incompressibility implies that the flow is approximately unidirectional, with the dominant streamwise velocity scaling as $\hat{u} \sim \hat{U}$, whereas the wall-normal and spanwise components scale as $\hat{v} \sim \epsilon\hat{U}$ and $\hat{w} \sim \epsilon\hat{U}$. The velocity scale is $\hat{U} = \hat{h}^2\hat{G}/\hat{\mu}$, with $\hat{\mu}$ the dynamic viscosity. At leading order in $\epsilon$, the Stokes equations for the flow are $\partial_{yy}u + \partial_{zz}u = \partial_{xx}p$ and $\partial_{yy}p = \partial_{zz}p = 0$ (SI Appendix, Flow Field Derivation), where $u(x,y,z) = \hat{u}/\hat{U}$ and $p(x) = \hat{p}/(\hat{G}\hat{L})$ are the dimensionless streamwise velocity and pressure. The unidirectional nature of this leading-order flow is a good approximation far from the downstream and upstream edges of the plastron, i.e., where $|x \pm \phi_x/2| \gg \epsilon$, with $\phi_x$ the streamwise gas fraction (Fig. 1B). Therefore, the asymptotic expansion in $\epsilon$ is singular as is common for thin-gap approximations (34). Since we consider slender gratings with $\epsilon \ll 1$, the regions of validity represent most of the domain, and useful approximations of both local and integrated quantities can be obtained.

No-slip boundary conditions $u = 0$ are imposed at solid walls and ridges. The interface imposes a Marangoni stress $\hat{y}$ determined by the local gradient of interfacial (adsorbed) surfactant. This stress is independent of transverse direction $z$ at leading order in $\epsilon$ (SI Appendix, Flow Field Derivation), and thus, the plastron boundary condition is $\partial_y u|_I = \gamma(x)$, where $\gamma(x) = \hat{y}/(\hat{\mu}\hat{U}/\epsilon\hat{L})$ and the subscript $I$ denotes conditions at the interface.

Note that, for finite gratings, the pressure gradient is not constant in $x$, and $p(x)$ must be determined from

two integral constraints. First, the volume flow rate $Q = \int_{-P/2}^{P/2} \int_{-1}^{1} u(x,y,z) \, dy \, dz$ (where $P = \hat{P}/\hat{h}$ is the normalized pitch; Fig. 1B) must be independent of $x$ to satisfy mass conservation. Second, the pressure drop across the whole unit cell must match the imposed mean pressure gradient, such that $\int_{-1/2}^{1/2} \partial_x p(x) \, dx = -1$. These two conditions lead to an expression for the flow field (as detailed in SI Appendix, Flow Field Derivation),

$$ u(x, y, z) = \begin{cases} \frac{[1 + q_d^\infty(1 - \phi_x)(\gamma)] u_P(y) + [1 - (\gamma)] u_d^\infty(y, z)}{1 + q_d^\infty(1 - \phi_x)} & -\frac{\gamma(x) - (\gamma)}{1 + q_d^\infty} [q_d^\infty u_P(y) - u_d^\infty(y, z)], \quad \text{if } |x| < \frac{\phi_x}{2} \\ \\ \frac{[1 + q_d^\infty(1 - \phi_x)(\gamma)]}{1 + q_d^\infty(1 - \phi_x)} u_P(y), & \quad \text{if } \frac{\phi_x}{2} < |x| \le \frac{1}{2}, \end{cases} \quad [1] $$

where $u_P(y) = (1 - y^2)/2$ is the Poiseuille profile and $u_d^\infty(y, z)$ is the deviation from $u_P(y)$ in the infinite-grating case, where the interface has no stagnation points and thus surfactant effects are absent. In other words, if $\phi_x = 1$ then $u(x, y, z) = u_P(y) + u_d^\infty(y, z)$ and $\gamma(x) = 0$, where $u_d^\infty(y, z)$ is known from previous studies (3, 35). In [1], $q_d^\infty = 3Q_d^\infty/(2P)$, where $Q_d^\infty = \int_{-P/2}^{P/2} \int_{-1}^{1} u_d^\infty(y, z) \, dy \, dz$, and $(\gamma) = \frac{1}{\phi_x} \int_{-\phi_x/2}^{\phi_x/2} \gamma(x) \, dx$ is the average Marangoni shear across the plastron, which varies between 0, for a clean interface, and 1 for a fully immobilized interface. Eq. 1 provides the velocity field as a linear combination of two known simpler solutions ($u_P$ and $u_d^\infty$), relying on parameters which are either prescribed ($\gamma(x)$) or known from the infinite-grating problem ($q_d^\infty$).

The flow field is linked to the surfactant dynamics via $\gamma$, which is found from the equations for soluble surfactant (SI Appendix, Governing Equations)

$$ u \frac{\partial c}{\partial x} + v \frac{\partial c}{\partial y} + w \frac{\partial c}{\partial z} = \frac{1}{\epsilon Pe} \left( \epsilon^2 \frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial y^2} + \frac{\partial^2 c}{\partial z^2} \right), \quad [2a] $$

$$ \frac{\partial (u\Gamma)}{\partial x} + \frac{\partial (w\Gamma)}{\partial z} = \frac{1}{\epsilon Pe_I} \left( \epsilon^2 \frac{\partial^2 \Gamma}{\partial x^2} + \frac{\partial^2 \Gamma}{\partial z^2} \right) + \frac{Bi}{\epsilon} (c_I - \Gamma), \quad [2b] $$

$$ \left. \frac{\partial c}{\partial y} \right|_I = Da(c_I - \Gamma) \quad \text{at the interface,} \quad [2c] $$

$$ \left. \frac{\partial u}{\partial y} \right|_I = \gamma(x) = \epsilon k Ma \frac{\partial \Gamma}{\partial x} \quad \text{at the interface.} \quad [2d] $$

Eqs. 2a and 2b describe the transport of the nondimensional bulk and interfacial surfactant concentrations $c = \hat{c}/\hat{c}_0$ and $\Gamma = \hat{\Gamma}/\hat{\Gamma}_0$, respectively, where $\hat{c}_0$ is the background bulk concentration and $\hat{\Gamma}_0$ is the equilibrium interfacial concentration (Fig. 1 C and D). The adsorption and desorption kinetics are modeled through 2c and the last term in 2b, whereas the Marangoni boundary condition 2d relates the shear stress to the gradient of surfactant concentration at the interface. Six dimensionless groups control the surfactant dynamics in Eq. 2. 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 diffusivities. The Marangoni number $Ma = n_s \hat{R} \hat{T} \hat{\Gamma}_m/(\hat{\mu}\hat{U})$
depends on the maximum interfacial packing concentration $\hat{\Gamma}_m$, the ideal gas constant $\hat{R}$, the temperature $\hat{T}$, and a parameter $n_s$ quantifying the effects of salinity. The Biot $Bi = \hat{h}\hat{k}_d/\hat{U}$ and Damköhler $Da = \hat{h}\hat{k}_a\hat{\Gamma}_m/\hat{D}$ numbers parameterize the effect of kinetics, with $\hat{k}_a$ and $\hat{k}_d$ being the adsorption and desorption rate constants. The normalized background concentration is $k = \hat{k}_a\hat{k}_0/\hat{k}_d$ and can be related to the interfacial concentration at kinetic equilibrium through $k = \hat{\Gamma}_0/\hat{\Gamma}_m$. These six dimensionless groups, in addition to four geometrical parameters $\phi_x$, $\phi_z$, P, and $g = \hat{g}/\hat{h} = \phi_x/\epsilon$, fully describe the flow and surfactant transport problem.

A scaling analysis of Eq. 2, similar to the one performed in ref. 15 for transverse gratings, leads to an expression for $\langle\gamma\rangle$. The derivation, which can be found in SI Appendix, Scaling Theory for Surfactant Transport, is based on the assumption of low normalized concentration ($k \ll 1$), which justifies the choice of Henry kinetics (19) in Eq. 2, and is the case in applications unless substantial amounts of surfactant are deliberately added (12–14). The stress is approximated as spatially uniform, i.e., $\gamma(x) \approx \langle\gamma\rangle$, as has been found to be the case in small-scale applications (15) such that the term on the second line of Eq. 1 is negligible.

To quantify slip and enable comparison with experiments, we use the centerline slip velocity $u_{Ic}$, defined as the velocity along the centerline of the interface, $u_{Ic} = u(x, y = -1, z = 0)$, as illustrated in Fig. 1C. We use $u_{Ic}$ because it can be measured with greater ease and accuracy than the local slip length $\lambda = u_I / \partial_y u_I$, which requires estimation of velocity gradients (12). From the combination of the flow field from Eq. 1 and the expression for $\langle\gamma\rangle$ obtained from the scaling of Eq. 2, we obtain

$$ \frac{u_{Ic}}{u_{Ic}^{\text{clean}}} = 1 - \frac{a_1 k M a u_{Ic}^{\text{clean}}}{a_1 k M a u_{Ic}^{\text{clean}} + a_2 \frac{B_i g^2}{1 + \delta Da} + \frac{1}{Pe_I}}. \quad [3] $$

Our model also yields an expression for the effective slip length $\lambda_e$, defined through the wall-averaged Navier slip boundary condition that would result in the same flow rate as alternating no-slip/slip boundary conditions on the bottom wall (1, 5, 15, 35),

$$ \frac{\lambda_e}{\lambda_e^{\text{clean}}} = 1 - \frac{a_1 k M a u_{Ic}^{\text{clean}} \left( 1 + \frac{\lambda_e^{\text{clean}}}{2} \right)}{a_1 k M a u_{Ic}^{\text{clean}} \left( 1 + \frac{\lambda_e^{\text{clean}}}{2} \right) + a_2 \frac{B_i g^2}{1 + \delta Da} + \frac{1}{Pe_I}}. \quad [4] $$

In Eqs. 3 and 4, $\delta = \hat{\delta}/\hat{h}$ is the concentration boundary layer thickness (Fig. 1C), modeled as $\delta(g, Pe) = a_3(1 + a_4 Pe/g)^{-1/3}$ following a canonical Lévêque scaling (34). Here, $u_{Ic}^{\text{clean}}$ is the centerline slip velocity for clean finite-length gratings, found setting $\gamma(x) = \langle\gamma\rangle = 0$ in Eq. 1. This leads to $u_{Ic}^{\text{clean}} = u_{Ic}^{\infty}/[1 + q_d^{\infty}(1 - \phi_x)]$, where $u_{Ic}^{\infty}$ is the centerline slip velocity for infinite gratings, known from previous studies (3, 35). The effective slip length for clean gratings is $\lambda_e^{\text{clean}} = 2\phi_x q_d^{\infty}/[3 + q_d^{\infty}(3 - 4\phi_x)]$. Although exact solutions of $u_{Ic}^{\infty}$ and $q_d^{\infty}$ require numerical calculations, useful approximations are $u_{Ic}^{\infty} \approx \{[(P/\pi) \cosh^{-1}(\sec(\pi\phi_z/2))]^{n_u} + 2^{n_u}\}^{1/n_u}$ and $q_d^{\infty} \approx \{[(3P/2\pi) \ln(\sec(\pi\phi_z/2))]^{n_q} + (3\phi_z)^{n_q}\}^{1/n_q}$, with $n_u = -1.46, n_q = -1.21$ (these become exact as $P \to 0$ or $P \to \infty$, SI Appendix, Flow Field Derivation).

Eqs. 3 and 4 link the loss of performance to the surfactant-induced stresses; these are enhanced by increasing the product of normalized concentration and Marangoni number $kMa$, which expresses the ratio of Marangoni effects and viscous shear, for a unit of surfactant gradient. Marangoni stresses can be reduced by decreasing the surfactant gradient, achieved by either i) increasing surfactant flux between the interface and the bulk (relative to advection), as captured by $B_i g^2/(1 + \delta Da)$, or ii) increasing interfacial diffusivity (again, compared to advection), expressed by $1/Pe_I$.

The scaling coefficients $a_1, a_2, a_3$, and $a_4$ in Eqs. 3 and 4 are estimated by performing 155 simulations of the full governing equations, spanning a wide range of values in the dimensionless groups to ensure proper coverage of the parameter space (Materials and Methods). Fig. 2A shows good agreement between the model Eq. 3 and simulations across four orders of magnitude in the slip velocity, for $a_1 \approx 0.345, a_2 \approx 0.275, a_3 \approx 5.581$, and $a_4 \approx 3.922$, which are values of order one as expected for scaling coefficients. The simulations also corroborate the assumption $\gamma(x) \approx \langle\gamma\rangle$, as discussed in SI Appendix, Finite-Element Simulations.

Equipped with a 3D theory and a set of numerical simulation results, we aim to identify realistic combinations of the 10 dimensionless parameters that maximize the drag reduction of SHSs. Fig. 2 B and C illustrate how $\lambda_e$ and $u_{Ic}$ change with the spanwise gas fraction $\phi_z$, for a fixed bulk concentration $\hat{c}_0 = 3 \cdot 10^{-4} \text{ mol m}^{-3}$ and $\hat{g} = 1 \text{ mm}$, and for several values of surfactant solubility, in a microchannel with $\hat{h} = 60 \ \mu\text{m}$. Slip is maximized for surfactant with higher solubility, corresponding to lower $\hat{k}_a/\hat{k}_d$. Incidentally, the results depend weakly on $\hat{k}_a$ or $\hat{k}_d$ individually; although these parameters appear separately in Eqs. 2 and 4 through $Bi$ and $Da$, $Da$ is large if $\hat{h}$ is larger than a

<figure><img src="image_4.png" /><figcaption>Fig. 2. (A) Comparison of the centerline slip velocity from 155 numerical simulations with our model prediction Eq. 3. (B) Effective slip length and (C) Centerline slip velocity as a function of the spanwise gas fraction $\phi_z$, for a fixed surfactant concentration $\hat{c}_0$ and varying surfactant solubility ($\hat{k}_a/\hat{k}_d$). Unless noted, the parameters are as in SI Appendix, Table SI.</figcaption></figure>
few micrometer, such that the term $Bi/(1 + \delta Da) \approx Bi/(\delta Da)$, which is a function of $\hat{\kappa}_a/\hat{\kappa}_d$ (SI Appendix, Discussion of the Mobilization Length).

Since the chemical properties of naturally occurring surfactants are virtually impossible to control in practice (32), we focus on geometrical parameters. We observe that $\phi_z$ has a negligible impact on surfactant effects, even as $\phi_z \to 1$. Mathematically, $\phi_z$ affects Eq. 3 only through the surfactant-independent term $u_{lc}^{\text{clean}}(P, \phi_z, \phi_x)$, which is at most of order one. Eqs. 3 and 4 reveal that slip is maximized by increasing the grating length $g$, which progressively overcomes surfactant effects, undergoing a transition of the form $u_{lc} \sim g^2$ and ultimately approaching the asymptotes $u_{lc} \to u_{lc}^{\text{clean}}$, $\lambda_e \to \lambda_e^{\text{clean}}$. This transition is challenging to simulate due to the large computational cost of long domains needed at large $g$.

## Experiments Demonstrate Effect of Grating Length

To acquire data at large $g$ and test the prediction of a slip transition, we build microfluidic devices using polydimethylsiloxane (PDMS), as shown in Fig. 3 A and B (Materials and Methods). The channel upper wall consists of streamwise gratings of pitch $\hat{P} = 60$ µm and spanwise gas fraction $\phi_z = 2/3$. The channel half-height is $\hat{h} = 60$ µm, and the depth of the grating trenches is $\hat{d} = 25$ µm, enough to ensure a stable plastron throughout each experiment. We test gratings with $\hat{g} = 15, 25, 35$ and $45$ mm, separated in the streamwise direction by solid ridges of length $\hat{L} - \hat{g} = 20$ µm. We employ a confocal microscope and microparticle image velocimetry (μ-PIV) in a setup similar to the one in ref. 13. A syringe pump provides a constant flow rate $\hat{Q}_{\text{TOT}} = 1.152 \, \mu\text{L min}^{-1}$. We use deionized water without any additives, since it has been established that the unavoidable amounts of surfactant naturally present in similar microfluidic settings are sufficient to induce significant stresses at the plastron (12–14). The μ-PIV beads used to seed the liquid (Fig. 3E) are thoroughly prewashed to remove their added surfactant (16), and we also follow a cleaning protocol for the syringes and tubing (Materials and Methods).

The flow velocity is measured over two adjacent gratings, at several distances $\Delta\hat{y}$ from the interface, as defined in Fig. 3C. Examples of velocity profiles are displayed in Fig. 3D, for $\hat{g} = 45$ mm. The flow over the solid ridges is consistent with the no-slip condition at the wall, whereas velocity increases noticeably over the gratings. These vertically spaced profiles around the grating centerline are extrapolated to obtain the slip velocity at the interface, shown by the black dashed lines in Fig. 3D; the centerline slip velocity is marked by the bull's eyes in the figure. The measured slip velocities $\hat{u}_{lc}$ are only 6 to 10% of the values $\hat{u}_{lc}^{\text{clean}}$ predicted by surfactant-free theories (Fig. 3F), consistently with prior experiments (10–14).

Comparing quantitatively these experimental measurements to the predictions from our model requires assumptions on the type and amount of surfactant present in the channel. Although some parameter values are known and others can be accurately estimated, the normalized surfactant concentration $k$ and the kinetic rate adsorption and desorption constants $\hat{\kappa}_a$ and $\hat{\kappa}_d$ can vary across a broad range. Nevertheless, it is possible to combine our model for the slip velocity Eq. 3 with previous experimental results (13) to obtain an estimate, as described in detail in SI Appendix, Estimate of Surfactant Parameters. We find approximate ranges for the normalized concentration

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** (A) Top view and (B) Cross-section of the microfluidic channels used in the experiments. (C) Cross-section of a unit cell, as defined in Fig. 1B (note the inverted microscope setup). The distance from the interface is $\Delta\hat{y}$. (D) Example of velocity profiles at different distances from the interface, for a grating length $\hat{g} = 45$ mm. The dashed line denotes the linearly extrapolated slip velocity, whereas the shadings show standard error (details in SI Appendix, Experimental Methods). The bull's eyes mark the centerline slip velocity. (E) Micrograph of the gratings shown in (D), with μ-PIV particles appearing in green. (F) Ratio between the actual centerline slip velocity and the surfactant-free value (i.e., "clean"), $u_{lc}/u_{lc}^{\text{clean}} = \hat{u}_{lc}/\hat{u}_{lc}^{\text{clean}}$, from our theory [3], simulations, and experiments, as well as prior experiments of refs. (13) and (14). Ref. 14 used an annular grating in a rheometer, with an effectively infinite groove length. The theoretical prediction from our model 3 and the simulations use a best estimate of the unknown surfactant parameters (SI Appendix, Estimate of Surfactant Parameters), with the shaded region denoting a range of plausible levels of contamination within our estimates. The uncertainty in the present experimental data is smaller than the size of the symbols.</figcaption></figure>

7.3 · 10⁻³ $\lesssim$ $k$ $\lesssim$ 1 · 10⁻¹ and for the ratio of constants 7.1 · 10¹ m³ mol⁻¹ $\lesssim$ $\hat{\kappa}_a/\hat{\kappa}_d$ $\lesssim$ 1.8 · 10³ m³ mol⁻¹.

Choosing the mid-range values $k = 3.6 · 10⁻²$ and $\hat{\kappa}_a/\hat{\kappa}_d = 1.2 · 10² m³ mol⁻¹$, our predictions of the slip velocity show
<figure><img src="image_6.png" /><figcaption>Fig. 4. (A) Lines (from our model) showing the normalized surfactant concentration k as a function of grating length g, that yields a slip length that is 50% of the ideal, "clean" value. Colors denote different surfactant solubilities, expressed by the ratio of the adsorption and desorption constants $\kappa_a/\kappa_d$. Shaded bands show the weak effect of changing the flow velocity $\dot{U}$ across two orders of magnitude. (B) Experimental and numerical data for the slip velocity $u_{lc}$, plotted against g, as one varies surfactant properties, flow velocity, and SHS geometry, together characterized by ten dimensionless groups (including g). (C) Normalizing g by the mobilization length $L_m$ approximately collapses the same data onto a single curve, governed by only one dimensionless group. Here, the factor $(u_{lc}^{\text{clean}})^{1/2}$ helps collapse data across a range of gas fractions, which give a wide range of values for $u_{lc}^{\text{clean}}$. In practice, the normalized clean interface velocity $u_{lc}^{\text{clean}}$ is usually of O(1), so this factor could be omitted in the normalization for g.</figcaption></figure>

good agreement with our experimental data and with previous studies (13, 14), as illustrated in Fig. 3F. At small g, measuring the small slip velocity with high precision is challenging; we performed finite-element simulations with the same surfactant properties as in the experiments, shown by the blue circles in Fig. 3F. These simulations are restricted to $g < 60$ (approximately 3.5 mm in practice), as computational cost increases with g. Simulations and experiments agree with Eq. 3, showing increased slip consistent with the theoretical prediction as g increases and confirming the key impact of this geometric parameter in controlling drag reduction. Furthermore, all our measurements are consistent with the range of k and $\hat{\kappa}_a/\hat{\kappa}_d$ that we estimated from previous experiments in a different laboratory (13).

To assess practical implications for general microfluidic flows, we plot in Fig. 4A the combinations of g and k at which the slip length reaches 50% of the surfactant-free value (i.e., $\lambda_e/\lambda_e^{\text{clean}} = 0.5$, plotted with solid lines), for surfactant solubility values $\hat{\kappa}_a/\hat{\kappa}_d$ (plotted with different colors) and a range of representative velocities $\hat{U}$ (shown by shaded colored bands around each line) found in small-scale applications. The gray horizontal band shows the range $7 \cdot 10^{-3} \lesssim k \lesssim 0.1$ estimated for microfluidic experiments. At larger values $k > 0.1$, the risk of plastron collapse increases significantly due to capillary effects. Increasing or reducing the velocity $\hat{U}$ by an order of magnitude has a relatively weak effect on the slip, as shown by the narrow bands around each line in Fig. 4A. In addition, Fig. 4A shows that varying $\hat{\kappa}_a/\hat{\kappa}_d$ leads to uniformly shifted contours, which remain approximately parallel. Together with the weak dependence on velocity, these results suggest that the 50% threshold for slip could be expressed through a simpler underlying criterion.

## Single Lengthscale Predicts Interface Mobility

Fig. 3F shows that surfactant impairment is strongly dependent on the grating length g. We therefore define the mobilization length $L_m$ as the value of g that gives significant slip, such that the

local slip length on the interface is $\lambda_I(L_m) \sim 1$ (or in dimensional terms, $\hat{\lambda}_I(\hat{L}_m/\hat{h}) \sim \hat{h}$). We seek a scaling for $\lambda_I$ in terms of g, starting from the definition of $\lambda_I$ and using Eq. 2d,

$$ \lambda_I = \frac{u_I}{\gamma} \sim \frac{u_I}{kMa \frac{\Delta\Gamma}{g}} \qquad [5] $$

Note that x is normalized by $\hat{L}$, whereas g is normalized by $\hat{h}$, such that gradients along the interface occur over an x length scale $\hat{g}/\hat{L} \sim \varepsilon g$. Here, $\Delta\Gamma$ is the order of magnitude of the change in $\Gamma$ along the interface, illustrated in Fig. 1D. More specifically, we assume that $c_I$ ranges from $1 - \Delta c_I$ to $1 + \Delta c_I$, and similarly, $\Gamma$ is between $1 - \Delta\Gamma$ and $1 + \Delta\Gamma$, such that at some point near the middle of the interface, $\Gamma \sim c \sim 1$. Scaling Eq. 2c,

$$ -\frac{\Delta c_I}{\delta} \sim Da(\Delta c_I - \Delta\Gamma), \quad \text{i.e.} \quad \Delta\Gamma \sim \Delta c_I \left( \frac{1}{\delta Da} + 1 \right). \qquad [6] $$

We find $\Delta c_I$ by scaling Eq. 2b, which represents a balance between advection, which acts to establish a surfactant gradient, and diffusion and kinetics, which act to extinguish the gradient. We consider changes along the x-direction, between the edge of the interface (where $u = 0$) and a point near the middle (where $u \sim u_I$ and $\Gamma \sim 1$). Then, $\partial_x(u\Gamma) \sim u_I/(\varepsilon g)$ and $\partial_{xx}\Gamma \sim \Delta\Gamma/(\varepsilon g)^2$. Writing the kinetic flux $Bi(c_I - \Gamma) = (Bi/Da) \partial_y c_I$ and using Eq. 6 to eliminate $\Delta\Gamma$, Eq. 2b yields

$$ \frac{u_I}{g} \sim \left( Bi + \frac{1 + \delta Da}{Pe_I g^2} \right) \frac{\Delta c_I}{\delta Da}. \qquad [7] $$

The second term in the parenthesis on the right-hand side is negligible if diffusion is weak compared to kinetics, such that $g^2 \gg (1 + \delta Da)/(Bi Pe_I)$. Note that $\delta \sim 1$ unless the gratings are very short, and recall that Da is large if $\hat{h}$ is larger than a few micrometers, such that $(1 + \delta Da) \sim Da$. In dimensional form, neglecting diffusion in Eq. 7 then requires $\hat{g}^2 \gg (\hat{L}_d^{\text{mod}})^2$, where
$$ \hat{L}_d^{\text{mod}} := \left( \hat{h} \hat{L}_d \frac{\hat{D}_I}{\hat{D}} \right)^{1/2}, \qquad [8] $$

and $\hat{L}_d = \hat{\kappa}_a \hat{\Gamma}_m / \hat{\kappa}_d$ is the depletion length, which arises commonly in models of soluble surfactant (19). We define $\hat{L}_d^{\text{mod}}$ as a modified depletion length, representing the grating length above which interface diffusion is small compared to kinetics. Omitting the diffusive term in [7], and combining with [5] and [6] to eliminate $\Delta\Gamma$ and $\Delta c_I$, we find

$$ \lambda_I \sim \frac{g^2 Bi}{k Ma Da}. \qquad [9] $$

This scaling shows concisely that $\lambda_I$ depends on $g^2$, consistently with Eq. 4 and with Fig. 3F. If $\lambda_I \sim 1$ when $g \sim L_m$, Eq. 9 defines the mobilization length; in dimensional form, we find

$$ \hat{L}_m := \hat{h} \frac{\hat{\kappa}_a \hat{\Gamma}_m}{\hat{\kappa}_d} \left( \frac{n_s \hat{R} \hat{T} \hat{c}_0}{\hat{\mu} \hat{D}} \right)^{1/2} . \qquad [10] $$

Marangoni stresses become negligible if $\hat{g}^2 \gg \hat{L}_m^2$ as transport between the bulk and the interface suppresses the surfactant gradient that would otherwise be established by advection. For our experiments (SI Appendix, Table SI), Eq. 10 predicts $\hat{L}_m \approx 4.3$ cm, indicating that gratings must be at least several centimeters long to minimize Marangoni stresses, consistently with Fig. 3F. Remarkably, $\hat{L}_m$ depends linearly on $\hat{\kappa}_a / \hat{\kappa}_d$ but only on the square root of $\hat{c}_0$, showing higher sensitivity to the type of surfactant than to its concentration.

For small-scale applications, and using properties of known surfactants, we find that $\hat{L}_d^{\text{mod}}$ is much smaller than $\hat{L}_m$; for example, $\hat{L}_d^{\text{mod}} \approx 170 \ \mu\text{m}$ in our experiments (SI Appendix, Table SI and SI Appendix, Discussion of the Mobilization Length). Therefore, the mobilization length $\hat{L}_m$ is the key length scale that determines the slip of a given SHS.

Our analysis suggests that it should be possible to write simplified expressions for $u_{Ic}$ and $\lambda_e$ in terms of $g/L_m$. With the approximation $(1 + \delta Da) \approx Da$ (as discussed above) and assuming $\hat{L}_m \gg \hat{L}_d^{\text{mod}}$ (such that the terms $1/Pe_I$ are negligible), Eqs. 3 and 4 simplify to the one-parameter curves

$$ \frac{u_{Ic}}{u_{Ic}^{\text{clean}}} = 1 - \frac{1}{1 + \frac{a_2}{a_1 a_3} \left( \frac{g}{L_m \sqrt{u_{Ic}^{\text{clean}}}} \right)^2}, \qquad [11] $$

$$ \frac{\lambda_e}{\lambda_e^{\text{clean}}} = 1 - \frac{1}{1 + \frac{a_2}{a_1 a_3} \left( \frac{g}{L_m \sqrt{u_{Ic}^{\text{clean}} \left(1 + \frac{\lambda_e^{\text{clean}}}{2}\right)}} \right)^2}. \qquad [12] $$

We reexamine the 155 simulations shown in Fig. 2A and set aside 24 cases that are not achievable in reality (e.g. involving unphysically small diffusivities). Fig. 4B shows the relative slip $u_{Ic}/u_{Ic}^{\text{clean}}$ versus $g$ for the remaining 131 simulations and for our experiments. These data span different surfactant properties, grating geometries and flow velocities, and exhibit large scatter, illustrating how $g$ alone is insufficient to predict slip. However, when $g$ is normalized by the mobilization length $L_m$, the same data collapse near the one-parameter curve given by Eq. 11, as shown in Fig. 4C.

## Discussion and Outlook

Regarding the surfactant type inherent to our experiments, Eqs. 10 and 11 suggest a surfactant with large $\hat{\kappa}_a / \hat{\kappa}_d$, implying low solubility. This is consistent with previous findings that PDMS used in microfluidic channels (including our experiments) releases uncrosslinked oligomer chains (36–38), which are surface active (31, 32), and which have also been detected in solution (30, 33). The mass fractions reported in ref. 33, in combination with the oligomer chain lengths found in ref. 30, lead to concentrations $\hat{c}_0 \sim O(10^{-4} - 10^{-2}) \text{ mol m}^{-3}$, compatible with our estimates. Incidentally, in other contexts, PDMS has sometimes been approximately modeled as insoluble (39, 40); under this assumption, our theory yields $u_{Ic} = u_{Ic}^{\text{clean}} / (1 + a_{\text{ins}} Ma_{\text{ins}} u_{Ic}^{\text{clean}})$ (SI Appendix, Scaling Theory for Surfactant Transport), where $a_{\text{ins}}$ is a scaling coefficient, $Ma_{\text{ins}} = n_s \hat{R} \hat{T} \hat{\Gamma}_0 \hat{h} / (\hat{\mu} \hat{D}_I)$ is a Marangoni number, and $\hat{\Gamma}_0$ is the average interfacial surfactant concentration. Note that this expression for $u_{Ic}$ does not depend on $g$, inconsistently with the experimental results in Fig. 3F. This highlights the importance of including solubility in models of surfactant dynamics on SHSs.

The results described here provide insights into the slip and drag of superhydrophobic surfaces in realistic conditions. Our theory for slender, finite gratings (which are widely used) enables comparisons with experiments, where inherent surfactants must be accounted for. The hydrodynamic component of the model from Eq. 1 can also quantify SHS performance in the presence of general, nonuniform shear stresses at the air-water interface, thus circumventing the need for computationally expensive simulations. In addition, we have shown that a single mobilization length scale arising from the theory can serve as a guide in the design of SHS textures that mitigate surfactant effects. This laminar theory is also a stepping stone toward predicting surfactant effects in turbulent flow. As a first approximation, the mobilization length could be estimated by replacing the channel half-height $\hat{h}$ in Eq. 10 with the thickness of the viscous sublayer. Finally, since $\hat{L}_m$ depends primarily on surfactant properties and on the shear length scale $\hat{h}$, we may also expect that the concept of mobilization length, derived here for streamwise gratings, could qualitatively apply to other SHS textures.

## Materials and Methods

**Finite-Element Simulations.** We solved the full governing equations and boundary conditions (in dimensional form, detailed in SI Appendix, Governing Equations) in three dimensions using COMSOL Multiphysics 5.5. We performed a total of 155 simulations using different grating geometries, flow velocities, and surfactant properties in order to span a large portion of the parameter space characterized by the ten dimensionless numbers of the problem. The domain was one half of the SHS unit cell depicted in Fig. 1B, with $\hat{z}$ between $\hat{z} = 0$ and $\hat{z} = \hat{P}/2$ due to the spanwise symmetry of the solution. The volume is meshed with tetrahedral elements, with the finest ones (with a minimum element size of $1.5 \cdot 10^{-9}$ m) around the upstream and downstream edges of the interface $\hat{x} = \pm \phi_x \hat{L}/2$ (SI Appendix, Finite-Element Simulations). We used the Creeping Flow module for the flow field, the Dilute Species Transport module for the transport of a bulk surfactant, and the transport of an interfacial surfactant was implemented through a General Form Boundary PDE. The Marangoni boundary conditions were enforced through a Weak Contribution constraint, as was the condition that fixed the mean bulk concentration to be $\hat{c}_0$. The system of nonlinear equations was solved through a Newton-type iterative method using the PARDISO direct solver for the linear system at each iteration. All simulations satisfied a relative tolerance for convergence of $1 \cdot 10^{-5}$. We used linear
elements for the pressure, bulk concentration and interfacial concentration, and either linear or quadratic elements for the velocity field, depending on the computational demands of each simulation.

**Microchannel Fabrication.** Microfluidic channels with an array of parallel SHS gratings on their ceiling (Fig. 3 A and B) were built by casting PDMS (Sylgard 184) over a mold obtained by two-layer photolithography. The photoresist was SU-8 (Microchem SU-8 3025 and Microchem SU-8 3050). The chips were bonded to 0.1 mm-thick glass coverslips (Bellco Glass 1916-25075) through untreated adhesion. Every coverslip was washed with isopropyl alcohol, then with 18 MΩ cm DI water, and finally dried with nitrogen before the microfluidic chip was attached. The static contact angle of water droplets was measured to be higher than 100° over samples of smooth PDMS, and to increase further over samples of textured PDMS. This is consistent with previous measurements for untreated PDMS (41) and demonstrates the superhydrophobicity of the substrate. The total width of the microfluidic channel was set to $\hat{W} = 2$ mm (Fig. 3A) to ensure an approximately periodic flow in $\hat{z}$ over the gratings far way from lateral walls (given that $\hat{W} \gg \hat{h}$).

**Experimental Setup.** A glass syringe (Hamilton Gastight) was filled with particle-seeded DI water, which was driven through the microchannels using a syringe pump (KD Legato 111). We used the barrel of a plastic syringe (BD Luer-Lok) as an outlet reservoir open to the room, to impose atmospheric pressure at the end of the circuit. The height of this reservoir was adjusted with a vertical translation stage (Thorlabs VAP10) at the beginning of each experiment to ensure that the plastron at each grating remained approximately flat, by controlling the average pressure in the microchannel. The microchannel was connected to the syringe and reservoir through plastic tubing (Tygon S3). All circuit elements were thoroughly prewashed with 18 MΩ cm DI water, following a protocol described in *SI Appendix, Experimental Methods*.

**Confocal Microscopy.** The tracer particles (ThermoFisher FluoSpheres carboxylate 0.5 µm diameter) were washed using a centrifuge (Eppendorf 5418) to separate them from the buffer solution, which was then discarded and replenished with 18 MΩ cm DI water. This process was repeated three times before each experiment to essentially eliminate surfactant contamination from the particle solution. The flow was observed with a confocal microscope (Leica SP8 Resonant Scanning), using a 40× water objective (as in Fig. 3B). The microfluidic device was enclosed in a stage top chamber (Okolab H101-K-FRAME) with a controlled temperature set to $\hat{T} = 296$ K. Using the bright field imaging of the

microscope (Fig. 3E), we focused on two adjacent gratings. We avoided imaging the five gratings closest to each lateral side wall of the channel to prevent effects related to the loss of periodicity. The fluorescence imaging of the microscope (superimposed on the image in Fig. 3E) revealed the positions of the tracer beads in each snapshot, which we obtained at a rate of between 20 and 28 frames per second. All the data were taken at the center of the grating in the streamwise direction (i.e. $\hat{x} \approx 0$) and at several distinct $\hat{y}$-planes close to the interface (Fig. 3D).

**Image Analysis and Micro-PIV.** The μ-PIV analysis was performed with the open-source MATLAB toolbox PIVlab (42), using an acquisition window of approximately 125 µm × 125 µm. The velocity field obtained for a given window was averaged in time and along the streamwise $\hat{x}$ direction to obtain the spanwise velocity profiles depicted in Fig. 3D, for different distances away from the interface. To extract the centerline slip velocity $\hat{u}_{Ic}$ (at $\hat{y} = -\hat{h}$ and $\hat{z} = 0$), we performed a linear least-squares fit using data from between three to five $\hat{y}$-planes. We only used data in a neighborhood of the grating center $\hat{z} = 0$ (Fig. 3D), since velocity profiles at $\hat{y} = -\hat{h}$ were not smooth due to the transitions between the interfaces and the solid ridges at $\hat{z} = \pm \phi_z \hat{P}/2$. The uncertainty for $\hat{u}_{Ic}$ was calculated by accounting for how uncertainties in the velocity measurements and in the $\hat{y}$-coordinate of the interface (±1 µm) propagated through the fitting procedure.

**Data, Materials, and Software Availability.** All study data are included in the article and/or *SI Appendix*.

**ACKNOWLEDGMENTS.** We thank David Bothman, Benjamin Lopez and Rachel Schoeppner for technical assistance. This work is supported by NSF CAREER 2048234, ARO MURI W911NF-17-1-0306, ONR MURI N00014-17-1-2676, the California NanoSystems Institute, and the European Union's Horizon 2020 research and innovation program under Marie Skłodowska-Curie grant agreement No. 798411. We acknowledge the use of the NRI-MCDB Microscopy Facility supported by NSF MRI grant DBI-1625770. A portion of this work was performed in the Microfluidics Laboratory within the California NanoSystems Institute, supported by the University of California Santa Barbara, and the University of California, Office of the President.

Author affiliations: <sup>a</sup>Department of Mechanical Engineering, University of California, Santa Barbara, CA 93106; <sup>b</sup>Institute of Environmental Engineering, Department of Civil, Environmental and Geomatic Engineering, Eidgenössische Technische Hochschule (ETH) Zürich, 8093 Zürich, Switzerland; and <sup>c</sup>Department of Mathematics, Alan Turing Building, University of Manchester, Manchester M13 9PL, United Kingdom

1. J. P. Rothstein, Slip on superhydrophobic surfaces. *Annu. Rev. Fluid Mech.* **42**, 89-109 (2010).

2. E. Lauga, M. Brenner, H. Stone, "Microfluidics: The no-slip boundary condition" in C. Tropea, A. L. Yarin, J. F. Foss, Eds. (Springer, Berlin, Heidelberg, 2007), pp. 1219-1240.

3. J. R. Philip, Flows satisfying mixed no-slip and no-shear conditions. *Z. Angew. Math. Phys.* **23**, 353-372 (1972).

4. J. R. Philip, Integral properties of flows satisfying mixed no-slip and no-shear conditions. *Z. Angew. Math. Phys.* **23**, 960-968 (1972).

5. E. Lauga, H. A. Stone, Effective slip in pressure-driven Stokes flow. *J. Fluid Mech.* **489**, 55-77 (2003).

6. J. Ou, B. Perot, J. P. Rothstein, Laminar drag reduction in microchannels using ultrahydrophobic surfaces. *Phys. Fluids* **16**, 4635-4643 (2004).

7. J. Ou, J. P. Rothstein, Direct velocity measurements of the flow past drag-reducing ultrahydrophobic surfaces. *Phys. Fluids* **17**, 103606 (2005).

8. R. Truesdell, A. Mammoli, P. Vorobieff, J. van Swol, C. J. Brinker, Drag reduction on a patterned superhydrophobic surface. *Phys. Rev. Lett.* **97**, 044504 (2006).

9. C. Lee, C.-H. Choi, C.-J. Kim, Structured surfaces for a giant liquid slip. *Phys. Rev. Lett.* **101**, 064501 (2008).

10. T. J. Kim, C. Hidrovo, Pressure and partial wetting effects on superhydrophobic friction reduction in microchannel flow. *Phys. Fluids* **24**, 112003 (2012).

11. G. Bolognesi, C. Cottin-Bizonne, C. Pirat, Evidence of slippage breakdown for a superhydrophobic microchannel. *Phys. Fluids* **26**, 082004 (2014).

12. D. Schäffel, K. Koynov, D. Vollmer, H. J. Butt, C. Schönecker, Local flow field and slip length of superhydrophobic surfaces. *Phys. Rev. Lett.* **116**, 134501 (2016).

13. F. J. Peaudecerf, J. R. Landel, R. E. Goldstein, P. Luzzatto-Fegiz, Traces of surfactants can severely limit the drag reduction of superhydrophobic surfaces. *Proc. Natl. Acad. Sci. U.S.A.* **114**, 7254-7259 (2017).

14. D. Song et al., Effect of a surface tension gradient on the slip flow along a superhydrophobic air-water interface. *Phys. Rev. Fluids* **3**, 033303 (2018).

15. J. R. Landel et al., A theory for the slip and drag of superhydrophobic surfaces with surfactant. *J. Fluid Mech.* **883**, A18 (2020).

16. H. Li et al., Three-dimensional backflow at liquid-gas interface induced by surfactant. *J. Fluid Mech.* **899**, A8 (2020).

17. T. Baier, S. Hardt, Influence of insoluble surfactants on shear flow over a surface in Cassie state at large péclet numbers. *J. Fluid Mech.* **907**, A3 (2021).

18. J. Sundin, S. Bagheri, Slip of submerged two-dimensional liquid-infused surfaces in the presence of surfactants. *J. Fluid Mech.* **950**, A35 (2022).

19. H. Manikantan, T. M. Squires, Surfactant dynamics: Hidden variables controlling fluid flows. *J. Fluid Mech.* **892**, P1 (2020).

20. W. Bond, D. A. Newton, Bubbles, drops and Stokes' law. *Phil. Mag.* **5**, 794-800 (1928).

21. A. N. Frumkin, V. G. Levich, Effect of surface-active substances on movements at the boundaries of liquid phases. *Zhuh. Fiz. Khim.* **21**, 1183-1204 (1947).

22. V. Levich, *Physicochemical Hydrodynamics* (Prentice Hall, 1962).

23. R. Palaparthi, D. T. Papageorgiou, C. Maldarelli, Theory and experiments on the stagnant cap regime in the motion of spherical surfactant-laden bubbles. *J. Fluid Mech.* **559**, 1-44 (2006).

24. O. Manor et al., Hydrodynamic boundary conditions and dynamic forces between bubbles and surfaces. *Phys. Rev. Lett.* **101**, 024501 (2008).

25. A. Maali et al., Viscoelastic drag forces and crossover from no-slip to slip boundary conditions for flow near air-water interfaces. *Phys. Rev. Lett.* **118**, 084501 (2017).

26. P. Tsai et al., Quantifying effective slip length over micropatterned hydrophobic surfaces. *Phys. Fluids* **21**, 112002 (2009).

27. C. J. Teo, B. C. Khoo, Analysis of Stokes flow in microchannels with superhydrophobic surfaces containing a periodic array of micro-grooves. *Microfluid. Nanofluid.* **7**, 353-382 (2009).

28. F. Feuillebois, M. Z. Bazant, O. I. Vinogradova, Effective slip over superhydrophobic surfaces in thin channels. *Phys. Rev. Lett.* **102**, 026001 (2009).

29. E. S. Asmolov, O. I. Vinogradova, Effective slip boundary conditions for arbitrary one-dimensional surfaces. *J. Fluid Mech.* **706**, 108-117 (2012).
30. K. J. Regehr et al., Biological implications of polydimethylsiloxane-based microfluidic cell culture. *Lab Chip* **9**, 2132–2139 (2009).

31. A. Hourlier-Fargette, A. Antkowiak, A. Chateauminois, S. Neukirch, Role of uncrosslinked chains in droplets dynamics on silicone elastomers. *Soft. Matter* **13**, 3484–3491 (2017).

32. A. Hourlier-Fargette, J. Dervaux, A. Antkowiak, S. Neukirch, Extraction of silicone uncrosslinked chains at air-water-polydimethylsiloxane triple lines. *Langmuir* **34**, 12244–12250 (2018).

33. S.-S.D. Carter et al., PDMS leaching and its implications for on-chip studies focusing on bone regeneration applications. *Organs-on-a-Chip* **2**, 100004 (2020).

34. L. G. Leal, *Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes. Cambridge Series in Chemical Engineering* (Cambridge University Press, 2007).

35. C. J. Teo, B. C. Khoo, Analysis of Stokes' flow in microchannels with superhydrophobic surfaces containing a periodic array of micro-grooves. *Microfluid. Nanofluid.* **7**, 353 (2008).

36. J. Kim, M. K. Chaudhury, M. J. Owen, Modeling hydrophobic recovery of electrically discharged polydimethylsiloxane elastomers. *J. Colloid Interface Sci.* **293**, 364–375 (2006).

37. D. T. Eddington, J. P. Puccinelli, D. J. Beebe, Thermal aging and reduced hydrophobic recovery of polydimethylsiloxane. *Sens. Actuators B Chem.* **114**, 170–172 (2006).

38. W. S. Y. Wong et al., Adaptive wetting of polydimethylsiloxane. *Langmuir* **36**, 7236–7245 (2020).

39. L. T. Lee, E. K. Mann, D. Langevin, B. Farnoux, Neutron reflectivity and ellipsometry studies of a polymer molecular layer spread on the water surface. *Langmuir* **7**, 3076–3080 (1991).

40. V. Bergeron, D. Langevin, Monolayer spreading of polydimethylsiloxane oil on surfactant solutions. *Phys. Rev. Lett.* **76**, 3152 (1996).

41. I. Miranda et al., Properties and applications of PDMS for biomedical engineering: A review. *J. Funct. Biomater.* **13**, 2 (2021).

42. W. Thielicke, E. J. Stamhuis, PIVlab - towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB. *J. Open Res. Softw.* **2**, e30 (2014).