Roughness on liquid-infused surfaces induced by capillary waves

Liquid-infused surfaces (LIS) are a promising technique for reducing friction, fouling and icing in both laminar and turbulent flows. Previous work has demonstrated that these surfaces are susceptible to shear-driven drainage. Here, we report a different failure mode using direct numerical simulations of a turbulent channel flow with liquid-infused longitudinal grooves. When the liquid-liquid surface tension is small and/or grooves are wide, we observe traveling-wave perturbations on the interface with amplitudes larger than the viscous sublayer of the turbulent flow. These capillary waves induce a roughness effect that increases drag. The generation mechanism of these waves is explained using the theory developed by Miles for gravity waves. Energy is transferred from the turbulent flow to the LIS provided that there is a negative curvature of the mean flow at the critical layer. Given the groove width, the Weber number and an estimate of the friction Reynolds number, we provide relations to determine whether a LIS behaves as a smooth or rough surface in a turbulent flow.


Introduction
A protective and functional surface coating can be created by lubricating a textured surface with an appropriate liquid. The drag-reducing properties of these liquid-infused surfaces (LIS) have been explored recently both numerically Cartagena et al. 2018;Arenas et al. 2019) and experimentally (Van Buren & Smits 2017;Fu et al. 2019). LIS can also prevent fouling (Epstein et al. 2012), corrosion (Wang et al. 2015) and ice formation (Kim et al. 2012).
The drag-reducing capabilities of LIS for liquid flows is often compared to those of superhydrophobic surfaces (SHS), where air is used as the infused medium. The low viscosity of air is beneficial for drag reduction, but the use of SHS in turbulent applications is restricted by mass diffusion (Ling et al. 2017) and instability of the gas pockets (Seo et al. 2018). For LIS, the mass diffusion is negligible if the liquids are immiscible, and LIS are not susceptible to failure due to hydrostatic pressure (Wong et al. 2011). However, also for LIS, the stability of the interface depends on the texture's geometry, the surface tension between the two liquids and the contact angle at the liquid-liquid-solid interface. In particular, these surfaces may experience shear-driven drainage of the infused liquid, but this can be mitigated, for example with chemical patterning (Wexler et al. 2015;Fu et al. 2019).
In this paper, we show that capillary motion of the liquid-liquid interface may drastically lower the drag-reducing performance of LIS. In the present study, the surface texture is fixed to longitudinal (streamwise-aligned) grooves. We use direct numerical simulations of a LIS in a liquid turbulent channel flow for frictional Reynolds numbers around Re τ ≈ 180. The employed volume-of-fluid (VOF) framework allows for large interface deformations (low surface tension) and moving contact lines. When the liquid-liquid surface tension is small and/or grooves are wide, we find travelling-wave perturbations on the interface with amplitudes larger than the viscous sublayer of the turbulent flow (a + ≈ 5 − 8). These capillary waves induce a roughness effect and increase friction drag.
The detrimental capillary waves develop for viscosity and density ratios of one, which excludes interface instability mechanisms driven by density and viscosity stratification (Boomkamp & Miesen 1996). Instead, we find that the linear instability can be described by the theory developed by  in the context of two-dimensional gravity waves. This inviscid instability is due to an energy transfer from the external flow to the waves that makes them grow in time at an exponential rate. The energy transfer can occur if (i) there is a negative curvature of the mean velocity profile where it equals the phase speed of a wave, i.e. at the critical layer and, (ii) the critical layer is not too far away from the surface, so that velocity fluctuations due to the wave (dispersive stresses) are non-zero.
The existence of energy transfer is not sufficient for failure of LIS, however. The interface fluctuations also need to grow sufficiently fast to reach large amplitudes that induce roughness effects. Figure 1a, which summarises our main contribution, shows three domains, namely rough, smooth and transitional (grey) in a plane spanned by a groove width (w + ) and a Weber number (We + ), both normalised with the viscous length scale. This design map is obtained from the critical-layer theory and provides a means to design LIS that can be predicted to achieve a balance between performance (large w + ) and stability (smooth domain).

