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

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.


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 π/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 π, 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 π/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 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 * x * W * and −H * y * 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 * x * 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).
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.

Model
We model the spreading of an insoluble surfactant at the surface of an incompressible liquid of dynamic viscosity µ * in a two-dimensional rectangular domain. The domain V spans from −W * x * W * and −H * y * 0 with x * and y * in the horizontal and vertical directions, respectively, with W * the half width of the cavity, H * its height, as shown in figure 1. Stars in superscript indicate dimensional variables. To model the flow induced by the surfactant spreading at the free surface, we use the Stokes equations to relate the velocity field u * (x * , t * ) = (u * , v * ) to the pressure p * (x * , t * ), with x * = (x * , y * ) ignoring gravity and inertia. We impose the no slip and no penetration conditions on the three solid boundaries found at x * = −W * and W * for −H * y * 0 for the side walls, and y * = −H * for −W * x * W * for the bottom wall. The free surface F , assumed flat to leading order, is located at y * = 0 for −W * x * W * . We balance the Cauchy stress σ * · n at F (with the normal unit vector n = (0, 1)) with the gradient of the surface tension γ * tangentially and a restoring force due to strong surface tension normally. The insoluble surfactant concentration Γ * 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 where subscript s denotes evaluation on F , γ * 0 is the reference surface tension and γ * c is its (lower) value when the surfactant is at a reference concentration Γ * c . The boundary conditions associated with (2.1) are with t = (1, 0) the tangential unit vector at the free surface. We have made the assumption that S * = γ * 0 − γ * c ≪ γ * 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 * /γ c ≪ 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 * ∇ * · n ≪ 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 F .
Using the length scale W * , velocity scale S * /µ * and pressure scale S * /W * , we relate dimensional starred variables to their dimensionless counterparts via After rescaling, the governing equations become, in the bulk, with, at the free surface, and u = 0, on x = ±1, and y = −H.
(2.4c) We introduce a stream function ψ(x, y, t) such that ψ y = u, and ψ x = −v, enforcing incompressibility. The problem reduces to the biharmonic equation in the bulk (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 Γ , given a surface velocity ψ y . Therefore writing Γ as the sum of endogenous and exogenous components, Γ = Γ 1 + Γ 2 say, the two components satisfy the same transport equation Γ it = −(Γ i ψ 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) = (±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.

Linearization
At large times, the surfactant relaxes to a uniform levelΓ 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 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ψ and αψ yy = −ψ xxyψ = 0, on y = 0. (2.10) We seek the set of decay rates α n , n = 1, 2, . . . , which are eigenvalues for which a solution exists with the corresponding eigenmodesψ n . We distinguish two families for the decay rates α n and eigenmodesψ n , such thatψ 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Γ n , the shear stress at the free surface τ n (x) = ∂û n (x, 0)/∂y, the vorticity fieldω n = −∇ 2ψ n , and the pressure fieldp n , obtained from the vorticity via ∇p n = (−∂ω n /∂y, ∂ω n /∂x).
This problem has two key global characteristics. An energy dissipation argument (see Appendix A) shows that all the decay rates satisfy 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 (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Γ 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).

Finite-difference numerical solution
We compute a numerical solution ψ n (where a wide tilde denotes a numerically computed variable) of the unknown perturbation modes of the stream function,ψ n , satisfying the biharmonic equation (2.5) and the boundary conditions (2.6b,c) and (2.10), using a finite-difference scheme. Using a row×column ordering convention, the domain described in figure 1 is discretized as an M ×N rectangular grid, in the y and x directions, respectively, with uniform spacing ∆y and ∆x, and supplemented with a set of ghostpoints around the periphery creating an (M + 2) × (N + 2) grid. We use a secondorder-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 ψ n (x j , y i ), with 3 i M , 3 j N . We impose ψ n (x j , y i ) = 0 at the boundaries of the domain j = 2 and j = N + 1, with 2 i M + 1, and i = 2 and i = M + 1, with 2 j 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) × (N − 2) grid of unknowns is expressed as a column vector ψ which is assembled by the vertical concatenation of the rows of unknowns of ψ n (x j , y i ). The numerical operator modelling the biharmonic operator on the grid can then be approximated by a matrix operating on ψ. 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 ×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 ψ satisfies,  Table 1: Decay rates predicted as eigenvalues α n computed numerically from (2.13) compared toα n , which are those computed from eigenmodes using (2.11) for H = 2, found using a 4000 × 4000 grid. The relative error provides a measure of global numerical error, suggesting that the values for α n are accurate up to three significant figures.
where B and C are sparse ) represents a generalised numerical eigenvalue problem for the eigenmodes ψ n and the associated numerically-calculated decay rates α n , from the free-surface boundary condition (2.10). We solved (2.13) using the MATLAB function eigs. From the solutions for each mode ψ n we compute numerical approximations of all other quantities of interest for each mode, such as the surfactant concentration profile Γ n (x), the surface stress τ n (x) = ψ n,yy (x, y = 0), the velocity field ( u n , v n )(x, y), the vorticity field ω n (x, y), and the pressure field p n (x, y). The solutions ψ n are normalised by requiring max( Γ n (x)) − min( Γ 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α n ) remains less than 4 × 10 −4 for all odd and even modes calculated (up to n = 4) for H = 2 using a 4000 × 4000 grid (see table 1). For the same refinement and the same modes, the concentration profiles Γ n are orthogonal with an absolute error less than 8 × 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) = (±1, 0).

