# Surfactant spreading in a two-dimensional cavity and emergent contact-line singularities

Richard Mcnair, Oliver E Jensen, Julien R Landel

► To cite this version:

Richard Mcnair, Oliver E Jensen, Julien R Landel. Surfactant spreading in a two-dimensional cavity and emergent contact-line singularities. Journal of Fluid Mechanics, 2021, 930, pp.A15. 10.1017/jfm.2021.911. hal-04651714

HAL Id: hal-04651714

https://hal.science/hal-04651714v1

Submitted on 17 Jul 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Surfactant spreading in a two-dimensional cavity and emergent contact-line singularities

Richard Mcnair¹,†, Oliver E. Jensen¹,† and Julien R. Landel¹,†

¹Department of Mathematics, University of Manchester, Oxford Road, Manchester M13 9PL, UK

(Received 26 May 2021; revised 9 September 2021; accepted 14 October 2021)

We model the advective Marangoni spreading of insoluble surfactant at the free surface of a viscous fluid that is confined within a two-dimensional rectangular cavity. Interfacial deflections are assumed small, with contact lines pinned to the walls of the cavity, and inertia is neglected. Linearising the surfactant transport equation about the equilibrium state allows a modal decomposition of the dynamics, with eigenvalues corresponding to decay rates of perturbations. Computation of the family of mutually orthogonal two-dimensional eigenfunctions reveals singular flow structures near each contact line, resulting in spatially oscillatory patterns of shear stress and a pressure field that diverges logarithmically. These singularities at a stationary contact line are associated with dynamic compression of the surfactant monolayer. We show how they can be regularised by weak surface diffusion. Their existence highlights the need for careful treatment in computations of unsteady advection-dominated surfactant transport in confined domains.

**Key words:** Marangoni convection, contact lines, convection in cavities

## 1. Introduction

Surfactants play a crucial role in a variety of natural and industrial flows (Kovalchuk & Simmons 2021), including cleaning and decontamination (Landel & Wilson 2021), transport in lung airways (Filoche, Tai & Grotberg 2015; Stetten et al. 2018) and microfluidic applications (Temprano-Coleto et al. 2018). Understanding how surfactants equilibrate near contact lines is also important for flows over superhydrophobic surfaces, and slippery-liquid-infused porous surfaces (Wang & Guo 2020). As recently shown, surfactants can significantly affect the drag-reducing properties of superhydrophobic surfaces when a flow concentrates them at stationary contact lines (Landel et al. 2020; Baier & Hardt 2021; Peaudecerf et al. 2017). Surfactants are amphiphilic materials that

†Email addresses for correspondence: richard.mcnair@postgrad.manchester.ac.uk, oliver.jensen@manchester.ac.uk, julien.landel@manchester.ac.uk

© The Author(s), 2021. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
accumulate at interfaces, where they lower surface tension. Surfactant concentration gradients induce surface tension gradients, leading to so-called Marangoni flows of adjacent liquids that transport the surfactant itself (Manikantan & Squires 2020). Here, we address a canonical surfactant transport problem in a two-dimensional confined geometry, showing how singular flow structures arise when spreading is impeded by a solid boundary.

The self-induced spreading of a surfactant monolayer over a gas-liquid interface generates a variety of striking flow features (Afsar-Siddiqui, Luckham & Matar 2003; Matar & Craster 2009; Liu, Peco & Dolbow 2019). If surface diffusion and solubility effects are weak, the leading edge of a localised surfactant monolayer spreading on an otherwise uncontaminated interface effectively rigidifies the interface locally. For a monolayer spreading on a thin viscous film, this induces a jump in film depth via mass conservation effects (Gaver & Grotberg 1990; Dussaud, Matar & Troian 2005), that is captured within lubrication theory as a kinematic shock (Borgas & Grotberg 1988). Lubrication theory also predicts a jump in surface shear stress, although a more refined analysis over shorter length scales reveals an integrable shear-stress singularity at the leading edge of the monolayer, which can be regularised by surface diffusion or the presence of low levels of pre-existing (endogenous) surfactant (Jensen & Halpern 1998). Out-of-plane displacement of the free surface may be suppressed by surface tension or gravity (Gaver & Grotberg 1992; Jensen & Grotberg 1992). Inertial effects (which may be important if the initial spreading flow is sufficiently rapid) can generate an interfacial deflection at the leading edge of the monolayer known as the Reynolds ridge, arising due to displacement effects in the viscous boundary layer beneath the spreading monolayer (Scott 1982; Jensen 1998). Spreading on thin films may also be accompanied by dramatic secondary fingering instabilities (Warner, Craster & Matar 2004; Jensen & Naire 2006; Liu et al. 2019), showing the richness of surfactant flow phenomena.

The present paper addresses a complementary spreading problem, namely surfactant spreading at the free surface of viscous fluid that is confined within a two-dimensional cavity. We make a number of simplifying assumptions to aid our analysis: the free surface remains (almost) flat as a result of a restoring force, provided for example by surface tension, with contact lines pinned to the lateral walls of the cavity; the surfactant is insoluble and has a linear equation of state (a weakness discussed by Swanson et al. 2015); inertial effects are neglected and the Stokes flow is two-dimensional; molecular diffusion of surfactant at the free surface is assumed negligible, except when analysing its impact on the regularisation of the singularities at the contact line; and the interface is pre-loaded with endogenous surfactant, to which exogenous surfactant is added. These modelling assumptions and their implications are discussed further in § 4.

The aim of this study is to describe the interfacial and bulk transient flows produced by self-induced Marangoni spreading of surfactant in a confined geometry. Exogenous surfactant added to the interface spreads, compressing the endogenous surfactant ahead of it (Grotberg, Halpern & Jensen 1995; Sauleda et al. 2021). Since the surfactant monolayer is laterally confined, the surfactant concentration rises at the pinned contact lines, while Marangoni stresses drive further surfactants towards the stationary boundary. Although there is no motion of the contact line with respect to the solid wall, and therefore no risk of generating a non-integrable stress singularity associated with contact-line motion (Huh & Scriven 1971), we nevertheless find that the unsteady Marangoni flow generates its own singular flow behaviour. We find two separate classes of singularities generated at the corner, one inducing an oscillatory shear-stress pattern and the other a logarithmic pressure singularity. In the main part of the paper, we focus our study to the case of a single-fluid flow with a free surface pinned at the contact line at an angle of $\pi/2$. This special case simplifies the numerical simulation and asymptotic analysis. In § 4 we relax
Surfactant cavity flow

these assumptions and discuss other wedge angles, between 0 and $\pi$, and two-fluid flows with arbitrary viscosities and where the surfactants lie at the interface. Our asymptotic analysis shows that the structure of the flow and the type of singularities identified for the single-fluid flow with $\pi/2$ contact angle extend to this broader class of problems.

These findings complement studies of the lid-driven and shear-driven cavity problems (Shankar & Deshpande 2000), benchmark problems for computational fluid dynamics, for which corner singularities are a recognised computational challenge (Kuhlmann & Romanò 2019). Rather than tackle the full unsteady nonlinear problem numerically, our approach is to address the problem of small surfactant gradients, allowing the advective surfactant transport equation to be linearised. When coupled to the linear Stokes flow in the cavity, we derive a problem that admits a decomposition into eigenmodes, with each eigenvalue representing the decay rate of a particular modal disturbance. This approach allows us to focus our numerical effort on capturing spatial structures. While the temporal dynamics resembles a purely diffusive process, with mutually orthogonal modes decaying exponentially in time, each eigenmode has a singular form near the contact lines in the absence of surface diffusion. We combine asymptotic and numerical approximations to obtain a full understanding of the flow structure at the interface and in the bulk.

The model and the methods used to solve this surfactant Marangoni-driven cavity flow problem are described in § 2, with results presented in § 3. Implications of the study are discussed in § 4.

## 2. Model

We model the spreading of an insoluble surfactant at the surface of an incompressible liquid of dynamic viscosity $\mu^*$ in a two-dimensional rectangular domain. The domain $V$ spans from $-W^* \le x^* \le W^*$ and $-H^* \le y^* \le 0$ with $x^*$ and $y^*$ in the horizontal and vertical directions, respectively, with $W^*$ the half-width of the cavity, $H^*$ its height, as shown in figure 1. Stars in superscript indicate dimensional variables. To model the flow induced by the surfactant spreading at the free surface, we use the Stokes equations to relate the velocity field $u^*(x^*, t^*) = (u^*, v^*)$ to the pressure $p^*(x^*, t^*)$, with $x^* = (x^*, y^*)$ ignoring gravity and inertia. We impose the no-slip and no-penetration conditions on the three solid boundaries found at $x^* = -W^*$ and $W^*$ for $-H^* \le y^* \le 0$ for the sidewalls, and $y^* = -H^*$ for $-W^* \le x^* \le W^*$ for the bottom wall. The free surface $\mathcal{F}$, assumed flat to leading order, is located at $y^* = 0$ for $-W^* \le x^* \le W^*$. We balance the Cauchy stress $\sigma^* \cdot n$ at $\mathcal{F}$ (with the normal unit vector $n = (0, 1)$) with the gradient of the surface tension $\gamma^*$ tangentially and a restoring force due to strong surface tension normally. The insoluble surfactant concentration $\Gamma^*$ at the surface is coupled to the flow via a time-dependent advective transport equation, and to the surface tension via an equation of state, assumed linear, which is valid for small variations of the concentration of surfactant from a reference concentration (Stone & Leal 1990). The governing equations are therefore

$$ \nabla^* p^* = \mu^* \nabla^{*2} u^*, \quad \nabla^* \cdot u^* = 0, \quad \Gamma_{t^*}^* = -\frac{\partial}{\partial x^*} (u_s^* \Gamma^*), \quad \gamma^* = \gamma_0^* - (\gamma_0^* - \gamma_c^*) \frac{\Gamma^*}{\Gamma_c^*}, \quad (2.1a-d) $$

