Confinement-induced drift in Marangoni-driven transport of surfactant: a Lagrangian perspective

Abstract

Successive drops of coloured ink mixed with surfactant are deposited onto a thin film of water to create marbling patterns in the Japanese art technique of Suminagashi. To understand the physics behind this and other applications where surfactant transports adsorbed passive matter at gas–liquid interfaces, we investigate the Lagrangian trajectories of material particles on the surface of a thin film of a confined viscous liquid under Marangoni-driven spreading by an insoluble surfactant. We study a model problem in which several deposits of exogenous surfactant simultaneously spread on a bounded rectangular surface containing a pre-existing endogenous surfactant. We derive Eulerian and Lagrangian formulations of the equations governing the Marangoni-driven surface flow. Both descriptions show how confinement can induce drift and flow reversal during spreading. The Lagrangian formulation captures trajectories without the need to calculate surfactant concentrations; however, concentrations can still be inferred from the Jacobian of the map from initial to current particle position. We explore a link between thin-film surfactant dynamics and optimal transport theory to find the approximate equilibrium locations of material particles for any given initial condition by solving a Monge–Ampère equation. We find that as the endogenous surfactant concentration $\delta$ vanishes, the equilibrium shapes of deposits using the Monge–Ampère approximation approach polygons with corners curving in a self-similar manner over lengths scaling as $\delta ^{1/2}$. We explore how Suminagashi patterns may be produced by using computationally efficient successive solutions of the Monge–Ampère equation.

1. Introduction

Successive drops of coloured inks mixed with surfactant create intricate patterns by Marangoni spreading in the Japanese art technique of Suminagashi (see figure 1a). A surfactant–ink drop is gently deposited at the surface of a thin layer of water, which may have a small initial concentration of endogenous surfactant due to normal environmental contamination. It then spreads outwards and equilibrates before reaching the edges of the container. Then successive drops deposited at different locations of the liquid surface form the intricate patterns. During pattern creation, the artist can blow on the surface with a straw after drop equilibrations to further deform the pattern. Eventually, the pattern is captured on pieces of paper placed onto the surface (Chambers 1991). Rouwet & Iorio (2017) noticed similar patterns occurring in volcanic crater lakes, and hypothesised that similar physics were responsible: thermal gradients in the lake create Marangoni flows, and wind action creates a blowing effect, resulting in marbling patterns of the adsorbed multicoloured sediments. In this study, we seek to understand the Lagrangian trajectories of material points and curves on a surface during the spreading of surfactants on a confined surface, and thus the dynamics of any adsorbed passive tracer, similar to the advection of ink by surfactant in the Suminagashi technique.

Figure 1. (a) Pictures showing the Japanese art technique of Suminagashi (Rouwet & Iorio 2017). Successive drops of a mixture of coloured ink and surfactant are deposited on the surface of a thin film of water to create a multicoloured pattern. Blowing on the surface then creates further intricate patterning. (b) Picture of a Suminagashi pattern of ink on water created by artist Bea Mahan (2011). (c) Schematic of the model problem. Circular deposits of insoluble exogenous surfactant (red) spread on the surface $\varOmega$ of a thin layer of viscous liquid (blue) of mean height $h$ confined in a rectangular region of dimensions $L_1$ and $L_2$, where the surface contains an initially uniform endogenous surfactant (green). We assume that the ratio of vertical to horizontal length scales is small enough, and that the Bond number (ratio of gravitational to surface tension forces) is large enough, for height deflections caused by spreading to be negligible, confining spreading to the flat plane of the surface $\varOmega$.

In addition to the cultural importance of Suminagashi, which has been part of Japanese art since the 12th century, and similar practices in China for even longer (Ishii & Muro 1989), understanding Marangoni-driven surface motion can help us to better understand various industrial and biological applications involving surfactants carrying passive, adsorbed material. For example, Deng et al. (2018) showed how small amounts of surfactant added to perovskite (a calcium titanium oxide mineral) can suppress the formation of islands during the drying phase of blade coating by creating Marangoni flows that keeps the solution coating even. Many other coating processes involve Marangoni flows induced by trace amounts of surfactants. Some methods of drug delivery in lungs mix pharmaceutical substances with exogenous surfactant (Haitsma, Lachmann & Lachmann 2001), so that the surfactant acts as a carrier to spread the drug through the airways. In particular, surfactant replacement therapies have been used successfully in lungs of neonates affected with respiratory distress syndrome (Avery & Mead 1959; Jobe 1993; Rodriguez 2003; Halliday 2008). The surfactant-driven spreading in the complex and confined tree-like geometry of the lungs acts against its natural endogenous surfactant (Espinosa et al. 1993; Jensen, Halpern & Grotberg 1994; Grotberg, Halpern & Jensen 1995; Halpern, Jensen & Grotberg 1998; Temprano-Coleto et al. 2018; Mcnair et al. 2023). These methods of delivery can help to overcome difficulties such as poor solubility of the pharmaceuticals (Hidalgo, Cruz & Pérez-Gil 2015).

Molecules and substances that act as surfactants are ubiquitous in the environment. They can cause unexpected fluid flows that have confounded scientists and engineers, as described by Manikantan & Squires (2020), who discussed the ‘hidden’ variables related to surfactant dynamics in many fluid flows. The present study addresses insoluble surfactant spreading into pre-existing, endogenous surfactant on a thin film of a bounded Newtonian viscous liquid, allowing us to use lubrication theory to approximate the Stokes flow in the liquid film. Lubrication theory for insoluble surfactant-driven flows has its origins with the work of Borgas & Grotberg (1988) who derived coupled partial differential equations (PDEs) describing the leading-order evolution of the liquid film height and surfactant concentration. The work was then extended theoretically and experimentally by Gaver & Grotberg (1990, 1992). Thess, Spirn & Jüttner (1997) and Jensen & Halpern (1998) showed that the coupled equations in the limit of large Bond number could be combined into a single nonlinear diffusion (or ‘porous medium’) equation governing surfactant concentration evolution as a function of space and time. The effect of gravity is to suppress deflections of the surface, removing the functional dependence of the spreading on the dynamic film height.

In this paper, we explore a link between the theory of surfactant spreading and the theory of optimal transport. This theory was initiated by Monge (1781), who was trying to find the optimal way to transport mounds of soil under some cost function. The theory was extended into its modern formulation by Kantorovich (1942, 2006). Most of its current uses are found in machine learning and image analysis (Kolouri et al. 2017). A powerful result, enabling significant simplification of optimal transport problems, occurs when the cost function takes a quadratic form, yielding the quadratic Monge–Kantorovich optimal transport problem (qMK). For such cost functions, solutions for the optimal map of material from initial to final location can be shown to be the gradient of a convex function that satisfies a so-called Monge–Ampère equation. A variety of approaches have been taken to find solutions of this nonlinear equation (Froese & Oberman 2011; Benamou, Froese & Oberman 2012, 2014). Otto (2001), building on work by Jordan, Kinderlehrer & Otto (1998) and Benamou & Brenier (2000), showed that porous medium equations (a class of equations to which the Jensen & Halpern (1998) surfactant equation that we use in this study belongs) have the variational structure of a gradient flow on a Riemannian manifold measured by the quadratic Wasserstein distance. The square of the Wasserstein distance, which is defined as the minimiser of a functional, doubles as qMK, which suggests that solutions to the surfactant-induced transport problem may be approximated by solving the Monge–Ampère equation under certain conditions. We explore these conditions in this paper, and consider whether the Monge–Ampère equation could be an efficient tool to determine equilibrium solutions to this complex confined transport problem.

The primary aim of this study is to understand the underlying physics behind surfactant-induced Marangoni dynamics in a confined environment when the surface contains an initial endogenous concentration of surfactant, which is the case for most environmental fluids. The spreading of multiple exogenous deposits is particularly considered; this was investigated experimentally and with COMSOL$^{\circledR}$ models recently by Iasella et al. (2024), showing how adjacent droplets interact and deform. A Lagrangian framework, which has been adopted in the analysis of other transport problems with nonlinear diffusive character (Meĭrmanov, Pukhnachev & Shmarev 1997), enables us to compute efficiently individual surface particle trajectories and equilibrium states as functions of initial distributions. Moreover, the Lagrangian framework reveals underpinning flow phenomena such as stretching, compression and rotational motion that govern the particle trajectories. While there have been limited investigations of Lagrangian surfactant dynamics in one spatial dimension (Grotberg et al. 1995), there is none (to our knowledge) in higher dimensions, despite the potential relevance to a variety of applications. Furthermore, while some authors have exploited the gradient flow structure of thin-film evolution equations (Thiele, Archer & Pismen 2016; Henkel, Snoeijer & Thiele 2021), we are not aware of prior studies linking thin-film flows to optimal transport. We show how we can exploit this link for practical purposes. In particular, we describe a procedure to reproduce the intricate patterns of Suminagashi art, through resolution of the Monge–Ampère equation associated with the surfactant transport model. These results appear to capture, at least qualitatively, the dominant physics behind Suminagashi art, suggesting a powerful tool for other applications where surface transport is dominated by surfactants in confined environments.