Numerical methods and configuration
We consider a fully developed turbulent open channel flow. The flow domain, shown schematically in fig. 1b, has the size (L x , L y , L z ) = (6.4h, h+k, 3.2h), where x, y and z correspond to the streamwise, wall-normal and spanwise directions, respectively, and h is the half-channel height. At the top boundary we impose a free-slip (symmetry) boundary condition (BC). Periodic boundary conditions are imposed in the streamwise and spanwise directions. The streamwise-aligned grooves at the bottom wall have a height k, a width w and square cross-section, k = w. The fluid-solid ratio is set to 0.5. The infused and external fluids have the same density ρ i = ρ ∞ = ρ, but different viscosities (µ i and µ ∞ ). We have used grooves of width w +0 = 18. Throughout this paper, +0 refers to normalisation using the friction velocity of a regular smooth wall (nominal wall units). A single superscript + refers to normalisation using the friction velocity of each individual case. The corresponding viscous length scale is δ ν = µ ∞ /(ρu τ ), where u τ is the friction velocity.
We impose a constant mass flow rate through a uniform pressure gradient over 0 < y < h, where y = 0 corresponds to the crest of the texture so that Re b = ρhU b /µ ∞ = 2820, giving Re τ = ρhu τ /µ ∞ ≈ 180 (with the pressure gradient implemented as a volume force). Here, Re b and Re τ are Reynolds numbers based on bulk velocity U b and friction velocity u τ , respectively.
Our simulations allow for a moving liquid-liquid-solid contact line with a dynamic contact angle different from the static value, which is θ = 45°with respect to the infused liquid. Our method also allows for interface deformation, which is typically quantified by the Weber number, defined as We = ρU 2 b h/γ, where γ is the surface tension. We have simulated LIS for We = 100, 150 and 200 and viscosity ratios µ i /µ ∞ = 0.5, 1 and 2. The Weber number in wall units is We + = ρu 2 τ δ ν /γ = µ ∞ u τ /γ. It can be noted that the Weber number based on the friction velocity and the width of the grooves is We + w + .
The numerical configuration described above corresponds to an infused liquid consisting of some alkane (with dynamic viscosities similar to that of water (Van Buren & Smits 2017)), a water channel with h = 0.5 cm and U b = 1 m/s. This results in Re b = 5000, which is close to the value in our simulations. A typical surface tension γ = 50 mN/m then results in We = 100.
The code used for the simulations is based on the PArallel, Robust, Interface Simulator (PARIS), which employs a VOF method for the multiphase description (Aniszewski et al. 2021;Arrufat et al. 2021). The cited papers also include additional test cases and validations. Height functions are used for curvature calculation for the surface tension. The interface is advected in the manner suggested by Weymouth & Yue (2010) at each substep. This advection scheme conserves the volume of both liquids to a high accuracy. At solid surfaces, a contact angle is imposed by using the height functions and a dynamic contact angle model for VOF based on hydrodynamic theory . Details are given in sec. S1 of the supplementary material (SM). Finally, the grid size is (N x , N y , N z ) = (256, 640, 1024), with constant grid spacing in each direction. The flow was converged after 600h/U b , and statistics were collected over a time of 500h/U b .

Results
3.1. Dependence of drag on Weber number Figure 2(a) shows an instantaneous snapshot of the liquid-liquid interface, viewed from the top, for We = 100 and µ i = µ ∞ . There are oscillations on the interface, due to the finite surface tension, but they remain small. The deformation of the interface increases with We, however. For We = 200, significantly larger waves develop on the interface, as shown in fig. 2 The consequences of the waves on the overlying flow can be quantified by the drag reduction where c f = 2τ w /(ρU 2 b ) is the friction coefficient and τ w is the total stress at the crest plane of the surface (computed from the pressure gradient). Here, c 0 f is the coefficient of a regular smooth wall at y = 0. For We = 100 and We = 200, we obtained DR = 0.09 and DR = −0.04, respectively. In other words, the capillary waves observed at We = 200 increase frictional drag compared to a smooth and homogeneous surface and have therefore induced failure of the LIS. Figure 3a shows DR for We = {100, 150, 200} and viscosity ratios µ i /µ ∞ = {0.5, 1, 2} as a function of the apparent slip length, b +0 . The slip length b is the distance at which the mean velocity would be zero if linearly extrapolated at the crests of the surface. It is largely unaffected by changes in We, but it increases with decreasing viscosity ratio. In fact, the slip lengths extracted from our numerical simulations with w + ≈ 18 are well approximated by the slip lengths obtained by solving the Stokes equations for a periodic array of grooves exposed to unit shear ) (black line in fig. 3b). A similar agreement was observed for smaller grooves (w + ≈ 9) in turbulent flows by Fu et al. (2017); Arenas et al.