Corner asymptotics
We anticipate from the outset the appearance of Moffatt vortices (Moffatt 1964) in the lower corners of the domain, at (x, y) = (±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, θ) centred on (x, y) = (−1, 0) with θ = 0 along y = 0 (and r increasing in the positive x direction) and θ = −π/2 along x = −1 (and r increasing in the negative y direction) and seeking separable solutions of the form (2.14) Here ℜ [·] indicates taking the real part, noting that the amplitudes A i and exponents Φ 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 Φ to represent exponents, and we will suppress the subscript n onψ and α associated with each eigenmode. To satisfy the governing biharmonic equation (2.5), ∇ 4ψ (r, θ) = 0, and the boundary conditions (2.6b,c),ψ(r, 0) =ψ(r, −π/2) =ψ θ (r, −π/2) = 0, starting from more general formulas for the θ-dependent functions (Moffatt 1964), we find that if Φ is an odd integer strictly greater than 1, if Φ is an even integer strictly greater than 2.
(2.15d ) The final boundary condition, the surfactant boundary condition at the free surface in (2.10): −αψ yy =ψ xxy at y = 0, is, in polar coordinates, Writing the first few terms in the expansion (2.14) forψ in the form The r a term is dominant as r → 0, imposing that This equation opens multiple possibilities. The first case a = 1 corresponds to a solution with non-integrable stress τ (r) =ψ θθ (θ = 0)/r 2 , and therefore unbounded surfactant concentration as r → 0, so we impose A 1 = 0. This leaves two other cases: a = 2 and f ′ a (0) = 0. The latter third case yields an infinite set of complex exponents of the type described by Moffatt (1964), each representing a homogeneous local solution, which we will analyse further below. The expansion (2.14) therefore constitutes multiple independent series.
In the second case, with a = 2, (2.18) becomes The dominant balance must be between the r b and r 3 terms since we imposed ℜ(a) < ℜ(b), hence b = 3. We then have Using (2.15b) and (2.15d) for f 2 and f 3 , respectively, (2.21) becomes αA 2 = 2A 3 . The next balance in (2.20) gives −αA 3 f ′′ 3 (0) = 6A 4 f ′ 4 (0), however from (2.15d) we can see that f ′′ 3 (0) = 0, which implies that A 4 = 0 and this series terminates. The first series contributing toψ is therefore simply In the third case, setting f ′ a (0) = 0 in (2.19) (with a = 1 and 2) gives cos 2 (aπ/2) − a + 1 sin 2 (aπ/2) − a a + (a − 2) = 0, (2.23) so that the complex roots satisfy sin 2 (aπ/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 , . . . with 0 < ℜ(a 1 ) < ℜ(a 2 ) < . . . . 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. Each complex root a generates its own asymptotic series of the form with the approximate value of a 1 given in table 2. From such relations for a 1 , b 1 , c 1 , . . . we can derive the associated contribution toψ, of the form 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 for integers i 1 and j 1 and with K i0 = 1. In summary, the full asymptotic series for the nth eigenmode in the neighbourhood of the contact line, which we originally stated in the form (2.14), iŝ (2.29) The coefficients K ij can be computed from the recurrence relation (2.28), whilst the coefficients A n2 and A nai must be determined from fitting to numerical solutions. Analysing (2.29) we can notice thatψ n approaches different values as r → 0 for different values of θ, 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ψ n enters into this expression through α n . We note for later reference that the surface shear stress τ n (x) has the local form which has a non-zero asymptotic value at the contact linesτ n (−1) = −4A n2 . Moreover, the complex roots a i imply that the value of the shear stress oscillates as x → −1. Near the contact line (x, y) = (−1, 0), the leading-order surfactant distribution has a finite non-zero gradient, computed from (2.8), (2.31) The leading-order vorticity near the contact line (x, y) = (−1, 0) has the form ω n ≈ 4A n2 1 + 4 π tan −1 y 1 + x + . . . , (2.32)  (a)) for the first two even and two odd modes for H = 2 using a solution for the streamfunction ψ n calculated numerically using 4000 × 4000 gridpoints.
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 isp n ≈ − 8A n2 π log(y 2 + (1 + x) 2 ) + . . . , (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) = (±1, 0), providing a complete understanding of the effect of the corner singularities on the relevant physical quantities in this problem.