In § 2.1, we use a two-dimensional extension of the model of Jensen & Halpern (1998) (derived in Appendix A) to describe transport of material particles on a surface. We outline a physical problem in Eulerian coordinates in a confined rectangular domain, implementing initial conditions that represent multiple deposits of exogenous surfactant spreading on a surface with an initially uniform endogenous surfactant concentration. We solve the particle-tracking problem using a finite-difference method by first solving for the evolution of surfactant concentration, and then interpolating the gradient of this solution onto a second Lagrangian grid where we integrate the surface velocity to find the trajectories of surface particles initially located at each grid point. In § 2.2, we reformulate the problem in Lagrangian coordinates, and show how it can be reduced from three to two scalar PDEs, enabling the same calculation without the intermediate step of finding the evolution of the surfactant concentration, and without the need to interpolate concentration gradients from an Eulerian to a Lagrangian grid. We solve the resulting scheme using a finite-element method. In § 2.3, we show how to approximate the equilibrium locations of surface particles as a function of their initial locations via a Monge–Ampère equation, without having to compute their intermediate trajectories. In § 3.1, we show consistency between the Eulerian and Lagrangian methods, and describe dynamical phenomena not normally associated with spreading surfactants, such as drift and flow reversals due to confinement. In § 3.2, we show that solutions of the Monge–Ampère equation approximate the equilibrium solution well when the endogenous and exogenous concentrations are of comparable magnitude, and also provide a credible approximation when the endogenous concentration is much smaller. We show how, in the limit of small endogenous concentration, the boundaries of the deposits become almost polygonal with self-similar structures at the corners, resembling a two-dimensional foam. We discuss subtle discrepancies between the Monge–Ampère solution, and a solution computed with the Eulerian particle-tracking method, indicating that surfactant transport can be considered almost, but not exactly, optimal. We analyse the two-dimensional mapping between the initial surfactant distribution and its equilibrium distribution, and discuss how the divergence and curl of the mapping can reveal regions of stretching, compression and rotational motion. Finally, we show that successive solutions of the Monge–Ampère equation, combined with divergence-free maps to mimic blowing, can be used to create a computational Suminagashi marbling pattern, illustrating the power of the optimal-transport approximation. Additional results are shown in supplementary material available at https://doi.org/10.1017/jfm.2024.334 to provide further evidence supporting the main findings and discussion presented in this paper.

2. Model and methods

2.1. The Eulerian particle-tracking problem

2.1.1. The problem and derivation of the model

We investigate the trajectories of particles on the surface of a viscous Newtonian liquid advected by surface tension gradients caused by a non-uniform concentration profile of insoluble surfactant, which is assumed to have negligible molecular diffusivity. Concentration gradients are caused by deposits of exogenous surfactant added to a uniform concentration field of endogenous surfactant. We assume that both species of surfactant have the same material properties, which combine to create a single concentration field that has a linear relationship with surface tension. A typical length scale is found from the initial size of an exogenous deposit, which is much greater than the initial height of the film. The thickness of the film is assumed to remain approximately uniform during the spreading, as we assume that any large vertical deflections are suppressed by gravity (in a large-Bond-number limit). The spreading takes place in a closed region with rectangular horizontal cross-section $\varOmega$, given in non-dimensional Cartesian coordinates as $0\leq x\leq L_1$, $0 \leq y\leq L_2$ confined by impermeable walls. Surfactant concentrations are scaled by the maximum initial concentration of one of the deposits. As explained in Appendix A, the surfactant is transported from its initial profile to its final equilibrium state via the nonlinear diffusion equation, which describes the evolution of the surfactant concentration as a function of space and time:

(2.1)\begin{equation} \varGamma_t = \tfrac{1}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\cdot} (\varGamma\, \boldsymbol{\nabla}_{\boldsymbol{x}}\varGamma), \end{equation}

where $\boldsymbol {\nabla }_{\boldsymbol {x}}$ is the gradient operator in the $\boldsymbol {x}= (x,y)$ plane of the Eulerian coordinates, and $\varGamma _t$ is the derivative of surfactant concentration with respect to non-dimensional time $t$; here, time is scaled by the ratio of liquid viscosity to maximum surface tension gradient (Appendix A). We impose a no-flux boundary condition at the periphery of the domain:

(2.2)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{x}}\varGamma \boldsymbol{\cdot} \boldsymbol{n}_b = 0 \quad \text{on } \partial \varOmega, \end{equation}

where $\partial \varOmega$ is the boundary of the domain $\varOmega$, and $\boldsymbol {n}_b$ is a unit normal vector to the boundary of the domain. Comparison of (2.1) with the non-dimensional surface transport equation $\varGamma _t+\boldsymbol {\nabla }_{\boldsymbol {x}}\boldsymbol {\cdot }(\boldsymbol {u}_s\varGamma )=0$ for a flat surface and non-diffusive surfactant shows that the surface velocity is

(2.3)\begin{equation} \boldsymbol{u}_s(\boldsymbol{x},t) ={-}\tfrac{1}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}} \varGamma. \end{equation}

Since we impose no flux of surfactant at the boundaries of $\varOmega$, as time goes to infinity, concentration gradients vanish to reach an equilibrium or steady state, so the initial concentration profile of surfactant $\varGamma _0(x,y)$ spreads to a uniform state with concentration $\bar {\varGamma }>\delta >0$, where $\delta$ is the initial endogenous concentration. We do not consider the singular limit $\delta =0$, which is beyond the scope of this study. In that case, spreading at the edges of the deposits would continue until the edges meet a solid boundary or the edges of another deposit. The final concentration relates to the initial concentration profile by

(2.4)\begin{equation} \bar{\varGamma}=\frac{\int_{\varOmega}\varGamma_0(x,y) \,\mathrm{d}\kern0.7pt x\, \mathrm{d}y}{\int_{\varOmega} \,\mathrm{d}\kern0.7pt x\, \mathrm{d}y} = \frac{1}{L_1L_2}\int_{\varOmega} \varGamma_0(x,y) \,\mathrm{d}\kern0.7pt x \,\mathrm{d}y . \end{equation}

Equation (2.1) represents a natural generalisation of the spatially one-dimensional nonlinear diffusion equation derived in Jensen & Halpern (1998), and aligns with the two-dimensional formulation of Thess et al. (1997). In stepping from one to two dimensions, an extra degree of freedom must be considered: any surface velocity field for which $\boldsymbol {u}_s\varGamma$ has zero divergence will not change surface concentrations but will nevertheless transport surface particles. This is illustrated in Appendix A by considering the influence of an imposed surface stress, as might arise from external blowing on the liquid film. For a monolayer close to equilibrium, the divergence of the stress field is area-changing; this is resisted by Marangoni effects (A10). However, the curl of the stress field in this simple model generates a flow that can redistribute surfactant (i.e. surface material elements carrying either endogenous or exogenous surfactant) without inducing surface tension gradients (A11). As well as being exploited by Suminagashi artists, this feature highlights a potential degeneracy in (2.1): namely, that the energetic cost of any flow that preserves concentrations of surface material elements is not captured by the evolution equation.

We now introduce a Lagrangian coordinate system $(x_0,y_0,\tau )$ to complement the Eulerian system $(x,y,t)$. We define a mapping $\boldsymbol {X}=(X,Y)$ between them, such that particles starting on the interface at $\boldsymbol {x}_0=(x_0,y_0)\in \varOmega$ at $t=0$ are advected at time $t=\tau$ to

(2.5a,b)\begin{equation} x = X(x_0,y_0,\tau), \quad y = Y(x_0,y_0,\tau). \end{equation}

Since surfactant transport is purely advective under (2.1), the mapping satisfies

(2.6)\begin{equation} \frac{\partial \boldsymbol{X}(x_0,y_0,\tau)}{\partial \tau} ={-}\frac{1}{4}\,\boldsymbol{\nabla}_{\boldsymbol{x}} \varGamma(\boldsymbol{X}(x_0,y_0,\tau),\tau). \end{equation}

The mapping function $\boldsymbol {X}(x_0,y_0,\tau )$ from initial to current particle location is the main quantity that we seek throughout this study. The initial conditions for each simulation that we perform in this study will be of the form

(2.7)\begin{equation} \varGamma_0(x_0,y_0) = \begin{cases} \delta + \mathcal{F}(x_0,y_0), & \text{in } \varOmega',\\ \delta, & \text{in } \varOmega -\varOmega', \end{cases} \end{equation}

where $\delta =\min {\varGamma _0(x_0,y_0)}$ for all $(x_0,y_0)$ in $\varOmega$ represents the initially uniform endogenous surfactant, and $\mathcal {F}$ is a function describing the initial distribution of exogenous surfactant deposited in $\varOmega '$, a region of $\varOmega$. In this study, we consider only non-overlapping depositions of exogenous surfactant that are axisymmetric about their own centre, and with a radially decreasing concentration profile. Although we have studied various initial distributions for the exogenous surfactant deposits (see the supplementary material), we focus on quadratic distributions, which we denote as

(2.8)\begin{equation} \mathcal{C}_q(\boldsymbol{x}_0;\boldsymbol{x}_c,r,\varGamma_{0,c}-\delta) = \begin{cases} (\varGamma_{0,c}-\delta)\left(1- \dfrac{|\boldsymbol{x}_0-\boldsymbol{x}_c|^2}{r^2}\right), & |\boldsymbol{x}_0-\boldsymbol{x}_c|\leq r,\\ 0, & |\boldsymbol{x}_0-\boldsymbol{x}_c|> r,\end{cases} \end{equation}

which is centred at $\boldsymbol {x}_c=(x_c,y_c)$, where the initial concentration has a local maximum $\varGamma _{0,c}$, with deposit radius $r$. The concentration profile $\varGamma _0$ is continuous when added to the endogenous field, and the Euclidean distance is given by $|\boldsymbol {x}_0-\boldsymbol {x}_c| \equiv \sqrt {(x_0-x_c)^2+(\kern0.7pt y_0-y_c)^2}$. The subscript $q$ in (2.8) refers to the quadratic nature of the initial concentration profile. In Appendix F and in § S5 of the supplementary material, we consider circular concentration profiles with other functional forms.

2.1.2. Scenarios studied