(2019).
When the interface is perfectly flat (We = 0), the drag reduction can be related to slip length as

2)
Roughness on liquid-infused surfaces due to capillary waves where U +0 b is the bulk velocity in nominal wall units. This relation -shown in 3a (black line) -can be obtained by neglecting changes in the Reynolds shear stress above a smooth wall (Rastegari & Akhavan 2015). We observe from fig. 3a that, for We = 100 (blue) and We = 150 (red), there is a drag reduction (DR > 0) for all three viscosity ratios. Moreover, the deviations from (3.2) are small, confirming that the drag reduction mechanism is indeed slippage. These small deviations are due to change of Reynolds shear stress. In contrast, the deviations from (3.2) are significant for We = 200 (yellow), where we observe a drag increase (DR < 0) for µ i /µ ∞ = 1 and 2 and a DR close to zero for µ i /µ ∞ = 0.5. The corresponding mean velocity and velocity fluctuations reflect the increase in drag, and these are described in the SM (sec. S2).
The waves formed on the interface at We = 200 are sufficiently large to cause roughness effects. Interface height profiles at different times are shown in fig. 4b, together with amplitudes of a wave in its initial stage ( fig. 4c, yellow). The wave amplitude is defined as the height of the local maximum of the wave. Wave amplitudes of a + > 5 are observed and these extend outside the viscous sublayer, indicating that the surface is transitionally rough. The amplitude grows initially at an exponential rate, before it levels off. In contrast, interface fluctuations for We = 100 have small amplitudes (a + < 1) (fig. 4a) and show a significantly smaller growth rate ( fig. 4c, blue). The exponential growth rate ( fig. 4c, dashed) is an indication of a linear instability. In the next section, we provide evidence of a critical-layer instability , where energy is transferred to the wave perturbation from the turbulent flow.

