2.1. Flow cell and LIS design

Figure 1(a) shows the flow channel used in this study. The height of the main chamber of the duct is delimited by the coordinates $y=0$ and $y=H$ , with $H= 770\,\mu \mathrm{m}$ , the width is $W=3.25\,\textrm {mm}$ and the length is $L=80\,\textrm {mm}$ . The top surface is patterned with rectangular grooves parallel to the flow direction, which constitute the solid substrate of the LIS. The inlet and outlet of the flow are located at the extremities of the bottom wall through holes of $1\,\textrm {mm}$ in diameter. The detailed description of the channel assembly is reported in Appendix A.1. Hexadecane (Sigma–Aldrich), with dynamic viscosity $\mu _{{l}}=3.1\,\mathrm{mPa\,s}^{-1}$ , is employed as a lubricant to impregnate the grooves, while the working fluid consists of an aqueous solution of glycerol–water mixed with milk at 20 % volume ratio (later referred to as simply water) with dynamic viscosities ranging from $\mu _{{w}}= 0.9 \,\mathrm{mPa\,s}^{-1}$ to $40.7\,\mathrm{mPa\,s}^{-1}$ depending on the glycerol ratio. The viscosity ratio between water and lubricant is defined as $N=\mu _{{w}}/\mu _{{l}}$ . Note that the lubricant infused in the grooves is located at coordinates $y\lt 0$ .

A total of five experimental cases with different geometrical features and viscosity ratio were designed (see table 1). Three grooved surfaces, differing by the groove aspect ratio, are fabricated in transparent resin using a casting technique (see Appendix A.2). An example of a groove profile is shown in figure 1( $b$ ). The grooves have depth $k$ of approximately $230\,\mu \mathrm{m}$ , while their widths $w$ vary from 223 to $376\,\mu \mathrm{m}$ , resulting in aspect ratios $A=k/w$ ranging from 1.06 to 0.59. The slipping area fraction, defined as $a=w/p$ , ranges between 0.45 and 0.61. The influence of the viscosity ratio on the slip length was tested using the LIS with aspect ratio $A=1.06$ . According to theoretical models, the effective slip length can vary significantly with both groove aspect ratio $A$ and the groove slipping fraction $a$ (Schönecker et al. Reference Schönecker, Baier and Hardt2014). However, for the LIS geometries considered here, the effective slip depends weakly on $A$ and $a$ , i.e. we do not investigate shallow grooves ( $A\lt 0.5)$ or very thin ( $a\lt 0.2$ ) or wide slits ( $a\gt 0.8$ ). We merely investigate different geometries to ensure that our observations are not specific to one particular LIS geometry. More specifically, as shown in Schönecker et al. (Reference Schönecker, Baier and Hardt2014), the effective slip length increases with $a$ since the interface relative to solid increases, allowing for larger slip. The effective slip length also increases with $A$ as the viscous dissipation in the lubricant is reduced with deeper grooves, until around $A=1$ , where the velocity gradients of the lubricant near the interface become independent of the flow deep in the groove. For our configurations: $A=0.59$ , $A=0.89$ and $A=1.06$ , with $a=0.61$ , $a=0.45$ and $a=0.45$ , respectively, the theoretically predicted effective slip lengths (reported in second column of table 2) show a modest variation.

Table 1. The parameters defining the LISs for the five configurations (see also Figure 1): groove width $w$ , groove depth $k$ , groove pitch $p$ , aspect ratio $A=k/w$ , slipping surface fraction $a=w/p$ , viscosity ratio $N=\mu _{w}/\mu _{{l}}$ and protrusion angle of the meniscus $\phi$ .

Table 2. Central slip length value predicted by Schönecker & Hardt (Reference Schönecker and Hardt2013), $\beta ^{{Sch}}(0)$ , measured average slip length over the lubricated groove, $\langle \beta ^{\textit{int}}\rangle$ , ratio $R=\beta ^{{Sch}}(0)/\langle \beta ^{\textit{int}}\rangle$ , estimated shear stress from conserved (i.e. closed cavity) lubricant $\tau _{{l}}$ , measured shear stress at the water–lubricant interface $\tau _{{w}}$ and their ratio.

The water solution was mixed with glycerol in mass ratios of $1:0.7$ and $1:3.5$ , increasing the dynamic viscosity of the water–glycerol solutions to $3.6 \,\textrm {mPa} \,\textrm {s}^{-1}$ and $40.7 \,\textrm {mPa} \,\textrm {s}^{-1}$ , respectively (computed using the empirical formulae of Cheng (Reference Cheng2008) and Volk & Kähler (Reference Volk and Kähler2018)). Accordingly, the viscosity ratio $N$ is increased from $0.3$ to $1.2$ and $13.2$ . As detailed in Appendix A.2, the LIS is prepared by first filling the duct with lubricant and then with the water solution at creeping flow. The lubricant remains trapped inside the grooves by capillary forces. After complete filling of the duct, the flow rate was gradually increased to ${0.4} \,\textrm {mL} \,\textrm {min}^{-1}$ , which provides a bulk velocity $U_b={2}\,\textrm {mm} \,\textrm {s}^{-1}$ and a Reynolds number $\textit{Re}\approx 1$ . The imposed velocity is a compromise between having a high D-OCT resolution (higher resolution at higher velocities) while maintaining the stability of the lubricant-filled groove.

2.2. Optical coherence tomography