We have investigated scenarios involving one, two or three distinct deposits (i.e. $\varOmega '$ is constituted of one, two or three disconnected regions in $\varOmega$). The different configurations studied for the one- and two-deposit problems are presented in the supplementary material (see table S1). These two problems are helpful to understand basic dynamical features and the impact of the relevant non-dimensional parameters, as we will discuss briefly in § 3. However, the one- and two-deposit problems miss topological features that appear only with three or more exogenous deposits, such as internal corners where the edges of the deposits meet away from the domain boundaries. As we will discuss in § 3, internal corners display self-similar features. For the sake of simplicity and to enable analytical progress, we focus mainly on the three-deposit problem for the rest of this paper. Nevertheless, we anticipate that many of the results found with the three-deposit problem will also apply to problems involving more deposits. Therefore, we devise a model problem where $\mathcal {F}(x_0,y_0)$ consists of three circular regions of different radii ($r_1=1$, $r_2$ and $r_3$ in non-dimensional variables; see figure 1c) containing exogenous surfactant with quadratic concentration profiles, with differing non-dimensional maximum values $1$, $\varGamma _2$ and $\varGamma _3$ in the different regions (the number in the subscript corresponds to the region). Deposit 1, the smallest, is centred at $(x_1,y_1)$; the second largest circular deposit is centred at $(x_2,y_2)$; the largest is centred at $(x_3,y_3)$. Therefore, using our notation for circular deposits (2.8), we have

(2.9)\begin{equation} \mathcal{F}(x_0,y_0) = \mathcal{C}_q(x_1,y_1,1,1-\delta) + \mathcal{C}_q(x_2,y_2,r_2,\varGamma_2-\delta) + \mathcal{C}_q(x_3,y_3,r_3,\varGamma_3-\delta). \end{equation}

For the three-deposit problem, we choose $r_2=2$, $r_3=3$, $\varGamma _2=1$ and $\varGamma _3=2$. For every problem tackled in this paper and in the supplementary material, we choose $L_1=13$ and $L_2=11$. The centres of the deposits are chosen to be $(x_1,y_1) = (6,2)$, $(x_2,y_2) = (10,5)$ and $(x_3,y_3)=(4,7)$ for most of the solutions presented, unless otherwise stated.

2.1.3. Numerical scheme for the Eulerian particle-tracking problem

A finite-difference approximation of (2.1) and (2.6) is calculated using two rectangular grids. The first grid is used to solve for an approximation of (2.1) subject to boundary conditions (2.2) and initial conditions (2.7) and (2.9) in an Eulerian reference frame, which is accomplished using a second-order central differencing system in space, and a first-order forward Euler method in time (choosing a sufficiently small time step to ensure stability). This is solved simultaneously with a forward Euler approximation of (2.6) for the dynamics of the particle paths on a second grid in the Lagrangian reference frame. At each time step, the concentration gradient is approximated on the Eulerian grid, and interpolated onto the Lagrangian grid at the current particle locations using a linear interpolation method, meaning that the method as a whole is first-order in space and time.

The simulation is computed from $t=0$ to a large time $t=t_f$ when the solution approximates the steady state. The value of $t_f$ is found by considering the analysis in Appendix B, which shows how to ensure that the map is within a small tolerance vector $[X_{tol},Y_{tol}]^{\rm T}$ of the steady state everywhere (we set $[X_{tol},Y_{tol}]^{\rm T}=[10^{-3},10^{-3}]^{\rm T}$).

2.2. The Lagrangian particle-tracking problem

2.2.1. Derivation of the Lagrangian method

Rather than solving the three scalar PDEs in (2.1) and (2.6) in an Eulerian framework, it is sufficient to solve only two PDEs by adopting a Lagrangian framework, as we now demonstrate, by calculating $\boldsymbol {X}(x_0,y_0,\tau )$ without the intermediate step of determining surfactant concentrations. We present a Lagrangian scheme reminiscent of that presented by Carrillo, Matthes & Wolfram (2021) for a general Wasserstein gradient flow. The chain rule combined with (2.5a,b) yields the material derivative $\partial /\partial \tau |_{x_0,y_0} = \partial /\partial t |_{x,y}+ \boldsymbol {u}_s\boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}}$, where $\boldsymbol {u}_s = \boldsymbol {X}_{\tau }$, with the $\tau$ subscript meaning the partial derivative with respect to $\tau$. It is also the case that

(2.10) \begin{equation} \begin{pmatrix} \dfrac{\partial}{\partial x_0} \\ \dfrac{\partial}{\partial y_0} \end{pmatrix} = \begin{pmatrix} X_{x_0} & Y_{x_0} \\ X_{y_0} & Y_{y_0} \end{pmatrix} \begin{pmatrix} \dfrac{\partial}{\partial x} \\ \dfrac{\partial}{\partial y} \end{pmatrix}, \quad \text{or} \quad \boldsymbol{\nabla}_{\boldsymbol{x_0}} = (\boldsymbol{\nabla}_{\boldsymbol{x_0}} \boldsymbol{X} )^{\rm T}\,\boldsymbol{\nabla}_{\boldsymbol{x}}. \end{equation}

We define tensor calculus operators as

(2.11)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{x_0}} \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} a_{1x_0} & a_{1y_0} \\ a_{2x_0} & a_{2y_0} \end{pmatrix}, \quad \boldsymbol{\nabla}_{\boldsymbol{x_0}} \boldsymbol{\cdot} \begin{pmatrix} a_1 & a_2 \\ a_3 & a_4 \end{pmatrix} = \begin{pmatrix} a_{1x_0} + a_{3y_0} & a_{2x_0} +a_{4y_0} \end{pmatrix}. \end{equation}

The Jacobian of the mapping (2.5a,b),

(2.12)\begin{equation} \alpha \equiv \det(\boldsymbol{\nabla}_{\boldsymbol{x_0}}\boldsymbol{X}) = X_{x_0}Y_{y_0}-X_{y_0}Y_{x_0}, \end{equation}

quantifies how area elements are deformed by the map between initial and current particle positions, such that area elements $\mathrm {d} A_{\boldsymbol {x_0}}$ and $\mathrm {d} A_{\boldsymbol {x}}$ are related by $\mathrm {d} A_{\boldsymbol {x}} = \alpha \,\mathrm {d} A_{\boldsymbol {x_0}}$. By conservation of mass, we can equate integrals of the surfactant concentration over the Lagrangian and Eulerian domains, respectively:

(2.13)\begin{equation} \int_{\boldsymbol{X}^{{-}1}(\Delta \varOmega)} \varGamma_0(x_0,y_0) \,\mathrm{d} A_{\boldsymbol{x_0}} = \int_{\Delta \varOmega} \varGamma(\boldsymbol{X},t) \, \mathrm{d} A_{\boldsymbol{x}} , \end{equation}

where $\boldsymbol {X}^{-1}(\Delta \varOmega )$ is the pre-image of any subset $\Delta \varOmega$ of the Eulerian domain $\varOmega$, and there is a one-to-one mapping between the domains. Using the Jacobian of the mapping, we can change variables on the right-hand side of (2.13) to give

(2.14)\begin{align} \int_{\boldsymbol{X}^{{-}1}(\Delta \varOmega)} \varGamma_0(x_0,y_0) \,\mathrm{d} A_{\boldsymbol{x_0}} = \int_{\boldsymbol{X}^{{-}1}(\Delta \varOmega)} \varGamma(\boldsymbol{X}(x_0,y_0,\tau),\tau)\, \alpha(\boldsymbol{X}(x_0,y_0,\tau),\tau) \,\mathrm{d} A_{\boldsymbol{x_0}} . \end{align}

We are now integrating over the same space with respect to the same variables, and as $\Delta \varOmega$ is arbitrary, the integrands must be equal, yielding

(2.15)\begin{equation} \varGamma(\boldsymbol{X}(x_0,y_0,\tau),\tau)\,\alpha(\boldsymbol{X}(x_0,y_0,\tau),\tau) = \varGamma_0(x_0,y_0). \end{equation}

This is the main statement of mass conservation in $\varOmega$, valid for any $\tau \geq 0$, and is key for our analysis in this subsection and the next.

The choice of Lagrangian coordinate system is arbitrary, and in the rest of this subsection, we choose a spatially non-uniform coordinate system $(\xi, \eta )$. This coordinate system, non-uniform in $\varOmega$, also defines a geometric transformation of the domain $\varOmega$, which is achieved by deforming $\varOmega$ such that $(\xi, \eta )$ become regularly spaced Cartesian coordinates. We call this new domain the deformed Lagrangian domain, with coordinates $(\xi,\eta )$ replacing $(x_0,y_0)$. In § 2.3 we will revert back to $(x_0,y_0)$, which there will refer to regular Cartesian coordinates in an undeformed copy of the Eulerian domain such that $(x,y)=(x_0,y_0)$ at $\tau =0$. (These two domains will be referred to as the deformed and undeformed Lagrangian domains, respectively.) For now, however, we choose a coordinate system $(\xi,\eta )$ such that the initial surfactant concentration is uniform in the deformed domain, with $\varGamma _0(\xi,\eta )=1$ everywhere. This new coordinate system $(\xi, \eta )$ defines a geometric transformation of the rectangular domain $\varOmega$, such that surface areas are stretched or compressed until the concentration per unit area in the deformed system is $1$ everywhere. To illustrate, if a region of unit area has an initial uniform concentration of $0.25$ in the undeformed domain, then in the deformed domain it would have an area of $0.25$ and therefore an initial uniform concentration of $1$. In the coordinate system of the deformed domain, (2.15) becomes

(2.16)\begin{equation} \alpha(\boldsymbol{X}(\xi,\eta,\tau),\tau)\,\varGamma(\boldsymbol{X}(\xi,\eta,\tau),\tau) = 1. \end{equation}

With this choice, and using (2.5a,b), it follows that $\boldsymbol {\nabla }_{\boldsymbol {x}} (\alpha \varGamma ) = \alpha \, \boldsymbol {\nabla }_{\boldsymbol {x}} \varGamma + \varGamma \, \boldsymbol {\nabla }_{\boldsymbol {x}}\alpha = 0$, and so