Conditions for phase speed and growth rate of capillary waves
We assume a small perturbation on the liquid-liquid interface of the form where z = 0 is located in the centre of the groove. As illustrated in fig. 1c, η is the height of the interface, A is the initial wave amplitude, c is a complex wave speed and t is time. Moreover, k x = 2π/λ x and k z = 2π/λ z are streamwise and spanwise wavenumbers, respectively. The spanwise wavelength can have a maximum value of λ z = 2w, due to the finite width of the grooves. This value can be seen to dominate in the snapshots of fig. 2, as most waves only have one crest or one trough in the spanwise direction. We also observe from fig. 2 that streamwise wavelengths are generally similar to, or larger than, λ z . This threedimensionality implies that both spanwise and streamwise curvatures contribute to the capillary pressure of a wave, Here, k = k 2 x + k 2 z and p + 0 (p − 0 ) is the pressure above (below) the interface. Next, we consider a wall-normal velocity disturbance v(x, y, z, t) on the turbulent mean flow U (y) with the same waveform as η. If we neglect viscous and nonlinear effects, the amplitudev(y) is governed by the Rayleigh equation (SM, sec. S3 A), where denotes a derivative with respect to y. The velocity perturbation must vanish at infinity and satisfy the kinematic condition at the interface, v/(U −c) = ik x η. The equation for the pressure amplitude,p(y), corresponding to eq. (3.5) iŝ Our aim is to find an approximate solution to the equations (3.4-3.6) in order to determine the phase speed (c) (real part) and growth rate (k x c) (imaginary part) of the interface perturbation (3.3).  formulated a similar set of equations for describing wind-induced water waves, where gravity -instead of capillarity -balances fluid pressure. He suggested the following approximate solution for v: This expression, which satisfies the boundary conditions at y → ∞, implies that 1/k is the relevant length scale over which v decreases. The assumption of exponential decay can also be used in the grooves: Roughness on liquid-infused surfaces due to capillary waves Here, we have assumed that the grooves are sufficiently deep such that the velocity perturbation is nearly zero at the bottom of the groove. With a depth w, kw > k z w (2π/(2w))w = π, and, since e −π 1, the assumption is valid for our configuration. We have also neglected U and its derivative inside the groove. Inserting (3.8) into (3.6) results in the following expression for the pressure immediately below the interface (y → 0 − ): Similarly, by inserting (3.7) into (3.6) the pressure just above the interface is where α and β are real constants and U 1 is an arbitrary reference velocity. It is shown in SM (sec. S3 C) that the parameter α can be decomposed into two parts, α = α 1 + α 2 , where α 1 corresponds to eq. (3.9) and α 2 incorporates the remaining contributions from the slip velocity and the shear. Using this decomposition and inserting (3.9) and (3.10) into (3.4), we obtain (SM, sec. S3 D), Here c w = γk 3 /(2ρk 2 x ) is the free phase speed, i.e. the speed of a capillary wave without forcing from the overlying flow.
The free phase speed is shown in fig. 5a (in nominal wall units) as a function of λ x /w for We = {100, 150, 200}. We note that two-dimensional capillary waves (k z = 0) have a phase speed that monotonically decreases with λ x . However, for LIS, there is a minimum phase speed due to the finite spanwise wavelength. This minimum is approximately (for the analytical expression see SM, sec. S3 D) For We = 200 (and µ i /µ ∞ = 1)¸c + w,min = 7.0, which is slightly lower than the phase speed of the wave shown in fig. 4b. One may use c + w,min as a lower bound of the actual phase speed of capillary waves on LIS.
We now turn our attention to the imaginary part of (3.11) to approximate the growth rate of the instability. As shown by  -and repeated in the SM (sec S3 E) -one may integrate the Rayleigh equation and evaluate the pressure equation (3.6) at the interface to find where the subscript c denotes values at the position of the critical layer y c . In order for an infinitesimal wave to have a positive growth rate, i.e. β > 0, a first requirement is that U c < 0, i.e. negative curvature at the critical layer. This is satisfied if the critical layer is outside of the viscous sublayer. A second requirement for β > 0 is that v c in eq. (3.13) is non-zero at the critical layer. The approximate solution of v in eq. (3.7) implies, however, that v is zero at the critical layer. As shown in SM (sec. 3 E), one may transform the condition for positive growth rate to an integral form to estimate v in the vicinity of the critical layer. By further assuming a logarithmic mean velocity profile and setting the reference velocity to U 1 = u τ /κ (where κ is the von Kármán constant), one may evaluate the expression for β (as a function of ky c ), and obtain what is shown in fig. 5b.
We observe from fig. 5b that, when ky c > 3/2, then β < 0.02, which results in very slow-growing waves, whereas when ky c < 1/2, we have β > 0.8, resulting in a factor 40 or more faster growth. The grey region in fig. 5b marks the range 1/2 < ky c < 3/2 where there is a transition from low to high growth rates. Now, since k > k z and the upper limit of λ z is 2w, we may formulate bounds for the position of the critical layer. When the growth rate can be expected to be significant, in contrast to when 15) for which the growth is negligible. Equations (3.12), (3.14) and (3.15) provide relationships between We + , w + and c + w,min , y + c,max that can be confirmed by our simulations. Equation (3.12) states that a large We + and/or w + give a small phase speed. This is observed qualitatively by following the travelling waves on the interface in fig. 4(a,b). More quantitatively, the space-time correlations of the interface height for We = 100 and We = 200 give c + = 15.1 and c + = 10.5, respectively (SM fig. S8).
Compared to We = 100, the lower phase speed for We = 200 results in a lower position of the critical layer. When the height of the critical layer approaches the interface and satisfies eq. (3.14), the growth rate coefficient β of the waves (eq. 3.13) is large. This is confirmed by our simulations, where we observe in fig. 4c that both the growth rate and interface amplitudes are larger for We = 200 compared to We = 100.
We use an inviscid model here to get a tractable analytical solution, and to illustrate the important physics involved. It has been shown that the effect of introducing viscosity on capillary waves with relevant wavenumbers would be a Roughness on liquid-infused surfaces due to capillary waves 9 slight damping (Jeng et al. 1998). However, it is possible that viscosity influences the velocity induced by the waves deep inside the grooves to a higher extent.

