1. Introduction
Surfactants play a crucial role in a variety of natural and industrial flows (Kovalchuk & Simmons Reference Kovalchuk and Simmons2021), including cleaning and decontamination (Landel & Wilson Reference Landel and Wilson2021), transport in lung airways (Filoche, Tai & Grotberg Reference Filoche, Tai and Grotberg2015; Stetten et al. Reference Stetten, Iasella, Corcoran, Garoff, Przybycien and Tilton2018) and microfluidic applications (Temprano-Coleto et al. Reference Temprano-Coleto, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2018). Understanding how surfactants equilibrate near contact lines is also important for flows over superhydrophobic surfaces, and slippery-liquid-infused porous surfaces (Wang & Guo Reference Wang and Guo2020). 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. Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020; Baier & Hardt Reference Baier and Hardt2021; Peaudecerf et al. Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017). 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 Reference Manikantan and Squires2020). 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 Reference Afsar-Siddiqui, Luckham and Matar2003; Matar & Craster Reference Matar and Craster2009; Liu, Peco & Dolbow Reference Liu, Peco and Dolbow2019). 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 Reference Gaver and Grotberg1990; Dussaud, Matar & Troian Reference Dussaud, Matar and Troian2005), that is captured within lubrication theory as a kinematic shock (Borgas & Grotberg Reference Borgas and Grotberg1988). 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 Reference Jensen and Halpern1998). Out-of-plane displacement of the free surface may be suppressed by surface tension or gravity (Gaver & Grotberg Reference Gaver and Grotberg1992; Jensen & Grotberg Reference Jensen and Grotberg1992). 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 Reference Scott1982; Jensen Reference Jensen1998). Spreading on thin films may also be accompanied by dramatic secondary fingering instabilities (Warner, Craster & Matar Reference Warner, Craster and Matar2004; Jensen & Naire Reference Jensen and Naire2006; Liu et al. Reference Liu, Peco and Dolbow2019), 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. Reference Swanson, Strickland, Shearer and Daniels2015); 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 Reference Grotberg, Halpern and Jensen1995; Sauleda et al. Reference Sauleda, Chu, Tilton and Garoff2021). 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 Reference Huh and Scriven1971), 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 
${\rm \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 
${\rm \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 
${\rm \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 Reference Shankar and Deshpande2000), benchmark problems for computational fluid dynamics, for which corner singularities are a recognised computational challenge (Kuhlmann & Romanò Reference Kuhlmann and 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^* \leqslant x^* \leqslant W^*$ and 
$-H^*\leqslant y^* \leqslant 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 
$\boldsymbol {u}^*(\boldsymbol {x}^*,t^*)=(u^*,v^*)$ to the pressure 
$p^*(\boldsymbol {x}^*,t^*)$, with 
$\boldsymbol {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^*\leqslant y^*\leqslant 0$ for the sidewalls, and 
$y^*=-H^*$ for 
$-W^*\leqslant x^*\leqslant W^*$ for the bottom wall. The free surface 
$\mathcal {F}$, assumed flat to leading order, is located at 
$y^*=0$ for 
$-W^*\leqslant x^*\leqslant W^*$. We balance the Cauchy stress 
$\boldsymbol {\sigma }^*\boldsymbol {\cdot } \boldsymbol {n}$ at 
$\mathcal {F}$ (with the normal unit vector 
$\boldsymbol {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 
$\varGamma ^*$ 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 Reference Stone and Leal1990). The governing equations are therefore
\begin{align} \boldsymbol{\nabla}^* p^* = \mu^*\nabla^{*2} \boldsymbol{u}^*, \quad \boldsymbol{\nabla}^*\boldsymbol{\cdot} \boldsymbol{u}^* = 0, \quad \varGamma^*_{t^*} ={-} \frac{\partial}{\partial x^*} \left(u_s^*\varGamma^* \right), \quad \gamma^* = \gamma_0^* - (\gamma_0^* - \gamma_c^*) \frac{\varGamma^*}{\varGamma_c^*}, \end{align}
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 
$\varGamma _c^*$. The boundary conditions associated with (2.1a–d) are
\begin{gather} \boldsymbol{u}^* = \boldsymbol{0} \quad \text{on } x^*={\pm} W^*, \quad y^*={-}H^*, \end{gather}
\begin{gather}\boldsymbol{u}^* \boldsymbol{\cdot} \boldsymbol{n} = 0 \quad \text{on } y^*=0, \end{gather}
\begin{gather}\boldsymbol{\sigma}^* \boldsymbol{\cdot} \boldsymbol{n} ={-}\gamma^* \boldsymbol{n} (\boldsymbol{\nabla}^*\boldsymbol{\cdot} \boldsymbol{n}) + \frac{\partial \gamma^*}{\partial x^*}\boldsymbol{t} \quad \text{on} \ \mathcal{F}, \end{gather}
with 
$\boldsymbol {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^*\boldsymbol {\nabla }^*\boldsymbol {\cdot } \boldsymbol {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}$.
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^*\leqslant x^* \leqslant W^*$ and 
$-H^*\leqslant y^* \leqslant 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^*\leqslant x^* \leqslant 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).
Using the length scale 
$W^*$, velocity scale 
$S^*/\mu ^*$ and pressure scale 
$S^*/W^*$, we relate dimensional starred variables to their dimensionless counterparts via
\begin{gather} \boldsymbol{x}= (x,y)= \left(\frac{x^*}{W^*},\frac{y^*}{W^*} \right), \quad H=\frac{H^*}{W^*}, \quad \varGamma = \frac{\varGamma^*}{\varGamma_c^*}, \quad \gamma = \frac{\gamma^*-\gamma_c^*}{S^*}, \end{gather}
\begin{gather}\boldsymbol{u}=(u,v)=\frac{\mu^* }{ S^*}(u^*,v^*), \quad t = \frac{S^* t^* }{W^* \mu^*}, \quad p = \frac{ W^* p^*}{S^*}. \end{gather}After rescaling, the governing equations become, in the bulk,
\begin{equation} \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 - \varGamma, \end{equation}with, at the free surface,
\begin{equation} {\varGamma}_{ t} ={-} (\varGamma u)_x , \quad u_y ={-} \varGamma_x, \quad v = 0, \quad \text{on } y = 0 \end{equation}and
\begin{equation} \boldsymbol{u} = \boldsymbol{0}, \quad \text{on } x=\pm1, \text{ and } y ={-}H. \end{equation}
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
\begin{equation} \nabla^4 \psi = 0, \end{equation}subject to
\begin{gather} \varGamma_t={-}(\varGamma \psi_y)_x,\quad \psi_{yy} ={-}\varGamma_{x}, \quad \psi=0 \quad \text{on} \ y=0, \end{gather}
\begin{gather}\psi_x({\pm} 1,y) = \psi({\pm} 1,y) = 0, \end{gather}
\begin{gather}\psi(x,-H)=\psi_y(x,-H) = 0. \end{gather}
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 
$\varGamma$, given a surface velocity 
$\psi _y$. Therefore writing 
$\varGamma$ as the sum of endogenous and exogenous components, 
$\varGamma =\varGamma _1+\varGamma _2$ say, the two components satisfy the same transport equation 
$\varGamma _{it}=-(\varGamma _i\psi _y)_x$ (
$i=1,2$), allowing the evolution of the components to be tracked individually if necessary (Grotberg et al. Reference Grotberg, Halpern and Jensen1995). 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 {\varGamma }$ 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
\begin{equation} \varGamma(x,t) = \bar{\varGamma} + \hat{\varGamma}(x)\textrm{e}^{-\lambda t}, \quad \mathrm{where}\ \int_{{-}1}^1 \hat{\varGamma}(x)\,\mathrm{d}\kern0.06em x=0, \end{equation}
assuming the same time dependence for other variables, for example 
$u(x,y,t) = \hat {u}(x,y)\textrm {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
\begin{equation} \alpha \hat{\varGamma} = \hat{u}_{x}, \quad \text{on}\ y=0, \end{equation}
where 
$\alpha = \lambda /\bar {\varGamma }$. Equation (2.8) is valid for small surfactant concentration perturbation, or at large times since we assume 
$\vert \hat {\varGamma }(x)\vert \textrm {e}^{-\lambda t}\ll \bar {\varGamma }$. Equation (2.8) may be combined with the stress condition 
$\varGamma _x = - u_y$ to give the homogeneous boundary condition
\begin{equation} {-}\alpha \hat{u}_{y} = \hat{u}_{xx}, \quad \text{on } y=0. \end{equation}
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
\begin{equation} \alpha \hat{\psi}_{yy} ={-}\hat{\psi}_{xxy} \quad \hat{\psi}=0, \quad \text{on } y=0. \end{equation}
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 {\varGamma }_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 
$\boldsymbol {\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
\begin{equation} \alpha_n=\frac{\displaystyle\int_V \hat{\omega}_n^2 \,\mathrm{d}A}{\displaystyle\int_{{-}1}^1 \hat{\varGamma}_n^2 \, \mathrm{d}\kern0.06em x}, \end{equation}
where 
$V$ is the rectangular flow domain (see figure 1). Application of the reciprocal theorem (Masoud & Stone Reference Masoud and Stone2019) (see Appendix B) yields the orthogonality condition
\begin{equation} \int_{{-}1}^1 \hat{\varGamma}_m \hat{\varGamma}_n \,\mathrm{d}\kern0.06em x=0, \quad \forall \ m\neq n. \end{equation}
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 {\varGamma }$ 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. Reference Grotberg, Halpern and Jensen1995).
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
$\times$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 Reference Fornberg1988), to approximate the biharmonic operator in (2.5) on the grid of unknowns. These unknowns correspond to the unknown values of the perturbation streamfunction of a given mode at a given grid point 
$\tilde {\psi }_n(x_j,y_i)$, with 
$3\leqslant i\leqslant M$, 
$3\leqslant j\leqslant 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\leqslant i\leqslant M+1$, and 
$i=2$ and 
$i=M+1$, with 
$2\leqslant j\leqslant 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 
$\boldsymbol {\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 
$\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 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 
$\boldsymbol {\psi }$ satisfies,
\begin{equation} \boldsymbol{\mathsf{B}} \boldsymbol{\psi}= \alpha\boldsymbol{\mathsf{C}}\boldsymbol{\psi}, \end{equation}
where 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\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 
$\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 {\varGamma }_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 {\varGamma }_n(x)) - \min (\tilde {\varGamma }_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 
$\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 
$\tilde {\varGamma }_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)$.
Decay rates predicted as eigenvalues 
$\tilde {\alpha }_n$ computed numerically from (2.13) compared with 
$\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 
$\tilde {\alpha }_n$ are accurate up to three significant figures.
2.3. Corner asymptotics
We anticipate from the outset the appearance of Moffatt vortices (Moffatt Reference Moffatt1964) 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 =-{\rm \pi} /2$ along 
$x=-1$ (and 
$r$ increasing in the negative 
$y$ direction) and seeking separable solutions of the form
\begin{equation} \hat{\psi}(r,\theta) \approx \textrm{Re}\left[ \sum_i A_i r^{\varPhi_i} f_{\varPhi_i}(\theta) \right], \quad \text{for} \ r\rightarrow 0. \end{equation}
Here, 
$\textrm {Re}[\boldsymbol {\cdot }]$ indicates taking the real part, noting that the amplitudes 
$A_i$ and exponents 
$\varPhi _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 
$\varPhi$ 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,-{\rm \pi} /2) = \hat {\psi }_\theta (r,-{\rm \pi} /2)=0$, starting from more general formulas for the 
$\theta$-dependent functions (Moffatt Reference Moffatt1964), we find that
\begin{gather} f_1(\theta) = \frac{\rm \pi}{2}\sin{\theta} - \frac{2}{\rm \pi}\theta \cos{\theta} +\theta \sin{\theta} \quad (\varPhi=1), \end{gather}
\begin{gather}f_2(\theta) = \cos{(2\theta)}-\frac{2}{\rm \pi}\sin{(2\theta)} - \frac{4 \theta}{\rm \pi} -1 \quad (\varPhi=2), \end{gather}
and generally for 
$\varPhi > 0$, 
$\varPhi \neq 1,2$,
\begin{align} f_{\varPhi}(\theta) &= \frac{\cos{(\varPhi{\rm \pi}/2)}\sin{(\varPhi {\rm \pi}/2)}}{\sin^2{(\varPhi{\rm \pi}/2)}-\varPhi }(\cos{(\varPhi\theta)}-\cos{((\varPhi-2)\theta)})\nonumber\\ &\quad + \frac{\cos^2{(\varPhi{\rm \pi}/2)}-\varPhi+1}{\sin^2{(\varPhi{\rm \pi}/2)}-\varPhi }\sin{(\varPhi\theta)} +\sin{((\varPhi-2)\theta)}, \end{align}
which in the special case of integer 
$\varPhi$ reduces to
\begin{align} f_{\varPhi}(\theta) = \begin{cases} \sin{(\varPhi\theta)}+\sin{((\varPhi-2)\theta)}, &\varPhi\ {\rm is\ an\ odd\ integer\ strictly\ greater\ than\ 1},\\ \dfrac{\varPhi-2}{\varPhi} \sin(\varPhi\theta)+\sin((\varPhi-2)\theta), & \varPhi\ {\rm is\ an\ even\ integer\ strictly\ greater\ than\ 2}.\end{cases} \end{align}
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,
\begin{equation} {-}\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. \end{equation}
Writing the first few terms in the expansion (2.14) for 
$\hat {\psi }$ in the form
\begin{equation} \hat{\psi} = \textrm{Re}\left[ A_a r^{a}f_a(\theta)+A_b r^bf_b(\theta)+ A_c r^cf_c(\theta)+\dots \right] \end{equation}
where 
$a$, 
$b$, 
$c\ldots$ are complex numbers (representing exponents 
$\varPhi _i$) such that 
${\textrm {Re}}(a)<{\textrm {Re}}(b)<{\textrm {Re}}(c)$ and 
$A_a\neq 0$, we obtain from (2.16) (after multiplying by 
$r^3$)
\begin{align} &{-}\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 )\nonumber\\ &\quad = \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)\nonumber\\ &\qquad + \left(c^2-3c +2\right)A_c r^{c}f_c^{\prime}(0) + \dots. \end{align}
The 
$r^a$ term is dominant as 
$r \to 0$, imposing that
\begin{equation} (a-1)(a-2)f_a^{\prime}(0)=0. \end{equation}
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 (Reference Moffatt1964), 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{align} &{-}\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 ) \nonumber\\ &\quad = \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. \end{align}
The dominant balance must be between the 
$r^b$ and 
$r^3$ terms since we imposed 
$\textrm {Re}(a)<\textrm {Re}(b)$, hence 
$b=3$. We then have
\begin{equation} {-}\alpha A_2 f_2^{\prime \prime}(0) = 2A_3 f_3^{\prime}(0). \end{equation}
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
\begin{align} &A_2r^2 \left(f_2(\theta) +\frac{\alpha}{2} r f_3{(\theta)} \right)\nonumber\\ &\quad = A_2r^2\left(\cos{(2\theta)}-\frac{2}{\rm \pi}\sin{(2\theta)} - \frac{4 \theta}{\rm \pi} - 1 +\frac{\alpha}{2} r(\sin{(3\theta)}+\sin{\theta})\right). \end{align} In the third case, setting 
$f_a^{\prime }(0)=0$ in (2.19) (with 
$a\neq 1$, 
$2$) gives
\begin{equation} \frac{\cos^2{(a{\rm \pi}/2)}-a+1}{\sin^2{(a{\rm \pi}/2)}-a }a +(a-2) = 0, \end{equation}so that the complex roots satisfy
\begin{equation} \sin^2{(a{\rm \pi}/2)} -a(2-a)=0. \end{equation}
We label the roots of (2.24) which have non-zero imaginary part as 
$a_1, a_2,\dots$ with 
$0<\textrm {Re} (a_1)<\textrm {Re} (a_2)<\dots$. The roots 
$a_i$ are independent of 
$H$, and correspond to the exponents in Moffatt's series (Moffatt Reference Moffatt1964) 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<\textrm {Re} (a_i)<\textrm {Re} (b_i)<\dots$, for 
$i=1,2,3,\dots$ For example, for 
$i=1$ (2.18) gives
\begin{align} &{-}\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 ) \nonumber\\ &\quad =\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, \end{align}
such that the dominant balance is 
$b_1=a_1+1$, requiring that
\begin{equation} {-}\alpha A_{a_1} f^{\prime \prime}_{a_1 }(0) = a_1(a_1-1)A_{b_1} f_{b_1 }^{\prime}(0), \end{equation}
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
\begin{equation} 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), \end{equation}
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
\begin{equation} 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)}, \end{equation}
for integers 
$i\geqslant 1$ and 
$j\geqslant 1$ and with 
$K_{i0}=1$.
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 Reference Moffatt1964) 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$.
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
\begin{align} \hat{\psi}_n(r,\theta)= A_{n2}r^2 \left[f_2(\theta) +\frac{\alpha_n}{2} r f_3{(\theta)} \right] + \textrm{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]. \end{align}
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 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
\begin{equation} \hat{\tau}_n(x) = \left.\frac{\partial\hat{u}_{n}}{\partial y}\right|_{y=0} \approx{-}4A_{n2} + \textrm{Re} \left\{ A_{na_1} (1+x)^{a_1-2} f_{a1}''(0) +\dots \right\} \quad\mathrm{as} \ x\rightarrow -1, \end{equation}
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),
\begin{equation} \hat{\varGamma}_n(x) \approx{-}\frac{8A_{n2}}{{\rm \pi}\alpha_n}+4A_{n2}(1+x)+\dots \quad\mathrm{as} \ x\rightarrow -1. \end{equation}
The leading-order vorticity near the contact line 
$(x,y)=(-1,0)$ has the form
\begin{equation} \hat{\omega}_n\approx 4A_{n2}\left(1+\frac{4}{\rm \pi}\tan^{{-}1}\frac{y}{1+x}\right)+\dots, \end{equation}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
\begin{equation} \hat{p}_n\approx{-}\frac{8A_{n2}}{\rm \pi}\log(y^2+(1+x)^2)+\dots, \end{equation}
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\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)
\begin{equation} \alpha_n \to \frac{n^2{\rm \pi}^2 H}{4},\ \text{(odd modes)}, \quad \alpha_n \to \frac{(2n-1)^2{\rm \pi}^2H}{16} \ \text{(even modes)}, \end{equation}
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) 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\times 4000$ grid points.
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.
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)$.
Figure 4(a) shows the surface shear stress, computed numerically from the viscous-Marangoni stress condition 
$\tilde {\tau }(x)=-\tilde {\varGamma }_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 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\rightarrow -1$ (as plotted with dashed lines when 
$1+x\leqslant 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.
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.175\textrm {e}^{0.157\textit{i}}$ 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\leqslant 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 
$\varGamma ^*_{t^*} = -(\varGamma ^* u^*)_{x^*} + D^*\varGamma ^*_{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
\begin{equation} \alpha\hat{\psi}_{yy} ={-}\hat{\psi}_{xxy} +D\hat{\psi}_{xxyy}, \end{equation}
where 
$D = D^*\mu ^*/(W^*S^*\bar {\varGamma }$). We then specify an additional boundary condition and impose no flux of surfactant, 
$\varGamma _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 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 
$\int _{-1}^1 \tilde {h}\,\mathrm {d} x=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.
(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\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.
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 
${\rm \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 
${\rm \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 
${\rm \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 
${\rm \pi} /2$ (see (2.22) for the streamfunction with 
${\rm \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 
${\rm \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}\ \textrm {m}^2\ \textrm {s}^{-1}$ for the common surfactant sodium dodecylsulphate (Chang & Franses Reference Chang and Franses1995) translates to a dimensionless diffusion coefficient 
$D=\mu ^* D^*/(S^* W^*)$ below 
$10^{-8}$ in magnitude, assuming a spreading coefficient 
$S^*=10^{-2}\ \textrm {kg}\ \textrm {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 
${\rm \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 Reference Sprittles and Shikhmurzaev2011; Game, Hodes & Papageorgiou Reference Game, Hodes and Papageorgiou2019), 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 Reference Meleshko1996, Reference Meleshko1997).
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 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. Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017), in maze solving by surfactant (Temprano-Coleto et al. Reference Temprano-Coleto, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2018) 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 {\varGamma }^*/\varGamma _c^*$ is the surface tension reduction due to the equilibrium surfactant distribution. Taking this as low as 
$10^{-3}\ \textrm {kg}\ \textrm {s}^{-2}$, say, for water in a narrow channel of width 1 mm (and comparable depth), the decay time scale 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.
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.
Appendix A. The energy budget
The Stokes equation can be written as 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\sigma } = \boldsymbol {0}$, where the non-dimensional Cauchy stress tensor is 
$\boldsymbol {\sigma } = -p\boldsymbol {I} + 2 \boldsymbol{\mathsf{e}}$ and the strain-rate tensor 
$\boldsymbol{\mathsf{e}} = (\boldsymbol {\nabla }\boldsymbol {u} + \boldsymbol {\nabla }\boldsymbol {u}^\textrm {T})/2$. Thus, for a Stokes flow, 
$\boldsymbol {\nabla }\boldsymbol {\cdot } (\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\sigma }) = \boldsymbol {\sigma }: \boldsymbol {\nabla } \boldsymbol {u} = \boldsymbol {\sigma }:\boldsymbol{\mathsf{e}}$, exploiting the symmetry of 
$\boldsymbol {\sigma }$. Given that 
$p\boldsymbol {I}:\boldsymbol {e} = \text {trace}(\boldsymbol{\mathsf{e}})p = 0$ by incompressibility, it follows that 
$\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\sigma }) = 2\boldsymbol{\mathsf{e}}: \boldsymbol{\mathsf{e}}$. Integrating this over the domain and applying the divergence theorem gives
\begin{equation} 2\int_V \boldsymbol{\mathsf{e}}\ \mathsf{:}\ \boldsymbol{\mathsf{e}} \,\mathrm{d}A = \int_{\partial V} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\sigma} \boldsymbol{\cdot} {\boldsymbol{n}} \,\mathrm{d}s, \end{equation}
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 (2.5), (2.6b,c) and (2.10) for each eigenmode, (A1) balancing work done by Marangoni forces with viscous dissipation becomes
\begin{equation} 2 \int_V \left(2\hat{\psi}_{yx}^2 + \tfrac{1}{2}\hat{\psi}_{yy}^2+\tfrac{1}{2}\hat{\psi}_{xx}^2 -\hat{\psi}_{xx}\hat{\psi}_{yy} \right) \mathrm{d}A ={-}\int_{{-}1}^{1} \left. \left(\hat{\varGamma}_x \hat{\psi}_y \right) \right|_{y=0} \,\mathrm{d}\kern0.06em x . \end{equation}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
\begin{align} 2 \int_V \left(2\hat{\psi}_{yx}^2 + \tfrac{1}{2}\hat{\psi}_{yy}^2+\tfrac{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, \end{align}
where 
$\hat {\omega } =-(\hat {\psi }_{yy} + \hat {\psi }_{xx})$ is the vorticity. Further integration by parts of the right-hand side of (A2) gives
\begin{equation} - \int_{{-}1}^{1} \left. \left(\hat{\varGamma}_x \hat{\psi}_y \right) \right|_{y=0} \, \mathrm{d}\kern0.06em x = \int_{{-}1}^{1} \left. \left(\hat{\varGamma} \hat{\psi}_{xy} \right) \right|_{y=0} \,\mathrm{d}\kern0.06em x . \end{equation}
Since 
$\hat {\psi }_{xy} = \alpha \hat {\varGamma }$ from (2.8), we obtain (2.11) for each decay rate 
$\alpha$ and corresponding eigenmode 
$\{\hat {\psi },\hat {\omega },\hat {\varGamma }\}$.
Appendix B. Orthogonality of modes
The reciprocal theorem (Masoud & Stone Reference Masoud and Stone2019) states that two Stokes flows, with velocity fields and Cauchy stress tensors 
$(\boldsymbol {u},\boldsymbol {\sigma })$ and 
$(\boldsymbol {u}',\boldsymbol {\sigma }')$, satisfy in the same region
\begin{equation} \int_{\partial V} \boldsymbol{u} \boldsymbol{\cdot}(\boldsymbol{\sigma}'\boldsymbol{\cdot} \boldsymbol{n})\, \mathrm{d}A = \int_{\partial V} \boldsymbol{u}' \boldsymbol{\cdot}(\boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n})\, \mathrm{d}A . \end{equation}
For the present perturbation problem, described by (2.5), (2.6b,c) and (2.10) for each eigenmode, we have 
$\boldsymbol {u}=\boldsymbol {0}$ on the three solid boundaries whilst at the free surface 
$u_y = - \varGamma _x$. Thus, for two distinct modes 
$m$ and 
$n$, (B1) implies
\begin{equation} \int_{{-}1}^1 \hat{u}_m \hat{\varGamma}_{n,x}\, \textrm{d} x = \int_{{-}1}^1 \hat{u}_n \hat{\varGamma}_{m,x}\, \textrm{d} x. \end{equation}
Integrating by parts and using 
$\alpha \hat {\varGamma } = -\hat {u}_x$ gives
\begin{equation} \left[\hat{u}_m\hat{\varGamma}_n\right]_{{-}1}^1 - \int_{{-}1}^1 \alpha_m\hat{\varGamma}_m\hat{\varGamma}_n \,\textrm{d} x = \left[\hat{u}_n\hat{\varGamma}_m\right]_{{-}1}^1 - \int_{{-}1}^1 \alpha_n\hat{\varGamma}_n\hat{\varGamma}_m \,\textrm{d} x. \end{equation}
As the surface velocity is zero at 
$x=\pm 1$, (B3) becomes
\begin{equation} (\alpha_m-\alpha_n)\int_{{-}1}^1 \hat{\varGamma}_m \hat{\varGamma}_n\, \textrm{d} x = 0, \end{equation}
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\leqslant m,n\leqslant 5$ gives
\begin{align} &\int_{{-}1}^1 \tilde{\varGamma}_m \tilde{\varGamma}_n \, \mathrm{d}\kern0.06em x \nonumber\\ &\quad = \boldsymbol{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} . \end{align}These results were computed using our numerical scheme presented in § 2.2 with 
$4000\times 4000$ grid points. The non-diagonal elements of the matrix 
$\boldsymbol {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\rightarrow 0$, the biharmonic equation (2.5) can be approximated as 
$\psi _{yyyy} = 0$, therefore the general solution in this limit is given by
\begin{equation} \psi = f_1(x)y^3 + f_2(x)y^2+ f_3(x)y + f_4(x). \end{equation}
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
\begin{equation} 0 ={-}f_1(x)H^3 + f_2(x)H^2 -f_3(x)H, \end{equation}and
\begin{equation} 0 = 3f_1(x)H^2 - 2 f_2(x)H +f_3(x). \end{equation}
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
\begin{equation} 4\alpha H f_1(x) ={-} f_1^{\prime \prime}(x)H^2. \end{equation}
The solution of this ordinary differential equation for 
$f_1(x)$ is
\begin{equation} 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)} , \end{equation}
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
\begin{equation} \alpha_n = \frac{(2n-1)^2 {\rm \pi}^2 H}{16}, \end{equation}
where 
$n$ can be any positive integer, 
$n\geqslant 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
\begin{equation} \alpha_n = \frac{n^2{\rm \pi}^2 H}{4}, \end{equation}
where 
$n$ can be any positive integer, 
$n\geqslant 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).
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 
${\rm \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 
$\varTheta$ to vary between 0 and 
${\rm \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 
$\varTheta _1$ and 
$\varTheta _2 = {\rm \pi}- \varTheta _1$, respectively.
(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 
$\varTheta _1$ and 
$\varTheta _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 
$\varTheta _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 
${\rm \pi} /2$ is generic.
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\geqslant 0,\theta =0)$, we also assume continuity of the tangential velocity, whilst the tangential dynamic boundary condition now becomes, in dimensional form,
\begin{equation} {-}{\boldsymbol{\tau}}\boldsymbol{\cdot} [\kern-1pt[ \boldsymbol{\sigma}^*]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n} = {\boldsymbol{\tau}} \boldsymbol{\cdot} {\boldsymbol{\nabla}}^*_s \gamma^*, \end{equation}
where the jump in tangential stress across the interface is balanced by the surfactant-induced Marangoni stress. The jump bracket is defined as 
$[\kern-1pt[ \boldsymbol {\sigma }^*]\kern-1pt] = \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 
${\boldsymbol {n}}=(0,1)$, in polar coordinates, as depicted in figure 6(a). The surface gradient operator is defined as 
${ \boldsymbol {\nabla }}^*_s={\boldsymbol {\nabla }}^* \boldsymbol {\cdot }(\boldsymbol{\mathsf{I}}-{\boldsymbol {n}} \otimes {\boldsymbol {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 
$\bar {\varGamma }$, and decompose the perturbation for each variable in the form 
$\hat {\varGamma }(r)\textrm {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{align} &{-}\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)\nonumber\\ &\quad =\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{align}
where 
$\alpha =\lambda /\bar {\varGamma }$ 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
\begin{equation} \left(a_1^2-3a_1 +2\right)f_{a_1}^{\prime}(0;\varTheta_1) = 0 \quad \text{or} \quad \left(a_2^2-3a_2 +2\right)f_{a_2}^{\prime}(0;-\varTheta_2) = 0. \end{equation}
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}^{\prime }(0;\varTheta _1) = 0$ or 
$f_{a_2}^{\prime }(0;-\varTheta _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 
${\rm \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}^{\prime }(0,\varTheta _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 
$\varTheta _1$ subject to zero velocity boundary conditions at the boundaries (and similarly for 
$a_2$) (Moffatt Reference Moffatt1964). This condition is
\begin{equation} (a_1-1)\sin{(\varTheta_1)} ={\pm} \sin{(\varTheta_1(a_1-1))} . \end{equation}
By inspection we can see that 
$a_1=0,1,2$ (and similarly for 
$a_2$) are solutions of (D4), for any value of 
$\varTheta _1$, and any integer is a solution when 
$\varTheta _1 = {\rm \pi}$. We show that a solution with 
$1 < \text {real}(a_1)<2$ cannot exist for 
$0\leqslant \varTheta _1 \leqslant {\rm \pi}$ in figure 6(b) where we have plotted the locations of the real parts of the admissible solutions of (D4) against 
$\varTheta _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.
Open access