(2.17)\begin{equation} \alpha\,\boldsymbol{\nabla}_{\boldsymbol{x}} \varGamma ={-}\varGamma\, \boldsymbol{\nabla}_{\boldsymbol{x}}\alpha ={-}\varGamma \left(\boldsymbol{\nabla}_{\boldsymbol{\xi}} \boldsymbol{X}\right)^{{-}{\rm T}} \boldsymbol{\nabla}_{\boldsymbol{\xi}}\alpha, \end{equation}

where $\alpha = \det (\boldsymbol {\nabla }_{\boldsymbol {\xi }}\boldsymbol {X})$, and $\boldsymbol {\nabla }_{\boldsymbol {\xi }} = [ \partial /\partial \xi, \partial /\partial \eta ]^{\rm T}$. The particle velocity (2.6) is given by $\boldsymbol {X}_{\tau } = - \boldsymbol {\nabla }_{\boldsymbol {x}} \varGamma /4$, so (2.12), (2.16) and (2.17) give

(2.18)\begin{equation} 4\alpha^2\left(\boldsymbol{\nabla}_{\boldsymbol{\xi}} \boldsymbol{X}\right)^{\rm T} \boldsymbol{X}_{\tau} = \boldsymbol{\nabla}_{\boldsymbol{\xi}}\alpha. \end{equation}

This expresses the time evolution of material particle locations in Eulerian coordinates as a function of the deformed Lagrangian coordinates. We can expand (2.18) as the system

(2.19a)$$\begin{gather} X_{\tau} = \frac{1}{4 \alpha^3} \left(\alpha_{\xi}Y_{\eta} - \alpha_{\eta}Y_{\xi}\right), \end{gather}$$

(2.19b)$$\begin{gather}Y_{\tau} = \frac{1}{4 \alpha^3} \left(\alpha_{\eta} X_{\xi} - \alpha_{\xi}X_{\eta}\right), \end{gather}$$

with $\alpha = X_{\xi } Y_{\eta }-X_{\eta }Y_{\xi }$. In turn, (2.19) can be rewritten as

(2.20)\begin{equation} \boldsymbol{X}^{\rm T}_{\tau} ={-}\frac{1}{8}\,\boldsymbol{\nabla}_{\boldsymbol{\xi}}\boldsymbol{\cdot} \left(\frac{1}{\alpha^2} \begin{pmatrix} Y_{\eta} & -X_{\eta} \\ -Y_{\xi} & X_{\xi} \end{pmatrix} \right), \end{equation}

which is an equation in divergence form that is easier to solve than (2.18) or (2.19) when using a finite-element method. Initial conditions are imposed via (2.16), so

(2.21)\begin{equation} \frac{1}{\varGamma(\boldsymbol{X}(\xi,\eta,0),0)} = X_{\xi}(\xi,\eta,0)\,Y_{\eta}(\xi,\eta,0) - X_{\eta}(\xi,\eta,0)\,Y_{\xi}(\xi,\eta,0) . \end{equation}

We choose $Y_{\eta }(\xi,\eta, 0 )=1$ and $Y_{\xi }(\xi,\eta, 0 )=0$, so that $X_{\xi }(\xi,\eta,0)= 1/\varGamma (\boldsymbol {X}(\xi,\eta,0),0)$. This yields a purely one-dimensional transformation, as illustrated in figure 2, from the undeformed to the deformed Lagrangian domain, simplifying the calculation of the deformed geometry. The initial conditions for $\xi$ are therefore obtained through

(2.22)\begin{equation} \int_0^x \varGamma_0(x',y) \,\mathrm{d}\kern0.7pt x' +C(\kern0.7pt y) = \xi(x,y,0), \quad \xi(0,y,0)=0. \end{equation}

Here, $C(\kern0.7pt y)$ is an arbitrary piecewise function chosen such that $\xi$ is continuous, and $1/X_{\xi } = \xi _x$ because we have fixed $Y$ and $t$. After finding the indefinite partial integral (2.22), we substitute $y=\eta$ and $x=X$, and invert (2.22) (numerically if needed) to find $X$ as an explicit function of $(\xi,\eta )$. Calling this solution $G(\xi,\eta )$, the initial conditions can be summarised as

(2.23a,b)\begin{equation} Y(\xi,\eta,0) = \eta, \quad X(\xi,\eta,0) = G(\xi,\eta). \end{equation}

Figure 2. (a) The Eulerian domain of the dynamic Lagrangian problem presented in § 2.2. This domain is broken into nine different regions (denoted R1 to R9) to compute the piecewise continuous definition of $\xi$, given in § S1 of the supplementary material. The red circles are the locations of the initial deposits of exogenous surfactant. (b) The deformed Lagrangian domain, calculated such that (2.16) holds for the Eulerian initial conditions (2.7) and (2.9) with the parameter choices taken in § 2.1.3. This is the domain in which we compute the numerical solution of (2.20) with boundary conditions (2.26).

2.2.2. The three-deposit problem

We illustrate the Lagrangian method introduced in § 2.2.1 by solving the model problem with the parameters outlined in § 2.1.2. For the initial conditions (2.7) and (2.9), the solution of (2.22) is

(2.24) \begin{align} \xi = \begin{cases} x- \left(\dfrac{(x-x_1)^3}{3}+x(\kern0.7pt y-y_1)^2\right)\left(1-\delta\right) +C_1(\kern0.7pt y), & |\boldsymbol{x}-\boldsymbol{x}_1| \leq 1,\\ \varGamma_2x- \dfrac{1}{r_2^2}\left(\dfrac{(x-x_2)^3}{3}+x(\kern0.7pt y-y_2)^2\right)\left(\varGamma_2-\delta\right) +C_2(\kern0.7pt y), & |\boldsymbol{x}-\boldsymbol{x}_2| \leq r_2,\\ \varGamma_3 x- \dfrac{1}{r_3^2}\left(\dfrac{(x-x_3)^3}{3}+x(\kern0.7pt y-y_3)^2\right)\left(\varGamma_3-\delta\right) +C_3(\kern0.7pt y), & |\boldsymbol{x}-\boldsymbol{x}_3| \leq r_3,\\ \delta x + C_4(\kern0.7pt y), & \text{everywhere else in } \varOmega. \end{cases} \end{align}

Here, $C_1(\kern0.7pt y)$, $C_2(\kern0.7pt y)$, $C_3(\kern0.7pt y)$ and $C_4(\kern0.7pt y)$ are determined for the choice $\varGamma _2=1$, $\varGamma _3=2$, $(x_1,y_1) = (6,2)$, $(x_2,y_2) = (10,5)$ and $(x_3,y_3)=(4,7)$ in § S1 of the supplementary material, along with the definition of the Lagrangian coordinates of the three circles. Finding this initial condition involves breaking the Lagrangian domain into nine regions, as shown in figure 2. By imposing $Y=\eta$, and imposing that the line $X=0$ corresponds to $\xi =0$, only the right-hand side of the Lagrangian domain, which we call $\partial \varOmega _R$ (defined for this problem in equation (S1.2) of the supplementary material), is not a straight line. We substitute $\eta =y$ into (2.24) and then invert (2.24) numerically to find the initial expression for $X$ as an explicit function of $\xi$ and $\eta$.

The boundary conditions for (2.20), and for the steady-state problem presented in the next subsection, are derived from the dynamic boundary condition (2.2). Analysis in Appendix C reveals that for corner angles less than ${\rm \pi}$, such as we have in the domain that we consider, a particle that begins on one of the four edges of the rectangle must stay on that edge for all time, and the appropriate boundary conditions accompanying (2.27) in the undeformed Lagrangian domain are the Dirichlet conditions

(2.25a,b)\begin{equation} X = 0, L_1 \text{ on } x_0=0,L_1 \quad \text{and} \quad Y=0, L_2 \text{ on } y_0=0,L_2, \end{equation}

which ensures that (2.2) is satisfied. This means that in the deformed Lagrangian domain,

(2.26)\begin{equation} X(0,\eta,\tau) = 0, \quad X(\xi,\eta,\tau) = L_1 \text{ on } \partial \varOmega_R, \quad Y(\xi,0,\tau) = 0, \quad Y(\xi,L_2,\tau) = L_2. \end{equation}

2.2.3. Numerical solution

Having inverted (2.24) numerically to find the initial conditions (2.23a,b), we use these initial conditions to solve (2.20) subject to boundary conditions (2.26), from $\tau =0$ to a final time taken to approximate the steady-state $\tau =t_f$, in the Lagrangian domain shown in figure 2(b), using COMSOL$^{\circledR}$. For reproducibility purposes we provide the details of the COMSOL$^{\circledR}$ settings chosen: we use the Mathematics suite, using the coefficient form PDE set-up that is designed to handle PDEs in divergence form such as (2.20). We discretise using standard COMSOL$^{\circledR}$ triangulation method, and we use quadratic Lagrange basis functions with $314\,198$ degrees of freedom plus $16\,578$ internal degrees of freedom, and set the relative tolerance to $10^{-9}$. We store the solution at every $2$ time units.

2.3. The steady-state problem

2.3.1. Formulation of the problem

We now consider the problem of approximating the equilibrium locations of surface particles (their locations as $t\to \infty$) as a function of their initial locations directly, i.e. without any intermediate calculation of surfactant concentrations, or intermediate calculation of particle trajectories. We return to the coordinate systems used in § 2.2; however, here we revert to calling the Lagrangian coordinates $(x_0,y_0)$ to indicate that the Lagrangian domain is a copy of the Eulerian domain, defined as $0\leq x_0\leq L_1$ and $0\leq y_0\leq L_2$, and $(X,Y) = (x_0,y_0)$ at $t=\tau =0$. Using these variables, at steady state, (2.15) becomes

(2.27)\begin{equation} X_{x_0}Y_{y_0} - X_{y_0}Y_{x_0} = \frac{\varGamma_0(x_0,y_0)}{\bar{\varGamma}}, \end{equation}

which is a PDE describing the mapping function $(X,Y)$ to the spatial coordinates $(x,y)$ for particles starting at $(x_0,y_0)$, in the limit $t\to \infty$. Equation (2.27) needs to be solved subject to boundary conditions (2.25a,b).

For one-dimensional problems, e.g. $\varGamma _0=\varGamma _0(x_0)$, we can impose $Y=y_0$ and (2.27) has a unique solution. However, in two dimensions, (2.27) constitutes only one equation for the two unknowns $(X,Y)$, and therefore does not have a unique solution, so we turn to the Helmholtz decomposition theorem to make progress. By this theorem, we know that we can write the map $\boldsymbol {X}=[X,Y]^{\rm T}$ in terms of two scalar potentials $\phi (x_0,y_0)$ and $\psi (x_0,y_0)$ such that

(2.28)\begin{equation} [X,Y]^{\rm T} = \boldsymbol{\nabla}_{\boldsymbol{x}_0}\phi + \boldsymbol{\nabla}_{\boldsymbol{x}_0}\times\boldsymbol{\psi}, \end{equation}

where $\boldsymbol {\psi }$ is a vector of magnitude $\psi$ pointing out of the plane (in the $z$-direction), with $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } \boldsymbol {X}=\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi$ and $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times \boldsymbol {X} = - \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\boldsymbol {\psi }$. To make this Helmholtz decomposition unique up to constants, we impose the boundary conditions