Optical coherence tomography is a low-coherence interferometry-based technique used to generate images that capture optical scattering within opaque samples on a millimetre scale, with a resolution of just a few micrometres. While OCT has primarily been used in the medical field (Yonetsu et al. Reference Yonetsu, Bouma, Kato, Fujimoto and Jang2013), recently it has been employed also in microfluidics and soft-matter research (Haavisto et al. Reference Haavisto, Salmela, Jäsberg, Saarinen, Karppinen and Koponen2015; Gowda et al. Reference Gowda, Brouzet, Lefranc, Söderberg and Lundell2019; Lauri et al. Reference Lauri, Haavisto, Salmela, Miettinen, Fabritius and Koponen2021; Huisman et al. Reference Huisman, Blankert, Horn, Wagner, Vrouwenvelder, Bucs and Fortunato2024; Gupta et al. Reference Gupta, Daneshi, Frigaard and Elfring2024; Lauga et al. Reference Lauga, Brenner and Stone2005; Amini et al. Reference Amini, Wittig, Saoncella, Tammisola, Lundell and Bagheri2025).

Figure 2. Schematic of velocity components obtained through D-OCT. The phase of the incoming reference beam (red) is compared with that of the backscattered beam (blue) from the moving particle to obtain the parallel component of the particle velocity.

The functioning principle of OCT is based on a Michelson interferometer. A broadband light source is split into two components using a beam splitter. One component is directed into a reference arm, while the other is transmitted through the sample. The light reflected from the sample is recombined at the beam splitter with the reference beam, producing an interference signal. When the optical-path-length difference is within the coherence length of the light source, interference occurs (Tomlins & Wang Reference Tomlins and Wang2005). The output spectrum is analysed by a spectrometer, generating a complete interference pattern that represents variations in intensity as a function of depth. Within this pattern, refractive index variations between layers in the sample manifest as intensity magnitude variations, providing detailed information about the sample’s depth.

When there is a flow, the OCT imaging technique can be combined with the acquisition of a Doppler frequency (D-OCT) to obtain velocimetry information (Chen et al. Reference Chen, Milner, Srinivas, Wang, Malekafzali, van Gemert and Nelson1997). The interference of light backscattered from a moving particle with the reference beam produces a beating at the Doppler frequency, which generates a corresponding Doppler phase shift. This phase shift is proportional to the projected component of the velocity measured along the OCT beam, $u_{\textit{OCT}}$ , as illustrated in figure 2. To reconstruct the velocity component in the $x$ (streamwise) direction, an angle $\alpha$ between the beam and the vertical direction is required, yielding $u=u_{\textit{OCT}}/\sin {\alpha }$ . The velocity resolution depends on the acquisition time and the scan angle. We used a Telesto II Spectral Domain OCT apparatus (Thorlabs Inc., NJ, USA), which has a central wavelength of 1310 nm and a nominal bandwidth of 270 nm, leading to a depth resolution of $2.58\,\mu \mathrm{m}$ in water.

The measurement of $u_{\textit{OCT}}$ is related to the doppler phase shift

(2.1) \begin{align} u_{\textit{OCT}} = \frac {\lambda f \Delta \varphi }{4\pi } \end{align}

where $\lambda = {1310}\,\textrm {nm}$ is the central wavelength of the infrared source, $f = {5.5}\,\textrm {kHz}$ is the scanning frequency and $\Delta \varphi \in [-\pi ,\pi ]$ (rad) is the phase difference between two consecutive acquisitions, corresponding to the Doppler shift. This leads to the maximal measurable velocity $u_{{OCT,max}} = {\lambda f}/{4} = {7.2}\, \textrm {mm} \,\textrm {s}^{-1}$ . For more details, see Amini et al. (Reference Amini, Wittig, Saoncella, Tammisola, Lundell and Bagheri2025).

In these experiments, commercial milk (lactose-free pasteurized semiskimmed cow milk: Arla Ekologisk Lactosfri Mellanmjölkdryck 1.5 % fett) is added to water at a volume ratio of $1:5$ as an infrared diffusing contrast agent for the OCT/D-OCT measurement. No significant change in the viscosity, interfacial tension or refractive index of water is measured at this concentration of milk. As milk mixes with water rather than with the lubricant, the contrast in the intensity signal allows the two fluid components to be clearly distinguished. Optical coherence tomography is used here to image the LIS inside the duct during flow and in Doppler mode for the acquisition of velocity profiles. The flow cell is positioned below the OCT objective at an angle $\alpha \approx {6}^{\circ }$ . The OCT beam enters the duct through the clear top wall.

2.2.1. Intensity measurements

The intensity signal acquired with OCT across the flow cell allowed us to monitor the distribution of lubricant within the grooves during the initial infusion and the subsequent measurements. As an example, figure 3 shows a cross-section of the flow cell during flow with the LIS mounted at the top wall (highlighted with a dashed line) and lubricant within the grooves. From this tomogram, it can be seen that the water–lubricant interface is curved with a meniscus bending towards the top. This is due to a pressure difference across the interface and to the hydrophobicity of the grooved wall. The shadows in the working fluid below the vertical walls of the grooves are due to the extreme grazing angle between the OCT beam and these vertical structures. These regions are excluded from the measurements in the following. The deflection $\phi$ of the liquid interfaces from the level $y=0$ was measured from the tomograms in streamwise positions corresponding to $4.5$ , $5.0$ , $5.5$ , $6.0$ and 6.5 cm from the duct inlet. For each case, the streamwise average of the deflection and the associated standard deviation were calculated and reported in table 1. In the following, it is assumed that the interface is of constant curvature $\kappa = 2 \tan \phi /w$ in the spanwise direction, and zero curvature in the streamwise (longitudinal) direction.

Figure 3. Cross-sectional scan of the duct acquired with OCT at streamwise position corresponding to 4.5 cm from the inlet. The grooves, marked as ‘solid’ are mounted on the top wall and their profile is highlighted with a red dashed line. The space between the grooves is filled with lubricant, marked with ‘lub’, whose profile bows at an angle $\phi$ to the top of the grooves. The volume of the duct is filled with a water–milk mixture which appears opaque. The insert shows an enlargement of the measurement positions marked with blue dots at the water–solid interface and with green dots at the water–lubricant interface.

2.2.2. Velocity measurements