where subscript $s$ denotes evaluation on $\mathcal{F}$, $\gamma_0^*$ is the reference surface tension and $\gamma_c^*$ is its (lower) value when the surfactant is at a reference concentration $\Gamma_c^*$. The boundary conditions associated with (2.1a-d) are

$$ u^* = 0 \quad \text{on } x^* = \pm W^*, \quad y^* = -H^*, \qquad (2.2a) $$

$$ u^* \cdot n = 0 \quad \text{on } y^* = 0, \qquad (2.2b) $$
R. Mcnair, O.E. Jensen and J.R. Landel

<figure><img src="image_3.png" /><figcaption>Figure 1. Diagram of the two-dimensional rectangular domain of the flow problem. The flow is confined within rigid walls (hashed lines) and a free surface, for $-W^* \le x^* \le W^*$ and $-H^* \le y^* \le 0$. The incompressible Stokes flow in the bulk has velocity $u^*$ in the $x^*$ coordinate direction, and $v^*$ in the $y^*$ direction. At the free surface located at $y^* = 0$ and $-W^* \le x^* \le W^*$, an arbitrary initial non-uniform concentration profile of surfactant leads to an unsteady Marangoni flow that drives the flow in the bulk. The arbitrary initial concentration profile can be formed by exogeneous surfactant (in red) deposited on the free surface at $t^* = 0$, and some pre-existing endogenous surfactant with uniform concentration (in blue).</figcaption></figure>

$$ \sigma^* \cdot n = -\gamma^* n (\nabla^* \cdot n) + \frac{\partial \gamma^*}{\partial x^*} t \quad \text{on } \mathcal{F}, \qquad (2.2c) $$

with $t = (1, 0)$ the tangential unit vector at the free surface. We have made the assumption that $S^* = \gamma_0^* - \gamma_c^* \ll \gamma_c^*$, so that the surface tension remains sufficiently large, in comparison with its reduction by surfactant $S^*$, to suppress deflections of the free surface from $y^* = 0$. Effectively, the capillary number in our problem is $S^*/\gamma_c \ll 1$, since $S^*$ is a characteristic viscosity–velocity scale. This small capillary number assumption is discussed further in §§ 3 and 4. Hence, the leading-order kinematic boundary condition reduces to (2.2b) and any surface curvature can be neglected, such that $W^*\nabla^* \cdot n \ll 1$. We will exploit the normal stress condition later to evaluate the small surface deflections induced by the flow, while imposing the tangential component of (2.2c) on $\mathcal{F}$.

Using the length scale $W^*$, velocity scale $S^*/\mu^*$ and pressure scale $S^*/W^*$, we relate dimensional starred variables to their dimensionless counterparts via

$$ x = (x, y) = \left( \frac{x^*}{W^*}, \frac{y^*}{W^*} \right), \quad H = \frac{H^*}{W^*}, \quad \Gamma = \frac{\Gamma^*}{\Gamma_c^*}, \quad \gamma = \frac{\gamma^* - \gamma_c^*}{S^*}, \qquad (2.3a) $$

$$ u = (u, v) = \frac{\mu^*}{S^*}(u^*, v^*), \quad t = \frac{S^*t^*}{W^*\mu^*}, \quad p = \frac{W^*p^*}{S^*}. \qquad (2.3b) $$
## Surfactant cavity flow

After rescaling, the governing equations become, in the bulk,

$$ \begin{pmatrix} p_x \\ p_y \end{pmatrix} = \begin{pmatrix} u_{xx} + u_{yy} \\ v_{xx} + w_{yy} \end{pmatrix}, \quad u_x + v_y = 0, \quad \gamma = 1 - \Gamma, \qquad (2.4a) $$

with, at the free surface,

$$ \Gamma_t = -(\Gamma u)_x, \quad u_y = -\Gamma_x, \quad v = 0, \quad \text{on } y=0 \qquad (2.4b) $$

and

$$ u = 0, \quad \text{on } x = \pm 1, \text{ and } y = -H. \qquad (2.4c) $$

We introduce a streamfunction $\psi(x, y, t)$ such that $\psi_y = u$, and $\psi_x = -v$, enforcing incompressibility. The problem reduces to the biharmonic equation in the bulk

$$ \nabla^4 \psi = 0, \qquad (2.5) $$

subject to

$$ \Gamma_t = -(\Gamma \psi_y)_x, \quad \psi_{yy} = -\Gamma_x, \quad \psi = 0 \quad \text{on } y=0, \qquad (2.6a) $$

$$ \psi_x(\pm 1, y) = \psi(\pm 1, y) = 0, \qquad (2.6b) $$

$$ \psi(x, -H) = \psi_y(x, -H) = 0. \qquad (2.6c) $$

The streamfunction vanishes at the four boundaries of the domain, in order to impose the no-flux boundary condition. The problem is closed by imposing an initial surfactant profile, representing the addition of exogenous surfactant to an endogenous monolayer initially present on the interface. We note that the transport equation (2.6a) is linear in $\Gamma$, given a surface velocity $\psi_y$. Therefore writing $\Gamma$ as the sum of endogenous and exogenous components, $\Gamma = \Gamma_1 + \Gamma_2$ say, the two components satisfy the same transport equation $\Gamma_{it} = -(\Gamma_i \psi_y)_x$ ($i=1, 2$), allowing the evolution of the components to be tracked individually if necessary (Grotberg et al. 1995). From (2.6a,b) we anticipate the presence of singularities at the contact lines $(x, y) = (\pm 1, 0)$, as the boundary conditions are discontinuous here. We discuss these in detail in § 2.3 below, to understand their impact on the surface and bulk velocity fields and the surfactant distribution.

### 2.1. Linearisation

At large times, the surfactant relaxes to a uniform level $\bar{\Gamma}$ and the velocity field decays to zero. We perturb the system around this state, noting that the resulting problem is homogeneous. We decompose the solution into individual eigenmodes, writing for one such mode

$$ \Gamma(x, t) = \bar{\Gamma} + \hat{\Gamma}(x)e^{-\lambda t}, \quad \text{where} \quad \int_{-1}^{1} \hat{\Gamma}(x) \, dx = 0, \qquad (2.7) $$

assuming the same time dependence for other variables, for example $u(x, y, t) = \hat{u}(x, y)e^{-\lambda t}$. The surfactant transport equation at the free surface (first equation in (2.6a)) is the only equation that changes under linearisation, becoming

$$ \alpha \hat{\Gamma} = \hat{u}_x, \quad \text{on } y=0, \qquad (2.8) $$

where $\alpha = \lambda/\bar{\Gamma}$. Equation (2.8) is valid for small surfactant concentration perturbation, or at large times since we assume $|\hat{\Gamma}(x)|e^{-\lambda t} \ll \bar{\Gamma}$. Equation (2.8) may be combined with
the stress condition $\Gamma_x = -u_y$ to give the homogeneous boundary condition

$$-\alpha \hat{u}_y = \hat{u}_{xx}, \quad \text{on } y = 0. \qquad (2.9)$$

The resulting eigenvalue problem for the perturbation streamfunction may then be stated as the biharmonic equation (2.5), under the boundary conditions (2.6b,c) for $\hat{\psi}$ and

$$\alpha \hat{\psi}_{yy} = -\hat{\psi}_{xxy} \quad \hat{\psi} = 0, \quad \text{on } y = 0. \qquad (2.10)$$

We seek the set of decay rates $\alpha_n$, $n = 1, 2, \dots$, which are eigenvalues for which a solution exists with the corresponding eigenmodes $\hat{\psi}_n$. We distinguish two families for the decay rates $\alpha_n$ and eigenmodes $\hat{\psi}_n$, such that $\hat{\psi}_n$ is either an even or odd function of $x$. For the rest of this study the subscript $n$ will refer to the odd modes unless otherwise specified. From each of these eigenmodes and eigenvalues we can derive the associated surfactant distribution $\hat{\Gamma}_n$, the shear stress at the free surface $\tau_n(x) = \partial \hat{u}_n(x, 0)/\partial y$, the vorticity field $\hat{\omega}_n = -\nabla^2 \hat{\psi}_n$ and the pressure field $\hat{p}_n$, obtained from the vorticity via $\nabla \hat{p}_n = (-\partial \hat{\omega}_n/\partial y, \partial \hat{\omega}_n/\partial x)$.

This problem has two key global characteristics. An energy dissipation argument (see Appendix A) shows that all the decay rates satisfy

$$\alpha_n = \frac{\int_V \hat{\omega}_n^2 dA}{\int_{-1}^1 \hat{\Gamma}_n^2 dx}, \qquad (2.11)$$

where $V$ is the rectangular flow domain (see figure 1). Application of the reciprocal theorem (Masoud & Stone 2019) (see Appendix B) yields the orthogonality condition

$$\int_{-1}^{1} \hat{\Gamma}_{m} \hat{\Gamma}_{n} dx = 0, \quad \forall m \neq n. \qquad (2.12)$$

When solving an initial value problem with time-dependent surfactant and velocity fields, the condition (2.12) can be used to project the initial condition for $\hat{\Gamma}$ onto its component modes. The initial Marangoni gradient that drives the flow can arise for example from the addition of exogenous surfactant to an otherwise uniform endogenous surfactant distribution. In this study, we do not track the interface between these distributions explicitly, and we assume that the exogenous and endogenous surfactant concentrations add together to form a single concentration field (Grotberg et al. 1995).

## 2.2. Finite-difference numerical solution

We compute a numerical solution $\tilde{\psi}_n$ (where a wide tilde denotes a numerically computed variable) of the unknown perturbation modes of the streamfunction, $\hat{\psi}_n$, satisfying the biharmonic equation (2.5) and the boundary conditions (2.6b,c) and (2.10), using a finite-difference scheme. Using a row × column ordering convention, the domain described in figure 1 is discretised as an $M \times N$ rectangular grid, in the $y$ and $x$ directions, respectively, with uniform spacing $\Delta y$ and $\Delta x$, and supplemented with a set of ghost-points around the periphery creating an $(M + 2) \times (N + 2)$ grid. We use a second-order accurate 13-point stencil, involving standard finite-difference approximations of derivatives (see e.g. Fornberg 1988), to approximate the biharmonic operator in (2.5) on the grid of unknowns. These unknowns correspond to the unknown values of the perturbation
## Surfactant cavity flow