Implications for the design of LIS
The conditions (3.12), (3.14) and (3.15) can be used as design criteria for LIS. One may expect a high-performing LIS by choosing a groove width and a surface tension of the infused liquid such that -for relevant friction Reynolds numbers -the design falls within the smooth region of fig. 1a. This region is defined by (We + , w + ), where y + c 1.5w + /π, and thus from eq. (3.15) very small growth rates of capillary waves are predicted. Conversely, the rough region in fig. 1a shows (We + , w + ), where y + c 0.5w + /π, and therefore waves will amplify rapidly. Here, we may expect either a very low-performing LIS or even a drag-increasing LIS due to roughness effects. In between the smooth and rough domains, we show in fig. 1a a transitional region (grey), which corresponds to 0.5w + /π y + c 1.5w + /π. Here, we cannot predict if the resulting waves induce roughness effects using our analytical approach. It should be mentioned that the boundaries of the transitional region in fig. 1a are determined in three steps: (i) given w + , determine y + c from (3.14) (lower boundary) or (3.15) (upper boundary); (ii) given y + c , determine c + (phase speed) from U + (y + c ) = c + , where U + (y) is a turbulent mean profile of a smooth wall; and finally (iii) given c + , assume c + ≈ c +0 w,min and determine We + from (3.12) (or the exact coefficient of (3.12) given in the SM). Figure 1a also shows scaling laws between smooth and rough regions. For small w + (and thus y + c ), we may assume that the critical-layer velocity is U + c ≈ y + c . This is acceptable right above the viscous sublayer where the mean flow has some curvature. Then eq. (3.12) gives that the height of the lowest critical layer is y + c ≈ π/(We + w + ). By assuming that y + c ∼ w + /π, we obtain 16) which is shown with dashed line in fig. 1a. It is observed that this asymptotic relation represent a reasonable scaling law for w + 20. For larger w + , away from the viscous sublayer, we assume U + = (1/κ) log(y + )+ B, where B is a constant. This gives a nonlinear relation 1/( We + w + ) ∼ (1/κ) log(w + /π). (3.17) This curve is shown in fig. 1a with a dashed-dotted line, and provides a scaling of the neutral curve for w + 30. The scaling laws illustrate that, when w + increases, there needs to be rapid decrease of We + to remain in the smooth region. For example, for w + ≈ 70, we need We + ≈ 3 · 10 −4 , which corresponds to We ≈ 10. This is relevant for drag reduction, since the width of the grooves should be maximised for a given surface tension to optimise DR (see fig. 3), but without entering the rough zone in fig. 1a. Note that for a fixed geometry, increasing the flow speed, and thereby u τ , increases both We + and w + , so that the design needs to be made for the largest flow speed to which the surface is exposed.
Finally, in fig. 1a, the values of our numerical simulations are shown with symbols. These also include more extreme Weber numbers, We = 50 and We = 400 (using µ i /µ ∞ = 1), which resulted in a drag reduction of 9.3% and −15%, respectively, confirming the trend of the other simulations. In addition to the  Figure 6: (a) If θ < φ < θ + 90°, the contact line remains pinned according to Gibbs' criterion (grey area). This is illustrated for θ = 45°. If φ is outside this range, the contact line depins, and moves in the direction indicated by the arrows. (b) The PDF of φ from simulations for the pinned cases with w +0 = 18, θ = 45°, µ i /µ ∞ = 0.5 (---), µ i /µ ∞ = 1 (--) and µ i /µ ∞ = 2 (-· -) and the Weber numbers We = 100 (blue) and We = 150 (red). The boundaries of the interval corresponding to a probability of 95% for the widest PDF are also shown (· · · · · ·).
simulations at w +0 = 18, we also show points (circles) for larger grooves of width w +0 = 36 (also using µ i /µ ∞ = 1). For these grooves, there was a drag reduction by 18% for We = 25 (purple circle) and 17% for We = 50 (green circle), whereas for We = 100 (blue circle), the drag reduction was lowered to 2% and we observed large waves. This implies that the growth rate rapidly increases between the last two cases as they fall in the transitional zone in fig. 1a.