The measurements of velocity profiles were conducted through the full depth of the duct at the groove closest to the axial centre of the duct, where the influence of sidewalls on the flow is minimal. The acquisition frequency was 5.5 kHz. The locations of the local measurements are shown in the insert of figure 3 where they are marked by blue-filled circles at the water–solid boundary and by green circles at the water–lubricant interface. The spanwise spacing between them is $30\,\mu \mathrm{m}$ . For each configuration, five spanwise series of velocity profiles straddling a pitch were recorded in the streamwise locations corresponding to 4.5, 5.0, 5.5, 6.0 and 6.5 cm from the duct inlet, where the flow was fully developed. The results presented are obtained from a single experiment performed for each configuration.

Figure 4. ( $a$ ) Contour plot of the velocity measurements in the central region of the duct for case A106N13. The arrows indicate the location of the measurements shown in ( $b$ ). ( $b$ ) Velocity profiles were acquired at the water–lubricant interface (green markers) and the water–solid interface (blue markers), for the same case. In the insert, the coloured area surrounding the markers represents the standard deviation of the velocity.

The velocities measured in the vicinity of the groove, normalized with the bulk velocity $U_b$ , are shown as a contour plot in figure 4( $a$ ) for case A106N13. Note that the y-axis is now directed upwards and the grooved surface is at the bottom. The black solid lines represent the solid ridges, while the blue and green arrows show the locations of the velocity profiles displayed in figure 4( $b$ ). The two velocity profiles in figure 4( $b$ ) are shown using the same colour code as the arrows in figure 4( $a$ ). The top and bottom walls of the flow cell have coordinates $y=0$ and $y=H$ , respectively. They both peak at a normalized velocity $u/U_b = 1.83$ which is compatible with the viscous flow in a rectangular channel of section $770 \,\mu \textrm {m}\, \times \,3.25\, \textrm {mm}$ where $u_{th}/U_b = 1.84$ is predicted Panton (Reference Panton1984) and Amini et al. (Reference Amini, Wittig, Saoncella, Tammisola, Lundell and Bagheri2025). The velocity profile measured on the solid ridge is analogous to a Poiseuille flow, exhibiting zero velocity at the walls. In contrast, the profile measured on the liquid–liquid interface deviates from the Poiseuille profile due to the presence of slip at the interface. The increased distance from the top wall of the green profile is due to the interface protruding towards the groove bottom. In the insert of figure 4( $b$ ), an enlargement of the velocity profiles shows the uncertainty intervals of the measurements calculated as the standard deviation which, near the boundary, equals to $\pm {5}\,\%$ of the measured value. In the following, the key parameter extracted from the velocity measurements is the streamwise-averaged interfacial slip velocity $u_{{s}}(z)$ , which depends on the $z$ position due to the interface bending.

3.1. Local slip length

Figure 5( $a$ ) shows the streamwise-averaged slip velocities $u_{{s}}(z)$ normalized by the bulk velocity $U_b$ , measured along both the liquid–solid and liquid–liquid interfaces within a single pitch. Figure 5( $a$ i–iii) correspond to cases with increasing aspect ratio and constant viscosity ratio (A059N03, A089N03 and A106N03). Figure 5( $a$ iv,v) represent configurations with the same geometry as figure 5(a iii) but with higher viscosity ratios (A106N1 and A106N13). The green background represents the physical width of the groove – containing the liquid–liquid interface – while the blue colour represents the area occupied by the liquid–solid interface. In all configurations, the slip velocity at the solid boundary is approximate ${0.2}\,\%$ of $U_b$ . Given the accuracy of our measurements, this is compatible with the expected no-slip condition there. At the liquid interface, the slip velocity averages ${1.3}\,\%$ of $U_b$ for the cases with low viscosity ratio $N=0.3$ , increasing to ${2}\,\%$ and ${3}\,\%$ with increasing viscosity of the water phase, with $N=1$ and $N=13$ , respectively. The variance associated with the slip velocities are calculated as the standard deviation of the repeated scans, one spanwise set per each of the five streamwise positions, whose average gives a single velocity measurement.

Figure 5. See table 1 for geometrical parameters: ( $a$ ) streamwise-averaged slip velocities as a function of the spanwise position $z$ . The error bars correspond to the standard deviations. The sketch in (i) indicates that the values are calculated at the interface. ( $b$ ) Local slip lengths obtained from measurements, $\beta ^{\textit{int}}(z)$ (markers) and from the analytical model of Schönecker & Hardt (Reference Schönecker and Hardt2013), $\beta ^{{Sch}}(z)$ (solid line). In (iv) and (v), the curve corresponding to the model exceeds the boundary of the plot. The background colour of the tiles identifies the region of the liquid–liquid interface (green) and the region of the liquid–solid interface (blue).

For each velocity profile, the ratio between the slip velocity and the velocity gradient $\partial u/\partial y$ is calculated and the local slip length at the interface, $\beta ^{\textit{int}}(z)$ , is evaluated as the average of the repeated measurements at $n$ streamwise positions, ( $x_j=4.5,5.0,5.5,6.0,{6.5}\,\textrm {cm}$ ),

(3.1) \begin{equation} \beta ^{\textit{int}}(z)=\dfrac {1}{n}\sum _{j=1}^n\frac {u_{{s}}(z,x_j)}{\frac {\partial u(z,x_j)}{\partial y}|_{\textit{interface}}}. \end{equation}

The superscript ‘int’ denotes variables evaluated at the curved liquid–liquid interface, corresponding to the red curve in the sketch in figure 5( $a$ i). The variation in slip velocity between the solid and liquid boundaries translates directly to the corresponding slip length values shown in figure 5( $b$ ). At the liquid–solid interface, the measured slip length remains consistently around $1\pm 1\,\mu \mathrm{m}$ (no-slip condition) across all cases. In contrast, the slip length at the liquid–liquid interface is $5\pm 1\,\mu \mathrm{m}$ for the first three cases, increasing to $7\pm 2\,\mu \mathrm{m}$ and $13\pm 2\,\mu \mathrm{m}$ for the fourth and fifth cases, respectively. The uncertainty associated with the local slip length is the standard deviation of the repeated measurements in the different streamwise positions (see Appendix B).