(2.29)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{x}_0}\phi \boldsymbol{\cdot} \boldsymbol{n}_b = [x_0,y_0]^{\rm T}\boldsymbol{\cdot} \boldsymbol{n}_b, \quad \boldsymbol{\nabla}_{\boldsymbol{x}_0}\times \boldsymbol{\psi}\boldsymbol{\cdot} \boldsymbol{n}_b = 0 \quad \text{on } \partial \varOmega , \end{equation}

which satisfies (2.25a,b).

The map at time $t$ is generated by (2.6), the right-hand-side of which is an Eulerian gradient of the instantaneous surfactant concentration. Thus the map remains irrotational with respect to the Eulerian coordinates. Now we investigate whether the map at time $t$ can be approximated by a map that is irrotational with respect to the Lagrangian coordinates, as this would allow us to remove the indeterminacy in (2.27), since $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times [X,Y]^{\rm T} =\boldsymbol {0}$ yields $\psi$ equal to a constant, reducing the problem (2.27) to finding a solution for a single scalar potential $\phi$. We summarise the statement that we want to test as that, for all time $t$,

(2.30)\begin{equation} |\boldsymbol{\nabla}_{\boldsymbol{x}_0}\phi| \gg |\boldsymbol{\nabla}_{\boldsymbol{x}_0}\times\boldsymbol{\psi}|. \end{equation}

In effect, we test the idea that because the Eulerian curl of $\boldsymbol {u}_s$ is zero, and material particles on boundaries are not allowed to traverse corners, (2.30) might hold for all time, at least when the rearrangement of the surface is small. We will test this hypothesis a posteriori in § 3.

Assuming that the map (2.28) is given by $[X,Y]^{\rm T} = \boldsymbol {\nabla }_{\boldsymbol {x}_0}\phi$, (2.27) and boundary conditions (2.25a,b) reduce to the Monge–Ampère equation

(2.31)\begin{equation} \phi_{x_0 x_0}\phi_{y_0 y_0} - \phi_{x_0 y_0}^2 = \frac{\varGamma_0(x_0, y_0)}{\bar{\varGamma}} \quad \text{on } 0\leq x_0 \leq L_1, 0\leq y_0 \leq L_2 , \end{equation}

subject to

(2.32)\begin{equation} \phi_{x_0} = x_0 \text{ on } x_0 = 0,L_1, \quad \phi_{y_0}=y_0 \text{ on } y_0=0,L_2, \quad \phi(0,0)=0. \end{equation}

The last boundary condition is necessary to close the problem, as $\phi$ is unique only up to a constant. The Monge–Ampère equation arises often in the theory of optimal transport, a connection that we will discuss further in § 4.

2.3.2. Numerical method

We solve (2.31) subject to the boundary conditions (2.32) for the initial concentration profile of surfactant (2.7) and (2.9) using an iterative Newton–Raphson scheme for a finite-difference approximation of the solution, the full details of which are in Appendix D. The Newton–Raphson scheme converges to the desired solution only if the initial guess is in the basin of attraction of the desired solution, which for a nonlinear problem such as (2.31) and (2.32) is difficult to determine a priori. We surmount this problem with the following continuation scheme. Using a parameter $\beta _j \in [0,1]$, we take advantage of the fact that the PDE

(2.33)\begin{equation} \phi_{x_0 x_0}\phi_{y_0 y_0} - \phi_{x_0 y_0}^2 - 1 + \beta_j \left(1 - \frac{\varGamma_0(x_0,y_0)}{\bar{\varGamma}} \right) = 0, \end{equation}

subject to boundary conditions (2.32), has a known solution when $\beta _j = 0$, namely $\phi = x_0^2/2 + y_0^2/2$; when $\beta _j=1$, we have the desired solution to (2.31) and (2.32). We step from $\beta _0=0$ to $\beta _J=1$, in steps of some fixed quantity $\Delta \beta =1/J$ (where $J$ is an integer), solving (2.33) and (2.32) each time. Starting from $\beta _0=0$ and $\phi _0=x_0^2/2+y_0^2/2$, we find $\phi _{j+1}$ by using $\phi _j$ as a guess solution for (2.33) and (2.32), where $\beta _{j+1}=\beta _j +\Delta \beta$. If we choose $\Delta \beta$ to be small enough, then we ensure that we stay inside the basin of attraction of solutions, finding the desired solution to (2.31) and (2.32) when $j=J$.

We use this process to solve (2.31) and (2.32) for intermediate and low values of endogenous surfactant, $\delta =0.25$ and $\delta = 0.002$. We solve for the larger value of $\delta$ using a grid with grid points spaced uniformly 0.05 units apart in MATLAB$^{\circledR}$, using the software's ‘sparse’ variable type to handle the large sparse matrices, and its efficient algorithms for finding solutions to linear systems such as (D3) with a direct LU factorisation scheme. This solution is obtained by using $\Delta \beta =0.1$. For the solution with the smaller value of $\delta$, we use grid points spaced evenly 0.05 units apart. We need $\Delta \beta =0.0025$ for this second solution, which means that the computational cost is increased. The convergence of the numerical scheme is presented in § D.2. In addition, we present a method for creating a computational Suminagashi picture in Appendix E.

To quantify how well the Monge–Ampère method approximates the solution found by the Eulerian particle-tracking method at $t=t_f$ (assumed to be an accurate solution of the steady state), we define metrics that characterise the difference between solutions found using the two methods for the same initial conditions. We define the Euclidean distance between final particle locations $X_{EU}$ and $X_{MA}$ predicted by both methods and normalised by the longest side of the domain,

(2.34)\begin{equation} \frac{1}{L_1}\,|{X}_{EU}-{X}_{MA}| \equiv \frac{1}{L_1}\,\sqrt{(X_{EU}-X_{MA})^2+(Y_{EU}-Y_{MA})^2}, \end{equation}

which we call the normalised absolute error between the two methods for a given initial particle location. Statistics of the error are then obtained by analysing distributions for a large number of the initial particle locations, particularly the median, the upper quartile, the 90th percentile and the maximum values of (2.34).

3. Results

Table 1 summarises all of the simulations and their parameters that are presented in the results section, with a key with which we refer to each simulation.

Table 1. A table presenting a summary of the simulations presented in § 3, together with parameters used, and a key with which we refer to each simulation. The methods used are the Eulerian particle-tracking method (2.1) and (2.6), the Lagrangian particle-tracking method (2.20), and the Monge–Ampère method (2.31). For all these simulations, we choose $r_2=2$, $r_3=3$, $\varGamma _2=1$ and $\varGamma _3=2$.

3.1. Particle-tracking solutions

The results for the Eulerian (Eul[0.25]) and Lagrangian (Lag[0.25]) particle-tracking methods (presented in §§ 2.1.3 and 2.2.3) with $\delta =0.25$ are shown in figures 3(a) and 3(b), respectively, and also as supplementary movies 1 and 2, where each thin coloured line represents a particle trajectory, terminating at a black dot at $t=t_f$. The trajectories shown in figure 3(a) represent 1 in every 225 trajectories calculated, selected such that their initial locations are evenly spaced. The data obtained by the solution for the Lagrangian method presented in figure 2(b) are spaced irregularly, with each data point corresponding to a node of the mesh used in COMSOL$^{\circledR}$ to discretise the deformed Lagrangian domain; we display 1 in every 50 particles from the data list obtained from the simulation, so the density of particles shown is not significant.

Figure 3. Solution to the example problem of three circular deposits of exogenous surfactant spreading together with $\delta =0.25$. (a) The results of Eul[0.25] (see supplementary movie 1). The initial boundaries of the exogenous surfactant circular deposits are the black circles, and the final locations are the thick red lines. The points described by the green curves map to the black circles at $t=t_f$. Individual particle trajectories are plotted using thin coloured lines terminating at black points. (b) The results of Lag[0.25] with the same colour scheme as in (a) (see supplementary movie 2). The particles represent $1/50$ of all the particle trajectories calculated, which are chosen at random, so the density of particles shown is not significant. (c) Graph showing the results of MA[0.25] overlaid onto Eul[0.25] and Lag[0.25]. The steady-state boundaries of the three deposits and the curves found by the inverse map (which spread from and to the black circles in the steady state, respectively) are given by the colour scheme shown in the figure legend.