streamfunction of a given mode at a given grid point $\tilde{\psi}_n(x_j, y_i)$, with $3 \le i \le M$, $3 \le j \le N$. We impose $\tilde{\psi}_n(x_j, y_i) = 0$ at the boundaries of the domain $j = 2$ and $j = N + 1$, with $2 \le i \le M + 1$, and $i = 2$ and $i = M + 1$, with $2 \le j \le N + 1$, to implement the boundary conditions in (2.6b,c) and (2.10) which state that the streamfunction vanishes at all the boundaries. The $(M - 2) \times (N - 2)$ grid of unknowns is expressed as a column vector $\psi$ which is assembled by the vertical concatenation of the rows of unknowns of $\tilde{\psi}_n(x_j, y_i)$. The numerical operator modelling the biharmonic operator on the grid can then be approximated by a matrix operating on $\psi$. The remaining no-slip boundary conditions in (2.6b,c) at the walls and the surfactant boundary condition in (2.10) at the free surface can be approximated by a finite-difference discretisation acting on the $M \times N$ grid and the ghost points. Values of the streamfunction at the ghost points are calculated as functions of the values in the interior points and added to the linear system such that the boundary conditions are satisfied. The system can then be rearranged so that $\psi$ satisfies,

$$B\psi = \alpha C\psi, \qquad (2.13)$$

where **B** and **C** are sparse $(N-2)(M-2) \times (N-2)(M-2)$ matrices. Equation (2.13) represents a generalised numerical eigenvalue problem for the eigenmodes $\tilde{\psi}_n$ and the associated numerically calculated decay rates $\tilde{\alpha}_n$, from the free-surface boundary condition (2.10). We solved (2.13) using the MATLAB function eigs. From the solutions for each mode $\tilde{\psi}_n$ we compute numerical approximations of all other quantities of interest for each mode, such as the surfactant concentration profile $\tilde{\Gamma}_n(x)$, the surface stress $\tilde{\tau}_n(x) = \tilde{\psi}_{n,yy}(x, y = 0)$, the velocity field $(\tilde{u}_n, \tilde{v}_n)(x, y)$, the vorticity field $\tilde{\omega}_n(x, y)$ and the pressure field $\tilde{p}_n(x, y)$. The solutions $\tilde{\psi}_n$ are normalised by requiring $\max(\tilde{\Gamma}_n(x)) - \min(\tilde{\Gamma}_n(x)) = 1$ for each $n$.

We use global integrated measures to estimate the accuracy of the computational scheme, such as the energy balance (2.11), see table 1, and mode orthogonality, see Appendix B. The relative error of the decay rates computed directly as eigenvalues from (2.13) and indirectly from the eigenmodes via (2.11) (denoted as $\tilde{\alpha}_n$) remains less than $4 \times 10^{-4}$ for all odd and even modes calculated (up to $n = 4$) for $H = 2$ using a $4000 \times 4000$ grid (see table 1). For the same refinement and the same modes, the concentration profiles $\tilde{\Gamma}_n$ are orthogonal with an absolute error less than $8 \times 10^{-6}$ [see (B5)]. We use asymptotic methods, described below, to assess and contain (via global grid refinement) the inevitable local numerical inaccuracies associated with corner singularities at the contact lines $(x, y) = (\pm 1, 0)$.

## 2.3. Corner asymptotics

We anticipate from the outset the appearance of Moffatt vortices (Moffatt 1964) in the lower corners of the domain, at $(x, y) = (\pm 1, -H)$, as we will show in our results presented in § 3. However, the structure of the flow at the top corners of the domain, where the surfactant-laden surface meets the wall, needs separate asymptotic treatment. We illustrate this at the top left corner, introducing polar coordinates $(r, \theta)$ centred on $(x, y) = (-1, 0)$ with $\theta = 0$ along $y = 0$ (and $r$ increasing in the positive $x$ direction) and $\theta = -\pi/2$ along $x = -1$ (and $r$ increasing in the negative $y$ direction) and seeking separable solutions of the form

$$\hat{\psi}(r, \theta) \approx \operatorname{Re} \left[ \sum_i A_i r^{\Phi_i} f_{\Phi_i}(\theta) \right], \quad \text{for } r \to 0. \qquad (2.14)$$
<table><caption>Table 1. Decay rates predicted as eigenvalues $\tilde{\alpha}_n$ computed numerically from (2.13) compared with $\tilde{\alpha}_n$, which are those computed from eigenmodes using (2.11) for $H = 2$, found using a 4000 × 4000 grid. The relative error provides a measure of global numerical error, suggesting that the values for $\tilde{\alpha}_n$ are accurate up to three significant figures.</caption><thead><tr><th>Even modes</th><th>Eigenvalue &alpha;&#x0303;<sub>n</sub></th><th>Eigenvalue &alpha;&#x0303;<sub>n</sub> from (2.11)</th><th>Relative error (&alpha;&#x0303;<sub>n</sub> - &alpha;&#x0303;<sub>n</sub>)/&alpha;&#x0303;<sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>0.57313</td><td>0.57290</td><td>0.00039</td></tr><tr><td>2</td><td>2.10723</td><td>2.10738</td><td>-0.00007</td></tr><tr><td>3</td><td>3.67670</td><td>3.67710</td><td>-0.00011</td></tr><tr><td>4</td><td>5.24715</td><td>5.24820</td><td>-0.00020</td></tr><tr><td>Odd modes</td><td></td><td></td><td></td></tr><tr><td>1</td><td>1.29728</td><td>1.29725</td><td>0.00002</td></tr><tr><td>2</td><td>2.88997</td><td>2.89018</td><td>-0.00007</td></tr><tr><td>3</td><td>4.46164</td><td>4.46232</td><td>-0.00015</td></tr><tr><td>4</td><td>6.03235</td><td>6.03390</td><td>-0.00026</td></tr></tbody></table>

Here, Re[·] indicates taking the real part, noting that the amplitudes $A_i$ and exponents $\Phi_i$ of local modes of the biharmonic equation may be complex. We will find that the sum in (2.14) is a sum of multiple countable series. For notational simplicity, we will use $\Phi$ to represent exponents, and we will suppress the subscript $n$ on $\hat{\psi}$ and $\alpha$ associated with each eigenmode. To satisfy the governing biharmonic equation (2.5), $\nabla^4 \hat{\psi}(r, \theta) = 0$, and the boundary conditions (2.6b,c), $\hat{\psi}(r, 0) = \hat{\psi}(r, -\pi/2) = \hat{\psi}_{\theta}(r, -\pi/2) = 0$, starting from more general formulas for the $\theta$-dependent functions (Moffatt 1964), we find that

$$f_1(\theta) = \frac{\pi}{2} \sin \theta - \frac{2}{\pi} \theta \cos \theta + \theta \sin \theta \quad (\Phi = 1), \quad (2.15a)$$

$$f_2(\theta) = \cos(2\theta) - \frac{2}{\pi} \sin(2\theta) - \frac{4\theta}{\pi} - 1 \quad (\Phi = 2), \quad (2.15b)$$

and generally for $\Phi > 0$, $\Phi \neq 1, 2$,

$$f_{\Phi}(\theta) = \frac{\cos(\Phi\pi/2)\sin(\Phi\pi/2)}{\sin^2(\Phi\pi/2) - \Phi}(\cos(\Phi\theta) - \cos((\Phi - 2)\theta)) \\ + \frac{\cos^2(\Phi\pi/2) - \Phi + 1}{\sin^2(\Phi\pi/2) - \Phi}\sin(\Phi\theta) + \sin((\Phi - 2)\theta), \quad (2.15c)$$

which in the special case of integer $\Phi$ reduces to

$$f_{\Phi}(\theta) = \begin{cases} \sin(\Phi\theta) + \sin((\Phi - 2)\theta), & \Phi \text{ is an odd integer strictly greater than 1,} \\ \frac{\Phi - 2}{\Phi} \sin(\Phi\theta) + \sin((\Phi - 2)\theta), & \Phi \text{ is an even integer strictly greater than 2.} \end{cases} \quad (2.15d)$$

The final boundary condition, the surfactant boundary condition at the free surface in (2.10): $-\alpha \hat{\psi}_{yy} = \hat{\psi}_{xxy}$ at $y = 0$, is, in polar coordinates,

$$-\alpha \left( \frac{1}{r^2} \hat{\psi}_{\theta\theta} + \frac{1}{r} \hat{\psi}_r \right) = \frac{1}{r} \hat{\psi}_{rr\theta} - \frac{2}{r^2} \hat{\psi}_{r\theta} + \frac{2}{r^3} \hat{\psi}_{\theta}, \quad \text{for } \theta = 0. \quad (2.16)$$
Surfactant cavity flow

Writing the first few terms in the expansion (2.14) for $\hat{\psi}$ in the form

$$ \hat{\psi} = \mathrm{Re} \left[ A_a r^a f_a(\theta) + A_b r^b f_b(\theta) + A_c r^c f_c(\theta) + \dots \right] \quad (2.17) $$

where $a, b, c...$ are complex numbers (representing exponents $\Phi_i$) such that $\mathrm{Re}(a) < \mathrm{Re}(b) < \mathrm{Re}(c)$ and $A_a \neq 0$, we obtain from (2.16) (after multiplying by $r^3$)

$$ \begin{aligned} & -\alpha(A_a r^{a+1} f_a''(0) + A_b r^{b+1} f_b''(0) + A_c r^{c+1} f_c''(0) + \dots) \\ & = (a^2 - 3a + 2) A_a r^a f_a'(0) + (b^2 - 3b + 2) A_b r^b f_b'(0) \\ & \quad + (c^2 - 3c + 2) A_c r^c f_c'(0) + \dots \end{aligned} \quad (2.18) $$