To understand how the slip length of a realistic LIS set-up differs from an ideal theoretical configuration, we compare our results to the local slip length predicted by the model of Schönecker & Hardt (Reference Schönecker and Hardt2013), here named $\beta ^{{Sch}}(z)$ . Their analytical equation characterizes the variation along the pitch of the local slip length of a shear flow over an infinitely long rectangular groove aligned with the flow direction. The cavity contains a second fluid with different viscosity, and the interface between the two fluids is flat. The local slip length is modelled as an elliptic function with a maximum value at the centre of the groove (where the flow is least affected by the sidewalls of the cavity itself), expressed by the dimensionless variable $D$ . These expression reads

(3.2) \begin{align} \beta ^{{Sch}}(z)&={\textrm{Re}}\Biggl \{-i2ND\sqrt {z^2-\frac {w^2}{4}}\Biggl \} \end{align}

(3.3) \begin{align} D&=d_0\times{\textit{erf}}\left ( \frac {\sqrt {\pi }}{2d_0}A\right ) \\[9pt] \nonumber \end{align}

where $\textrm{Re}$ stands for the real part, $d_0 = 0.347$ and erf is the error function. The modelling of $D$ is discussed in depth by Schönecker & Hardt (Reference Schönecker and Hardt2013) and is based on previous investigations of lid-driven cavity flows. Note that, in their consecutive work (Schönecker et al. Reference Schönecker, Baier and Hardt2014), the analytical expression of the local slip length distribution was derived for an array of parallel lubricant-filled grooves, which accounts for a modification of the velocity field when the distance between the grooves is small ( $a\to 1$ ). In our case, however, the simpler expression of single groove (3.2) is sufficiently accurate since the slip fraction is small (the difference between the expressions is smaller than 2 %).

Figure 5(b) compares the experimental results with the analytical solution (3.2). The measured values deviate from the elliptical profile predicted by the model; we observe instead a flat shape with a nearly constant local slip along the interface. To best of our knowledge, this constant spanwise slip velocity profile is not predicted by any existing model. One may speculate that it is due more complex Marangoni stresses than what is predicted by existing models (Landel et al. Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2019; Sundin et al. Reference Sundin, Zaleski and Bagheri2021; Baier Reference Baier2023).

Next, we compute the spanwise-average slip length measured at the interface, $\langle \beta ^{\textit{int}}\rangle$ , and use it as a measure of the local interfacial slip. Note that the effective slip length, which is widely used to characterize surface slip, is related to the resistance the fluid experiences as it flows through the channel, and has thus a different physical interpretation than the average slip length. Moreover, for our experimental set-up and for non-flat interfaces, it is not straightforward to compute the effective slip (see more details in Appendix F). The average slip length on the other hand is experimentally accessible and allows us to quantify the deviation of measured interfacial slip with predictions from theoretical and numerical simulations. We observe that $\langle \beta ^{\textit{int}}\rangle$ , is significantly lower than the central value of the model distribution, $\beta ^{{Sch}}(0)$ . Specifically, as reported in table 2, the ratio $R$ between those two values varies from $5$ to $7$ for cases with small $N$ and increases to $80$ in the case with the highest viscosity ratio.

These observations may be explained by considering how our set-up deviates from the ideal model. The three most evident differences are: (i) the finite boundaries of the duct which limit the number and the length of the grooves; (ii) the curvature of the water–lubricant interface; (iii) the presence of the scattering medium in the water phase that can affect the dynamics of the interface.

Regarding the first point, we can deduce from the velocity measurements that in the central region of the duct (see figure 4 $a$ ), the spanwise velocity variations are small, as the velocity field is nearly uniform in this region. On the other hand, the finite length of the grooves can cause a backflow of the lubricant opposing the water flow, contributing to a reduction of slip velocity with respect to the case of infinitely long grooves. The contribution of the shear stress arising from this backflow is discussed in the next section.

To investigate the impact of the meniscus at the liquid–liquid interface on the flow, numerical simulations were conducted. The computational set-up replicated the LIS designs and enforced the interface curvature observed in the experiments (table 1). A detailed description of the numerical set-up and methods is provided in § 4. Figure 6 shows the results obtained for each of the LIS configurations. The slip lengths obtained from simulations (dashed lines) at a curved interface are slightly smaller than the theoretical predictions of Schönecker & Hardt (Reference Schönecker and Hardt2013). This indicates that the presence of a concave interface does not significantly modify the slip length distribution as observed in figure 5( $a$ ). The third point, i.e. the effect of contrast medium in the working fluid, is treated in the following section.

Figure 6. Effect of the curved interface on the local slip length. Continuous black lines correspond to $\beta ^{{Sch}}(z)$ , red dashed lines to the results of direct numerical simulations (DNS) (with meniscus curvature according to table 1).

3.2. Surfactant contamination

As described in § 2, to measure the flow velocity using D-OCT, the fluid must be light-scattering in the near-infrared. For this purpose, commercial milk (as previously described in § 2.2) is mixed with water, as it provides a solid Doppler signal. However, milk might affect the interfacial properties between water and the lubricant infusing the LIS. Bovine milk contains two major proteins: casein and $\beta$ -lactoglobulin, which account for approximately 80 % and 20 % of the total protein content, respectively (O’Mahony & Fox Reference O’Mahony and Fox2014). Both these proteins have a hydrophilic and an oleophilic part to their structure, as detailed in McSweeney & Fox (Reference McSweeney and Fox2013). We may therefore regard them as soluble surfactants, which can accumulate at the water–lubricant interface, as well as being in solution. In a surfactant-laden fluid, the surfactant molecules adsorb at the interface (gas–liquid or liquid–liquid) and, advected by the flow, may accumulate at stagnation points (Manikantan & Squires Reference Manikantan and Squires2020). This effect leads to a gradient of surfactant concentration which generates Marangoni stresses opposing the flow.