Contact line depinning
The capillary waves modify the contact angle between the interface and the wall, and may potentially result in a depinning of the interface from the corners of the ridges. According to Gibbs' criterion, which is a purely geometrical criterion, the interface remains pinned if θ < φ < 90°+ θ, where φ is the angle the interface makes to the inner wall of the groove. The lower limit is the limit for when the contact line moves into the groove, while the upper is the limit for when it moves on top of the ridge (Gibbs 1906). This is illustrated in fig. 6a. In contrast to the cases We = 100 or 150 ( fig. 2), we observed that for We = 200 the interface depinned occasionally due to the waves on the interface.
The measured probability density functions (PDF) of φ for We = 100 and 150 with θ = 45°are plotted in fig. 6b for all three viscosity ratios. Since the contact line was observed to remain pinned for these parameters, the PDF are independent of θ and can be used to predict limits for the contact angle. It can be noted that the PDF for all parameters shown are centred around φ = 90°and that φ is unlikely to reach below 70°or above 110°. The interval between these values corresponds to a probability of more than 95% for the widest PDF. The standard deviation of φ decreases with µ i /µ ∞ and increases with We, as is indicated by the width of the PDF. This is to be expected, since the dissipation rate increases with the viscosity and the restoring force of surface tension becomes weaker with increasing We.
A restoring force for the contact line also comes from mass conservation. If the contact line occasionally does depin into the groove on one position, it will be raised elsewhere. Based on this observation and the statistics in fig. 6b, depinning is not the main failure mode of LIS for the geometry chosen in this study. Depinning is, however, expected to be important for a LIS with grooves of finite length.

Conclusions
We have explored the behaviour of LIS in a turbulent channel flow with square longitudinal grooves for Re τ ≈ 180. By allowing the interface and the contact line to move, we could investigate the unconstrained motion of the interface. For a fixed groove width, we found a rapid increase in drag of LIS above a certain Weber number due to the appearance of large capillary waves. The generation mechanism of these waves was elucidated using the theory developed by . The limit for when these waves act as roughness is set by the width of the grooves w + and the Weber number We + , as illustrated in fig. 1a. It should also be noted that these non-dimensional parameters depend on the flow speed. Using an analytical analysis, we have provided scaling laws and design criteria for robust drag-reducing LIS. Specifically, the relations show how to achieve a balance between large groove widths (enhancing drag reduction) and high surface tension of the infused liquid (enhancing stability) for different flow speeds.