The $r^a$ term is dominant as $r \to 0$, imposing that

$$ (a-1)(a-2)f_a'(0) = 0. \quad (2.19) $$

This equation opens multiple possibilities. The first case $a=1$ corresponds to a solution with non-integrable stress $\tau(r) = \hat{\psi}_{\theta\theta}(\theta=0)/r^2$, and therefore unbounded surfactant concentration as $r \to 0$, so we impose $A_1=0$. This leaves two other cases: $a=2$ and $f_a'(0)=0$. The latter third case yields an infinite set of complex exponents of the type described by Moffatt (1964), each representing a homogeneous local solution, which we will analyse further below. The expansion (2.14) therefore constitutes multiple independent series.

In the second case, with $a=2$, (2.18) becomes

$$ \begin{aligned} & -\alpha(A_2 r^3 f_2''(0) + A_b r^{b+1} f_b''(0) + A_c r^{c+1} f_c''(0) + \dots) \\ & = (b^2 - 3b + 2) A_b r^b f_b'(0) + (c^2 - 3c + 2) A_c r^c f_c'(0) + \dots \end{aligned} \quad (2.20) $$

The dominant balance must be between the $r^b$ and $r^3$ terms since we imposed $\mathrm{Re}(a) < \mathrm{Re}(b)$, hence $b=3$. We then have

$$ -\alpha A_2 f_2''(0) = 2A_3 f_3'(0). \quad (2.21) $$

Using (2.15b) and (2.15d) for $f_2$ and $f_3$, respectively, (2.21) becomes $\alpha A_2 = 2A_3$. The next balance in (2.20) gives $-\alpha A_3 f_3''(0) = 6A_4 f_4'(0)$, however, from (2.15d) we can see that $f_3''(0) = 0$, which implies that $A_4=0$ and this series terminates. The first series contributing to $\hat{\psi}$ is therefore simply

$$ \begin{aligned} & A_2 r^2 \left( f_2(\theta) + \frac{\alpha}{2} r f_3(\theta) \right) \\ & = A_2 r^2 \left( \cos(2\theta) - \frac{2}{\pi} \sin(2\theta) - \frac{4\theta}{\pi} - 1 + \frac{\alpha}{2} r(\sin(3\theta) + \sin\theta) \right). \end{aligned} \quad (2.22) $$

In the third case, setting $f_a'(0) = 0$ in (2.19) (with $a \neq 1, 2$) gives

$$ \frac{\cos^2(a\pi/2) - a + 1}{\sin^2(a\pi/2) - a} + (a-2) = 0, \quad (2.23) $$

so that the complex roots satisfy

$$ \sin^2(a\pi/2) - a(2-a) = 0. \quad (2.24) $$

We label the roots of (2.24) which have non-zero imaginary part as $a_1, a_2, ...$ with $0 < \mathrm{Re}(a_1) < \mathrm{Re}(a_2) < ...$. The roots $a_i$ are independent of $H$, and correspond to
<table><caption>Table 2. Left: approximate numerical values of the first five complex roots $a_i$ of (2.24). These roots are the exponents in Moffatt's series for anti-symmetric Stokes flow (Moffatt 1964) subject to a disturbance at a large distance in the case of right-angle corners. Right: approximate values of the first five coefficients $K_{ij}$ for the dominant root $i = 1$ in the series expansion (2.27) associated with $\hat{\psi}$ and computed using (2.28). All of these quantities are independent of the global parameter $H$.</caption><tr><td>a<sub>1</sub> &approx; 3.7396 + 1.1190i</td><td>K<sub>11</sub> &approx; -0.036666 + 0.076161i</td></tr><tr><td>a<sub>2</sub> &approx; 5.8083 + 1.4639i</td><td>K<sub>12</sub> &approx; 0.046960 + 0.002377i</td></tr><tr><td>a<sub>3</sub> &approx; 7.8451 + 1.6816i</td><td>K<sub>13</sub> &approx; -0.001195 + 0.000349i</td></tr><tr><td>a<sub>4</sub> &approx; 9.8688 + 1.8424i</td><td>K<sub>14</sub> &approx; 0.000198 + 0.000122i</td></tr><tr><td>a<sub>5</sub> &approx; 11.886 + 1.9702i</td><td>K<sub>15</sub> &approx; -0.000002 - 0.000002i</td></tr></table>

the exponents in Moffatt's series (Moffatt 1964) for anti-symmetric Stokes flow in a right-angle corner subject to arbitrary disturbance at a large distance. The first five complex roots are shown in table 2. Each complex root $a$ generates its own asymptotic series of the form (2.17) with $a = a_i$, $b = b_i$, $c = c_i$, etc. where $0 < \text{Re}(a_i) < \text{Re}(b_i) < \dots$, for $i = 1, 2, 3, \dots$. For example, for $i = 1$ (2.18) gives

$$ -\alpha(A_{a_1}r^{a_1+1}f''_{a_1}(0) + A_{b_1}r^{b_1+1}f''_{b_1}(0) + A_{c_1}r^{c_1+1}f''_{c_1}(0) + \dots) \\ = (b_1^2 - 3b_1 + 2)A_{b_1}r^{b_1}f'_{b_1}(0) + (c_1^2 - 3c_1 + 2)A_{c_1}r^{c_1}f'_{c_1}(0) + \dots, \quad (2.25) $$

such that the dominant balance is $b_1 = a_1 + 1$, requiring that

$$ -\alpha A_{a_1} f''_{a_1}(0) = a_1(a_1 - 1)A_{b_1}f'_{b_1}(0), \quad (2.26) $$

with the approximate value of $a_1$ given in table 2. From such relations for $a_1, b_1, c_1, \dots$ we can derive the associated contribution to $\hat{\psi}$, of the form

$$ A_{a_1} \left( r^{a_1} f_{a_1}(\theta) + \alpha K_{11} r^{a_1+1} f_{a_1+1}(\theta) + \alpha^2 K_{12} r^{a_1+2} f_{a_1+2}(\theta) + \dots \right), \quad (2.27) $$

where the first five coefficients $K_{ij}$ related to $a_1$ in (2.27) are shown in table 2. These coefficients can be computed through the recurrence relation