The observation that the slip velocity at the interface is significantly lower than predicted by the analytical model can be explained by taking into account the Marangoni stress into the tangential stress balance at the interface. At the water–lubricant interface the shear stress of the working fluid, $\tau _{{w}}$ , is balanced by two contributions. The first is given by the viscous resistance of the lubricant inside the groove, $\tau _{{l}}$ , while the second contribution is given by the Marangoni stress, $\tau _{{\textit{Ma}}}$ . Thus, we can write

(3.4) \begin{equation} \tau _{{w}} = \tau _{{l}} + \tau _{{\textit{Ma}}}. \end{equation}

The shear stress of the bulk fluid is calculated from the velocity measurements as $\tau _{{w}}=\mu _{{w}} {\partial u}/{\partial y}$ at the interface. The shear stress imposed by the lubricant on the interface can be estimated by assuming a velocity profile in the cavity. Neglecting interface curvature and confinement effects from the sidewalls of the duct, two limiting cases can be distinguished. In one case, we may assume that the lubricant is not conserved, i.e. lubricant is allowed to leave the grooves at the very downstream position, and if no new lubricant is entering upstream, the groove will eventually be drained (albeit very slowly). In this case, we expect a Couette profile in the groove and $\tau _{{l}} = \mu _{{l}} {u_{{s}}}/{k}$ , where $u_{{s}}$ is the slip velocity. In the second case, the lubricant is conserved because of a physical barrier downstream that creates a recirculating flow in the groove. We then have a pressure gradient that opposes the Couette flow to enforce lubricant conservation across the cavity section and $\tau _{{l}} = 4 \mu _{{l}} {u_{{s}}}/{k}$ (Busse et al. Reference Busse, Sandham, McHale and Newton2013). We estimate the lubricant shear stress using the expression for the conserved lubricant, since the grooves at the very end of channel encounter a wall. The lubricant shear stress calculated with this hypothesis is reported in table 2.

An estimation of the Marangoni stress can be obtained by computing the ratio of $\tau _{{l}}$ to $\tau _{{w}}$ . A stress ratio near unity means that $\tau _{{\textit{Ma}}}\ll \tau _{{l}}$ , whereas a ratio near zero indicates $\tau _{{\textit{Ma}}}\gg \tau _{{l}}$ . As shown in table 2, where $\tau _l$ is estimated assuming lubricant conservation in the groove, Marangoni stresses dominate over the lubricant shear stress as $\tau _l/\tau _w \leqslant 0.3$ for all the configurations. In particular, we find that $\tau _{{\textit{Ma}}}\approx 0.70\tau _{{w}}$ for the cases at low viscosity ratios and increases to $0.88\tau _{{w}}$ and $0.98\tau _{{w}}$ for $N=1$ and $N=13$ , respectively.

A related configuration was treated by Sundin & Bagheri (Reference Sundin and Bagheri2022), who modelled a two-dimensional flow carrying surfactants over transverse lubricated grooves. Their model is adapted to the present longitudinal configuration. The main difference lies in the fact that the Marangoni gradient builds along a much longer length (along $L = 8$ cm instead of the width of the groove $w \simeq 250\,\mu \textrm {m}$ ). This calls for redefining some of the non-dimensional numbers of the problem from Sundin & Bagheri (Reference Sundin and Bagheri2022). We derive the longitudinal model as detailed in Appendix C. The interpretation of the expressions is essentially the same for longitudinal and transverse configurations. We obtain the following expression for the effective slip length in the presence of surfactants in the bulk flow:

(3.5) \begin{equation} \beta _l^{\textit{Sun}} = \beta _{\textit{SHS}}^{0,l} b_{\textit{LIS}}^l \left (1- \frac {1}{1+\alpha _d^l+\alpha _s^l} \right ). \end{equation}

Here, $\beta _{\textit{SHS}}^{0,l} b_{\textit{LIS}}^l$ is the effective slip length for a clean LIS. The surfactant diffusivity number

(3.6) \begin{align} \alpha _{{d}}^l =\frac {1}{\textit{Pe}_{{s}} {\boldsymbol{Ma}} c_k} \end{align}

contains interfacial Péclet number

(3.7) \begin{align} \textit{Pe}_{{s}} = \frac {\hat {u}_{\textit{SHS}}^{0,l} b_{\textit{LIS}}^lw}{D_s} \end{align}

where $D_s$ ( $\mathrm{mol\,m}^{-2}$ ) is the surface diffusion coefficient of the surfactant and $\hat {u}_{\textit{SHS}}^{0,l} b_{\textit{LIS}}^l$ ( $\mathrm{m\,s}^{-1}$ ) is the estimated velocity at the interface of a clean LIS (depending only on the geometry and the viscosities of the working fluid and the lubricant). It also contains the Marangoni number

(3.8) \begin{align} {\textit{Ma}} = \dfrac {\textit{nRT} \varGamma _m}{\mu _{{w}} U}, \end{align}

which corresponds to the ratio of Marangoni forces ( $\textit{nRT} \varGamma _m$ ) to viscous forces $\mu _w U$ . Here $\varGamma _{{m}}$ ( $\mathrm{mol\,m}^{-2}$ ) is the saturating surface concentration of surfactants, while $n, R, T$ are provided in table 3. We thus note the effects of surfactants on the effective slip is negligible if $\alpha _{d}^l\gg 1$ , which corresponds to either a very small interfacial Péclet number (surface diffusion smears out all gradients) or a very small Marangoni number (Marangoni effect too small). Note that $\alpha _{{d}}^l$ also contains the non-dimensional bulk concentration of surfactants $c_k=\kappa _a c_0/\kappa _d$ , where $\kappa _a$ ( $\mathrm{m^{3}\,mol^{-1}s^{-1}}$ ), $\kappa _d$ ( $\mathrm{s}^{-1}$ ) and $c_0$ ( $\mathrm{mol\,m}^{-3}$ ) are, respectively, the adsorption coefficient, desorption coefficient the reference bulk concentration of surfactants.

