A single parameter can predict surfactant impairment of superhydrophobic drag reduction

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.


Problem parameters
The parameters appearing in Eqs. S1a-S3k are defined in Tables SI and SII, together with the values that they take in our experiments. Since, as explained in the main text, the surfactant type and concentration in the liquid are unknown, only an estimate can be obtained in some cases (see Section Estimate of surfactant parameters for details). We choose g =ĝ/ĥ, P =P /ĥ, ϕx and ϕz as the four independent geometric parameters of the problem, noting that ε can then be obtained as ε = ϕx/g. Tables SI and SII also include the values of the length scalesL mod d andLm, defined in Eqs. 10 and 9 of the main text.

Flow field derivation
Assumption of a spanwise constant interface shear stress. Note that, although ε ≪ 1 and k ≪ 1 in the conditions considered in our study (see Section Scaling theory for surfactant transport), the product εkM a appearing in [S3h] and [S3i] is typically not small, since the Marangoni number is expected to be large, i.e. M a ≫ 1 (see estimates in Table SII) and the term N (Γ) ≈ 1 as long as k and k|A| remain small. In fact, Eq. S3h implies that only when εkM a ≳ 1 the Marangoni stresses at the interface are non-negligible, as it is observed experimentally (4-7). Since ε ≪ 1, it is possible to assume that εkM a ≳ 1 ≫ ε 2 , and in that case it follows from [S3i] that ∂zΓ ≈ 0 at leading order in ε. As detailed in the main text, the asymptotic expansion leading to [S3h] and [S3i] is singular, and thus the approximation ∂zΓ ≈ 0 is valid only in regions far from the upstream and downstream stagnation points, i.e. for |x ± ϕx/2| ≫ ε. Indeed, our finite-element simulations of the full problem confirm that this approximation remains valid in all the regimes considered (as illustrated in Fig. S1b, showing the contours of Γ). The Marangoni shear γ(x) = ∂yu| I is thus also assumed to be independent of z and only dependent on x, following [S3h]. Channel half height (see Fig. 1)ĥ Pitch (see Fig. 1  Interface Péclet number P e I =ĥÛ /D I 1.7 · 10 −1 -6 · 10 2 2 · 10 1
The function u2 satisfies [S4a] with the no-slip boundary conditions u2 = 0 at y = ±1, and the solution is given by u2 = − [∂xp(x)] uP (y). Consequently, the following linear combination of uP (y) and u ∞ d (y, z) solves [S4]: To determine the pressure gradient term in [S8], we first pose a piecewise pressure field Integrating the continuity equation where 2P/3 and Q ∞ d (ϕz, P ) are the flow rates given by uP (y) and u ∞ d (y, z), respectively. Since Q2 must be constant in x, the pressure gradient ∂xp2 is necessarily independent of x as well. Equating Q1 = Q2 yields a relationship between the two pressure gradients, [S11] The last condition that must be satisfied by the solution is the fixed pressure drop across the domain given by [S1o]. The nondimensional version of this equation, taking α = 1 and β = 0 in [S1o], leads to p(x) = p(x + 1) + 1. This equation can be made specific to x = −1/2 and recast into an integral equation for the gradient [S12] Substituting [S11] into [S12], we arrive at which, after defining q ∞ d = 3Q ∞ d /(2P ), can finally be substituted into [S8] to produce the closed form solution for the flow field Eq. 1 in the main text. The term ⟨γ⟩ in [S13] represents the average value of γ(x) at the interface, i.e. ⟨γ⟩ = 1 Once the leading-order velocity field [S8] is fully determined from known parameters, the relevant quantities characterizing the performance of the SHS can be readily obtained. The local centerline slip velocity uIc = u(x, y = −1, z = 0) is . With the additional assumption of a uniform shear stress γ(x) = ⟨γ⟩, justified in Section Scaling theory for surfactant transport, [S14] further simplifies to where we define u clean Ic (ϕx, ϕz, P ) as the centerline slip velocity for the finite-grating clean case (i.e. γ(x) = ⟨γ⟩ = 0). We show in Section Scaling theory for surfactant transport that we can model the average shear stress ⟨γ⟩ using a scaling analysis of the surfactant transport equations, leading to Eq. S28. Introducing [S28] into [S15], we arrive at Eq. 3 of the main text. Another common, global measure of SHS performance is the increase in flow rate with respect to that of a Poiseuille flow. Our theory predicts where we again introduce Q clean d (ϕx, ϕz, P ) as the increase in flow rate for the finite-grating, clean problem. Perhaps the most common global quantity sought in theoretical SHS studies is the effective slip length λe. Here λe is defined as the quantity that yields the same increase Q d in flow rate if the mixed boundary conditions on y = −1 are replaced with a uniform Navier-slip condition u = λe∂yu. Such a flow yields a solution u λe (y) = uP (y) + λe(1 − y)/(2 + λe) and thus an increase in flow rate of 2P λe/(2 + λe) which, when equated to Q d , yields an expression for the slip length [S17] Infinite-grating problem. The calculation of Q ∞ d and u ∞ Ic for any values of P and ϕz require either the numerical solution of a dual trigonometric series (9) or the solution of nonlinear algebraic equations involving elliptic integrals and elliptic functions (8). We will not provide such level of detail here and instead refer the reader to (8)(9)(10). However, for completeness we provide the asymptotic limits Since the values of Q ∞ d and u ∞ Ic are monotonically increasing with P , the above limit [S19] reveals that u ∞ Ic < 2. In addition, the limits [S18], [S19] can be combined to obtain useful approximations for all P , which are within 12% of the exact values if 0.5 < ϕz < 0.95, and are of course most accurate at large or small P : Scaling theory for surfactant transport Full problem. The analysis of the surfactant transport equations is similar to that in (11), which we recapitulate here in order to clarify the differences between the two-dimensional and three-dimensional cases. The first assumption of our model for [S3e]-[S3k] is that the concentration of surfactant is low enough to ensure a dilute regime, that is k ≪ 1. We expect this assumption to be the case for most situations in which surfactants are not artificially added, for instance, when unwanted contaminants are naturally present in water (as discussed in (11)). Additionally, since the interaction parameter A is typically not large in absolute value, with |A| ≲ 20 (2), it is possible to assume that k|A| ≪ 1 as well.
where the spanwise average across the plastron of a function f (z) is defined as ⟨f ⟩z = 1 ϕz P ϕz P/2 −ϕz P/2 f (z) dz, and where the terms in [S3f] associated with derivatives in z vanish due to the no-slip (w = 0) and no-flux (∂zΓ = 0) boundary conditions at the edges z = ±ϕzP/2 of the plastron. If [S21] is further integrated from x = −ϕx/2 to x = ϕx/2 and boundary conditions u = 0 and ∂xΓ = 0 are applied at x = ±ϕx/2, we have that and thus, by virtue of the mean value theorem, an equilibrium condition ⟨cI ⟩z = ⟨Γ⟩z must occur at some coordinate along the interface, which we call x0. Downstream from x0, the flow advection promotes the accumulation of interfacial surfactant, which in turn triggers a net desorption flux and an increase in bulk surfactant with respect to the background level. Upstream from x0, the situation is the opposite, with a deficit of Γ and cI with respect to the equilibrium values and a net adsorption flux. Figure 1C, D depicts this physical scenario with the two distinct regions along the interface. The second main assumption is to consider the interfacial concentration Γ as approximately linear. In this case, [S22] implies that the equilibrium point must be approximately at the center of the interface (i.e. x0 ≈ 0), and thus the bulk concentration at x0 is approximately the background concentration and we have ⟨c⟩z(x0) = ⟨Γ⟩z(x0) ≈ 1. Consequently, this assumption allows to scale the concentrations at both ends of the interface x = ±ϕx/2 as with ∆c and ∆Γ the characteristic variation of the concentrations (see Fig. 1C, D). Additionally, note that an approximately linear Γ also implies, from [S3h], that the Marangoni shear at the interface is approximately constant (i.e. γ(x) ≈ ⟨γ⟩). This assumption is expected to hold as long as the flow is not in the so-called stagnant cap regime (11), characterized by a strongly nonuniform interfacial concentration. The stagnant cap regime is reached when advection at the interface overcomes both diffusion and kinetic effects (12), that is, when εP eI ≫ 1 and either Bi/ε ≪ 1 or Da ≫ 1 (11). Given the typical parameter values in small-scale flows like the ones considered in this study (see Section Estimate of surfactant parameters and Table SII), we conclude that for long gratings we have εP eI ≲ 1, which justifies the assumption of an approximately linear Γ. Furthermore, we verify a posteriori that our simulation results show an approximately linear profile of Γ (see Section Finite-element simulations). Using these two assumptions, it is possible to use scaling arguments on [S3] to obtain an expression for ⟨γ⟩ as a function of the nondimensional groups of the problem. We start by scaling the terms in [S3h] as ∂yu| I ∼ ⟨γ⟩ and ∂xΓ ∼ ∆Γ/ϕx, leading to ∆Γ ∼ ϕx⟨γ⟩ εkM a . [S24] Next, we evaluate the terms in [S3j] at the interface ends x = ±ϕx/2. We take ∂yc| I ∼ [1 − (1 ± ∆cI )]/δ ∼ ∓∆cI /δ and (cI − Γ) ∼ [1 ± ∆cI − (1 ± ∆Γ)] ∼ ±(∆cI − ∆Γ), where δ =δ/ĥ is the characteristic boundary layer thickness of the bulk concentration (Fig. 1C). We arrive at Eq. S21 is integrated from x = −ϕx/2 to x = x0, leading to whose terms we scale as ⟨uI Γ⟩z(x0) ∼ uIc, ∂x⟨Γ⟩z(x0) ∼ ∆Γ/ϕx, and . [S28] Equation S28 can now be introduced in [S14] to obtain the formula for the slip velocity [3] in the main text. Similarly, combining [S28] with [S16] and [S17], expressions for the increase in flow rate and effective slip length can be reached.
The only yet undetermined part of the model is an expression for the boundary layer thickness δ, which we seek through scaling of the conservation law for the bulk surfactant [S3e]. In situations with εP e ≫ 1, streamwise advection must balance wall-normal diffusion u∂xc ∼ 1 εP e ∂yyc, which is only possible if c varies over a small length scale δ ≪ 1 (12). We take ∂xc ∼ ∆cI /ϕx, ∂yyc ∼ ∆cI /δ 2 and the velocity inside the boundary layer as u ∼ uIc + ⟨γ⟩δ. In the case of an interface close to immobilization, i.e. uIc ∼ 0 and ⟨γ⟩ ∼ 1, these scalings indicate that δ ∼ (P e/g) −1/3 when εP e ≫ 1. In the opposite case of εP e ≪ 1, Eq. S3e is dominated by diffusion, and thus the characteristic length scale of variation of c in the wall-normal direction is the whole half height of the domain, implying δ ∼ 1. We choose to satisfy these two extremes, with a3 and a4 empirical parameters. It is also possible to obtain a similar expression with an exponent of −1/2 instead, by assuming that the boundary layer is essentially shear-free (i.e. uIc ∼ 1 and ⟨γ⟩ ∼ 0). In practice, the overall value of quantities like uIc are only weakly dependent on the specific functional form of δ, so we only consider the expression [S29]. Additionally, in the case of interest of long gratings ε ≪ 1 in small-scale flows we typically have P e/g ≲ 1 (Section Estimate of surfactant parameters) and thus the boundary layer thickness is approximately independent of P e or g. with M ains = k M a P eI = nsRTΓ0ĥ/(μDI ), and ains another empirical parameter.

Finite-element simulations
We solve numerically the full governing equations [S1a]-[S1t] of the problem in dimensional form, performing a total of 155 simulations. The objectives are to (i) determine the values of the empirical parameters a1, a2, a3 and a4 in our model (Section Scaling theory for surfactant transport), (ii) confirm the modeling assumptions, for the interfacial concentration Γ, of an approximately constant profile in the spanwise direction, and of a linear profile in the streamwise direction (Section Scaling theory for surfactant transport), and (iii) compare the theory to simulations of realistic microchannels in conditions representative of our experiments (Section Experimental methods). We implemented the three-dimensional simulations using the finite-element software COMSOL Multiphysics 5.5 ® . The simulation domain is one half of the SHS unit cell depicted in Fig. 1, withẑ spanning only betweenẑ = 0 andẑ =P /2 due to the spanwise symmetry of the solution. The volume is meshed with tetrahedral elements, concentrating the finest regions around the upstream and downstream edges of the interfacex = ±ϕxL/2 since it is in those areas where the most abrupt variations of the solution occur (see Fig. S1). Across all the simulations, the minimum element size (understood as the diameter of a sphere circumscribing the smallest element) is set to 1.5 · 10 −9 m.
The solution of the governing equations is achieved with a combination of the Creeping Flow module for the flow equations [S1a]-[S1d] and the Dilute Species Transport module for the transport of bulk surfactant [S1e]. The conservation law for the interfacial surfactant [S1f] is implemented through a General Form Boundary PDE, using [S1j] as the source term. The Marangoni boundary conditions [S1h] and [S3i] are enforced through a Weak Contribution constraint, as is the condition that fixes the mean bulk concentration [S1l].
The system of nonlinear equations is solved through a Newton-type iterative method using the PARDISO direct solver for the linear system at each iteration. We set a relative tolerance of 10 −5 as a convergence criterion for the solution, which is  satisfied by all of our simulations. The pressure, bulk concentration and interfacial concentration are discretized using linear elements, and the velocity field uses either quadratic or linear elements, depending on the computational demands of each simulation.
We vary the problem parameters to ensure that each of the distinct terms that are pre-multiplied by an empirical coefficient in [S28] changes over a few orders of magnitude. The ranges of variation of each dimensional quantity in the simulations, as well as of the corresponding nondimensional numbers, is indicated in Tables SI and SII. A small number of simulations were chosen with the same parameters as those estimated in the experiments, in order to achieve a direct comparison (see Fig. 3 in the main text). However, due to constraints in computational power, the value of the grating lengthĝ was smaller than that of the microfluidic devices.
The parameters a1, a2, a3 and a4 are obtained through least-squares fitting using the MATLAB function lsqnonlin. We define the error as the total sum of squares of the difference between the centerline slip velocities computed in the simulations and those predicted by the theory [3], i.e. ERR = (u theory Ic − u sim Ic ) 2 . We find a1 = 0.345, a2 = 0.275, a3 = 5.581, and a4 = 3.922. As illustrated in Fig. 2A, the agreement between simulations and theory is excellent over more than four orders of magnitude in the slip velocity.

Experimental methods
The experimental setup is centered around the custom-built microfluidic device depicted in Fig. 3A-C. The chips are made by casting PDMS (Sylgard 184) with a 1:10 mass ratio of elastomer-to-curing agent, using a master mold fabricated by two-layer photolithography. The photoresist used for the mold is SU-8 (Microchem SU-8 3050 for the first layer, corresponding to the main channel, and Microchem SU-3025 for the second layer, corresponding to the gratings). The flat PDMS substrate is determined to be hydrophobic, with static contact angles of water droplets measured to be higher than 90°. The textured PDMS exhibits a further increase of the static contact angle, with air pockets within the texture that are observable under the microscope, demonstrating that the PDMS textured substrate is in a Cassie-Baxter superhydrophobic state. Measured values of static contact angles are consistent with previous measurements for untreated PDMS (13).
The chip is bonded to a 0.1 mm-thick glass coverslip (Bellco Glass 1916-25075) through untreated adhesion, and a 40X water objective is used to image the interior of the channels through the coverslip using a confocal microscope (Leica SP8 Resonant Scanning). The device is placed inside a stage top chamber (Okolab H101-K-FRAME) that ensures precise control of temperature, which we set toT = 296 K. The fluid is initially contained in a glass syringe (Hamilton Gastight), and driven by a syringe pump (KD Legato 111) at a constant flow rate through plastic tubing (Tygon S3) into and out of the microfluidic channel. We use the barrel of a plastic syringe (BD Luer-Lok) as an outlet reservoir open to the room. To ensure that the air-water interface in the observed channels remains flat and that plastron curvature effects can be safely neglected, the magnitude of the pressure inside the channel is adjusted by varying the height of the outlet reservoir, which is mounted on a vertical translation stage (Thorlabs VAP10). The maximum deflection of the interface at the centerline, relative to the edges, is estimated to be less than ±1 µm.
Due to the extreme difficulty of removing all traces of surface-active contaminants, even in controlled experimental conditions (4), we do not attempt an exhaustive cleaning protocol with that aim. Nevertheless, we follow standard cleaning procedures on all syringes and tubing, ensuring that they are rinsed with 18 MΩ cm DI water with at least twice their volume before they are used. In addition, we follow a specific cleaning protocol for the µ-PIV particles (ThermoFisher FluoSpheres carboxylate 0.5-µm diameter), since they typically contain surfactants to prevent particle aggregation (14). We use a centrifuge (Eppendorf 5418) to separate the beads from the buffer solution, which is discarded and replenished with clean 18 MΩ cm DI water, and we repeat the process three times. These cleaning procedures ensure that the traces of surfactants responsible for the non-negligible Marangoni stresses that we observe in the experiments are the result of contamination that would naturally occur in typical small-scale flows through microfluidic devices, and not as a byproduct of the specific experimental methods used in this study.
The µ-PIV analysis is performed using the open-source MATLAB toolbox PIVlab. The acquisition window has an approximate size of 125 µm × 125 µm, which is sufficient to cover the span of two pitches (see Fig. 3A,B,C), which are imaged around the center of the grating in the streamwise direction (i.e. x = 0). We image the motion of the µ-PIV particles across time intervals between 20 s and 60 s, at different distances from the interface, with frame rates between 20 fps and 28 fps.The velocity field is averaged in time and in the streamwise (x) direction to obtain the velocity profiles depicted in Fig. 3D. To calculate the value ofûIc, we perform a linear least-squares fit, typically using between three and five velocity profiles to obtain an extrapolated slip velocity, from which we extract its value at z = 0. This linear fit is performed in MATLAB with a script that takes into account the propagation of uncertainties in the distance ∆ŷ from the interface, as well as the uncertainty inû inherent to the measurement and the averaging in the x direction.

Estimate of surfactant parameters
The main challenge in comparing the experimental measurements ofûIc to the predictions from our model is the absence of information regarding the type and amount of surfactant present in the channels. Some parameters in the problem are known from the experimental conditions, and hence we fix those as ns = 2 (3),R = 8.314 J mol −1 K −1 ,T = 296 K and µ = 8.9 · 10 −4 kg m −1 s −1 (15). Others can be accurately estimated, since most surfactants have diffusivities (bothD and DI ) bounded between 5 · 10 −10 and 5 · 10 −9 m 2 s −1 , and have values ofΓm between 1 · 10 −6 and 5 · 10 −5 mol m −2 (2). We thus use as a reference surfactant the well-studied sodium dodecyl sulfate (SDS), settingD =DI = 7 · 10 −10 m 2 s −1 and Γm = 3.9 · 10 −6 mol m −2 . However, the values of the rate constantsκa andκ d can vary significantly across surfactants (2), and they cannot be assumed a priori. In addition, we must estimate the background bulk concentrationĉ0.
Here we describe a strategy to estimate these quantities. The value of k =κaĉ0/κa has an upper bound, since we expect k ≪ 1 not only because this is the case when surfactants are not artificially added, but also because values of k ≳ 1 would significantly decrease the mean surface tension, leading to a rapid plastron collapse, which is not observed in experiments. We therefore choose the bound k < kmax = 1 · 10 −1 , which ensures that k remains at least one order of magnitude smaller than 1 and that the absolute surface tension decrease ∆σ is small compared to the clean surface tension valueσ0 = 7.2 · 10 −2 N m −1 . Indeed, an estimation using an equation of state derived from the Langmuir isotherm (11) yields ∆σ/σ0 = nsRTΓm ln(1 + kmax)/σ0 ≈ 0.025.

Discussion of the mobilization length
Figures S2A, B display the same data collapse as Figs. 4B, C, but include the data from all 155 simulations. The yellow triangles indicate simulations with a thin boundary layer δ < 1, which can only occur if P e/g ≫ 1. In practical applications with long gratings g ≫ 1, such a regime would require extremely large Péclet numbers P e ≫ g ≫ 1 that are outside the scope of this study. The green diamonds denote cases in which the grating length is below the modified depletion length g < L mod d , which are of little practical relevance. This is due to the fact that, for typical surfactants, we find Lm ≫ L mod d (see Section Problem parameters), and thus Lm is ultimately the lengthscale that g needs to overcome for a SHS grating to display significant slip. Figure S2C shows the same data collapse when the parameter (u clean Ic ) 1/2 is omitted from the normalization of g. Note that u clean Ic depends only on the SHS texture geometry. Since, for SHS capable of drag reduction, one needs u clean Ic of order one (see Section Infinite-grating problem), the collapse in Fig. S2C is only marginally worse than in Fig. S2B. Furthermore, a good approximation to the curve shown in Fig. S2B is found by using Eq. 11 with u clean Ic = 1/2, as shown in Fig. S2C. For convenience, here we provide a detailed discussion of the assumptions under which the surfactant problem is only a function of g/Lm: Eq. 11 with Eq. 11  simulations with a thin boundary layer δ < 1, which is not representative of long gratings, whereas the green diamonds denote cases where the grating length is below the modified depletion length g < L mod d , which are of little practical relevance, as explained in the text. (C) Data collapse when u clean Ic is ignored in the normalization of g. The curve shows Eq. 11 with a constant value of u clean Ic = 1/2.
• Dilute surfactant regime, k ≪ 1: This is consistent with environmental traces of surfactant, intrinsic to the microfluidic system, without artificially added surfactant. Furthermore, for values of k ∼ O(1) the plastron would risk collapse due the decrease in mean surface tension. Based on estimates from experimental data (4), we take bounds 1 · 10 −4 ≲ k ≲ 1 · 10 −1 (see Section Estimate of surfactant parameters), consistently also with values found in the literature for engineered (16) and natural settings (17).
• Kinetics are fast compared to diffusion, Da ≳ 1: For microchannel half-heights and diffusivities as in the points above, and finding the approximate boundsκaΓm ∈ [10 −5 , 10 −3 ] m s −1 in the literature (satisfied for all surfactants except for two in Table 3 of (2)), one has that Da =κaΓmĥ/D ∈ [1, 10 3 ]. Only forĥ considerably smaller than tens of microns one would have very small values of Da.
• Uniform shear stress, γ(x) ≈ ⟨γ⟩: Deviations from this assumption are associated with the 'stagnant cap regime', which requires a dominant interface advection term, with negligible kinetics and diffusion (12). The stagnant cap regime is difficult to achieve for long gratings, since Bi/ε ≳ 1 will be easily satisfied, thereby making interface kinetic fluxes significant in the interfacial transport equation (Eq. 2b in the main text). We assume that the channel half-height is in the rangeĥ ∈ [10 −4 , 10 −3 ] m, and that the flow velocity is in the rangeÛ ∈ [10 −5 , 10 −3 ] m s −1 . Based on the review of physicochemical properties of surfactants in Ref.
(2) (see their Table 3), we assumeκ d is in the rangeκ d ∈ [10 −1 , 10 3 ] s −1 , which covers all but one of the surfactants that they describe. With these bounds, we find Bi =κ dĥ /Û ∈ [10 −2 , 10 5 ]. Even with this wide range of parameters, the Biot number remains high enough such that Bi/ε will be at least of order one for gratings that are practical for drag reduction, where ε =ĥ/L ≈ĥ/ĝ ∼ 10 −2 or smaller. The uniform shear stress assumption is also supported by our simulations, which span a wide range of parameters (see Table SII).