$$ K_{ij} = - \frac{K_{i(j-1)} f''_{a_i + j - 1}(0)}{((a_i + j)^2 - 3(a_i + j) + 2) f'_{a_i + j}(0)}, \quad (2.28) $$

for integers $i \ge 1$ and $j \ge 1$ and with $K_{i0} = 1$.

In summary, the full asymptotic series for the nth eigenmode in the neighbourhood of the contact line, which we originally stated in the form (2.14), is

$$ \hat{\psi}_n(r, \theta) = A_{n2}r^2 \left[ f_2(\theta) + \frac{\alpha_n}{2}rf_3(\theta) \right] + \text{Re} \left[ \sum_{i=1}^{\infty} A_{na_i} \left( \sum_{j=0}^{\infty} \alpha_n^j K_{ij} r^{a_i+j} f_{a_i+j}(\theta) \right) \right]. \quad (2.29) $$

The coefficients $K_{ij}$ can be computed from the recurrence relation (2.28), whilst the coefficients $A_{n2}$ and $A_{na_i}$ must be determined from fitting to numerical solutions. Analysing (2.29) we can notice that $\hat{\psi}_n$ approaches different values as $r \to 0$ for different values of $\theta$, capturing the streamfunction's singular behaviour in this limit. Away from the corner, the complex exponents involved in the second term of (2.29) produce oscillatory
Surfactant cavity flow

behaviour as a function of $r$, such that for the corresponding parameters (2.29) is capable of producing Moffat-type eddies at the top corners, which are observed in the numerical solution for higher-order modes (shown in § 3). This is a local expansion at the top corner, however, global information about the aspect ratio of the domain and the global symmetry or anti-symmetry of $\hat{\psi}_n$ enters into this expression through $\alpha_n$. We note for later reference that the surface shear stress $\tau_n(x)$ has the local form

$$ \hat{\tau}_n(x) = \left. \frac{\partial \hat{u}_n}{\partial y} \right|_{y=0} \approx -4A_{n2} + \text{Re} \left\{ A_{na_1} (1+x)^{a_1-2} f''_{a_1}(0) + \dots \right\} \quad \text{as } x \to -1, \quad (2.30) $$

which has a non-zero asymptotic value at the contact lines $\hat{\tau}_n(-1) = -4A_{n2}$. Moreover, the complex roots $a_i$ imply that the value of the shear stress oscillates as $x \to -1$. Near the contact line $(x, y) = (-1, 0)$, the leading-order surfactant distribution has a finite non-zero gradient, computed from (2.8),

$$ \hat{\Gamma}_n(x) \approx -\frac{8A_{n2}}{\pi\alpha_n} + 4A_{n2}(1+x) + \dots \quad \text{as } x \to -1. \quad (2.31) $$

The leading-order vorticity near the contact line $(x, y) = (-1, 0)$ has the form

$$ \hat{\omega}_n \approx 4A_{n2} \left( 1 + \frac{4}{\pi} \tan^{-1} \frac{y}{1+x} \right) + \dots, \quad (2.32) $$

which shows that the vorticity is multi-valued at the corner owing to the singularity and depending on the angle of approach. From the vorticity, we find that the leading-order pressure locally is

$$ \hat{p}_n \approx -\frac{8A_{n2}}{\pi} \log(y^2 + (1+x)^2) + \dots, \quad (2.33) $$

which shows that the pressure diverges in a logarithmic fashion at the corner. The above expressions give the local expansion at the top left corner of the domain. A similar expansion is then trivially obtained at the top right corner by symmetry. These asymptotic results complement the numerical solutions, less accurate near the contact lines $(x, y) = (\pm 1, 0)$, providing a complete understanding of the effect of the corner singularities on the relevant physical quantities in this problem.

### 3. Results

Figure 2(a) shows how the decay rates $\tilde{\alpha}_n$ of the odd and even modes, computed numerically from the eigenvalue problem (2.13), vary with the depth of the domain $H$. The decay rates become independent of the cavity depth for $H \gtrsim 2$. We can find an exact asymptotic expression for the odd and even decay rates in the limit $H \to 0$ using a lubrication approximation (see Appendix C)

$$ \alpha_n \to \frac{n^2 \pi^2 H}{4}, \quad (\text{odd modes}), \quad \alpha_n \to \frac{(2n-1)^2 \pi^2 H}{16}, \quad (\text{even modes}), \quad (3.1a,b) $$

as shown with dashed lines in figure 2(a) for the first two odd and even modes. The surfactant concentration profiles of the corresponding modes are shown in figure 2(b) (using the same colours as in figure 2a). All the surfactant profiles have a non-zero finite value and slope at the boundaries $x = \pm 1$, as anticipated from (2.31). We compute the dominant singularity strength $-4A_{n2}$ from the slope of the surfactant profile at the boundaries $x = \pm 1$ according to (2.31).
<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Decay rates computed numerically by solving (2.13) (solid lines) as functions of H, and compared with the analytical predictions obtained from lubrication theory (3.1a,b) in the limit $H \to 0$ (dashed lines). (b) Plots of the surfactant concentration profiles (using the same colour code as in (a)) for the first two even and two odd modes for $H = 2$ using a solution for the streamfunction $\tilde{\psi}_n$ calculated numerically using 4000 × 4000 grid points.</figcaption></figure>

The contour plots computed numerically in figure 3(a,b) show the streamfunction and vorticity of the dominant mode (first odd eigenmode $n=1$) for $H=2$. Weak Moffatt eddies can be seen in the lower corners of the cavity $(x, y) = (\pm 1, -2)$. Vorticity contours converge at the upper corners, indicating that $\hat{\omega}$ is multi-valued there, consistent with (2.32). The numerical results of the contour plot of the streamfunction in a deep channel with $H=8$, as shown in figure 3(c), reveals a sequence of recirculations. The strength of the streamfunction decreases rapidly with increasing depth, by approximately three orders of magnitude for the amplitude of the streamfunction between successive cores. In a shallow channel with $H=0.2$ (figure 3d), elongated eddies appear, consistent with predictions of lubrication theory.

Figure 4(a) shows the surface shear stress, computed numerically from the viscous-Marangoni stress condition $\tilde{\tau}(x) = -\tilde{\Gamma}_x(x)$, for the dominant mode (first odd mode, $n=1$) for $H=2$, revealing an oscillatory structure near the contact lines as $x \to \pm 1$, consistent with (2.30). The log-log plot in figure 4(b) reveals in more detail how the stress calculated from the numerical solution matches against the stress found using the asymptotic approximation (2.30). The finite-difference approximation, with a finite grid size, necessarily fails to capture the increasingly short-wavelength oscillations as $x \to -1$ (as plotted with dashed lines when $1+x \le 100\Delta x$), and the asymptotic solution can be expected to fail as $1+x$ becomes too large. However, there is an overlap region (indicated by solid lines for the numerical results), which grows in size with increasing grid resolution, over which the agreement is sufficiently strong to provide confidence in the numerical predictions throughout the rest of the domain. Thus, at the maximum grid resolution ($\Delta x = 1/8000$) the numerical results are close to the asymptotic results for $\log(1+x) \approx -1.8$. We note that the asymptotic results could be made more accurate by including more terms in the series (2.30), which would for instance show variations in the wavelength of the shear stress oscillations as $x \to -1$. However, the dominant odd mode $n=1$ clearly captures the oscillatory behaviour of the shear stress near the corner.

### 3.1. Regularising the corner singularities

The corner singularities can be regularised by adding a small amount of surface diffusion in the problem. In this case the transport equation for the surfactant, the first equation
<figure><img src="image_5.png" /><figcaption>Figure 3. Contour plots of the dominant mode, the first odd eigenmode ($n = 1$), computed from numerical simulations in a square domain ($H = 2$), showing (a) the streamfunction and (b) the vorticity. Similarly, (c) shows numerical results of the contour plots of the streamfunction for the dominant mode (odd mode $n = 1$) in a deep domain ($H = 8$) and (d) a shallow domain ($H = 0.2$).</figcaption></figure>

in (2.6a), modifies to $\Gamma_{t*}^* = -(\Gamma^*u^*)_{x*} + D^*\Gamma_{x*x*}^*$, in dimensional form with $D^*$ the surface diffusivity of the surfactant. Hence, the stress boundary condition in the eigenvalue problem for the perturbation streamfunction (first equation in (2.10)) becomes

$$ \alpha \hat{\psi}_{yy} = -\hat{\psi}_{xxy} + D \hat{\psi}_{xxyy}, \quad (3.2) $$

where $D = D^*\mu^*/(W^*S^*\bar{\Gamma})$. We then specify an additional boundary condition and impose no flux of surfactant, $\Gamma_x = 0$, at the contact lines $(x, y) = (\pm 1, 0)$. Thus, in the
R. Mcnair, O.E. Jensen and J.R. Landel

<figure><img src="image_6.png" /><figcaption>Figure 4. Distribution of the surface shear stress computed from the dominant mode, the first odd mode $n = 1$. (a) Numerical results (red) vs asymptotic results including two terms in (2.30) (blue) with coefficients $4A_{12} = 1.216$ and $A_{1a_1} \approx -0.175e^{0.157i}$ which were locally optimised to fit the numerical solution. (b) Logarithmic plot showing similar results as in (a) vs $\log(1+x)$. This reveals overlap between asymptotics and numerics for different grid resolutions; dashed lines for the numerical results are used for $1+x \le 100\Delta x$, where numerical errors increase. (c) Distributions of the shear stress, calculated numerically and when surface diffusion is included following (3.2). (d) Semi-log plots of the profiles in (c), with the horizontal coordinate scaled by $D$ in the inset. The dashed line in the main graph is the asymptotic value of the shear stress at $x = -1$, i.e. $\tau = -4A_{12} = -1.216$.</figcaption></figure>

presence of weak surface diffusion, surface stress falls abruptly to zero in small boundary layers at the wall for $D \ll 1$ (figure 4c). Increasing $D$ causes the surface stress of the first odd mode $n = 1$ to take a smoother more sinusoidal profile, revealing the impact of surface diffusion on the free surface. The profile of the shear-stress distribution near the contact line is shown in greater detail in figure 4(d) (using a logarithmic spatial scale), demonstrating its adjustment from the constant value $-4A_{n2}$ (value of the shear stress at the corner for $D = 0$, plotted with a dashed line) to zero over a very short length scale. Collapse of these data when $x$ is rescaled by $D$ (inset) provides evidence that weak surface diffusion regularises the singularity over a boundary layer of characteristic length scale $D$.

We assumed initially that a restoring force is present that is sufficiently strong to suppress interfacial deflections by imposing a flat interface and ignoring the normal stress condition. Given the singularity in the pressure at the contact lines (2.33), it is prudent to revisit this boundary condition. Linearising the normal stress condition around $y = 0$, we can state this as $\tilde{p}(x, 0) - p_{ext} - 2\tilde{\psi}_{xy}(x, 0) = \tilde{h}_{xx}$, where the dimensional deflection is $(W*S^*/\gamma_0^*)\tilde{h}$ and $p_{ext}$ is a reference pressure (recall that, for the small surfactant concentration variations considered here, we can assume $S^* \ll \gamma_0^*$). The displacement is computed by imposing $\tilde{h} = 0$ at each contact line, and imposing a volume constraint
Surfactant cavity flow

<figure><img src="image_7.png" /><figcaption>Figure 5. (a) Surface pressure profiles computed numerically for the different values of the diffusion constant *D* given using the colours indicated in figure 4(c). The inset shows a scaled semi-log graph showing the pressure profiles collapsing on to each other in a diffusive boundary layer located for *x* ∈ [-1, 0] as *D* decreases. (b) Non-dimensional interfacial deflections (relative to the length scale W*S*/γ₀*), computed numerically as the leading-order correction to the flat state for the first two odd and even modes.</figcaption></figure>

$\int_{-1}^{1} \tilde{h} \, dx = 0$. The pressure field at the free surface (computed numerically with gauge $\tilde{p}(0, 0) = 0$) is shown in figure 5(a) for the same values of $D$ as used in figure 4(c,d). Collapse of the data (inset in figure 5a) demonstrates that the pressure, like the stress, is regularised over a length scale in $x$ of $O(D)$ at the contact lines, reaching a maximum value of $O(\log(1/D))$. The corresponding interfacial displacement (figure 5b), computed numerically with $D=0$, shows that, despite the strong local forcing, displacements remain bounded. For the first odd mode, for example, there is weak upwelling near each contact line compensated by lowering of the free surface in the middle of the domain.

4. Discussion

Confined gas-liquid interfaces are commonly contaminated by surfactant, deliberately or by trace amounts naturally present in the environment. Here, we have addressed Marangoni-driven spreading of insoluble surfactant, towards an equilibrium state with uniform concentration, taking place across the width of an interface in a channel, when viscosity dominates the flow. Many features of this spreading are diffusive in character, particularly the decomposition of the flow into a set of mutually orthogonal eigenmodes. The modes decay exponentially in time at different rates, with the longest-wave modes being the most long lasting. Despite this benign temporal structure, the dynamic compression of surfactant near each lateral boundary gives eigenmodes more exotic spatial features. The logarithmic pressure singularity near each contact line (2.33) has the potential to generate a measurable surface deflection (figure 5), while the shear stress has a pronounced oscillatory structure (figure 4b).

So far, we have focussed our study to the special case of a single-fluid flow with a pinned contact line forming a $\pi/2$ contact angle with the cavity walls. However, we can relax these assumptions to show that our results extend to a broader class of problems.
In Appendix D, we present a local asymptotic analysis for the case of a two-fluid incompressible Stokes flow, with fluids of arbitrary viscosities and with an arbitrary contact angle between 0 and $\pi$ for the interface which is still assumed locally flat and pinned to the flat walls of the cavity. We show that the structure of the (two-fluid) flow near the contact line and the type of singularities generated by the time-dependent spreading of the surfactants lying at the interface between the two fluids is similar to what we found for the single-fluid flow with $\pi/2$ contact angle. Indeed, the singularity is always associated with a real exponent of 2 in the $r$-dependence of the streamfunction, which leads to the logarithmic pressure singularity and the multi-valuedness of the vorticity near each contact line. The main difference is that the part of the series of the streamfunction associated with real exponents does not terminate at $r^3$ for contact angles different from $\pi/2$ (see (2.22) for the streamfunction with $\pi/2$ contact angle). The oscillatory structure of the shear stress, associated with complex exponents in the $r$-dependence of the streamfunction, depends on the wedge angle, but not on the viscosity ratio between the two fluids. The viscosity ratio only affects the coefficients of higher-order terms in the series for each eigenmode. We note that a local analysis is sufficient to establish the local flow structure and type of singularities near the contact line, which can be generated by an arbitrary far-field disturbance of the surfactant distribution. Such fundamental surfactant flows are found across the range of applications mentioned in the introduction.

We have also shown how, in the case of a single-fluid flow with $\pi/2$ contact angle, the introduction of a small amount of surface diffusion regularises the contact line singularities, leading to prominent changes in stress distributions (figure 4c,d). Surface diffusivity of $7 \times 10^{-10} \text{ m}^2 \text{ s}^{-1}$ for the common surfactant sodium dodecylsulphate (Chang & Franses 1995) translates to a dimensionless diffusion coefficient $D = \mu^*D^*/(S^*W^*)$ below $10^{-8}$ in magnitude, assuming a spreading coefficient $S^* = 10^{-2} \text{ kg s}^{-2}$ over water in a channel of width 1 cm. At such low levels, the impact of the singularities may still be visible, with the pressure maximum near each contact line, proportional to $\log(1/D)$, generating surface deflections resembling that shown in figure 5. We expect that diffusion will act in the same way for the broader class of problems involving two-fluid viscous flows and arbitrary contact angles between 0 and $\pi$.

The singular flow structures near contact lines that we have identified have the potential to contaminate dynamic computations that do not take proper account of these small-scale local structures. In our calculations of spatial eigenmodes, we chose to combine asymptotic analysis with a dense numerical grid to capture the dominant spatial features of the flow. There are a range of alternative strategies that could be deployed, notably singularity removal (Sprittles & Shikhmurzaev 2011; Game, Hodes & Papageorgiou 2019), although (at least) two distinct singularities would require removal in the present problem. As indicated above, singularity regularisation via the introduction of surface diffusion can operate over extremely small length scales when using realistic parameter values, itself presenting a computational challenge. Singularities can be expected to become a particular difficulty in time-dependent studies, when artificial diffusion associated with a computational scheme may generate spurious disturbances propagating outward from the contact lines. Corner singularities can also present convergence difficulties for numerical schemes that represent solutions using (Fadle-Papkovich) eigenfunctions that assume separability of spatial variables (Meleshko 1996, 1997).

The present model rests on numerous assumptions. We have restricted attention to the near-equilibrium state, ignoring nonlinearities associated with large concentration gradients. One benefit of this assumption is that small concentration changes support our assumption of a linearised equation of state, and the assumption that the Marangoni flow is sufficiently weak for restoring forces to maintain a nearly flat interface. We have
Surfactant cavity flow

assumed that the surfactant is insoluble, but anticipate that desorption near the contact line may contribute to regularisation of the singularity. While the planar flow studied here could readily be extended to an axisymmetric geometry, a potentially more interesting avenue, with regard to three-dimensionality, will be to examine how the transverse flow studied here interacts with axial flows along the channel, as may occur in plastrons used in superhydrophobic drag reduction (Peaudecerf et al. 2017), in maze solving by surfactant (Temprano-Coleto et al. 2018) or in microfluidic applications. Furthermore, although we have considered a purely viscous regime, Marangoni spreading can be very rapid. The dimensional decay rates predicted here are 1/α times $W^*\mu^*/\Delta\gamma^*$, where α is an O(1) modal decay rate and $\Delta\gamma^* = (\gamma_0^* - \gamma_c^*)\bar{\Gamma}^*/\Gamma_c^*$ is the surface tension reduction due to the equilibrium surfactant distribution. Taking this as low as $10^{-3}$ kg s$^{-2}$, say, for water in a narrow channel of width 1 mm (and comparable depth), the decay time scale is approximately 1 ms/α, with a Reynolds number of order unity. Despite decaying exponentially with respect to time, the structure of the flow in wider and deeper channels can therefore be expected to be influenced by inertia, at least initially.

In summary, this study shows how the unsteady spreading of a surfactant monolayer along a liquid-liquid or liquid-gas interface, confined by a lateral rigid boundary, can generate singular flow structures near stationary contact lines, including a logarithmically divergent pressure field and an oscillatory shear-stress distribution. Careful treatment of these structures is needed in computational simulations involving dynamic surfactant transport in confined domains.

**Funding.** The authors are grateful to Professor D. Silvester for enlightening discussions relating to the numerical methods used in this study. J.R.L. and O.E.J. acknowledge financial support from EPSRC (grant: EP/T030739/1).

**Declaration of interests.** The authors report no conflict of interest.

Author ORCIDs.

Richard Mcnair https://orcid.org/0000-0002-6539-8203;

Oliver E. Jensen https://orcid.org/0000-0003-0172-6578;

Julien R. Landel https://orcid.org/0000-0003-3159-8749.

Appendix A. The energy budget

The Stokes equation can be written as $\nabla \cdot \boldsymbol{\sigma} = 0$, where the non-dimensional Cauchy stress tensor is $\boldsymbol{\sigma} = -p\boldsymbol{I} + 2\boldsymbol{e}$ and the strain-rate tensor $\boldsymbol{e} = (\nabla\boldsymbol{u} + \nabla\boldsymbol{u}^T)/2$. Thus, for a Stokes flow, $\nabla \cdot (\boldsymbol{u} \cdot \boldsymbol{\sigma}) = \boldsymbol{\sigma} : \nabla\boldsymbol{u} = \boldsymbol{\sigma} : \boldsymbol{e}$, exploiting the symmetry of $\boldsymbol{\sigma}$. Given that $p\boldsymbol{I}$: $\boldsymbol{e}$: trace($\boldsymbol{e}$) $p = 0$ by incompressibility, it follows that $\nabla \cdot (\boldsymbol{u} \cdot \boldsymbol{\sigma}) = 2\boldsymbol{e} : \boldsymbol{e}$. Integrating this over the domain and applying the divergence theorem gives

$$2 \int_V \boldsymbol{e} : \boldsymbol{e} \, dA = \int_{\partial V} \boldsymbol{u} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{n} \, ds, \qquad (A1)$$

where $\partial V$ represents the boundary of V and $ds$ is the curvilinear element along the boundary. For the present perturbation problem, described by (2.5), (2.6b,c) and (2.10) for each eigenmode, (A1) balancing work done by Marangoni forces with viscous dissipation becomes

$$2 \int_V \left( 2\hat{\psi}_{yx}^2 + \frac{1}{2}\hat{\psi}_{yy}^2 + \frac{1}{2}\hat{\psi}_{xx}^2 - \hat{\psi}_{xx}\hat{\psi}_{yy} \right) dA = - \int_{-1}^{1} \left. \left( \hat{\Gamma}_x \hat{\psi}_y \right) \right|_{y=0} \, dx. \qquad (A2)$$

Integration by parts of the left-hand side of (A2) and the use of the no-slip boundary condition on vertical boundaries and no penetration boundary condition on horizontal
boundaries gives

$$2 \int_V (2\hat{\psi}_{yx}^2 + \frac{1}{2}\hat{\psi}_{yy}^2 + \frac{1}{2}\hat{\psi}_{xx}^2 - \hat{\psi}_{xx}\hat{\psi}_{yy}) \, dA = \int_V (\hat{\psi}_{yy}^2 + \hat{\psi}_{xx}^2 + 2\hat{\psi}_{xx}\hat{\psi}_{yy}) \, dA = \int_V \hat{\omega}^2 \, dA, \quad (A3)$$

where $\hat{\omega} = -(\hat{\psi}_{yy} + \hat{\psi}_{xx})$ is the vorticity. Further integration by parts of the right-hand side of (A2) gives

$$-\int_{-1}^{1} (\hat{\Gamma}_{x} \hat{\psi}_{y})|_{y=0} \, dx = \int_{-1}^{1} (\hat{\Gamma} \hat{\psi}_{xy})|_{y=0} \, dx. \quad (A4)$$

Since $\hat{\psi}_{xy} = \alpha \hat{\Gamma}$ from (2.8), we obtain (2.11) for each decay rate $\alpha$ and corresponding eigenmode $\{\hat{\psi}, \hat{\omega}, \hat{\Gamma}\}$.

## Appendix B. Orthogonality of modes

The reciprocal theorem (Masoud & Stone 2019) states that two Stokes flows, with velocity fields and Cauchy stress tensors ($u, \sigma$) and ($u', \sigma'$), satisfy in the same region

$$\int_{\partial V} u \cdot (\sigma' \cdot n) \, dA = \int_{\partial V} u' \cdot (\sigma \cdot n) \, dA. \quad (B1)$$

For the present perturbation problem, described by (2.5), (2.6b,c) and (2.10) for each eigenmode, we have $u = 0$ on the three solid boundaries whilst at the free surface $u_y = -\hat{\Gamma}_x$. Thus, for two distinct modes $m$ and $n$, (B1) implies

$$\int_{-1}^{1} \hat{u}_m \hat{\Gamma}_{n,x} \, dx = \int_{-1}^{1} \hat{u}_n \hat{\Gamma}_{m,x} \, dx. \quad (B2)$$

Integrating by parts and using $\alpha \hat{\Gamma} = -\hat{u}_x$ gives

$$[\hat{u}_m \hat{\Gamma}_n]_{-1}^{1} - \int_{-1}^{1} \alpha_m \hat{\Gamma}_m \hat{\Gamma}_n \, dx = [\hat{u}_n \hat{\Gamma}_m]_{-1}^{1} - \int_{-1}^{1} \alpha_n \hat{\Gamma}_n \hat{\Gamma}_m \, dx. \quad (B3)$$

As the surface velocity is zero at $x = \pm 1$, (B3) becomes

$$(\alpha_m - \alpha_n) \int_{-1}^{1} \hat{\Gamma}_m \hat{\Gamma}_n \, dx = 0, \quad (B4)$$

which shows that the surfactant concentration profiles corresponding to different modes are orthogonal since $\alpha_m \neq \alpha_n$, as stated in (2.12). Equivalently, we note that this result can be derived by exploiting the self-adjointness of (2.5), (2.6b,c) and (2.10).

Numerical computation of the scalar product of the surfactant concentration modes for $1 \le m, n \le 5$ gives

$$\int_{-1}^{1} \tilde{\Gamma}_m \tilde{\Gamma}_n \, dx = A_{m,n} = \begin{pmatrix} 0.20243519 & -0.00000524 & -0.00000370 & -0.00000345 & -0.00000358 \\ -0.00000524 & 0.19832581 & -0.00000394 & -0.00000457 & -0.00000508 \\ -0.00000370 & -0.00000394 & 0.19136202 & -0.00000570 & -0.00000635 \\ -0.00000345 & -0.00000457 & -0.00000570 & 0.19112895 & -0.00000744 \\ -0.00000358 & -0.00000508 & -0.00000635 & -0.00000744 & 0.19096022 \end{pmatrix}. \quad (B5)$$
These results were computed using our numerical scheme presented in § 2.2 with 4000 × 4000 grid points. The non-diagonal elements of the matrix $A_{m,n}$ in (B5) are very small, showing the orthogonality of the modes calculated numerically and the global accuracy of our numerical scheme.

## Appendix C. The thin-film limit

As $H \to 0$, the biharmonic equation (2.5) can be approximated as $\psi_{yyyy} = 0$, therefore the general solution in this limit is given by

$$ \psi = f_1(x)y^3 + f_2(x)y^2 + f_3(x)y + f_4(x). \qquad (C1) $$

The boundary condition $\psi(x, y = 0) = 0$ means $f_4(x) = 0$. The boundary conditions $\psi(x, y = -H) = 0$ and $\psi_y(x, y = -H) = 0$ give

$$ 0 = -f_1(x)H^3 + f_2(x)H^2 - f_3(x)H, \qquad (C2) $$

and

$$ 0 = 3f_1(x)H^2 - 2f_2(x)H + f_3(x). \qquad (C3) $$

Eliminating $f_3$ gives $f_2(x) = 2f_1(x)H$, implying $f_3(x) = f_1(x)H^2$. Hence, $\psi = f_1(x)(y^3 + 2Hy^2 + H^2y)$. The stress boundary condition at the surface $y = 0$ is $\alpha\psi_{yy} = -\psi_{xxy}$, which imposes

$$ 4\alpha H f_1(x) = -f_1''(x)H^2. \qquad (C4) $$

The solution of this ordinary differential equation for $f_1(x)$ is

$$ f_1(x) = C_1 \cos \left( \sqrt{\frac{4\alpha}{H}} x \right) + C_2 \sin \left( \sqrt{\frac{4\alpha}{H}} x \right), \qquad (C5) $$

for some arbitrary integration constants $C_1$ and $C_2$. The odd and even modes are given by setting either of these constants to zero. Let $C_2 = 0$, then applying the boundary conditions $\psi(x = \pm 1) = 0$ at the sidewalls gives $\cos(\sqrt{4\alpha/H}) = 0$, which implies

$$ \alpha_n = \frac{(2n - 1)^2 \pi^2 H}{16}, \qquad (C6) $$

where $n$ can be any positive integer, $n \ge 1$. Equation (C6) is an approximation for the decay rates for the modes even of the streamfunction in $x$ as $H$ becomes small. Similarly, letting $C_1 = 0$ in (C5), and then applying the boundary conditions $\psi(x = \pm 1) = 0$ at the sidewalls gives $\sin(\sqrt{4\alpha/H}) = 0$, which implies

$$ \alpha_n = \frac{n^2 \pi^2 H}{4}, \qquad (C7) $$

where $n$ can be any positive integer, $n \ge 1$. Equation (C7) is an approximation for the decay rates for the modes even of the streamfunction in $x$ as $H$ becomes small. We note finally that the sinusoidal modes (C5) in this long-wave theory predict zero slope of the surfactant concentration at the contact lines, so fail to capture the finite slope predicted by (2.31) and demonstrated in figure 2(b).
R. Mcnair, O.E. Jensen and J.R. Landel

<figure><img src="image_8.png" /><figcaption>Figure 6. (a) Schematic of the local problem near the contact line pinned at O, with surfactant on the interface ($\theta = 0$) between two incompressible fluids with viscosities $\mu_1$ (bottom fluid) and $\mu_2$ (top fluid) and contact angles $\Theta_1$ and $\Theta_2$, respectively. (b) Plot of the real part of admissible exponents for the radial dependence of the streamfunction, calculated from (D4), against the contact angle $\Theta_1$. This gives the exponents in the asymptotic series capturing the behaviour of the fluid as $r \to 0$. Importantly, this shows no admissible exponents with $1 < \text{real}(a_1) < 2$, which means that for any contact angle and viscosity ratio, the nature of the dominant singularity presented in the main text for a single-fluid flow with contact angle $\pi/2$ is generic.</figcaption></figure>

## Appendix D. Asymptotics with arbitrary contact angle and two fluids of arbitrary viscosities

We consider the generalisation of the local analysis of the single-fluid flow with $\pi/2$ contact angle made in §2.3. As depicted in figure 6(a), we now assume a two-fluid incompressible Stokes flow with arbitrary viscosities. We assume that the interface is locally flat near the contact lines and remains pinned to the flat walls of the cavity at the contact-line location at point O, which corresponds to the coordinates $(x, y) = (-1, 0)$ in the original problem described in figure 1. We allow the contact angle $\Theta$ to vary between 0 and $\pi$, thereby relaxing the assumption made in the main part of the text. We find local approximations to the streamfunctions in each fluid, $\psi_{(1)}$ in fluid 1 (bottom fluid), and $\psi_{(2)}$ in fluid 2 (top fluid). These fluids have viscosities $\mu_1$ and $\mu_2$, and contact angles $\Theta_1$ and $\Theta_2 = \pi - \Theta_1$, respectively.

We use a plane polar coordinate system with $r$ being the radial direction from the origin, and $\theta$ the angular coordinate, with the interface along the line $\theta = 0$ (see figure 6a). Similar to the model presented in §2, we formulate the flow problem with the streamfunction, which follows the biharmonic equation (2.5) in each fluid, with no flux and no penetration at the cavity wall, and no penetration at the interface, which is assumed fixed and locally flat. At the interface ($r \ge 0$, $\theta = 0$), we also assume continuity of the tangential velocity, whilst the tangential dynamic boundary condition now becomes, in dimensional form,

$$-\tau \cdot [[\sigma^*]] \cdot n = \tau \cdot \nabla_s^* \gamma^*, \quad (D1)$$

where the jump in tangential stress across the interface is balanced by the surfactant-induced Marangoni stress. The jump bracket is defined as $[[\sigma^*]] = \sigma_2^* - \sigma_1^*$. The stress tensor $\sigma^*$ is assumed Newtonian for both fluids. The unit tangential and normal vectors at the interface follow $\tau = (1, 0)$ and $n = (0, 1)$, in polar coordinates, as depicted in figure 6(a). The surface gradient operator is defined as $\nabla_s^* = \nabla^* \cdot (l - n \otimes n)$, with $\otimes$ the outer product. Similar to before, we non-dimensionalise this problem using (2.3), taking fluid 1 as the reference fluid, then we linearise the surfactant distribution around
Surfactant cavity flow

$\bar{\Gamma}$, and decompose the perturbation for each variable in the form $\hat{\Gamma}(r)e^{-\lambda t}$ for example. Relating all variables to the streamfunction, as done previously in § 2.1, the tangential dynamic boundary condition (D1) becomes for each mode, in polar coordinates,

$$ 
\begin{aligned}
& -\alpha \left( \frac{1}{r^2} \psi_{(1)\theta\theta} + \frac{1}{r} \psi_{(1)r} - \frac{\mu_r}{r^2} \psi_{(2)\theta\theta} - \frac{\mu_r}{r} \psi_{(2)r} \right) \\
& = \frac{1}{r} \psi_{(1)rr\theta} - \frac{2}{r^2} \psi_{(1)r\theta} + \frac{2}{r^3} \psi_{(1)\theta} = \frac{1}{r} \psi_{(2)rr\theta} - \frac{2}{r^2} \psi_{(2)r\theta} + \frac{2}{r^3} \psi_{(2)\theta}, \quad \text{on } \theta = 0,
\end{aligned}
\quad (D2) $$

where $\alpha = \lambda/\bar{\Gamma}$ is the decay rate of each mode, and $\mu_r = \mu_2/\mu_1$ is the viscosity ratio between the two fluids. The first term incorporates the difference in shear stress between the lower and the upper fluid; the middle and final terms describe stretching of the interface. Continuity of the tangential velocity field along the interface requires both fluids to stretch at an equal rate.

We seek separable solutions to the above problem such that the leading-order terms in each series is $\psi_{(1)} = r^{a_1}f_{a_1}(\theta)$ and $\psi_{(2)} = r^{a_2}f_{a_2}(\theta)$ where the functions $f$ are given by (2.15a) to (2.15d). Applying the boundary conditions we find that the constants in the functions $f_{a_1}$ and $f_{a_2}$ depend on the contact angle, whilst the exponents $a_1$ and $a_2$ satisfy

$$ (a_1^2 - 3a_1 + 2) f'_{a_1}(0; \Theta_1) = 0 \quad \text{or} \quad (a_2^2 - 3a_2 + 2) f'_{a_2}(0; -\Theta_2) = 0. \quad (D3) $$

The conditions in (D3) can be satisfied with $a_1 = a_2 = 2$, based on the first bracket in each condition, such that the streamfunction solutions must involve series with exponent equal to 2. However, we ask the question whether taking either $f'_{a_1}(0; \Theta_1) = 0$ or $f'_{a_2}(0; -\Theta_2) = 0$ in (D3) can give us an exponent with real part between 1 and 2. This is important as an exponent less than 2 would give us a stronger corner singularity than that discussed in the main text for the case of a single-fluid flow with a contact angle of $\pi/2$. Exponents with real part less than or equal to 1, from the conditions (D3), are rejected on physical grounds to avoid the radial velocity to diverge or be non-zero as $r \to 0$.

When we impose the condition $f'_{a_1}(0; \Theta_1) = 0$, we find that the exponent must obey the condition which is the same as found in Moffatt's problem for a flow near a corner of angle $\Theta_1$ subject to zero velocity boundary conditions at the boundaries (and similarly for $a_2$) (Moffatt 1964). This condition is

$$ (a_1 - 1) \sin(\Theta_1) = \pm \sin(\Theta_1(a_1 - 1)). \quad (D4) $$

By inspection we can see that $a_1 = 0, 1, 2$ (and similarly for $a_2$) are solutions of (D4), for any value of $\Theta_1$, and any integer is a solution when $\Theta_1 = \pi$. We show that a solution with $1 < \text{real}(a_1) < 2$ cannot exist for $0 \le \Theta_1 \le \pi$ in figure 6(b) where we have plotted the locations of the real parts of the admissible solutions of (D4) against $\Theta_1$. We also note that none of the exponents in the expansions for $\psi_{(1)}$ or $\psi_{(2)}$ depend on the viscosity ratio. Hence, in this problem the exponent in the dominant term in both streamfunctions will be 2 for any contact angle and viscosity ratio, giving the same type of singularities as presented in the main body of the paper for this broader class of problems.

REFERENCES

AFSAR-SIDDIQUI, A.B., LUCKHAM, P.F. & MATAR, O.K. 2003 The spreading of surfactant solutions on thin liquid films. *Adv. Colloid Interface Sci.* **106** (1–3), 183–236.

BAIER, T. & HARDT, S. 2021 Influence of insoluble surfactants on shear flow over a surface in Cassie state at large Péclet numbers. *J. Fluid Mech.* **907**, A3.
BORGAS, M.S. & GROTBERG, J.B. 1988 Monolayer flow on a thin film. *J. Fluid Mech.* **193**, 151–170.

CHANG, C.-H. & FRANSES, E.I. 1995 Adsorption dynamics of surfactants at the air/water interface: a critical review of mathematical models, data, and mechanisms. *Colloids Surf.* (A) **100**, 1–45.

DUSSAUD, A.D., MATAR, O.K. & TROIAN, S.M. 2005 Spreading characteristics of an insoluble surfactant film on a thin liquid layer: comparison between theory and experiment. *J. Fluid Mech.* **544**, 23–51.

FILOCHE, M., TAI, C.-F. & GROTBERG, J.B. 2015 Three-dimensional model of surfactant replacement therapy. *Proc. Natl Acad. Sci. USA* **112** (30), 9287–9292.

FORNBERG, B. 1988 Generation of finite difference formulas on arbitrarily spaced grids. *Maths Comput.* **51** (184), 699–706.

GAME, S.E., HODES, M. & PAPAGEORGIOU, D.T. 2019 Effects of slowly varying meniscus curvature on internal flows in the Cassie state. *J. Fluid Mech.* **872**, 272–307.

GAVER, D.P. III & GROTBERG, J.B. 1992 Droplet spreading on a thin viscous film. *J. Fluid Mech.* **235**, 399–414.

GAVER, D.P. III & GROTBERG, J.B. 1990 The dynamics of a localized surfactant on a thin film. *J. Fluid Mech.* **213**, 127–148.

GROTBERG, J.B., HALPERN, D. & JENSEN, O.E. 1995 Interaction of exogenous and endogenous surfactant: spreading-rate effects. *J. Appl. Physiol.* **78** (2), 750–756.

HUH, C. & SCRIVEN, L.E. 1971 Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. *J. Colloid Interface Sci.* **35** (1), 85–101.

JENSEN, O.E. 1998 The stress singularity in surfactant-driven thin-film flows. Part 2. Inertial effects. *J. Fluid Mech.* **372**, 301–322.

JENSEN, O.E. & GROTBERG, J.B. 1992 Insoluble surfactant spreading on a thin viscous film: shock evolution and film rupture. *J. Fluid Mech.* **240**, 259–288.

JENSEN, O.E. & HALPERN, D. 1998 The stress singularity in surfactant-driven thin-film flows. Part 1. Viscous effects. *J. Fluid Mech.* **372**, 273–300.

JENSEN, O.E. & NAIRE, S. 2006 The spreading and stability of a surfactant-laden drop on a prewetted substrate. *J. Fluid Mech.* **554**, 5–24.

KOVALCHUK, N.M. & SIMMONS, M.J.H. 2021 Surfactant-mediated wetting and spreading: recent advances and applications. *Curr. Opin. Colloid Interface Sci.* **51**, 101375.

KUHLMANN, H.C. & ROMANÒ, F. 2019 The lid-driven cavity. In *Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics* (ed. A. Gelfgat), pp. 233–309. Springer.

LANDEL, J.R., PEAUDECERF, F.J., TEMPRANO-COLETO, F., GIBOU, F., GOLDSTEIN, R.E. & LUZZATTO-FEGIZ, P. 2020 A theory for the slip and drag of superhydrophobic surfaces with surfactant. *J. Fluid Mech.* **883**, A18.

LANDEL, J.R. & WILSON, D.I. 2021 The fluid mechanics of cleaning and decontamination of surfaces. *Annu. Rev. Fluid Mech.* **53**, 147–171.

LIU, Y., PECO, C. & DOLBOW, J. 2019 A fully coupled mixed finite element method for surfactants spreading on thin liquid films. *Comput. Meth. Appl. Mech. Engng* **345**, 429–453.

MANIKANTAN, H. & SQUIRES, T.M. 2020 Surfactant dynamics: hidden variables controlling fluid flows. *J. Fluid Mech.* **892**, P1.

MASOUD, H. & STONE, H.A. 2019 The reciprocal theorem in fluid dynamics and transport phenomena. *J. Fluid Mech.* **879**, P1.

MATAR, O.K. & CRASTER, R.V. 2009 Dynamics of surfactant-assisted spreading. *Soft Matt.* **5** (20), 3801–3809.

MELESHKO, V.V. 1996 Steady stokes flow in a rectangular cavity. *Proc. R. Soc. Lond. A* **452**, 1999–2022.

MELESHKO, V.V. 1997 Biharmonic problem in a rectangle. *Appl. Sci. Res.* **58**, 217–249.

MOFFATT, H.K. 1964 Viscous and resistive eddies near a sharp corner. *J. Fluid Mech.* **18** (1), 1–18.

PEAUDECERF, F.J., LANDEL, J.R., GOLDSTEIN, R.E. & LUZZATTO-FEGIZ, P. 2017 Traces of surfactants can severely limit the drag reduction of superhydrophobic surfaces. *Proc. Natl Acad. Sci. USA* **114** (28), 7254–7259.

SAULEDA, M.L., CHU, H.C.W., TILTON, R.D. & GAROFF, S. 2021 Surfactant driven Marangoni spreading in the presence of predeposited insoluble surfactant monolayers. *Langmuir* **37** (11), 3309–3320.

SCOTT, J.C. 1982 Flow beneath a stagnant film on water: the Reynolds ridge. *J. Fluid Mech.* **116**, 283–296.

SHANKAR, P.N. & DESHPANDE, M.D. 2000 Fluid mechanics in the driven cavity. *Annu. Rev. Fluid Mech.* **32** (1), 93–136.

SPRITTLES, J.E. & SHIKHMURZAEV, Y.D. 2011 Viscous flow in domains with corners: numerical artifacts, their origin and removal. *Comput. Meth. Appl. Mech. Engng* **200** (9–12), 1087–1099.
STETTEN, A.Z., IASELLA, S.V., CORCORAN, T.E., GAROFF, S., PRZYBYCIEN, T.M. & TILTON, R.D. 2018 Surfactant-induced Marangoni transport of lipids and therapeutics within the lung. *Curr. Opin. Colloid Interface Sci.* **36**, 58–69.

STONE, H.A. & LEAL, L.G. 1990 The effects of surfactants on drop deformation and breakup. *J. Fluid Mech.* **220**, 161–186.

SWANSON, E.R., STRICKLAND, S.L., SHEARER, M. & DANIELS, K.E. 2015 Surfactant spreading on a thin liquid film: reconciling models and experiments. *J. Engng Maths* **94** (1), 63–79.

TEMPRANO-COLETO, F., PEAUDECERF, F.J., LANDEL, J.R., GIBOU, F. & LUZZATTO-FEGIZ, P. 2018 Soap opera in the maze: geometry matters in Marangoni flows. *Phys. Rev. Fluids* **3** (10), 100507.

WANG, C. & GUO, Z. 2020 A comparison between superhydrophobic surfaces (SHS) and slippery liquid-infused porous surfaces (SLIPS) in application. *Nanoscale* **12** (44), 22398–22424.

WARNER, M.R.E., CRASTER, R.V. & MATAR, O.K. 2004 Fingering phenomena associated with insoluble surfactant spreading on thin liquid films. *J. Fluid Mech.* **510**, 169–200.