In figures 3(a) and 3(b), the largest deposit spreads out through Marangoni stresses, and compresses the other two deposits. Flow reversals (sharp turns of particle trajectories of more than $90^{\circ }$) arise in several areas for two reasons. First, reversals in the top left-hand corner are due to confinement. Early outward spreading is into a region containing endogenous surfactant at low concentration $\delta$; later reversals arise once the surfactant concentration in this region is much larger due to non-local compression of the endogenous material. Second, points that begin on the edges of the smaller two deposits nearest the centre of the domain first spread into the centre, but soon the effect of the largest deposit spreading is felt, and these points reverse their trajectories. The final shapes of the smallest deposits are non-trivial oval shapes, the centres of which are shifted away from their initial locations. Some particles to the top left of the centre of the largest deposit traverse distances close to $1$ unit in length and then move approximately the same distance back, close to where the particles started. Particles compress into the top and bottom right-hand corners. A variety of trajectories are evident: for example, particles in the top and left have trajectories that involve straight lines and sharp turns, whereas particles towards the bottom right describe gentle arcs. In figure S1 in § S2 of the supplementary material, we present an overlay of the contour plots of $X(x_0,y_0)$ and $Y(x_0,y_0)$ for the solutions at $t=t_f$ obtained from the Eulerian and Lagrangian particle-tracking methods, respectively. The methods find the same particle locations to within a distance of $0.05$ almost everywhere, apart from the locations of small oscillations in the Lagrangian solution that appear to be an artefact of the domain deformation as discussed in § S2, and much closer than that in most places. Some of the small discrepancies that do exist can be explained partly by the fact that small errors arise by interpolating the Lagrangian solution onto a regular, rectangular grid to make the comparison, and errors occur in the Eulerian solution by the interpolation of the gradient of the evolving concentration shown in § S3 (figure S2) of the supplementary material at every time step in that solution.

3.2. The steady-state solution

We investigated the steady-state solution of a variety of configurations involving one and two deposits (see table S1 in the supplementary material). We compare the results obtained between Eulerian (EU) particle tracking and the Monge–Ampère (MA) method in § S5 of the supplementary material. We tested how the final equilibrium shape of the deposits is influenced by the proximity of the domain boundaries. In the case of a single initial deposit, we find only small differences in the discrepancy between the EU and MA results, quantified using (2.34) for deposit locations at various distances from the domain boundaries. The median normalised absolute error is approximately $10^{-4}$, and the maximum error is bounded by $2\times 10^{-3}$. In the case of two initial deposits, the median normalised absolute error is approximately $5\times 10^{-4}$, and the maximum error is bounded by $5\times 10^{-3}$. The median discrepancy between the two methods tends to be inversely correlated with the symmetry of the initial configuration, whereas the upper quartile, 90th percentile and maximum discrepancy are much noisier for both the one-deposit and two-deposit problems studied. Discrepancies between the EU and MA methods increase with an increase in the number of deposits, and with a decrease in $\delta$ (the normalised initial endogenous surfactant concentration), as shown in Appendix F. As stated previously, we choose to focus on the three-deposit case.

3.2.1. The three-deposit problem with $\delta =0.25$

The solution for the approximation of the edges of the three deposits in the steady state (MA[0.25]) found by the Monge–Ampère method (outlined in § 2.3.2) is presented in figure 3(c) for $\delta = 0.25$, with the approximations of the steady state found from the Eulerian (Eul[0.25]) and Lagrangian (Lag[0.25]) particle-tracking solutions overlaid. The final edges of the deposits predicted by MA[0.25] are almost indistinguishable except in a few places, which supports assumption (2.30). Predictions of Eul[0.25] and Lag[0.25] are indistinguishable to the naked eye in figure S1 of the supplementary material, providing a reliable benchmark against which to test the prediction of MA[0.25]. We also use the inverse maps to calculate the contours that map to the initial drop boundaries under the spreading in figure 3; again, only very small discrepancies between MA[0.25] and Eul[0.25] are evident.

A comparison of the global behaviour of the Monge–Ampère approximation of the map from initial to final particle configuration (MA[0.25]) with the map calculated from the particle-tracking solution (Eul[0.25]) is given by contour plots in figure 4 (see also a colour map of the absolute error between the two predictions for the final particle location in § S4, figure S3, of the supplementary material). Dense contours in figures 4(a,b) indicate that surface areas starting at these locations are stretched by the mapping, and similarly large gaps between contours indicate that the map compresses the surface. Conversely, in figures 4(c,d), dense contours of the inverse maps indicate that surface areas finishing at these locations have been compressed by the spreading, and large gaps between contours indicate that the spreading has stretched the surface. The Monge–Ampère approximation agrees with the particle-tracking solution in most places, although some noticeable discrepancies exist, such as in the left half of the smallest deposit most notably. The median error across the solution is approximately 0.25 % of the domain length, and the error for every particle is within $1.5$% of the domain length, as shown in figure S7 of the supplementary material.

Figure 4. Contour plots of solutions for the map from initial configuration to steady state found from MA[0.25] and Eul[0.25]. (a) 25 evenly spaced contours of $X_{MA}$ taken from MA[0.25] (red) overlaid with the same-valued contours of $X_{EU}$ taken from Eul[0.25] (blue). (b) 25 contours of $Y_{MA}$ (red) overlaid with $Y_{EU}$ (blue). (c) The inverse map $X_{MA}^{-1}$ (red) overlaid with $X_{EU}^{-1}$ (blue), with the same contour scheme. (d) The same for $Y_{MA}^{-1}$ (red) overlaid with $Y_{EU}^{-1}$ (blue).

To illustrate where and how the discrepancies arise, figure 5 shows the divergence and curl of the map found from Eul[0.25], shown in both Lagrangian and Eulerian coordinate systems. From the Helmholtz decomposition (2.28), figures 5(a,c) show $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)=\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$, and figures 5(b,d) show $(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$, where $(\cdot )_{\perp }$ means the $z$-component, perpendicular to the plane of the solution. The vector field $\boldsymbol {X}-\boldsymbol {x}_0$, which points from initial to final particle locations, is an easier quantity to interpret physically than $\boldsymbol {X}$ itself. Given boundary conditions (2.29) (the boundary condition for $\psi$ can be taken to be equivalent to the Dirichlet condition $\psi =0$ on all four boundaries), the fact that $|\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi |$ is an order of magnitude greater than $| \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi |$ almost everywhere is further evidence justifying our assumption (2.30) (the ratio $\nabla ^2_{x_0}\psi /\nabla ^2_{x_0}\phi$ is plotted in § S4, figure S4, of the supplementary material), although it is certainly not the case that the curl of the map vanishes. In figures 5(a,c), $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)<0$ represents surface areas with net compression, and $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)>0$ represents areas with net expansion by the map. Areas within the deposits expand, as do area elements connecting the largest deposit with the smaller deposits, and the corner regions compress. Saddle-like area elements directly between each of the initial deposit locations, and between each deposit and the nearest boundary, are compressed and expanded in orthogonal directions. In figures 5(b,d), positive values of $(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp }$ refer to anticlockwise net local rotation (twist) by the map, and negative values for clockwise twist. It is notable that weak twisting motions arise where interfaces spread towards a nearby boundary, or near another drop interface, with regions of oppositely oriented twisting typically appearing in pairs. The most intense twisting appears to be confined to regions immediately outside the boundaries of the exogenous surfactant drops.

Figure 5. Contour plots showing the divergence and curl of the vector field from initial to final particle location taken from Eul[0.25]. (a) Plot of $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0)= \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Lagrangian coordinates. (b) Plot of $(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Lagrangian coordinates. (c) Plot of $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0)= \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Eulerian coordinates. (d) Plot of $(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Eulerian coordinates.

3.2.2. Weak endogenous surfactant

Many fluids in the environment have low levels of contaminant surfactant, and in some cases (especially in controlled laboratory conditions) vanishingly small endogenous surfactant concentrations, so we would like to understand the Lagrangian motion of surface particles in the limit of very small $\delta$. The Lagrangian dynamic method (§ 2.2.3) is not capable of handling small values of $\delta$, because the deformed Lagrangian domain becomes extremely stretched. The Eulerian particle-tracking method (§ 2.1.3), and the Monge–Ampère approximation (§ 2.3.2), however, can both find well-behaved results with extremely small $\delta$.

We present the solution to the Monge–Ampère approximation MA[0.002] overlaid with the solution for the Eulerian particle-tracking method Eul[0.002] in figure 6. Figures 6(a,b) show the overlaid contour plots of the approximations for the steady-state solutions for $X$ and $Y$ in the Lagrangian coordinates, and figures 6(c,d) show the inverse mapping in the Eulerian coordinates. Figure 6(e) shows an overlay of the prediction for the final boundary locations for the three deposits using both methods.

Figure 6. Plots showing a comparison between MA[0.002] and Eul[0.002]. (a) Contour plots of $X_{MA}$ taken from MA[0.002] (red) overlaid with the same-valued contours from $X_{EU}$ (blue) taken from Eul[0.002]. (b) The same for $Y_{MA}$ (red) and $Y_{EU}$ (blue). (c) The same scheme for the inverse maps $X_{MA}^{-1}$ (red) and $X_{EU}^{-1}$ (blue). (d) Similarly for $Y_{MA}^{-1}$ (red) and $Y_{EU}^{-1}$ (blue). (e) An overlay of the final deposit boundaries from MA[0.002] (red) and Eul[0.002] (blue). ( f) Particle trajectories, each given by a thin coloured line terminating at a black dot, from Eul[0.002] (see supplementary movie 3). The three red dashed ellipses each contain a complete particle trajectory that involves two sharp changes of direction.

Figure 6 shows that when $\delta$ is small, the Monge–Ampère method does less well in approximating the solution, which is expected as particles spread further, leading to large discrepancies between the Lagrangian and Eulerian curls of the maps; however, it is still credible as a first-order approximation. In particular, for particles that are initially located within the three deposits, the approximation is still accurate, with the largest discrepancy occurring for particles of endogenous surfactant. This is shown clearly when comparing figures 6(c,d), which show the inverse map in the Eulerian coordinates, to figures 6(a,b), which show the map in the Lagrangian coordinates. (This is also shown clearly in the plot of the absolute and relative errors between the predictions for final particle locations in § S4, figures S3(b,d), of the supplementary material.) Contours coincide over an appreciably greater region of the final configuration in comparison to the initial configuration. In Appendix F, we present box-and-whisker plots of the absolute error between final particle locations predicted by the Monge–Ampère and Eulerian particle-tracking methods for several values of $\delta$, which show how the maximum discrepancy grows as $\delta \to 0$ (10 % of the domain length for $\delta =0.002$), but the median discrepancy remains an order of magnitude or more smaller than the maximum discrepancy for $\delta =0.002$.

The final locations of particles beginning on the region occupied by endogenous surfactant are compressed into effectively three lines, as shown by the blue curves in figure 6(e) for Eul[0.002]. The Monge–Ampère approximation (MA[0.002]) for the location of the lines (in red) is relatively accurate in most places. Figure 6( f) and supplementary movie 3 show the particle trajectories from the Eulerian particle-tracking method. The variety of trajectories is remarkable, with many particles having two sharp changes in direction during the spreading (see, for example, the three trajectories highlighted by red dashed ellipses). This likely reflects the fact that particle trajectories can be influenced by different deposits at different times.

Figure 7 shows the divergence and curl of the map computed from Eul[0.002]. The fact that the Monge–Ampère approximation is less accurate for small $\delta$ is reflected in the fact that $|\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi |$ and $|\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi |$ are shown to be of the same order of magnitude in certain places within the initial configuration (the ratio $\nabla ^2_{x_0}\psi /\nabla ^2_{x_0}\phi$ is plotted in figure S4 in § S4 of the supplementary material). However, figures 7(c,d), plotted in the final configuration, once again show that, for particles that start within the deposits, it is still the case that (2.30) holds. Once again, $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)<0$ represents area elements with net compression, and $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot } (\boldsymbol {X}-\boldsymbol {x}_0)>0$ represents area elements with net expansion by the map. Areas within the deposits expand, as do area elements connecting the largest deposit with the smaller deposits, and we also now see expansion for areas connecting the deposits with the boundaries. In figure 7(b), the patterns created by $(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp }$ are similar to the case where $\delta =0.25$, although the modulus of twist is greater in magnitude.

