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

Richard Mcnair1<EMAIL_ADDRESS>Oliver E. Jensen1
<EMAIL_ADDRESS>Julien R. Landel 1
<EMAIL_ADDRESS>1Department of Mathematics, University of
Manchester, Oxford Road, Manchester M13 9PL, UK

###### Abstract

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. Linearizing
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 wall 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 regularized by weak surface
diffusion. Their existence highlights the need for careful treatment in
computations of unsteady advection-dominated surfactant transport in confined
domains.

###### keywords:

##  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 et al., 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 (SLIPS) (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 (Baier & Hardt,
2021; Landel et al., 2020; Peaudecerf et al., 2017). Surfactants are
amphiphilic materials that 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 et
al., 2003; Matar & Craster, 2009; Liu et al., 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 (Dussaud
et al., 2005; Gaver III & Grotberg, 1990), 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 lengthscales 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 et al.,
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 et al., 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 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 linearized. 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

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^{*}\leq x^{*}\leq W^{*}$ and $-H^{*}\leq y^{*}\leq 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^{*}\leq x^{*}\leq 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).

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^{*}\leq x^{*}\leq W^{*}$ and
$-H^{*}\leq y^{*}\leq 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
$\mathbf{u}^{*}(\mathbf{x}^{*},t^{*})=(u^{*},v^{*})$ to the pressure
$p^{*}(\mathbf{x}^{*},t^{*})$, with $\mathbf{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^{*}\leq
y^{*}\leq 0$ for the side walls, and $y^{*}=-H^{*}$ for $-W^{*}\leq x^{*}\leq
W^{*}$ for the bottom wall. The free surface $\mathcal{F}$, assumed flat to
leading order, is located at $y^{*}=0$ for $-W^{*}\leq x^{*}\leq W^{*}$. We
balance the Cauchy stress $\boldsymbol{\sigma}^{*}\cdot\mathbf{n}$ at
$\mathcal{F}$ (with the normal unit vector $\mathbf{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

$\displaystyle\boldsymbol{\nabla}^{*}p^{*}$
$\displaystyle=\mu^{*}\nabla^{*2}\mathbf{u}^{*},\quad\boldsymbol{\nabla}^{*}\cdot\mathbf{u}^{*}=0,\quad\Gamma^{*}_{t^{*}}=-\frac{\partial}{\partial
x^{*}}\left(u_{s}^{*}\Gamma^{*}\right),\quad\gamma^{*}=\gamma_{0}^{*}-(\gamma_{0}^{*}-\gamma_{c}^{*})\frac{\Gamma^{*}}{\Gamma_{c}^{*}},$
(2.1)

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.1) are

$\displaystyle\mathbf{u}^{*}$ $\displaystyle=\mathbf{0}\quad\text{ on
}x^{*}=\pm W^{*},\quad y^{*}=-H^{*},$ (2.2a)
$\displaystyle\mathbf{u}^{*}\cdot\mathbf{n}$ $\displaystyle=0\quad\text{ on
}y^{*}=0,$ (2.2b) $\displaystyle\boldsymbol{\sigma}^{*}\cdot\mathbf{n}$
$\displaystyle=-\gamma^{*}\mathbf{n}(\boldsymbol{\nabla}^{*}\cdot\mathbf{n})+\frac{\partial\gamma^{*}}{\partial
x^{*}}\mathbf{t}\quad\text{on}\quad\mathcal{F},$ (2.2c)

with $\mathbf{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 to 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^{*}\boldsymbol{\nabla}^{*}\cdot\mathbf{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

$\mathbf{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^{*}},$
(2.3a) $\mathbf{u}=(u,v)=\frac{\mu^{*}}{S^{*}}(u^{*},v^{*}),\quad
t=\frac{S^{*}t^{*}}{W^{*}\mu^{*}},\quad p=\frac{W^{*}p^{*}}{S^{*}}.$ (2.3b)

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,$ (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$ (2.4b) and
$\mathbf{u}=\mathbf{0},\quad\text{ on }x=\pm 1,\text{ and }y=-H.$ (2.4c)

We introduce a stream function $\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$ (2.5)

subject to

$\displaystyle\Gamma_{t}=-(\Gamma\psi_{y})_{x},\quad\psi_{yy}=-\Gamma_{x},\quad\psi=0\quad\text{on}\quad
y=0,$ (2.6a) $\displaystyle\psi_{x}(\pm 1,y)=\psi(\pm 1,y)=0,$ (2.6b)
$\displaystyle\psi(x,-H)=\psi_{y}(x,-H)=0.$ (2.6c)

The stream function 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 Linearization

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\mathrm{where}\quad\int_{-1}^{1}\hat{\Gamma}(x)\,\mathrm{d}x=0,$ (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 linearization, becoming

$\alpha\hat{\Gamma}=\hat{u}_{x},\quad\text{on }y=0,$ (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.$ (2.9)

The resulting eigenvalue problem for the perturbation stream function 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}\qquad\hat{\psi}=0,\quad\text{on
}y=0.$ (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}\,\mathrm{d}A}{\int_{-1}^{1}\hat{\Gamma}_{n}^{2}\,\mathrm{d}x},$
(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}\,\mathrm{d}x=0,\quad\forall\
m\neq n.$ (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 $\widetilde{\psi}_{n}$ (where a wide tilde
denotes a numerically computed variable) of the unknown perturbation modes of
the stream function, $\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$\times$column ordering convention, the domain
described in figure 1 is discretized 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 stream function of a given mode at a given grid point
$\widetilde{\psi}_{n}(x_{j},y_{i})$, with $3\leq i\leq M$, $3\leq j\leq N$. We
impose $\widetilde{\psi}_{n}(x_{j},y_{i})=0$ at the boundaries of the domain
$j=2$ and $j=N+1$, with $2\leq i\leq M+1$, and $i=2$ and $i=M+1$, with $2\leq
j\leq N+1$, to implement the boundary conditions in (2.6b,c) and (2.10) which
state that the stream function vanishes at all the boundaries. The
$(M-2)\times(N-2)$ grid of unknowns is expressed as a column vector
$\boldsymbol{\psi}$ which is assembled by the vertical concatenation of the
rows of unknowns of $\widetilde{\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 $\boldsymbol{\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 stream function 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 $\boldsymbol{\psi}$ satisfies,

$\mathsf{B}\boldsymbol{\psi}=\alpha\mathsf{C}\boldsymbol{\psi},$ (2.13)

where $\mathsf{B}$ and $\mathsf{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 $\widetilde{\psi}_{n}$ and the associated
numerically-calculated decay rates $\widetilde{\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 $\widetilde{\psi}_{n}$ we compute
numerical approximations of all other quantities of interest for each mode,
such as the surfactant concentration profile $\widetilde{\Gamma}_{n}(x)$, the
surface stress $\widetilde{\tau}_{n}(x)=\widetilde{\psi}_{n,yy}(x,y=0)$, the
velocity field $(\widetilde{u}_{n},\widetilde{v}_{n})(x,y)$, the vorticity
field $\widetilde{\omega}_{n}(x,y)$, and the pressure field
$\widetilde{p}_{n}(x,y)$. The solutions $\widetilde{\psi}_{n}$ are normalised
by requiring
$\max(\widetilde{\Gamma}_{n}(x))-\min(\widetilde{\Gamma}_{n}(x))=1$ for each
$n$.

Even modes | Eigenvalue $\widetilde{\alpha}_{n}$ | Eigenvalue $\breve{\alpha}_{n}$ from (2.11) | Relative error $(\widetilde{\alpha}_{n}-\breve{\alpha}_{n})/\widetilde{\alpha}_{n}$
---|---|---|---
1 | 0.57313 | 0.57290 | 0.00039
2 | 2.10723 | 2.10738 | -0.00007
3 | 3.67670 | 3.67710 | -0.00011
4 | 5.24715 | 5.24820 | -0.00020
Odd modes |  |  |
1 | 1.29728 | 1.29725 | 0.00002
2 | 2.88997 | 2.89018 | -0.00007
3 | 4.46164 | 4.46232 | -0.00015
4 | 6.03235 | 6.03390 | -0.00026
Table 1: Decay rates predicted as eigenvalues $\widetilde{\alpha}_{n}$
computed numerically from (2.13) compared to $\breve{\alpha}_{n}$, which are
those computed from eigenmodes using (2.11) for $H=2$, found using a
$4000\times 4000$ grid. The relative error provides a measure of global
numerical error, suggesting that the values for $\widetilde{\alpha}_{n}$ are
accurate up to three significant figures.

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 $\breve{\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 $\widetilde{\Gamma}_{n}$ are orthogonal with
an absolute error less than $8\times 10^{-6}$ (see (B.5)). 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\Re\left[\sum_{i}A_{i}r^{\Phi_{i}}f_{\Phi_{i}}(\theta)\right],\quad\text{for}\
r\rightarrow 0.$ (2.14)

Here $\Re\left[\cdot\right]$ 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

$\displaystyle f_{1}(\theta)$
$\displaystyle=\frac{\pi}{2}\sin{\theta}-\frac{2}{\pi}\theta\cos{\theta}+\theta\sin{\theta}$
$\displaystyle(\Phi=1),$ (2.15a) $\displaystyle f_{2}(\theta)$
$\displaystyle=\cos{(2\theta)}-\frac{2}{\pi}\sin{(2\theta)}-\frac{4\theta}{\pi}-1$
$\displaystyle(\Phi=2),$ (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)},$
(2.15c) which in the special case of integer $\Phi$ reduces to
$f_{\Phi}(\theta)=\begin{cases}\sin{(\Phi\theta)}+\sin{((\Phi-2)\theta)},&\quad$if
$\Phi$ is an odd integer strictly greater than $1$,$\\\
\frac{\Phi-2}{\Phi}\sin(\Phi\theta)+\sin((\Phi-2)\theta),&\quad$if $\Phi$ is
an even integer strictly greater than $2$.$\end{cases}$ (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.$ (2.16)

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

$\hat{\psi}=\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]$
(2.17)

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

$-\alpha(A_{a}r^{a+1}f^{\prime\prime}_{a}(0)+A_{b}r^{b+1}f^{\prime\prime}_{b}(0)+A_{c}r^{c+1}f^{\prime\prime}_{c}(0)+\dots)\\\
=\left(a^{2}-3a+2\right)A_{a}r^{a}f_{a}^{\prime}(0)+\left(b^{2}-3b+2\right)A_{b}r^{b}f_{b}^{\prime}(0)+\left(c^{2}-3c+2\right)A_{c}r^{c}f_{c}^{\prime}(0)+\dots{}$
(2.18)

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

$(a-1)(a-2)f_{a}^{\prime}(0)=0.$ (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\rightarrow 0$, so we impose $A_{1}=0$. This
leaves two other cases: $a=2$ and $f^{\prime}_{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

$-\alpha(A_{2}r^{3}f^{\prime\prime}_{2}(0)+A_{b}r^{b+1}f^{\prime\prime}_{b}(0)+A_{c}r^{c+1}f^{\prime\prime}_{c}(0)+\dots)=\left(b^{2}-3b+2\right)A_{b}r^{b}f_{b}^{\prime}(0)\\\
+\left(c^{2}-3c+2\right)A_{c}r^{c}f_{c}^{\prime}(0)+\dots{}$ (2.20)

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

$-\alpha A_{2}f_{2}^{\prime\prime}(0)=2A_{3}f_{3}^{\prime}(0).$ (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}^{\prime\prime}(0)=6A_{4}f_{4}^{\prime}(0)$, however from (2.15d) we
can see that $f_{3}^{\prime\prime}(0)=0$, which implies that $A_{4}=0$ and
this series terminates. The first series contributing to $\hat{\psi}$ is
therefore simply

$A_{2}r^{2}\left(f_{2}(\theta)+\frac{\alpha}{2}rf_{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).$
(2.22)

In the third case, setting $f_{a}^{\prime}(0)=0$ in (2.19) (with $a\neq 1$ and
$2$) gives

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

so that the complex roots satisfy

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

We label the roots of (2.24) which have non-zero imaginary part as
$a_{1},a_{2},\dots$ with $0<\Re(a_{1})<\Re(a_{2})<\dots$. The roots $a_{i}$
are independent of $H$, and correspond to 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.

$a_{1}\approx 3.7396+1.1190i\quad$ | $\quad K_{11}\approx-0.036666+0.076161i$
---|---
$a_{2}\approx 5.8083+1.4639i$ | $\quad K_{12}\approx 0.046960+0.002377i$
$a_{3}\approx 7.8451+1.6816i$ | $\quad K_{13}\approx-0.001195+0.000349i$
$a_{4}\approx 9.8688+1.8424i$ | $\quad K_{14}\approx 0.000198+0.000122i$
$a_{5}\approx 11.886+1.9702i$ | $\quad K_{15}\approx-0.000002-0.000002i$
Table 2: Left: Approximate numerical values of the first five complex roots
$a_{i}$ of equation (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$.

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<\Re(a_{i})<\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^{\prime\prime}_{a_{1}}(0)+A_{b_{1}}r^{b_{1}+1}f^{\prime\prime}_{b_{1}}(0)+A_{c_{1}}r^{c_{1}+1}f^{\prime\prime}_{c_{1}}(0)+\dots)=\\\
\left(b_{1}^{2}-3b_{1}+2\right)A_{b_{1}}r^{b_{1}}f_{b_{1}}^{\prime}(0)+\left(c_{1}^{2}-3c_{1}+2\right)A_{c_{1}}r^{c_{1}}f_{c_{1}}^{\prime}(0)+\dots{}$
(2.25)

such that the dominant balance is $b_{1}=a_{1}+1$, requiring that

$-\alpha
A_{a_{1}}f^{\prime\prime}_{a_{1}}(0)=a_{1}(a_{1}-1)A_{b_{1}}f_{b_{1}}^{\prime}(0),$
(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),$
(2.27)

where the first five coefficients $K_{1j}$ 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^{\prime\prime}_{a_{i}+j-1}(0)}{((a_{i}+j)^{2}-3(a_{i}+j)+2)f^{\prime}_{a_{i}+j}(0)},$
(2.28)

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

In summary, the full asymptotic series for the $n$th 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]+\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].$
(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 stream function’s singular behaviour in this limit.
Away from the corner, the complex exponents involved in the second term of
(2.29) produce oscillatory 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}+\Re\left\\{A_{na_{1}}(1+x)^{a_{1}-2}f_{a1}^{\prime\prime}(0)+\dots\right\\}\quad\mathrm{as}\quad
x\rightarrow-1,$ (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\mathrm{as}\quad
x\rightarrow-1.$ (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,$ (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,$ (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 $\widetilde{\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\gtrapprox 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},\ \text{(odd
modes)},\quad\alpha_{n}\to\frac{(2n-1)^{2}\pi^{2}H}{16}\ \text{(even modes)},$
(3.1)

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).

(a)

(b)

Figure 2: (a) Decay rates computed numerically by solving (2.13) (solid lines)
as functions of $H$, and compared to the analytical predictions obtained from
lubrication theory (3.1) 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 $\widetilde{\psi}_{n}$ calculated numerically using $4000\times
4000$ gridpoints.

The contour plots computed numerically in figure 4(a,b) show the stream
function 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 multivalued there, consistent with (2.32).
The numerical results of the contour plot of the stream function in a deep
channel with $H=8$, as shown in figure 4(c), reveals a sequence of
recirculations. The strength of the stream function decreases rapidly with
increasing depth, by approximately three orders of magnitude for the amplitude
of the stream function between successive cores. In a shallow channel with
$H=0.2$ (figure 4d), elongated eddies appear, consistent with predictions of
lubrication theory.

(a)

(b)

(c)

(d)

Figure 4: 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 stream function and (b) the vorticity. Similarly, (c) shows numerical
results of the contour plots of the stream function for the dominant mode (odd
mode $n=1$) in a deep domain $(H=8)$ and (d) a shallow domain $(H=0.2)$.

Figure 5(a) shows the surface shear stress, computed numerically from the
viscous-Marangoni stress condition
$\widetilde{\tau}(x)=-\widetilde{\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\rightarrow\pm 1$, consistent with (2.30). The log–log plot
in figure 5(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\rightarrow-1$ (as plotted with dashed lines when $1+x\leq
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.

(a)

(b)

(c)

(d)

Figure 5: Distribution of the surface shear stress computed from the dominant
mode, the first odd mode $n=1$. (a) Numerical results (red) versus 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)
versus $\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\leq 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$

###  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 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 stream function (first equation in (2.10)) becomes

$\alpha\hat{\psi}_{yy}=-\hat{\psi}_{xxy}+D\hat{\psi}_{xxyy},$ (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 presence of weak surface
diffusion, surface stress falls abruptly to zero in small boundary layers at
the wall for ${D}\ll 1$ (figure 5c). 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 5(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 this 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}$.

(a)

(b)

Figure 6: (a) Surface pressure profiles computed numerically for the different
values of the diffusion constant $D$ given using the colours indicated in
figure 5(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\in[-1,0]$ as $D$ decreases. (b) Non-dimensional interfacial deflections
(relative to the length-scale $W^{*}S^{*}/\gamma_{0}^{*}$), computed
numerically as the leading-order correction to the flat state for the first
two odd and even modes.

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.
Linearizing the normal stress condition around $y=0$, we can state this as
$\widetilde{p}(x,0)-p_{ext}-2\widetilde{\psi}_{xy}(x,0)=\widetilde{h}_{xx}$,
where the dimensional deflection is $(W^{*}S^{*}/\gamma_{0}^{*})\widetilde{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
$\widetilde{h}=0$ at each contact line, and imposing a volume constraint
$\int_{-1}^{1}\widetilde{h}\,\mathrm{d}x=0$. The pressure field at the free
surface (computed numerically with gauge $\widetilde{p}(0,0)=0$) is shown in
figure 6(a) for the same values of D as used in figure 5(c,d). Collapse of the
data (inset in figure 6a) 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 6b), 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 6), while the shear stress has a pronounced oscillatory
structure (figure 5b).

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 5c,d). Surface diffusivity of $7\times
10^{-10}\,$m2/s for the common surfactant sodium dodecylsulfate (SDS) (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}\,$kg/s2 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 6.
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 et al., 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 lengthscales 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 linearized equation of
state, and the assumption that the Marangoni flow is sufficiently weak for
restoring forces to maintain a nearly flat interface. We have 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 superhydrophic 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/\alpha$ times
$W^{*}\mu^{*}/\overline{\Delta\gamma}^{*}$, where $\alpha$ is an $O(1)$ modal
decay rate and
$\overline{\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/s2, say, for water in a narrow
channel of width 1 mm (and comparable depth), the decay timescale is
approximately 1 ms$/\alpha$, 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.

## Acknowledgements

The authors are grateful to Prof. David Silvester for enlightening discussions
relating to the numerical methods used in this study.

Declaration of Interests. The authors report no conflict of interest.

##  Appendix A The energy budget

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

$2\int_{V}\mathsf{e:e}\ \mathrm{d}A=\int_{\partial
V}\mathbf{u}\cdot\boldsymbol{\sigma}\cdot{\mathbf{n}}\ \mathrm{d}s,$ (A.1)

where $\partial V$ represents the boundary of $V$ and $\mathrm{d}s$ is the
curvilinear element along the boundary. For the present perturbation problem,
described by equations (2.5), (2.6b,c) and (2.10) for each eigenmode, equation
(A.1) 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)\
\mathrm{d}A=-\int_{-1}^{1}\left.\left(\hat{\Gamma}_{x}\hat{\psi}_{y}\right)\right|_{y=0}\
\mathrm{d}x.$ (A.2)

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

$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)\
\mathrm{d}A=\int_{V}\left(\hat{\psi}_{yy}^{2}+\hat{\psi}_{xx}^{2}+2\hat{\psi}_{xx}\hat{\psi}_{yy}\right)\
\mathrm{d}A=\int_{V}\hat{\omega}^{2}\ \mathrm{d}A,$ (A.3)

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

$-\int_{-1}^{1}\left.\left(\hat{\Gamma}_{x}\hat{\psi}_{y}\right)\right|_{y=0}\
\mathrm{d}x=\int_{-1}^{1}\left.\left(\hat{\Gamma}\hat{\psi}_{xy}\right)\right|_{y=0}\
\mathrm{d}x.$ (A.4)

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
$(\mathbf{u},\boldsymbol{\sigma})$ and
$(\mathbf{u}^{\prime},\boldsymbol{\sigma}^{\prime})$, satisfy in the same
region

$\int_{\partial
V}\mathbf{u}\cdot(\boldsymbol{\sigma}^{\prime}\cdot\mathbf{n})\,\mathrm{d}A=\int_{\partial
V}\mathbf{u}^{\prime}\cdot(\boldsymbol{\sigma}\cdot\mathbf{n})\,\mathrm{d}A$
(B.1)

For the present perturbation problem, described by equations (2.5), (2.6b,c)
and (2.10) for each eigenmode, we have $\mathbf{u}=\mathbf{0}$ on the three
solid boundaries whilst at the free surface $u_{y}=-\Gamma_{x}$. Thus, for two
distinct modes $m$ and $n$, (B.1) implies

$\int_{-1}^{1}\hat{u}_{m}\hat{\Gamma}_{n,x}\,\mathrm{d}x=\int_{-1}^{1}\hat{u}_{n}\hat{\Gamma}_{m,x}\,\mathrm{d}x.$
(B.2)

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

$\left[\hat{u}_{m}\hat{\Gamma}_{n}\right]_{-1}^{1}-\int_{-1}^{1}\alpha_{m}\hat{\Gamma}_{m}\hat{\Gamma}_{n}\mathrm{d}x=\left[\hat{u}_{n}\hat{\Gamma}_{m}\right]_{-1}^{1}-\int_{-1}^{1}\alpha_{n}\hat{\Gamma}_{n}\hat{\Gamma}_{m}\mathrm{d}x.$
(B.3)

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

$(\alpha_{m}-\alpha_{n})\int_{-1}^{1}\hat{\Gamma}_{m}\hat{\Gamma}_{n}\,\mathrm{d}x=0,$
(B.4)

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\leq m,n\leq 5$ gives

$\int_{-1}^{1}\widetilde{\Gamma}_{m}\widetilde{\Gamma}_{n}\
\mathrm{d}x=\mathbf{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}.$
(B.5)

These results were computed using our numerical scheme presented in section
2.2 with $4000\times 4000$ gridpoints. The non-diagonal elements of the matrix
$\mathbf{A}_{m,n}$ in (B.5) 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\rightarrow 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).$ (C.1)

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,$ (C.2)

and

$0=3f_{1}(x)H^{2}-2f_{2}(x)H+f_{3}(x).$ (C.3)

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^{2}y)$. The
stress boundary condition at the surface $y=0$ is
$\alpha\psi_{yy}=-\psi_{xxy}$, which imposes

$4\alpha Hf_{1}(x)=-f_{1}^{\prime\prime}(x)H^{2}.$ (C.4)

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)},$
(C.5)

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 side walls
gives $\cos{\left(\sqrt{{4\alpha}/{H}}\right)}=0$, which implies

$\alpha_{n}=\frac{(2n-1)^{2}\pi^{2}H}{16},$ (C.6)

where $n$ can be any positive integer, $n\geq 1$. Equation (C.6) is an
approximation for the decay rates for the modes even of the stream function in
$x$ as $H$ becomes small. Similarly, letting $C_{1}=0$ in (C.5), and then
applying the boundary conditions $\psi(x=\pm 1)=0$ at the side walls gives
$\sin{\left(\sqrt{\frac{4\alpha}{H}}\right)}=0$, which implies

$\alpha_{n}=\frac{n^{2}\pi^{2}H}{4},$ (C.7)

where $n$ can be any positive integer, $n\geq 1$. Equation (C.7) is an
approximation for the decay rates for the modes even of the stream function in
$x$ as $H$ becomes small. We note finally that the sinusoidal modes (C.5) 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).

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

(a)

(b)

Figure 7: (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 stream
function, calculated from (D.4), 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.

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 7(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 stream functions 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 7a). Similar to the model presented in §2, we
formulate the flow problem with the stream function, 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\geq 0,\theta=0)$, we also assume
continuity of the tangential velocity, whilst the tangential dynamic boundary
condition now becomes, in dimensional form,

$-{\boldsymbol{\tau}}\cdot\llbracket\boldsymbol{\sigma}^{*}\rrbracket\cdot\mathbf{n}={\boldsymbol{\tau}}\cdot{\bf\nabla}^{*}_{s}\gamma^{*}$
(D.1)

where the jump in tangential stress across the interface is balanced by the
surfactant-induced Marangoni stress. The jump bracket is defined as
$\llbracket\boldsymbol{\sigma}^{*}\rrbracket=\boldsymbol{\sigma}_{2}^{*}-\boldsymbol{\sigma}_{1}^{*}$.
The stress tensor $\boldsymbol{\sigma}^{*}$ is assumed Newtonian for both
fluids. The unit tangential and normal vectors at the interface follow
${\boldsymbol{\tau}}=(1,0)$ and ${\bf n}=(0,1)$, in polar coordinates, as
depicted in figure 7(a). The surface gradient operator is defined as
${\bf\nabla}^{*}_{s}={\bf\nabla}^{*}\cdot(\mathsf{I}-{\bf n}\otimes{\bf n})$,
with $\otimes$ the outer product. Similar to before, we non-dimensionalize
this problem using (2.3b), taking fluid 1 as the reference fluid, then we
linearize the surfactant distribution around $\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 (D.1) becomes for each mode,
in polar coordinates,

$-\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},\qquad\text{on
}\theta=0.$ (D.2)

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

$\left(a_{1}^{2}-3a_{1}+2\right)f_{a_{1}}^{\prime}(0;\Theta_{1})=0\qquad\text{or}\qquad\left(a_{2}^{2}-3a_{2}+2\right)f_{a_{2}}^{\prime}(0;-\Theta_{2})=0.$
(D.3)

The conditions in (D.3) can be satisfied with $a_{1}=a_{2}=2$, based on the
first bracket in each condition, such that the stream function solutions must
involve series with exponent equal to $2$. However, we ask the question
whether taking either $f_{a_{1}}^{\prime}(0;\Theta_{1})=0$ or
$f_{a_{2}}^{\prime}(0;-\Theta_{2})=0$ in (D.3) 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 (D.3),
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}}^{\prime}(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))}.$ (D.4)

By inspection we can see that $a_{1}=0,1,2$ (and similarly for $a_{2}$) are
solutions of (D.4), 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\leq\Theta_{1}\leq\pi$ in figure
7(b) where we have plotted the locations of the real parts of the admissible
solutions of (D.4) 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 stream functions 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 et al. (2003) 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 & Hardt (2021) 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.
  * Borgas & Grotberg (1988) Borgas, M S & Grotberg, J B 1988 Monolayer flow on a thin film. J. Fluid Mech. 193, 151–170.
  * Chang & Franses (1995) 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: Physicochem. Eng. Aspects 100, 1–45.
  * Dussaud et al. (2005) 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 et al. (2015) Filoche, M, Tai, C-F & Grotberg, J B 2015 Three-dimensional model of surfactant replacement therapy. Proc. Nat. Acad. Sci. 112 (30), 9287–9292.
  * Fornberg (1988) Fornberg, Bengt 1988 Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of computation 51 (184), 699–706.
  * Game et al. (2019) 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 & Grotberg (1992) Gaver, D P III & Grotberg, J B 1992 Droplet spreading on a thin viscous film. J. Fluid Mech. 235, 399–414.
  * Gaver III & Grotberg (1990) Gaver III, D P & Grotberg, J B 1990 The dynamics of a localized surfactant on a thin film. J. Fluid Mech. 213, 127–148.
  * Grotberg et al. (1995) 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 & Scriven (1971) 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 (1998) Jensen, O E 1998 The stress singularity in surfactant-driven thin-film flows. part 2. inertial effects. J. Fluid Mech. 372, 301–322.
  * Jensen & Grotberg (1992) 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 & Halpern (1998) 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 & Naire (2006) 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 & Simmons (2021) Kovalchuk, N M & Simmons, M J H 2021 Surfactant-mediated wetting and spreading: recent advances and applications. Curr. Op. Colloid Interface Sci. 51, 101375\.
  * Kuhlmann & Romanò (2019) Kuhlmann, H C & Romanò, F 2019 The lid-driven cavity. Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics pp. 233–309.
  * Landel et al. (2020) 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.
  * Landel & Wilson (2021) Landel, J R & Wilson, D I 2021 The fluid mechanics of cleaning and decontamination of surfaces. Ann. Rev. Fluid Mech. 53.
  * Liu et al. (2019) Liu, Y, Peco, C & Dolbow, J 2019 A fully coupled mixed finite element method for surfactants spreading on thin liquid films. Comp. Meth. Appl. Mech. Eng. 345, 429–453.
  * Manikantan & Squires (2020) Manikantan, Harishankar & Squires, Todd M 2020 Surfactant dynamics: hidden variables controlling fluid flows. J. Fluid Mech. 892.
  * Masoud & Stone (2019) Masoud, H. & Stone, H. A. 2019 The reciprocal theorem in fluid dynamics and transport phenomena. J. Fluid Mech. 879, P1.
  * Matar & Craster (2009) Matar, OK & Craster, RV 2009 Dynamics of surfactant-assisted spreading. Soft Matter 5 (20), 3801–3809.
  * Meleshko (1996) Meleshko, V. V. 1996 Steady stokes flow in a rectangular cavity. Proc. Roy. Soc. A: Math. Phys. Eng. Sci. 452, 1999–2022.
  * Meleshko (1997) Meleshko, V. V. 1997 Biharmonic problem in a rectangle. Appl. Sci. Res. 58, 217–249.
  * Moffatt (1964) Moffatt, H. K. 1964 Viscous and resistive eddies near a sharp corner. J. Fluid Mech. 18 (1), 1–18.
  * Peaudecerf et al. (2017) 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. Nat. Acad. Sci. 114 (28), 7254–7259.
  * Sauleda et al. (2021) 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 (1982) Scott, J C 1982 Flow beneath a stagnant film on water: the Reynolds ridge. J. Fluid Mech. 116, 283–296.
  * Shankar & Deshpande (2000) Shankar, P. N. & Deshpande, M. D. 2000 Fluid mechanics in the driven cavity. Ann. Rev. Fluid Mech. 32 (1), 93–136.
  * Sprittles & Shikhmurzaev (2011) Sprittles, JE & Shikhmurzaev, YD 2011 Viscous flow in domains with corners: Numerical artifacts, their origin and removal. Comp. Meth. Appl. Mech. Eng. 200 (9-12), 1087–1099.
  * Stetten et al. (2018) 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. Op. Colloid Interface Sci. 36, 58–69.
  * Stone & Leal (1990) Stone, H A & Leal, L G 1990 The effects of surfactants on drop deformation and breakup. J. Fluid Mech. 220, 161–186.
  * Swanson et al. (2015) Swanson, E R, Strickland, S L, Shearer, M & Daniels, K E 2015 Surfactant spreading on a thin liquid film: reconciling models and experiments. J. Eng. Math. 94 (1), 63–79.
  * Temprano-Coleto et al. (2018) 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 & Guo (2020) 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 et al. (2004) 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.