Results
Figure 2(a) shows how the decay rates α 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 2. We can find an exact asymptotic expression for the odd and even decay rates in the limit H → 0 using a lubrication approximation (see Appendix C) α n → n 2 π 2 H 4 , (odd modes), α n → (2n − 1) 2 π 2 H 16 (even modes), The contour plots computed numerically in figure 3(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) = (±1, −2). Vorticity contours converge at the upper corners, indicating thatω 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 3(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 3d ), elongated eddies appear, consistent with predictions of lubrication theory. Figure 4(a) shows the surface shear stress, computed numerically from the viscous-Marangoni stress condition τ (x) = − Γ x (x), for the dominant mode (first odd mode, n = 1) for H = 2, revealing an oscillatory structure near the contact lines as x → ±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 → −1 (as plotted with dashed lines when 1 + x 100∆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 (∆x = 1/8000) the numerical results are close to the asymptotic results for log(1 + x) ≈ −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 → −1. However, the dominant odd mode n = 1 clearly captures the oscillatory behaviour of the shear stress near the corner.

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 Γ * t * = −(Γ * u * ) x * + D * Γ * 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 where D = D * µ * /(W * S * Γ ). We then specify an additional boundary condition and impose no-flux of surfactant, Γ x = 0, at the contact lines (x, y) = (±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 ≪ 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 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.
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 p(x, 0) − p ext − 2 ψ xy (x, 0) = h xx , where the dimensional deflection is (W * S * /γ * 0 ) h and p ext is a reference pressure (recall that, for the small surfactant concentration variations considered here, we can assume S * ≪ γ * 0 ). The displacement is computed by imposing h = 0 at each contact line, and imposing a volume constraint

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 longestwave 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 π/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 π 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 singlefluid flow with π/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 π/2 (see (2.22) for the streamfunction with π/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 π/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× 10 −10 m 2 /s for the common surfactant sodium dodecylsulfate (SDS) (Chang & Franses 1995) translates to a dimensionless diffusion coefficient D = µ * D * /(S * W * ) below 10 −8 in magnitude, assuming a spreading coefficient S * = 10 −2 kg/s 2 over water in a channel of width 1 cm. At such low levels, the impact of the singularities may still be visible, with the pressure maximum near each contact line, proportional to log(1/D), generating surface deflections resembling that shown in figure 5. We expect that diffusion will act in the same way for the broader class of problems involving two-fluid viscous flows and arbitrary contact angles between 0 and π.
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(Meleshko , 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/α times W * µ * /∆γ * , where α is an O(1) modal decay rate and ∆γ * = (γ * 0 − γ * c )Γ * /Γ * c is the surface tension reduction due to the equilibrium surfactant distribution. Taking this as low as 10 −3 kg/s 2 , say, for water in a narrow channel of width 1 mm (and comparable depth), the decay timescale is approximately 1 ms/α, with a Reynolds number of order unity. Despite decaying exponentially with respect to time, the structure of the flow in wider and deeper channels can therefore be expected to be influenced by inertia, at least initially.
In summary, this study shows how the unsteady spreading of a surfactant monolayer along a liquid-liquid or liquid-gas interface, confined by a lateral rigid boundary, can generate singular flow structures near stationary contact lines, including a logarithmically divergent pressure field and an oscillatory shear stress distribution. Careful treatment of these structures is needed in computational simulations involving dynamic surfactant transport in confined domains. Sinceψ xy = αΓ from (2.8), we obtain (2.11) for each decay rate α and corresponding eigenmode {ψ,ω,Γ }.

Appendix B. Orthogonality of modes
The reciprocal theorem (Masoud & Stone 2019) states that two Stokes flows, with velocity fields and Cauchy stress tensors (u, σ) and (u ′ , σ ′ ), satisfy in the same region For the present perturbation problem, described by equations (2.5), (2.6b,c) and (2.10) for each eigenmode, we have u = 0 on the three solid boundaries whilst at the free surface u y = −Γ x . Thus, for two distinct modes m and n, (B 1) implies Integrating by parts and using αΓ = −û x gives û mΓn As the surface velocity is zero at which shows that the surfactant concentration profiles corresponding to different modes are orthogonal since α m = α 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 m, n 5 gives  Figure 6: (a) Schematic of the local problem near the contact line pinned at O, with surfactant on the interface (θ = 0) between two incompressible fluids with viscosities µ 1 (bottom fluid) and µ 2 (top fluid) and contact angles Θ 1 and Θ 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 Θ 1 . This gives the exponents in the asymptotic series capturing the behaviour of the fluid as r → 0. Importantly, this shows no admissible exponents with 1 < 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 π/2 is generic.

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 π/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 Θ to vary between 0 and π, thereby relaxing the assumption made in the main part of the text. We find local approximations to the stream functions in each fluid, ψ (1) in fluid 1 (bottom fluid), and ψ (2) in fluid 2 (top fluid). These fluids have viscosities µ 1 and µ 2 , and contact angles Θ 1 and Θ 2 = π − Θ 1 , respectively.
We use a plane polar coordinate system with r being the radial direction from the origin, and θ the angular coordinate, with the interface along the line θ = 0 (see figure  6a). 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 nopenetration at the cavity wall, and no penetration at the interface, which is assumed fixed and locally flat. At the interface (r 0, θ = 0), we also assume continuity of the tangential velocity, whilst the tangential dynamic boundary condition now becomes, in dimensional form, − τ · σ * · n = τ · ∇ * s γ * (D 1) where the jump in tangential stress across the interface is balanced by the surfactantinduced Marangoni stress. The jump bracket is defined as σ * = σ * 2 − σ * 1 . The stress tensor σ * is assumed Newtonian for both fluids. The unit tangential and normal vectors at the interface follow τ = (1, 0) and n = (0, 1), in polar coordinates, as depicted in figure 6(a). The surface gradient operator is defined as ∇ * s = ∇ * · (I − n ⊗ n), with ⊗ the outer product. Similar to before, we non-dimensionalize this problem using (2.3), taking fluid 1 as the reference fluid, then we linearize the surfactant distribution aroundΓ , and decompose the perturbation for each variable in the formΓ (r)e −λ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, − α 1 r 2 ψ (1)θθ + 1 r ψ (1)r − µ r r 2 ψ (2)θθ − µ r r ψ (2)r = 1 r ψ (1)rrθ − 2 r 2 ψ (1)rθ + 2 r 3 ψ (1)θ = 1 r ψ (2)rrθ − 2 r 2 ψ (2)rθ + 2 r 3 ψ (2)θ , on θ = 0.
( D 2) where α = λ/Γ is the decay rate of each mode, and µ r = µ 2 /µ 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 ψ (1) = r a1 f a1 (θ) and ψ (2) = r a2 f a2 (θ) 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 a1 and f a2 depend on the contact angle, whilst the exponents a 1 and a 2 satisfy a 2 1 − 3a 1 + 2 f ′ a1 (0; Θ 1 ) = 0 or a 2 2 − 3a 2 + 2 f ′ a2 (0; −Θ 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 ′ a1 (0; Θ 1 ) = 0 or f ′ a2 (0; −Θ 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 π/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 → 0.
When we impose the condition f ′ a1 (0, Θ 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 Θ 1 subject to zero velocity boundary conditions at the boundaries (and similarly for a 2 ) (Moffatt 1964). This condition is (a 1 − 1) sin (Θ 1 ) = ± sin (Θ 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 Θ 1 , and any integer is a solution when Θ 1 = π. We show that a solution with 1 < real(a 1 ) < 2 cannot exist for 0 Θ 1 π in figure 6(b) where we have plotted the locations of the real parts of the admissible solutions of (D 4) against Θ 1 . We also note that none of the exponents in the expansions for ψ (1) or ψ (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.