Table 3. Physicochemical parameters of $\beta$ -lactoglobulin used for the calculation of $\beta _l^{{Sun}}$ (from Wahlgren & Elofsson (Reference Wahlgren and Elofsson1997) and Rabe et al. (Reference Rabe, Verdes, Rankl, Artus and Seeger2007)). The concentration $c_0$ was calculated as a 20 % diluted solution of the standard concentration $0.3\,\mathrm{g\,l}^{-1}$ of $\beta$ -lactoglobulin for commercial cow milk (Walstra, Wouters & Geurts (Reference Walstra, Wouters and Geurts2005)).

The interfacial solubility number

(3.9) \begin{align} \alpha _{{s}}^l =\frac {1}{8}\frac {{Bi}'}{1+{\textit{Da}}_\delta } \frac {1}{{\textit{Ma}}c_k} \end{align}

contains two additional non-dimensional numbers, namely, the Biot number

(3.10) \begin{align} {Bi}' = \frac {\kappa _{{d}} L^2}{w\hat {u}^{0,l}_{{SHS}}b^l_{{LIS}}} \end{align}

and the Damköhler number

(3.11) \begin{align} {\textit{Da}}_\delta =\dfrac {\kappa _{{a}}\delta _{{L}}\varGamma _m}{D_s}, \end{align}

where $\delta _{{L}}$ (m) is the typical length of depletion/excess in the bulk due to surface adsorption/desorption in the longitudinal case. Here, the effects of surfactants are small when $\alpha _{{s}}^l\gg 1$ , which can be achieved for a very large Biot number (desorption kinetics dominates). The above explanations of $\alpha _{{s}}^l$ and $\alpha _{{d}}^l$ are simplified here, and more details of the different transport regimes can be found in Sundin et al. (Reference Sundin, Zaleski and Bagheri2021).

The expression for $\beta ^{\textit{Sun}}_l(c_0)$ depends on the reference bulk concentration of surfactants $c_0$ through $\alpha _{s}^l$ and $\alpha _{{d}}^l$ , and allows us to calculate the expected value at the assumed concentration of $\beta$ -lactoglobuline of our solution. The physicochemical parameters of $\beta$ -lactoglobulin used to apply this model to the present case are provided by Wahlgren & Elofsson (Reference Wahlgren and Elofsson1997) and Rabe et al. (Reference Rabe, Verdes, Rankl, Artus and Seeger2007) and reported in table 3. They were taken assuming the early times ( $t\lt {3000}\,\textrm {s}$ ) kinetics from Rabe et al. (Reference Rabe, Verdes, Rankl, Artus and Seeger2007) (figure 7a in their article), where the adsorption/desorption rates are known and slower than the interfacial dynamics of the protein, making it behaving effectively as a simple poorly soluble surfactant.

Figure 7. Effective slip lengths $\beta ^{\textit{Sun}}_l$ predicted by Sundin’s model as a function of the concentration $c_0$ in $\beta$ -lactoglobuline (lines) and spanwise-averaged slip lengths measured at the interface $\langle \beta ^{\textit{int}}\rangle$ at the present estimated concentration (markers). The colours associated with the cases are: blue for A106N03, A106N1, A106N13; orange for A089N03; yellow for A059N03. Solid lines represent cases at $N=0.3$ , the dashed line represents $N=1.2$ and the dash–dotted line $N=13.2$ .

Figure 7 shows $\beta ^{\textit{Sun}}_l$ as a function of the bulk concentration of $\beta$ -lactoglobulin (lines) compared with the spanwise-averaged slip length measured at the interface (symbols). For all our geometries, the model predicts $\beta ^{\textit{Sun}}_l\approx 0$ for bulk concentrations above $c_0 \sim 1 \times 10^{-6}\, \textrm {mol} \,\textrm {m}^{-3}$ . This threshold value is considerably lower than the estimated concentration of $\beta$ -lactoglobulin for the working fluid, $c_{0,w}={0.03}\,\textrm {mol} \,\textrm {m}^{-3}$ . The model thus indicates a total immobilization of the interface. In other words, the model predicts that the Marangoni stress at the interface completely dominates over the lubricant stress, i.e. $\tau _{{w}} = \tau _{{l}} + \tau _{{\textit{Ma}}} \simeq \tau _{{\textit{Ma}}}$ .

This prediction of surfactant effects using Sundin’s model and the simplified milk chemistry should be interpreted as a minimal estimate for interface immobilization, as other surfactants – such as casein, additional lactoglobulins or even potential contaminants from the channel walls – may also contribute. Including these components would likely lead to even lower critical concentrations of $c_0$ for the interface to remain mobile in figure 7. Our prediction therefore serves to highlight how far the system operates from the mobile regime defined by the parameters in the Sundin and Landel models. It also does not imply that $\beta$ -lactoglobulin alone is responsible for the observed immobilization, but rather that, if it were the only contributing factor, it would already be sufficient to account for the strong Marangoni gradients measured.

Note, however, that we measure non-zero slip lengths (figure 5 b), which means that the interface retains more mobility than predicted by Sundin’s model adapted to longitudinal grooves. It is not clear how the flow of surfactants at the interface behaves in this system and how it relates to adsorption and desorption dynamics. In fact, the experimental surfactant flow field (including surfactants at interfaces and in the bulk) can be more complex than what is assumed in the simplified geometry of the model (Baier Reference Baier2023). A deeper understanding would require further investigations under controlled chemical conditions as well as accurate observation of the stagnation points. The model by Sundin & Bagheri (Reference Sundin and Bagheri2022) nevertheless indicates that the reduced slip length, relative to the idealized LIS model, may result from surfactants in the system. To test this hypothesis, numerical simulations are performed, applying boundary conditions for either a clean mobile interface or a surfactant-laden immobile interface for the experimental configurations.

4.1. Configuration and method

We consider the two-phase configuration shown in figure 8. This domain represents a unit cell of the experimental configuration shown in figures 1 and 3. The two-dimensional domain is a square whose side (and therefore also the height) is the pitch $p$ . Moreover, $k$ is defined as in figure 1 and $l = p - k$ is the distance between the groove crest and the top boundary. Finally, $\theta$ is the angle that the local normal to the interface $\boldsymbol{n}$ forms with the vertical axis and $\phi$ is the deflection angle of the meniscus (also shown in figure 3). Periodic boundary conditions are applied in the spanwise direction, effectively simulating an infinite array of grooves. Moreover, no-slip boundary conditions are imposed on the solid surface. To drive the flow, a constant shear $\tau _\infty$ is applied at the top boundary. The variables are non-dimensionalized using the surface tension $\sigma$ , fluid density $\rho$ , and characteristic length scale $l$ , matching the order of magnitude of the dimensionless groups from the experiments: $\textit{Re} = \rho U l/\mu _{{w}}$ , ${We} = \rho U^2 l/\sigma$ and ${\textit{Ca}} = U \mu _{{w}}/\sigma$ , where $U$ is the velocity at $y = l$ . Geometric properties ( $a$ , $A$ ) and the viscosity ratio $N$ are also matched with the experiments. The two fluids, water and lubricant, have the same densities. Since the capillary and Weber numbers are small ( ${\textit{Ca}} \approx 10^{-4}$ and ${We} \approx 10^{-4}$ ), the interface is essentially static under the flow. It is initialized with a shape extracted from the OCT measurements (i.e. figure 3).

Numerical simulations were conducted using Basilisk (Popinet Reference Popinet2003, Reference Popinet2009, Reference Popinet2015), an open-source framework for solving partial differential equations on adaptive Cartesian meshes. The interface was tracked through a geometric volume-of-fluid (VOF) method (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999). Despite the low Reynolds number, we solve the time-dependent equations, as it provides a more stable algorithm. The equations are discretized using a time-staggered approximate projection method on a Cartesian grid. The time discretization of the viscous diffusion term is treated via a second-order Crank–Nicolson fully implicit scheme, while spatial discretization employs second-order finite volume on an octree grid. The solid wall is modelled using an immersed boundary (cut-cell) method (Johansen & Colella Reference Johansen and Colella1998; Mittal & Iaccarino Reference Mittal and Iaccarino2005; Schwartz et al. Reference Schwartz, Barad, Colella and Ligocki2006). The computational domain is a cubic volume with a side length equal to the pitch $p$ of the investigated configuration, with periodic boundary conditions in the spanwise and streamwise directions, no-slip condition at the bottom boundary and the solid wall, and a Neumann condition at the top one. A two-dimensional slice of the domain is shown in figure 8 and represents the physical domain under consideration. The cubic computational domain was chosen to satisfy the constraints imposed by the VOF solver. Specifically, in this approach, the interface is reconstructed using a piecewise linear approximation and is advected based on the velocity at the cell faces. To accurately capture the interface dynamics, it is therefore essential to resolve the streamwise velocity components at both the front and back faces of the computational cells. Note that we impose a shear at the top boundary, rather than reproducing the full experimental channel geometry in the simulations. The latter is computationally too expensive.

The mesh refinement is based on a tree-grid discretization, with cell dimensions calculated as the domain length $p$ divided by $2^{{n}_{\textit{ref}}}$ , where ${n}_{\textit{ref}}$ denotes the refinement level. The mesh refinement ranges from level 5 in the far field to level 8 near the interface and solid wall, ensuring that the cell size is at least 10 times smaller than the maximum local slip length. This refinement strategy is validated through a mesh independence study (see Appendix D).

The numerical method is validated by simulating the idealized case of a longitudinal groove with a flat interface, as considered in the models of Schönecker & Hardt (Reference Schönecker and Hardt2013) and Schönecker et al. (Reference Schönecker, Baier and Hardt2014). The effective slip lengths obtained numerically are compared with the analytical results and reported in Appendix D.

Figure 8. Schematic representation of the numerical set-up.

4.2. Comparison between experiments and simulations

To quantify the interfacial mobility of the experiments, we make use the numerical simulations to model the experimental configurations, which provides the full velocity fields above and within the cavity. The boundary conditions at the liquid–liquid interface in the simulations are tuned to represent two scenarios: (i) immobilized surfactant-laden interfaces and (ii) surfactant-free interfaces. In scenario (i), the immobilized regime is modelled by imposing zero velocity at the curved liquid–liquid interface. In scenario (ii), the continuity of velocity and tangential stress at the interface is enforced implicitly through the VOF to allow for local slip.

From the velocity fields, the local slip length distribution is computed by taking the ratio between the local normal shear rate and the velocity along a reference plane (as detailed in Appendix F). Since case (i) has, by definition, zero velocity at the liquid–liquid interface, we choose the reference plane to coincide with the crest plane $y=0$ (shown in figure 9 a). The slip lengths obtained from a simulation in scenarios (i) and (ii) are, respectively, denoted by $\beta ^0_{\textit{NS}}(z)$ and $\beta ^0_{{S}}(z)$ , where the subscripts ‘NS’ and ‘S’ refer to ‘no-slip’ and ‘slip’. The local slip length distributions from the simulations are compared with the experimental slip length obtained from (3.1). Note that experimental slip lengths are also evaluated at reference plane $y=0$ . The experimental slip length evaluated there is denoted by $\beta ^0(z)$ .

Figure 9 compares the slip lengths extracted from the simulations with those from the experiments. The five frames correspond to the five experimental configurations listed in table 1. The experimental results are represented with black symbols, while the simulated results are shown with solid lines (orange for case (i) and purple for case (ii)). For all cases, the experimental values $\beta ^0(z)$ lie between the two predictions, but much closer to the immobilized case (i). The difference between maximum values of $\beta ^0(z)$ and $\beta ^0_{\textit{NS}}(z)$ is smaller than 20 % for all cases, despite the simplifications inherent in the numerical model. In contrast, the slip lengths simulated with a fully mobile interface, $\langle \beta ^0_{{S}}\rangle$ are much larger than the experimental values. Furthermore, the difference between the numerical distributions with slip and no-slip conditions is significantly larger at higher viscosity ratios (cases A106N1 and A106N13). This is expected, given that the slip length in the absence of surfactants increases with $N$ (see figure 11 of Appendix D).

Figure 9. Experimental slip lengths $\beta ^0(z)$ and the DNS results at the reference plane ( $x,y=0$ ) with no slip (orange line) and slip (purple line) boundary condition at the interface. The sketch in ( $a$ ) represents the reference plane.

To gain further insight into how geometrical features affect the slip length, the average slip length is calculated and related to the aspect ratio of the grooves. The average slip length is defined as the arithmetic mean of the local slip length over a pitch (see Appendices E and F), i.e.

(4.1) \begin{equation} \langle \beta \rangle = \frac {1}{n}\sum _n\beta (z_n) \qquad \text{with} \qquad -\frac {p}{2} \lt |z_n| \lt \frac {p}{2} .\end{equation}

Table 4 reports the average slip lengths of the LIS computed over the pitch (that is taking into account both the liquid–liquid and liquid–solid interfaces) and we observe again a good match between the experimental values $\langle \beta ^0\rangle$ and $\langle \beta ^0_{\textit{NS}}\rangle$ .

Table 4. Average slip lengths derived from experiments $\langle \beta ^0\rangle$ compared with the numerical results in case of immobile $\langle \beta ^0_{\textit{NS}}\rangle$ and slipping $\langle \beta ^0_{{S}}\rangle$ interface.

At fixed groove geometry, increasing the viscosity ratio (more viscous working fluid and/or less viscous lubricant) tends to get the system closer to the immobilized case (i). Increasing $N$ by nearly a factor 50 ( $ N$ goes from 0.3 to 13) results in a 30 % increase of the average slip length in the presence of surfactants, whereas in the absence of contaminants, the increase is 88 %.

4.3. Discussion

For case (i), which models a fully immobilized interface, Marangoni stress dominates $\tau _{{\textit{Ma}}}\approx \tau _{{w}}$ , whereas for the surfactant-free case (ii), $\tau _{{\textit{Ma}}}=0$ . These two limiting scenarios can be used to estimate the Marangoni stresses present in the experiments using the model by Sundin & Bagheri (Reference Sundin and Bagheri2022). In their model, the effective slip is linearly dependent on the Marangoni stress. Assuming that this is also valid for the average slip length, we may fit a linear function using the two limiting values and estimate the stress ratio $ \tau _{{\textit{Ma}}}/\tau _{{w}}$ given measurements of the average slip length.

The last column of table 4 reports the stress ratio for the five experimental cases. It can be observed that, the stress ratio $\tau _{{\textit{Ma}}}/\tau _{{w}}$ ranges between $0.7$ to unity. From the tangential stress balance at the interface,

(4.2) \begin{equation} \frac {\tau _{{\textit{Ma}}}}{\tau _{{w}}}=1-\frac {\tau _{{l}}}{\tau _{{w}}}. \end{equation}

The estimates confirm that for all configurations, the Marangoni stress is the main source of resistance at the interface, i.e. it dominates over the resistance imposed by the lubricant stress. The ratio between the stress exerted by the water and that exerted by the lubricant is proportional to the viscosity ratio, i.e. ${\tau _{{l}}}/{\tau _{{w}}}\propto {\mu _{{l}}}/{\mu _{{w}}}={1}/{N}$ . Therefore from (4.2), we expect a larger Marangoni stress with increasing $N$ , which is also confirmed in table 4. Note that in § 3.2 we estimated the Marangoni stress by estimating the lubricant shear stress as $\tau _{{l}} \approx 4 \mu _{{l}} u_{{s}}/k$ . The estimates using this approach (last column in table 2) are in good agreement with the estimate made here through a calibration with limiting values of average slip length obtained from numerical simulations (last column in table 4).

It is important to discuss the measured average slip length $\langle \beta ^0\rangle$ in table 4, which, despite the minimum mobility of the liquid–liquid interface, amounts to 10 to $20\,\mu \mathrm{m}$ . The reported values represent the average slip at the crest plane of the ridges, i.e. the slip velocity is evaluated in the water phase rather than at the interface of the two fluids. This approach was previously used by Crowdy (Reference Crowdy2017) who treated a mathematically similar problem to ours. Crowdy studied a laminar flow over longitudinal grooves where the interface is immobile and the interface has a meniscus deflected by an angle $\phi$ from the plane $y=0$ . To derive the explicit expression of the effective slip length, Crowdy used a conformal mapping method to map the curved interface into the plane of the ridges. The simulated slip lengths of a completely immobile interface are in good agreement with Crowdy’s predictions (see figure 13 $a$ in Appendix D). Another analytical approach accounting for the presence of surfactant is the model of Baier (Reference Baier2023), where the curved surfactant-laden interfaces are not immobilized but rather incompressible. In this case, the velocity profile at the interface is non-zero with a zero average across the groove, thus predicting the presence of interfacial back-flows. This model explores other ways to add surfactants as a physical ingredient in LIS in a more refined way than immobilizing the interface. The average slip-length prediction remains the same as in the immobilized case of Crowdy (Reference Crowdy2017) (that is due to the presence of the meniscus). The predicted recirculations are not observed in our experiments.