Figure 7. Contour plots showing the divergence and curl of the vector field from initial to final particle location taken from Eul[0.002]. (a) Plot of $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0)= \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Lagrangian coordinates. (b) Plot of $(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Lagrangian coordinates. (c) Plot of $\boldsymbol {\nabla }_{\boldsymbol {x}_0}\boldsymbol {\cdot }(\boldsymbol {X}-\boldsymbol {x}_0) = \boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\phi -2$ in Eulerian coordinates. (d) Plot of $(\boldsymbol {\nabla }_{\boldsymbol {x}_0}\times (\boldsymbol {X}-\boldsymbol {x}_0))_{\perp } = -\boldsymbol {\nabla }_{\boldsymbol {x}_0}^2\psi$ in Eulerian coordinates.

3.2.3. Self-similarity of corner regions

Computationally, the Monge–Ampère solution is significantly cheaper to solve compared to either of the particle-tracking solutions. We take advantage of this to analyse the shape of the corners of deposits as $\delta \to 0$ using finely discretised Monge–Ampère approximations. We note that the errors found between the Eulerian and Monge–Ampère predictions are small near the corner regions, even with $\delta =0.002$ (see figure S3 in the supplementary material). This suggests that predictions from the Monge–Ampère solution remain accurate in the corner regions in the limit $\delta \to 0$. Examples of these corner regions are numbered in green in figure 8(a), which shows the results of MA[0.04]. To compute the deposit edges as sufficiently smooth curves to accurately analyse curvature of the edges was not possible from the particle-tracking solutions with available computational power, but is possible with the Monge–Ampère method. During spreading, endogenous surfactant at initial concentration $\delta$, occupying an $O(1)$ area, is compressed into narrow threads (between pairs of drops) and into what resemble seven Plateau borders (at corners between drops, and at the domain boundary), at final concentration of order unity. Assuming that the bulk of the endogenous surfactant is driven into these seven regions, and noting that each has an $O(1)$ aspect ratio, we can expect each drop corner to curve over a length scale of $O(\delta ^{1/2})$. This scaling is motivated by the observation that an area with an $O(1)$ length scale and with an initial concentration of endogenous surfactant $\delta$ (giving a total mass of $O(\delta )$) is ultimately compressed into an area with final concentration $\bar {\varGamma }\sim 1$ over an $O(\delta ^{1/2})$ length scale (preserving the total mass of $O(\delta )$).

Figure 8. Scaled curvature plots for a selection of the corners of the boundaries of the three circular deposits at the steady state for small values of $\delta$. (a) A graph of the solution from MA[0.04] that shows in green our numbering system for corners inside each deposit. (b–f) Plots of the natural logarithm of the curvature scaled by $\delta ^{1/2}$ against the arc length scaled by $\delta ^{-1/2}$, where $s=0$ identifies the vertex of the corner in each case.

Figure 8 shows how the equilibrium shapes of the deposits calculated by Monge–Ampère become more and more polygonal in the limit $\delta \to 0$, with corner regions adopting a self-similar form. For several values of small $\delta$, we calculate the curvature $K$ of the equilibrium boundary of each deposit as a function of arc length $s$. To show how each of a selection of corners sharpens as $\delta \to 0$, we set $s=0$ at the local curvature maximum, and plot $\log {(\delta ^{1/2}K)}$ against $s\delta ^{1/2}$. Figures 8(b–f) show collapse of the data, illustrating how, as $\delta \to 0$, the curvature of corners of the polygons (numbered in green in figure 8a) is proportional to $\delta ^{-1/2}$, with this curvature occurring over arc lengths proportional to $\delta ^{1/2}$, and the edges of the deposits becoming effectively straight lines away from the corners. We see that in small regions away from $s=0$, the quantity $\ln {(\delta ^{1/2}K)}$ becomes linear, indicating a functional dependence such as $\delta ^{1/2}K \sim \exp (-\lambda \delta ^{-1/2}|s|)$ for some constant $\lambda$. In practice, the boundaries between the deposits have a small curvature in the Eulerian particle-tracking solution (figure 6e), so the deposits instead resemble slightly distorted polygons, with corner regions slightly distorted from the Monge–Ampère prediction. As discussed above, the self-similarity observed for the internal corners, and the appearance of a characteristic length scale of order $\delta ^{1/2}$, are based on the fundamental physical principle of mass conservation. Since this principle is independent of the choice of method, we expect that the results observed using the Monge–Ampère method will remain valid with the Eulerian and Lagrangian methods. This provides further evidence that the self-similarity behaviour observed in figure 8 is physically grounded.

Figure 9 shows the Monge–Ampère approximation of the three-deposit problem with $\delta = 0.005$ (MA[0.005]Alt1–Alt6), using multiple different centres of the deposits (given in the caption), revealing a variety of near-polygonal structures. Each deposit approaches a triangle, a quadrilateral, or in the case of the largest deposit (figures 9b,d), a hexagon. The final configuration approaches a structure that can be characterised by up to four coordinates, i.e. the locations of corners shared by multiple deposits, which we call the characteristic coordinates of a given initial configuration. As in figure 8, we expect the Monge–Ampère approximation to provide a good first-order estimate of equilibrium drop shape predicted by the full dynamic problem. However, at the present time we are unable to offer any simple strategy for determining the characteristic coordinates directly from the initial conditions.

Figure 9. Plots showing the steady-state mapping and inverse mapping for $\delta =0.005$ with varying initial locations for some of the circular deposits 1 and 2 (deposit 3, as shown in figure 8a, remains fixed). Initial locations of the boundaries of the circular deposits are in blue, with the final locations of those boundaries in red. The green curves map to the blue circles under the same map. (a–f) The results for MA[0.005]Alt1 to MA[0.005]Alt6, respectively, from the key shown in table 1.

4. Discussion

In this study, we have evaluated the trajectories of passive surface particles confined in a rectangular domain under the action of surfactant spreading on a thin film in a large-gravity and large-Péclet-number limit using two separate methods (Eulerian particle tracking, outlined in § 2.1, and Lagrangian particle tracking, outlined in § 2.2), the results of which corroborate each other. We have also identified a direct method (using the Monge–Ampère equation outlined in § 2.3) for approximating the equilibrium configuration as $t\to \infty$. The solutions to this problem show how confinement and drop–drop interactions can lead to drift of surface particles (in addition to spreading), and how drop–boundary interactions lead to transient flow reversals. It is striking the way that some particles move a significant distance from their initial location, before moving back close to where they started, and that in figure 6( f) some particles sharply change directions twice, which is reminiscent of the multiple regimes identified for multiple surfactant sources noted by Iasella et al. (2024). Predicting the equilibrium location of the interface between exogenous and endogenous surfactants is straightforward for spreading in one spatial dimension (relying on mass conservation arguments). However, in two-dimensions (as here) it becomes a non-trivial task, even for simple scenarios with one or two initial deposits, particularly when symmetry is broken. With three initial deposits, new topological features appear in the form of internal corners between deposits. Nevertheless, equilibrium shapes possess asymptotic structures, such as the self-similar geometry of the internal corners in the limit of small initial endogenous surfactant concentration (figures 8 and 9). Such observations might offer a route to explicit predictions.

In order to compute the particle trajectories involved in the surfactant-driven spreading, we constructed a finite-difference numerical method on a regular grid. We used an interpolated concentration gradient to determine the local velocity of the particles. The surfactant concentration gradients were calculated using a finite-difference solution from a second grid. In addition, we reformulated the Eulerian surfactant transport equation (2.1) and the Lagrangian particle transport equation (2.6) into a single Lagrangian vector equation (2.20), which describes the dynamics of material particles. This was achieved by choosing a Lagrangian coordinate system for which the surfactant concentration is initially uniform in the deformed Lagrangian domain (i.e. having uniform mass per unit Lagrangian area), and choosing this domain such that the initial conditions are simple to compute from (2.21). Although the calculation of surfactant concentration is bypassed by this approach, the concentration can still be recovered at any time from the Jacobian via (2.16).

From a computational point of view, the Eulerian and Lagrangian numerical methods used to compute the particle trajectories each have some advantages and drawbacks. They are similar in computational expense; however, the Eulerian method introduces additional errors at every time step by first approximating $\boldsymbol {\nabla }_{\boldsymbol {x}}\varGamma$ on the Eulerian grid, and then interpolating this quantity onto the current position of surface particles on the Lagrangian grid. The Eulerian method also requires storage of two solutions (Eulerian surfactant concentration and Lagrangian particle position) on two separate grids, whereas the Lagrangian method needs storage of only a single solution. Some of the drawbacks of the Lagrangian method are that: the deformed domain needs to be calculated beforehand (although this needs to be done only once); the Lagrangian domain has an irregular shape, which makes numerical resolution of (2.20) difficult; and the domain can become so deformed for small $\delta$ that numerical resolution of (2.20) is not practical. For moderate $\delta$, however, the Lagrangian formulation is a mathematically elegant approach that finds the solution without needing to introduce interpolation errors, and the close agreement of the results between the two methods (figures 4 and S1) provides corroboration of the particle-tracking method, lending support to our calculations of Eulerian solutions for small $\delta$.

Since equilibrium drop shapes, obtained in the limit of large time, can be computationally expensive to calculate with the Eulerian and Lagrangian methods, we have also explored how a Monge–Ampère equation can be used to find approximations of the equilibrium configurations. This was achieved by approximating the map between initial and final configuration as the gradient of a scalar potential, and neglecting the rotational part of the map. The computational benefits of this approach are beyond dispute, as only a single, scalar, time-independent PDE needs to be solved. For example, the Monge–Ampère solution shown in figure 4 can be computed for approximately $2.3\times 10^7$ grid points in a shorter time than a solution for approximately $3\times 10^5$ grid points using the Eulerian particle-tracking method from § 2.2.3 for the same parameters. Beyond computational issues, the Monge–Ampère method reveals a counterintuitive feature of surfactant-spreading dynamics in two dimensions: despite always being transported by an instantaneously irrotational flow (2.3), in the Eulerian sense, the mapping of particles from their initial to their current state can accumulate a weak rotational component with respect to Lagrangian coordinates, as illustrated in figures 5(b,d) and 7(b,d). This small Lagrangian rotational component can lead to weak distortions of final equilibrium configurations between the Monge–Ampère predictions and the Eulerian and Lagrangian predictions (figure 4). In the supplementary material, we explored a variety of alternate exogenous configurations to analyse further the accuracy of the Monge–Ampère predictions. These show how symmetry of the initial conditions reduces the median error between the Monge–Ampère approximation and the Eulerian solution. In general, the error between the Monge–Ampère approximation and the Eulerian solution reduces with larger initial endogenous surfactant concentration $\delta$, and with having fewer initial deposits. It is worth noting that when spreading is one-dimensional, the Monge–Ampère equation gives the exact solution. This is due to the fact that topologically, no rearrangement of the particles can be performed in one dimension without an energetic cost captured by the evolution equation. In contrast, in two dimensions, particles can rearrange through a divergence-free stress field (such as blowing by a Suminagashi artist after reaching equilibrium), which has no energetic cost for the evolution equation (2.1).

Furthermore, the Monge–Ampère formulation makes a connection between the theory of surfactant-driven transport and the theory of optimal transport. Monge–Ampère equations arise for optimal transport problems where the transport satisfies the so-called quadratic Monge–Kantorovich optimal transport problem (qMK). If surfactant were transported optimally according to qMK, then the map $\boldsymbol {X}$ would satisfy

(4.1)\begin{equation} \min_{\boldsymbol{X}:\varGamma_0\to \bar{\varGamma}} \int_{\varOmega}|\boldsymbol{X}-\boldsymbol{x}_0|^2\varGamma_0 \, \mathrm{d} A_{\boldsymbol{x}_0}, \end{equation}

where $\min _{\boldsymbol {X}:\varGamma _0\to \bar {\varGamma }}$ means that we choose the minimum over all possible maps $\boldsymbol {X}$ that transport the initial surfactant concentration profile $\varGamma _0$ to the final uniform profile $\bar {\varGamma }$ (equivalently, from all possible maps $\boldsymbol {X}$ that satisfy (2.27)). It is known that such maps $\boldsymbol {X}$, which satisfy (4.1), are the gradients of (convex) scalar potentials that satisfy the Monge–Ampère equation (Rockafellar 1970; Caffarelli & McCann 2010). Unfortunately, the surfactant problem (2.1) does not satisfy (4.1) precisely, as we can see by the discrepancies in figure 6, for example, which shows a different prediction for the final configuration compared to the Eulerian particle-tracking method. However, Otto (2001) showed that solutions to equations of porous medium type (a class of equations to which (2.1) belongs) are gradient flows on the function space $\mathcal {M}$ of possible solutions; this function space can be shown to satisfy the definition of a Riemannian manifold when distances on the manifold are measured by the Wasserstein-2 distance, defined variationally as

(4.2)\begin{equation} W_2(\varGamma_0,\bar{\varGamma}) = \sqrt{ \min_{\boldsymbol{X}:\varGamma_0\to \bar{\varGamma}} \int_{\varOmega}|\boldsymbol{X}-\boldsymbol{x}_0|^2\varGamma_0 \, \mathrm{d} A_{\boldsymbol{x}_0} }. \end{equation}

The Wasserstein-2 distance (4.2) is simply the square root of qMK in (4.1). Hence optimal maps that satisfy (4.1) must follow the shortest path between $\varGamma _0$ and $\bar {\varGamma }$ on $\mathcal {M}$ (a geodesic). The question of how closely the Monge–Ampère method approximates the correct steady-state solution for surfactant spreading is therefore equivalent to the question of how far the gradient flow deviates from the geodesic on $\mathcal {M}$. This leads to the possibility that a quantification of how well the Monge–Ampère method approximates the correct solution could be made with variational analysis. We do not pursue this further here, except to recall the degeneracy in the evolution equation revealed in Appendix A, namely that (2.1) does not account for the dissipation associated with any flow that preserves surface concentration; this raises the possibility that such flows may account for differences between equilibria predicted by the Monge–Ampère and Eulerian descriptions. Similarly, we leave as open the question of whether the solutions as $\delta \to 0$ can be used to infer the behaviour of the final state with an initially clean interface ($\delta =0$).

In the limit of vanishing endogenous surfactant concentration $\delta \to 0$, using the Monge–Ampère approximation, we find that exogenous deposits approach equilibrium shapes with polygonal boundaries (figure 9). Moreover, the internal corners between deposits present self-similar geometries, with a typical spatial extent scaling as $\delta ^{1/2}$ (figure 8). This self-similar behaviour can be explained by mass conservation. In the limit of $\delta \to 0$, we also observe that the topology at equilibrium is reduced to a small finite number of characteristic coordinates that are the intersections between the polygonal final shapes of the deposits. Predicting these coordinates a priori remains an open problem.

In figure 10 and supplementary movie 4, we explore how the Monge–Ampère method could be used to create a pattern resembling Suminagashi art. We have used divergence-free maps to mimic the blowing process used by Suminagashi artists. (The details for the creation of this picture are in Appendix E.) Such a computation, which involves 20 consecutive calculations of the Monge–Ampère equation, would be extremely expensive to run using a full dynamic solution. The choice of a divergence-free map found by time-stepping (E7) for the blowing step used by artists is credible, as gentle blowing on a surfactant-laden surface would likely deform the surface in such a way as to not create concentration gradients. Indeed, at equilibrium, the surfactant prevents the formation of concentration gradients (Manikantan & Squires 2020), because the Marangoni force opposes them (Appendix A). Moreover, we exploited the facts that gravitational forces in the Suminagashi process are sufficiently strong to suppress deformations of the gas–liquid interface, and that surface diffusion is sufficiently weak compared to advection to avoid smearing of the edge deposits. In other contexts, Marangoni flows can induce surface deformations. Despite the fact that the coupled evolution equations between the bulk and the surface for such problems can be written in a gradient flow formulation (Thiele et al. 2018), the link between such flows and optimal transport is not clear; the bulk flow that transports fluid is distinct from the surface flow that transports surfactant, such that distinct maps may be needed to characterise such problems. The numerical results from the Monge–Ampère method presented in figure 10 are only suggestive in their reproduction of Suminagashi art; nevertheless, we hope that they will encourage experimentalists to explore this link further.

Figure 10. Successive solutions of the Monge–Ampère equation are shown in (a–e), as detailed in Appendix E, showing a final pattern that is reminiscent of a Suminagashi pattern in ( f) (see also supplementary movie 4). The image was created by 20 solutions of the Monge–Ampère equation with a stirring or blowing step after the release of every four deposits. (a–c) Creation of the pattern after four depositions, after the first stirring step, and after $12$ depositions, respectively. (d–f) The solution just before the final stirring step after 20 depositions, the final result where a monochrome colour scheme is added, and a Suminagashi pattern made by Bea Mahan (2011) for qualitative comparison.

Many applications involving surfactants are concerned with the transport of passive solutes by the Marangoni effect in confined geometries, ranging from pharmaceutical delivery to the human lung to the creation of artistic patterns using Suminagashi techniques. We have presented methods for finding the dynamic behaviour of material surface particles, and for finding an approximation of the equilibrium location of material particles without having to solve for transient dynamics. The equilibrium approximation was achieved by showing a connection between surfactant dynamics and the theory of optimal transport, a research area at the forefront of modern mathematics, through a Monge–Ampère equation associated with the surfactant-driven transport. We hope that this connection, and the methods presented here, will spark the imaginations of researchers interested in both fundamental understanding and practical work related to surfactant-induced Marangoni flows carrying solutes.