S1. DETAILS OF NUMERICAL METHODS
In this section, we describe the schemes of the numerical code that was used to solve the two-component flow over and inside the textured surfaces. For momentum convection, a second-order central scheme was used, whereas a third-order Runge-Kutta scheme was used for time integration [1]. The equation for the pressure was solved with a fast Fourier transform (FFT) solver (FFTW) and handling domain decomposition with the 2DECOMP&FFT library. To describe solids, we used a grid-aligned immersed-boundary method (IBM) [2]. With this method, the edge of the solid is located at the edge of the cells. For a staggered grid, this means that the top of the solid ridges are at the location of the wall-normal velocity nodes.
To impose a contact angle at the liquid-liquid-solid contact line, we specified the height function of the first ghost layer [3]. Since the IBM is grid-aligned, the implementation of the contact angle is the same on an the immersed boundaries and the domain boundaries. We used a dynamic model of the contact angle based on hydrodynamic theory [4] -adapted to VOF [5] -, together with a no-slip velocity condition. As shown in [5], the dynamic contact angle also improves grid independence. The method is explained further in sec. S1 A. It should be noted that the ridge corners' ghost cells were set to always be interface cells, and that the contact angle was imposed when the interface moved to adjacent cells. Depinning occured for the higher Weber numbers.
The turbulent flow was validated by comparing the mean flow and the r.m.s. velocities for a full channel with smooth walls to data from ref. [6]. The results are shown in fig. S1, having for the mean velocity a deviation of 0.8% at the channel center. Shown are also statistics of a smooth open channel, with only small deviations for y + < 100. The description of the streamwise liquid filled grooves was validated against the analytical expressions of Schönecker et al. [7] for the slip length in a laminar flow. These expressions have recently also been used as a reference for turbulent data of LIS simulations [8]. For these tests, we use a computational box corresponding to a unit cell of the surface, with height three times the height of the groove. On the top boundary, a constant shear was applied in the streamwise direction. A similar setup was recently used to investigate the robustness of LIS with spanwise grooves [9]. Having 50 cells in the grooves, the errors were less than 4% for µ i /µ ∞ = 0.5, 1 and 2, see fig. S2.
A grid refinement study was also performed for µ i /µ ∞ = 1, We = 100 and θ = 45°, where another grid with (N x , N y , N z ) = (384, 960, 1536) was used, increasing the total number of grid cells by a factor of 3.4. This gave a change of the friction coefficient by 0.76%. Plots are shown in fig. S3. A. Dynamic contact angle Following ref. [5], the contact angle can be found from the equation where θ num is the numerically imposed contact angle, θ stat is the static angle, Ca = U cl √ µ i µ ∞ /γ, with contact line velocity U cl , λ is a microscopic length scale corresponding to an effective slip length of the contact line, here taken to be constant, and ∆ is the wall-normal cell height. The angle θ stat is in this manuscript denoted only θ for simplicity. The contact line velocity is measured half a cell above the wall. The function G * (θ) is a monotonically increasing function. It is defined as G * (θ) = √ qG(θ), where q = µ ∞ /µ i and G(θ) is the original function derived by ref. [10], with the extra factor for increased symmetry, with When written in this form, it is apparent that f * (θ, q) = f * (π − θ, q −1 ), which is a necessary requirement, since eq. (S1) should be independent of which of the two components we consider. This is identical to the hydrodynamic theory of an apparent angle model suggested in ref. [10], where θ num is the apparent angle. Legendre and Maglio [5]  has suggested that λ should be on the order of the the physical slip length (typically on the order of 1 nm), and may or may not be used in combination with a numerically applied slip length. We use a no-slip condition, with a numerical value of λ = 2 · 10 −7 h, that with h ≈ 0.5 cm attains a realistic value. The full equation (S1) was solved by the Newton's method, and f * −1 (θ, N ) was integrated with the trapezoidal rule. For numerical reasons, the value has been limited to 30°≤ θ num ≤ 150°. We evaluated the grid-independence of the contact line model by looking at droplet spreading for three different cell sizes. A three dimensional half-sphere droplet of radius R init = 0.5 was placed in a box of size (L x , L y , L z ) = (2,1,2), with periodic boundary conditions in the x-and z-directions, and no-slip and shear-free boundary conditions in the negative and positive y-direction, respectively. The initial radius sets the length scale of the problem. The contact angle was initially 90°, when the static contact was set to θ stat = 60°, so that the droplet started to spread. The droplet had a density ρ = 1, viscosity µ = 0.25 and surface tension γ = 7.5. Using the density and the viscosity we can define a time scale τ = ρR 2 init /µ. The surrounding fluid had the same density, and the viscosity ratio was changed by varying the viscosity of the surrounding fluid. Viscosity ratios of both 0.5 and 2 were tested, to span the range of this study. The number of cells of the wall parallel directions were N x = 64 and N z = 64, and the number of cells of the wall-normal direction was varied between N y = 32, N y = 64 and N y = 128. The microscopic length scale was set to λ = 2 · 10 −7 . Resulting spreading radii are shown in fig. S4a. The differences between the curves are a few percentage, considered small enough for this study. Corresponding droplet contraction tests, with contraction from 90°to 120°, are shown in fig. S4b. The tests using the intermediate resolution were also repeated on an immersed boundary by extending the domain in the y-direction below the droplet by ∆L y = 0.5 and adding a solid slab there. The curves are hardly distinguishable between the two setups, as shown in fig. S4. The fluid parameters and the dynamic contact angle model are the same as those used by Legendre & Maglio [5] (Dyn2) for the spreading case, except for the viscosity ratio (there equal to 1) and a slightly different value of λ. Their results agree well with those shown here.

S2. TURBULENCE STATISTICS
In this section, we show the turbulent statistics corresponding to the simulations We = {100, 150, 200} and viscosity ratios µ i /µ ∞ = {0.5, 1, 2}. Mean velocity profiles are plotted in fig. S5a in wall units and fig. S5b in outer units. As can be seen from the plot in wall units, the centerline velocity is increased by the LIS for low We. This is related to the drag reduction achieved at these We, whereas the opposite occurs for the cases of drag increase [11]. A reduction of the centerline velocity is followed by a decrease in the streamwise r.m.s. component, and an increase of the wall-normal and the spanwise fluctuations, shown in fig. S5c. Over all, drag increase results in an increased isotropy of the velocity fluctuations. Increased wall-normal velocity fluctuations at the surface has a strong connection to the increase in drag [12,13]. It implies stronger flow ejections, caused by the roughness the waves impose.
The pressure r.m.s. profiles are shown in fig. S5d. Outside the grooves, the profiles for We = 100 and We = 150 are similar to the smooth wall profile, but are higher inside the grooves. For We = 200, when there are large waves, the pressure fluctuations are increased also outside the grooves. This could be due to the increased roughness, but also due to formation of droplets.
In the streamwise r.m.s. component, fig. S5c, there is a peak at y = 0, but not for the spanwise nor for the S4 wall-normal. This can be related to the dispersive stresses from the solid structures, which are the r.m.s. of the roughness-coherent velocity, u RC . Due to the symmetry in the streamwise direction, the roughness-coherent velocity is the velocity field averaged in the streamwise direction and time. A peak is then created because the streamwise velocity component on average is larger over the interface at the grooves than over the ridges. The peak height differs between the viscosity ratios, but not so much between the different Weber numbers. This reflects the behaviour of the mean velocity at the interface, fig. S5a for small y + , as well as the slip length, fig. 3b. It was seen that both the spanwise and wall-normal roughness-coherent components were close to zero. The streamwise part is shown in fig. S6 for We = 100.  fig. 3 (We = 100blue, We = 150-red and We = 200-yellow), and the different lines represent different viscosity ratios, µ i /µ ∞ = 0.5 (---), µ i /µ ∞ = 1 (--) and µ i /µ ∞ = 2 (-· -). Also shown are statistics for a smooth wall (· · · · · ·). The spanwise velocity r.m.s. value is represented by w + rms for convenience, but elsewhere w refers to the groove width. In this section, we look at the governing equations for the Miles instability. We start from the equations for a infinitesimal disturbance on a baseflow U = U (y) [14], and ∂u ∂x + ∂v ∂y where denotes a derivative in the y-direction. The spanwise velocity component is here represented by w for convenience, which elsewhere refer to the groove width. Eliminating p gives Rayleigh's equation for v (inviscid Orr-Sommerfeld equation), with p given by 1 ρ We further assume an interface shape of a wave travelling in the streamwise direction, η = Ae ikx(x−ct) cos(k z z) = Ae kxcit e ikx(x−crt) cos(k z z) = a(t)e ikx(x−crt) cos(k z z), where η is the location of the interface, A is the wave amplitude, k x is the streamwise wavenumber, k z is the spanwise wavenumber and c = c r + ic i . Henceforth, we'll use the notation k = k 2 x + k 2 z . The curvature of the interface gives rise to a pressure difference over the surface of (going from negative to positive y), The perturbation must die out at infinity, v → 0 as y → ∞.
In addition to this, fluid parcels on the interface must remain on the interface, (so that the interface remains a streamline) [15], Returning now to Rayleigh's equation, eq. (S8), with (x, z, t)-dependence of v as η, From eq. (S9), with (x, z, t)-dependence of p as η,