coefficients $\sigma_{mn}$ are given by

$$ \sigma_{mn} = \begin{cases} \frac{1}{L_1 L_2} \int_{x=0}^{L_1} \int_{y=0}^{L_2} \Gamma(x, y, t_f) \, dy \, dx = \bar{\Gamma}, & \text{for } m, n = 0, \\ \frac{2}{L_1 L_2} \int_{x=0}^{L_1} \int_{y=0}^{L_2} \cos\left(\frac{n\pi y}{L_2}\right) \Gamma(x, y, t_f) \, dy \, dx, & \text{for } n > 0, m = 0, \\ \frac{2}{L_1 L_2} \int_{x=0}^{L_1} \int_{y=0}^{L_2} \cos\left(\frac{m\pi x}{L_1}\right) \Gamma(x, y, t_f) \, dy \, dx, & \text{for } m > 0, n = 0, \\ \frac{4}{L_1 L_2} \int_{x=0}^{L_1} \int_{y=0}^{L_2} \cos\left(\frac{m\pi x}{L_1}\right) \cos\left(\frac{n\pi y}{L_2}\right) \Gamma(x, y, t_f) \, dy \, dx, & \text{for } m, n > 0. \end{cases} \quad (B6) $$

It is relevant to our arguments for using the approximation (2.30) to point out that the late-time mapping approximation (B5) is exactly the gradient of the scalar potential

$$ \phi_{cr} = \left( \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \frac{\sigma_{mn}}{4\omega_{mn}} \cos(m\pi x/L_1) \cos(n\pi y/L_2) \right) - \frac{\sigma_{00}}{4\omega_{00}}. \quad (B7) $$

At a large time $t = t_f$, we expect the three leading order terms for $\Gamma(x, y, t_f)$ to be

$$ \Gamma(x, y, t_f) \approx \bar{\Gamma} + \sigma_{01} \cos\left(\frac{\pi y}{L_2}\right) + \sigma_{10} \cos\left(\frac{\pi x}{L_1}\right), \quad (B8) $$

as higher-order modes are subject to exponential decay over a smaller time scale through (B3). (Taking $t_f = 1047.8$ with $\delta = 0.25$, for example, we found that the next largest coefficient $\sigma_{mn}$ was a factor of $O(10^{-3})$ smaller than either $\sigma_{01}$ or $\sigma_{10}$, justifying (B8).) If we write $\Gamma_a$ for the concentration at $(0, L_2)$ and $\Gamma_b$ for the concentration at $(L_1, L_2)$, then

$$ \Gamma_a \approx \bar{\Gamma} + \sigma_{10} - \sigma_{01}, \quad \Gamma_b \approx \bar{\Gamma} - \sigma_{10} - \sigma_{01}, \quad (B9a,b) $$

which means

$$ \sigma_{10} \approx \frac{1}{2}(\Gamma_a - \Gamma_b), \quad \sigma_{01} \approx \frac{1}{2}(2\bar{\Gamma} - \Gamma_a - \Gamma_b). \quad (B10a,b) $$

The leading-order terms for the correction to the mapping $X_{cr}$ that maps particles from $t_f$ to $t \to \infty$ can therefore be approximated as

$$ X_{cr} \approx \left( \frac{\frac{L_1(\Gamma_a - \Gamma_b)}{2\pi\bar{\Gamma}} \sin(\pi x/L_1)}{\frac{L_2(2\bar{\Gamma} - \Gamma_a - \Gamma_b)}{2\pi\bar{\Gamma}} \sin(\pi y/L_2)} \right). \quad (B11) $$

We run a given Eulerian particle-tracking simulation to a time $t_f$ where the leading-order amplitude of the correction is below a pair of chosen tolerance levels $[X_{tol}, Y_{tol}]$, or in other words, until

$$ \left| \frac{L_1(\Gamma_a - \Gamma_b)}{2\pi \bar{\Gamma}} \right| < X_{tol}, \quad \left| \frac{L_2(2\bar{\Gamma} - \Gamma_a - \Gamma_b)}{2\pi \bar{\Gamma}} \right| < Y_{tol}. \quad (B12a,b) $$
## Appendix C. Behaviour of solutions to the nonlinear diffusion equation (2.1) near sharp corners

We wish to understand the behaviour of the two-dimensional nonlinear diffusion equation (2.1) subject to boundary conditions (2.2), near a sharp corner of a wedge-shaped domain. We introduce a polar coordinate system $(r, \theta)$ with the origin at the corner, with one boundary located at $\theta = 0$, and the other at $\theta = \Phi$. The surfactant concentration $\Gamma$ must neither diverge nor go to zero at the corner, so we look for expansions of the form

$$ \Gamma_{(m)} = A_{m,0}(\theta, t) + A_{m,1}(\theta, t) r^{a_{m,1}} + A_{m,2}(\theta, t) r^{a_{m,2}} + \dots, \quad (C1) $$

where $0 < a_{m,1} < a_{m,2} < \dots$, and $m$ is an index as we anticipate multiple expansions, each indexed by a different value of $m$. The result will be a sum of asymptotic series, a technique used to derive a corner expansion for surfactant-induced flow in Mcnair, Jensen & Landel (2022). Substituting an expansion for a given $m$ into (2.1), we obtain

$$ \begin{aligned} & 4 \left( \frac{\partial A_{m,0}}{\partial t} + r^{a_{m,1}} \frac{\partial A_{m,1}}{\partial t} + r^{a_{m,2}} \frac{\partial A_{m,2}}{\partial t} + \dots \right) \\ & = \frac{1}{2r^2} \frac{\partial^2 A_{m,0}^2}{\partial \theta^2} + a_{m,1}^2 A_{m,0} A_{m,1} r^{a_{m,1}-2} + r^{a_{m,1}-2} \frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{m,1}) \\ & \quad + \frac{1}{2} r^{2a_{m,1}-2} \frac{\partial^2 A_{m,1}^2}{\partial \theta^2} + A_{m,1}^2 (2a_{m,1})^2 r^{2a_{m,1}-2} + r^{a_{m,2}-2} \frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{m,2}) \\ & \quad + a_{m,2}^2 A_{m,0} A_{m,2} r^{a_{m,2}-2} + \dots. \end{aligned} \quad (C2) $$

As $a_{m,1} > 0$, the leading-order equation is $\partial^2 A_{m,0}^2 / \partial \theta^2 = 0$, which, when solved subject to boundary conditions $\partial A_{m,0} / \partial \theta = 0$ (from (2.2)), gives $A_{m,0} = A_{m,0}(t)$. The balance at the next order is

$$ 4 \frac{\partial A_{m,0}}{\partial t} = a_{m,1}^2 A_{m,0} A_{m,1} r^{a_{m,1}-2} + r^{a_{m,1}-2} \frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{m,1}). \quad (C3) $$

One possible exponent, which we index with $m = 0$, in the expansion is $a_{0,1} = 2$, yielding an expansion driven by $\partial A_{0,0} / \partial t$,

$$ \Gamma_{(0)} = A_{0,0}(t) + \frac{\partial A_{0,0}/\partial t}{A_{0,0}} r^2 + O(r^4), \quad (C4) $$

representing a purely radial flow that drives surfactant into or out of the corner, leading to changes in the corner concentration $A_{0,0}(t)$. The fact that the series goes in powers of $r^{2n}$ can be obtained by examining (C2) at the next order.

We now look for other possible expansions indexed by $m = 1, 2, 3, \dots$, and to avoid duplication of the primary flow contribution (C4), we set $A_{1,0} = A_{2,0} = \dots = 0$. (When $\Phi = \pi/2$, again assuming $a_{0,1} = 2$, (C3) also possesses a homogeneous solution $A_{0,1} = f_{0,1}(t) \cos(2\theta)$, representing a stagnation point flow in the corner, with strength $f_{0,1}(t)$ determined by conditions far from the corner.) More generally, (C3) possesses
## Appendix C. Behaviour of solutions to the nonlinear diffusion equation (2.1) near sharp corners

We wish to understand the behaviour of the two-dimensional nonlinear diffusion equation (2.1) subject to boundary conditions (2.2), near a sharp corner of a wedge-shaped domain. We introduce a polar coordinate system $(r, \theta)$ with the origin at the corner, with one boundary located at $\theta = 0$, and the other at $\theta = \Phi$. The surfactant concentration $\Gamma$ must neither diverge nor go to zero at the corner, so we look for expansions of the form

$$ \Gamma_{(m)} = A_{m,0}(\theta, t) + A_{m,1}(\theta, t) r^{a_{m,1}} + A_{m,2}(\theta, t) r^{a_{m,2}} + \dots, \quad (C1) $$

where $0 < a_{m,1} < a_{m,2} < \dots$, and $m$ is an index as we anticipate multiple expansions, each indexed by a different value of $m$. The result will be a sum of asymptotic series, a technique used to derive a corner expansion for surfactant-induced flow in Mcnair, Jensen & Landel (2022). Substituting an expansion for a given $m$ into (2.1), we obtain

$$ \begin{aligned} & 4 \left( \frac{\partial A_{m,0}}{\partial t} + r^{a_{m,1}} \frac{\partial A_{m,1}}{\partial t} + r^{a_{m,2}} \frac{\partial A_{m,2}}{\partial t} + \dots \right) \\ & = \frac{1}{2r^2} \frac{\partial^2 A_{m,0}^2}{\partial \theta^2} + a_{m,1}^2 A_{m,0} A_{m,1} r^{a_{m,1}-2} + r^{a_{m,1}-2} \frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{m,1}) \\ & \quad + \frac{1}{2} r^{2a_{m,1}-2} \frac{\partial^2 A_{m,1}^2}{\partial \theta^2} + A_{m,1}^2 (2a_{m,1})^2 r^{2a_{m,1}-2} + r^{a_{m,2}-2} \frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{m,2}) \\ & \quad + a_{m,2}^2 A_{m,0} A_{m,2} r^{a_{m,2}-2} + \dots. \end{aligned} \quad (C2) $$

As $a_{m,1} > 0$, the leading-order equation is $\partial^2 A_{m,0}^2 / \partial \theta^2 = 0$, which, when solved subject to boundary conditions $\partial A_{m,0} / \partial \theta = 0$ (from (2.2)), gives $A_{m,0} = A_{m,0}(t)$. The balance at the next order is

$$ 4 \frac{\partial A_{m,0}}{\partial t} = a_{m,1}^2 A_{m,0} A_{m,1} r^{a_{m,1}-2} + r^{a_{m,1}-2} \frac{\partial^2}{\partial \theta^2} (A_{m,0} A_{m,1}). \quad (C3) $$

One possible exponent, which we index with $m = 0$, in the expansion is $a_{0,1} = 2$, yielding an expansion driven by $\partial A_{0,0} / \partial t$,

$$ \Gamma_{(0)} = A_{0,0}(t) + \frac{\partial A_{0,0}/\partial t}{A_{0,0}} r^2 + O(r^4), \quad (C4) $$

representing a purely radial flow that drives surfactant into or out of the corner, leading to changes in the corner concentration $A_{0,0}(t)$. The fact that the series goes in powers of $r^{2n}$ can be obtained by examining (C2) at the next order.

We now look for other possible expansions indexed by $m = 1, 2, 3, \dots$, and to avoid duplication of the primary flow contribution (C4), we set $A_{1,0} = A_{2,0} = \dots = 0$. (When $\Phi = \pi/2$, again assuming $a_{0,1} = 2$, (C3) also possesses a homogeneous solution $A_{0,1} = f_{0,1}(t) \cos(2\theta)$, representing a stagnation point flow in the corner, with strength $f_{0,1}(t)$ determined by conditions far from the corner.) More generally, (C3) possesses
homogeneous solutions satisfying, for $m = 1, 2, \dots$,

$$a_{m,1}^2 A_{m,1} + \frac{\partial^2}{\partial \theta^2} (A_{m,1}) = 0. \qquad (C5)$$

This means that $A_{m,1} = f_{m,1}(t) \cos(a_{m,1}\theta)$, after applying $\partial A_{m,1}/\partial \theta = 0$ at $\theta = 0$. Applying $\partial A_{m,1}/\partial \theta = 0$ at $\theta = \Phi$ means that

$$a_{m,1} = \frac{m\pi}{\Phi}, \qquad (C6)$$

for any integer $m$, and this provides an infinite set of possible exponents. (When $\Phi = \pi/2$, (C6) with $m = 1$ recovers the stagnation point flow cited above with amplitude $f_{0,1}$.) By examining the next-order balance in (C2), it must also be the case that $a_{m,2} = m\pi/\Phi + 2$, and so on. Therefore, the full solution is given by the sum of asymptotic series

$$\Gamma = \Gamma(0) + \Gamma(1) + \Gamma(2) + \Gamma(3) + \dots, \qquad (C7)$$

where $\Gamma(0)$ is given by (C4), and

$$\Gamma(1) = f_{1,1}(t) \cos(\pi\theta/\Phi) r^{\pi/\Phi} + A_{1,2}(\theta, t) r^{\pi/\Phi+2} + A_{1,3}(\theta, t) r^{\pi/\Phi+4} + \dots, \qquad (C8)$$

$$\Gamma(2) = f_{2,1}(t) \cos(2\pi\theta/\Phi) r^{2\pi/\Phi} + A_{2,2}(\theta, t) r^{2\pi/\Phi+2} + A_{2,3}(\theta, t) r^{2\pi/\Phi+4} + \dots, \qquad (C9)$$

$$\Gamma(3) = f_{3,1}(t) \cos(3\pi\theta/\Phi) r^{3\pi/\Phi} + A_{3,2}(\theta, t) r^{3\pi/\Phi+2} + A_{3,3}(\theta, t) r^{3\pi/\Phi+4} + \dots, \qquad (C10)$$

and so on. The series can be summarised as

$$\Gamma(r, \theta, t) = \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} A_{m,n}(\theta, t) r^{m\pi/\Phi+2n}. \qquad (C11)$$

The velocity field is $u_s = -\nabla_x \Gamma/4$, so as $r \to 0$, it is a combination of a radial and a stagnation point flow:

$$u_s \approx \hat{r} \left( 2 \frac{\partial A_{0,0}}{\partial t} r + \frac{\pi}{4\Phi} f_{1,1}(t) \cos\left(\frac{\pi\theta}{\Phi}\right) r^{(\pi/\Phi)-1} \right) - \hat{\theta} f_{1,1}(t) \frac{\pi}{4\Phi} \sin\left(\frac{\pi\theta}{\Phi}\right) r^{(\pi/\Phi)-1}. \qquad (C12)$$

For internal angles in a convex domain (when $\Phi < \pi$), the velocity is proportional to $r$ to a positive power, so it goes to zero at the corner. In particular, when $\Phi = \pi/2$, a particle on the boundary $\theta = 0$ or $\theta = \pi/2$ has an inward radial velocity bounded above by $Fr$ for some finite $F > 0$. A particle starting at $r = r_0$ at $t = 0$ will then lie in $r > r_0 \exp(-Ft)$ for $t > 0$, never reaching the corner in finite time, supporting the use of (2.25a,b). For wedge angles $\Phi < \pi/2$, the radial flow is dominant as $r \to 0$, and (2.25a,b) again applies. For wedge angles satisfying $\Phi \in (\pi/2, \pi)$, the stagnation point flow dominates as $r \to 0$. In this case, although the velocity field along a boundary vanishes as $r \to 0$, a sustained inward flow $dr/dt = -Fr(\pi/\Phi)^{-1}$ for constant $F > 0$ has the potential to drive particles to the corner in finite time, with $r = [F(t_0 - t)(2 - (\pi/\Phi))]^{1/(2-(\pi/\Phi))}$ as $t \to t_0$ for some $t_0$, calling into question the validity of (2.25a,b) in this case. For even larger wedge angles ($\Phi > \pi$), regularity of the velocity field demands $f_{1,1} = 0$. This highlights the influence of any smoothing over a length scale of order $\delta \ll 1$ that might be applied to computations of flows around such corners in non-convex domains, where velocities of magnitude $1/\delta^{1-(\pi/\Phi)}$ can be anticipated.
homogeneous solutions satisfying, for $m = 1, 2, \dots$,

$$a_{m,1}^2 A_{m,1} + \frac{\partial^2}{\partial \theta^2} (A_{m,1}) = 0. \qquad (C5)$$

This means that $A_{m,1} = f_{m,1}(t) \cos(a_{m,1}\theta)$, after applying $\partial A_{m,1}/\partial\theta = 0$ at $\theta = 0$. Applying $\partial A_{m,1}/\partial\theta = 0$ at $\theta = \Phi$ means that

$$a_{m,1} = \frac{m\pi}{\Phi}, \qquad (C6)$$

for any integer $m$, and this provides an infinite set of possible exponents. (When $\Phi = \pi/2$, (C6) with $m = 1$ recovers the stagnation point flow cited above with amplitude $f_{0,1}$.) By examining the next-order balance in (C2), it must also be the case that $a_{m,2} = m\pi/\Phi + 2$, and so on. Therefore, the full solution is given by the sum of asymptotic series

$$\Gamma = \Gamma(0) + \Gamma(1) + \Gamma(2) + \Gamma(3) + \dots, \qquad (C7)$$

where $\Gamma(0)$ is given by (C4), and

$$\Gamma(1) = f_{1,1}(t) \cos(\pi\theta/\Phi) r^{\pi/\Phi} + A_{1,2}(\theta, t) r^{\pi/\Phi+2} + A_{1,3}(\theta, t) r^{\pi/\Phi+4} + \dots, \qquad (C8)$$

$$\Gamma(2) = f_{2,1}(t) \cos(2\pi\theta/\Phi) r^{2\pi/\Phi} + A_{2,2}(\theta, t) r^{2\pi/\Phi+2} + A_{2,3}(\theta, t) r^{2\pi/\Phi+4} + \dots, \qquad (C9)$$

$$\Gamma(3) = f_{3,1}(t) \cos(3\pi\theta/\Phi) r^{3\pi/\Phi} + A_{3,2}(\theta, t) r^{3\pi/\Phi+2} + A_{3,3}(\theta, t) r^{3\pi/\Phi+4} + \dots, \qquad (C10)$$

and so on. The series can be summarised as

$$\Gamma(r, \theta, t) = \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} A_{m,n}(\theta, t) r^{m\pi/\Phi+2n}. \qquad (C11)$$

The velocity field is $u_s = -\nabla_x \Gamma/4$, so as $r \to 0$, it is a combination of a radial and a stagnation point flow:

$$u_s \approx \hat{r} \left( 2 \frac{\partial A_{0,0}}{\partial t} r + \frac{\pi}{4\Phi} f_{1,1}(t) \cos\left(\frac{\pi\theta}{\Phi}\right) r^{(\pi/\Phi)-1} \right) - \hat{\theta} f_{1,1}(t) \frac{\pi}{4\Phi} \sin\left(\frac{\pi\theta}{\Phi}\right) r^{(\pi/\Phi)-1}. \qquad (C12)$$

For internal angles in a convex domain (when $\Phi < \pi$), the velocity is proportional to $r$ to a positive power, so it goes to zero at the corner. In particular, when $\Phi = \pi/2$, a particle on the boundary $\theta = 0$ or $\theta = \pi/2$ has an inward radial velocity bounded above by $Fr$ for some finite $F > 0$. A particle starting at $r = r_0$ at $t = 0$ will then lie in $r > r_0 \exp(-Ft)$ for $t > 0$, never reaching the corner in finite time, supporting the use of (2.25a,b). For wedge angles $\Phi < \pi/2$, the radial flow is dominant as $r \to 0$, and (2.25a,b) again applies. For wedge angles satisfying $\Phi \in (\pi/2, \pi)$, the stagnation point flow dominates as $r \to 0$. In this case, although the velocity field along a boundary vanishes as $r \to 0$, a sustained inward flow $dr/dt = -Fr(\pi/\Phi)^{-1}$ for constant $F > 0$ has the potential to drive particles to the corner in finite time, with $r = [F(t_0 - t)(2 - (\pi/\Phi))]^{1/(2-(\pi/\Phi))}$ as $t \to t_0$ for some $t_0$, calling into question the validity of (2.25a,b) in this case. For even larger wedge angles ($\Phi > \pi$), regularity of the velocity field demands $f_{1,1} = 0$. This highlights the influence of any smoothing over a length scale of order $\delta \ll 1$ that might be applied to computations of flows around such corners in non-convex domains, where velocities of magnitude $1/\delta^{1-(\pi/\Phi)}$ can be anticipated.
Appendix D. The numerical scheme to solve the Monge-Ampère equation

D.1. *Outline of the scheme*

We approximate (2.33) and (2.32) on an ($M + 2) \times (N + 2)$ rectangular grid. The points on the four boundaries we consider as known, as they can be expressed as functions of the interior points by the boundary conditions using one-sided second-order differences. This leaves us with an $M \times N$ grid of unknowns. We create an $MN \times 1$ vector of the values of the solution $\phi$ at these grid points by stacking the rows of the grid of unknowns into a column vector $\phi$, starting from the top and proceeding downwards. So, for example, $\phi_{i+N}$ refers to the grid point directly below $\phi_i$, and $\phi_{i+1}$ refers to the grid point directly to the right of $\phi_i$ unless $i = pN$ for some integer $p$.

We approximate (2.33) using second-order differences at the $i$th grid point as

$$f_i(\phi) = \frac{(\phi_{i+1} - 2\phi_i + \phi_{i-1})(\phi_{i+N} - 2\phi_i + \phi_{i-N})}{\Delta x^2 \Delta y^2} - \left( \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{4 \Delta x \Delta y} \right)^2 - G_i = 0, \quad (D1)$$

which we define to be the $i$th component of the vector $f(\phi)$. Here, $\Delta x$ and $\Delta y$ refer to the grid spacing in the $x_0$ and $y_0$ coordinate directions, respectively, and we define

$$G = 1 - \beta_j \left( 1 - \frac{\Gamma_0(x_0, y_0)}{\bar{\Gamma}} \right). \quad (D2)$$

The $k$th iteration of the Newton scheme for the vector $\phi$ of the values of the function $\phi$ approximated at the grid points is

$$\phi^k = \phi^{k-1} - (\nabla_{\phi} f(\phi^{k-1}))^{-1} f(\phi^{k-1}), \quad (D3)$$

where the Jacobian $\nabla_{\phi} f(\phi)$ is the sparse matrix of derivatives of the components $f_i$ with respect to solution at each grid point $\phi_j$, so that the element $\partial f_i / \partial \phi_j$ is zero, except for nine diagonals given by

$$\frac{\partial f_i}{\partial \phi_i} = \frac{2(4\phi_i - \phi_{i+1} - \phi_{i-1} - \phi_{i+N} - \phi_{i-N})}{\Delta x^2 \Delta y^2}, \quad (D4a)$$

$$\frac{\partial f_i}{\partial \phi_{i+1}} = \frac{\phi_{i+N} - 2\phi_i + \phi_{i-N}}{\Delta y^2 \Delta x^2}, \quad (D4b)$$

$$\frac{\partial f_i}{\partial \phi_{i-1}} = \frac{\phi_{i+N} - 2\phi_i + \phi_{i-N}}{\Delta y^2 \Delta x^2}, \quad (D4c)$$

$$\frac{\partial f_i}{\partial \phi_{i+N}} = \frac{\phi_{i+1} - 2\phi_i + \phi_{i-1}}{\Delta y^2 \Delta x^2}, \quad (D4d)$$

$$\frac{\partial f_i}{\partial \phi_{i-N}} = \frac{\phi_{i+N} - 2\phi_i + \phi_{i-N}}{\Delta y^2 \Delta x^2}, \quad (D4e)$$

$$\frac{\partial f_i}{\partial \phi_{i+N+1}} = \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}, \quad (D4f)$$
Appendix D. The numerical scheme to solve the Monge-Ampère equation

D.1. *Outline of the scheme*

We approximate (2.33) and (2.32) on an ($M + 2) \times (N + 2)$ rectangular grid. The points on the four boundaries we consider as known, as they can be expressed as functions of the interior points by the boundary conditions using one-sided second-order differences. This leaves us with an $M \times N$ grid of unknowns. We create an $MN \times 1$ vector of the values of the solution $\phi$ at these grid points by stacking the rows of the grid of unknowns into a column vector $\phi$, starting from the top and proceeding downwards. So, for example, $\phi_{i+N}$ refers to the grid point directly below $\phi_i$, and $\phi_{i+1}$ refers to the grid point directly to the right of $\phi_i$ unless $i = pN$ for some integer $p$.

We approximate (2.33) using second-order differences at the $i$th grid point as

$$f_i(\phi) = \frac{(\phi_{i+1} - 2\phi_i + \phi_{i-1})(\phi_{i+N} - 2\phi_i + \phi_{i-N})}{\Delta x^2 \Delta y^2} - \left( \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{4 \Delta x \Delta y} \right)^2 - G_i = 0, \quad (D1)$$

which we define to be the $i$th component of the vector $f(\phi)$. Here, $\Delta x$ and $\Delta y$ refer to the grid spacing in the $x_0$ and $y_0$ coordinate directions, respectively, and we define

$$G = 1 - \beta_j \left( 1 - \frac{\Gamma_0(x_0, y_0)}{\bar{\Gamma}} \right). \quad (D2)$$

The $k$th iteration of the Newton scheme for the vector $\phi$ of the values of the function $\phi$ approximated at the grid points is

$$\phi^k = \phi^{k-1} - (\nabla_{\phi} f(\phi^{k-1}))^{-1} f(\phi^{k-1}), \quad (D3)$$

where the Jacobian $\nabla_{\phi} f(\phi)$ is the sparse matrix of derivatives of the components $f_i$ with respect to solution at each grid point $\phi_j$, so that the element $\partial f_i / \partial \phi_j$ is zero, except for nine diagonals given by

$$\frac{\partial f_i}{\partial \phi_i} = \frac{2(4\phi_i - \phi_{i+1} - \phi_{i-1} - \phi_{i+N} - \phi_{i-N})}{\Delta x^2 \Delta y^2}, \quad (D4a)$$

$$\frac{\partial f_i}{\partial \phi_{i+1}} = \frac{\phi_{i+N} - 2\phi_i + \phi_{i-N}}{\Delta y^2 \Delta x^2}, \quad (D4b)$$

$$\frac{\partial f_i}{\partial \phi_{i-1}} = \frac{\phi_{i+N} - 2\phi_i + \phi_{i-N}}{\Delta y^2 \Delta x^2}, \quad (D4c)$$

$$\frac{\partial f_i}{\partial \phi_{i+N}} = \frac{\phi_{i+1} - 2\phi_i + \phi_{i-1}}{\Delta y^2 \Delta x^2}, \quad (D4d)$$

$$\frac{\partial f_i}{\partial \phi_{i-N}} = \frac{\phi_{i+N} - 2\phi_i + \phi_{i-N}}{\Delta y^2 \Delta x^2}, \quad (D4e)$$

$$\frac{\partial f_i}{\partial \phi_{i+N+1}} = \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}, \quad (D4f)$$
Lagrangian surfactant dynamics in confined domains

$$ \frac{\partial f_i}{\partial \phi_{i-N-1}} = \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}, \qquad (D4g) $$

$$ \frac{\partial f_i}{\partial \phi_{i+N-1}} = - \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}, \qquad (D4h) $$

$$ \frac{\partial f_i}{\partial \phi_{i-N+1}} = - \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}. \qquad (D4i) $$

The definitions (D1) and (D4) are valid for every interior point of the $M \times N$ grid of unknowns; however, for the values $i = 1$ to $N$, $i = (M - 1)N + 1$ to $MN$, $(p - 1)N + 1$ and $pn$ for every integer $1 < p < M$, the elements of $f$ and the Jacobian are slightly different to this as they incorporate the boundary conditions, but we omit them here for conciseness. We iterate (D3) until the sum of absolute differences of the components of $\phi$ from one iteration to the next falls below a small tolerance value, which we choose to be $10^{-6}$.

We use this numerical method to perform the calculation using the method presented in § 2.3.2 for a relatively coarse grid until $\beta_j = 1$, and then increase the refinement by interpolating the solution onto a more refined grid, and re-solving (2.33) and (2.32) with $\beta_j = 1$ (a multi-grid method). We repeat this second process until we have a sufficiently refined solution.

## D.2. Convergence of the numerical scheme

Figure 11 demonstrates the convergence of the numerical scheme presented in detail in § D.1 for the solution to the problem set out in § 2.3.2 with $\delta = 0.25, 0.05$ and $0.002$. For the purpose of illustration, we have picked a point $(x_0, y_0) = (5, 5)$ and displayed the error between an approximation for the correct solution and solutions for $\phi(5, 5)$, $X = \phi_{x_0}(5, 5)$ and $Y = \phi_{y_0}(5, 5)$ at various values for a uniformly spaced grid with $h = \Delta x = \Delta y$. The approximation for the correct solution is found by selecting $h \equiv h_{end}$ smaller than all of the other values of $h$ used, and using this as the grid spacing to solve for $\phi(5, 5) \equiv \phi_{end}(5, 5)$, $X(5, 5) = \phi_{x_0}(5, 5) \equiv X_{end}(5, 5)$ and $Y(5, 5) = \phi_{y_0}(5, 5) \equiv Y_{end}(5, 5)$, for each value of $\delta$.

The solutions for all three values in figure 11 converge approximately with slope 2 on the log scale indicated by a dotted black line, which is expected for our second-order finite-difference scheme. The graphs are only approximately showing this convergence, as there are numerous sources of noise, e.g. the approximation that we have used for the correct solution. Also, the code stops iterating once the sum of the absolute difference of each grid point in a solution between one iteration and the next falls below a certain small tolerance value. The exact locations of the edges of the deposits do not coincide precisely with grid points, and this is a further source of error as the initial conditions (and therefore the right-hand side of (2.31)) for differently discretised grids are effectively slightly different. The solutions for $X$ and $Y$ also involve calculating a further derivative using a second-order finite-difference approximation, introducing further error. However, the dominant error is clearly the square of the discretisation parameter $h$ used in the approximation (2.31). Solutions with $\delta = 0.002$ are significantly more expensive computationally, which is why the curves terminate at larger values of $h$ than the other curves. The fact that we have used a larger value of $h$ to calculate $\phi_{end}(5, 5)$ and its derivatives for $\delta = 0.002$ might also partly explain why some of these curves are noisier than for larger values of $\delta$.
Lagrangian surfactant dynamics in confined domains

$$ \frac{\partial f_i}{\partial \phi_{i-N-1}} = \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}, \quad (D4g) $$

$$ \frac{\partial f_i}{\partial \phi_{i+N-1}} = - \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}, \quad (D4h) $$

$$ \frac{\partial f_i}{\partial \phi_{i-N+1}} = - \frac{-\phi_{i+N+1} - \phi_{i-N-1} + \phi_{i+N-1} + \phi_{i-N+1}}{8 \Delta x^2 \Delta y^2}. \quad (D4i) $$

The definitions (D1) and (D4) are valid for every interior point of the $M \times N$ grid of unknowns; however, for the values $i = 1$ to $N$, $i = (M - 1)N + 1$ to $MN$, $(p - 1)N + 1$ and $pn$ for every integer $1 < p < M$, the elements of $f$ and the Jacobian are slightly different to this as they incorporate the boundary conditions, but we omit them here for conciseness. We iterate (D3) until the sum of absolute differences of the components of $\phi$ from one iteration to the next falls below a small tolerance value, which we choose to be $10^{-6}$.

We use this numerical method to perform the calculation using the method presented in § 2.3.2 for a relatively coarse grid until $\beta_j = 1$, and then increase the refinement by interpolating the solution onto a more refined grid, and re-solving (2.33) and (2.32) with $\beta_j = 1$ (a multi-grid method). We repeat this second process until we have a sufficiently refined solution.

## D.2. Convergence of the numerical scheme

Figure 11 demonstrates the convergence of the numerical scheme presented in detail in § D.1 for the solution to the problem set out in § 2.3.2 with $\delta = 0.25, 0.05$ and $0.002$. For the purpose of illustration, we have picked a point $(x_0, y_0) = (5, 5)$ and displayed the error between an approximation for the correct solution and solutions for $\phi(5, 5)$, $X = \phi_{x_0}(5, 5)$ and $Y = \phi_{y_0}(5, 5)$ at various values for a uniformly spaced grid with $h = \Delta x = \Delta y$. The approximation for the correct solution is found by selecting $h \equiv h_{end}$ smaller than all of the other values of $h$ used, and using this as the grid spacing to solve for $\phi(5, 5) \equiv \phi_{end}(5, 5)$, $X(5, 5) = \phi_{x_0}(5, 5) \equiv X_{end}(5, 5)$ and $Y(5, 5) = \phi_{y_0}(5, 5) \equiv Y_{end}(5, 5)$, for each value of $\delta$.

The solutions for all three values in figure 11 converge approximately with slope 2 on the log scale indicated by a dotted black line, which is expected for our second-order finite-difference scheme. The graphs are only approximately showing this convergence, as there are numerous sources of noise, e.g. the approximation that we have used for the correct solution. Also, the code stops iterating once the sum of the absolute difference of each grid point in a solution between one iteration and the next falls below a certain small tolerance value. The exact locations of the edges of the deposits do not coincide precisely with grid points, and this is a further source of error as the initial conditions (and therefore the right-hand side of (2.31)) for differently discretised grids are effectively slightly different. The solutions for $X$ and $Y$ also involve calculating a further derivative using a second-order finite-difference approximation, introducing further error. However, the dominant error is clearly the square of the discretisation parameter $h$ used in the approximation (2.31). Solutions with $\delta = 0.002$ are significantly more expensive computationally, which is why the curves terminate at larger values of $h$ than the other curves. The fact that we have used a larger value of $h$ to calculate $\phi_{end}(5, 5)$ and its derivatives for $\delta = 0.002$ might also partly explain why some of these curves are noisier than for larger values of $\delta$.
R. Mcnair, O.E. Jensen and J.R. Landel

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

Figure 11. Graph showing the convergence of the finite-difference approximation to the solution of the Monge–Ampère equation (2.31). To illustrate convergence, we choose a point on the Lagrangian domain $(x_0, y_0) = (5, 5)$ and calculate the solution for the three-deposit problem outlined in § 2.3.2 with $\delta = 0.25$, 0.05 and 0.002 (each $\delta$ is assigned a different line style – see legend) for multiple values of a uniform grid spacing $h$ (data points are shown as circles). A set of values for each $\delta$, which we call $\phi_{end}(5, 5)$, $X_{end}(5, 5)$ and $Y_{end}(5, 5)$, is calculated for a discretisation parameter $h_{end}$ that is smaller than the rest of the values of $h$ used. We use this solution as our approximation to the correct solution. The curves show the difference between $X(5, 5)$ (blue), $Y(5, 5)$ (orange) and $\phi(5, 5)$ (purple), and the approximated correct solution, showing convergence at a rate of order $h^2$ (dotted black line).

## Appendix E. Simulating Suminagashi patterns

We create a Suminagashi pattern using successive solutions of the Monge–Ampère equation. As described in the Introduction (§ 1), the creation of Suminagashi art patterns consists of repeated deposition of an ink–surfactant mixture onto a liquid surface, followed by the artist blowing on the surface between ink depositions. We simulate this procedure using the Monge–Ampère equation (2.31) to simulate the final locations of ink depositions in the spreading steps, and a divergence-free map to create the effects of the blowing steps. We show in (A11) how the rotational component of a surface shear stress can redistribute surfactant without changing its surface concentration.

The physical assumptions that we consider in this study for the Monge–Ampère method are appropriate for Suminagashi. Exploiting conditions (A8) and (A9): for a typical 1 cm layer of water with a surface tension reduction of 1–10 g s⁻², we find that the Bond number is large and of the order of 10²–10³; surface diffusivity for most surfactants is of the order of $D^* = 10^{-10}–10^{-9} \text{ m}^2 \text{ s}^{-1}$ (Chang & Franses 1995), which leads to a very large surface Péclet number of the order of 10⁷–10⁹; furthermore, Suminagashi patterns can be transferred to paper within a minute, during which any boundary thickens under molecular (surface) diffusivity by less than 0.3 mm.

### E.1. Spreading steps

Each spreading step starts from initial conditions (2.7) with $F(x_0, y_0) = C_q(x_1, y_1, 0.5, 1 – \delta)$, and by imposing $L_1 = 13, L_2 = 11$, so that this is now a single deposit of exogenous
R. Mcnair, O.E. Jensen and J.R. Landel

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

Figure 11. Graph showing the convergence of the finite-difference approximation to the solution of the Monge–Ampère equation (2.31). To illustrate convergence, we choose a point on the Lagrangian domain $(x_0, y_0) = (5, 5)$ and calculate the solution for the three-deposit problem outlined in § 2.3.2 with $\delta = 0.25$, 0.05 and 0.002 (each $\delta$ is assigned a different line style – see legend) for multiple values of a uniform grid spacing $h$ (data points are shown as circles). A set of values for each $\delta$, which we call $\phi_{end}(5, 5)$, $X_{end}(5, 5)$ and $Y_{end}(5, 5)$, is calculated for a discretisation parameter $h_{end}$ that is smaller than the rest of the values of $h$ used. We use this solution as our approximation to the correct solution. The curves show the difference between $X(5, 5)$ (blue), $Y(5, 5)$ (orange) and $\phi(5, 5)$ (purple), and the approximated correct solution, showing convergence at a rate of order $h^2$ (dotted black line).

## Appendix E. Simulating Suminagashi patterns

We create a Suminagashi pattern using successive solutions of the Monge–Ampère equation. As described in the Introduction (§ 1), the creation of Suminagashi art patterns consists of repeated deposition of an ink–surfactant mixture onto a liquid surface, followed by the artist blowing on the surface between ink depositions. We simulate this procedure using the Monge–Ampère equation (2.31) to simulate the final locations of ink depositions in the spreading steps, and a divergence-free map to create the effects of the blowing steps. We show in (A11) how the rotational component of a surface shear stress can redistribute surfactant without changing its surface concentration.

The physical assumptions that we consider in this study for the Monge–Ampère method are appropriate for Suminagashi. Exploiting conditions (A8) and (A9): for a typical 1 cm layer of water with a surface tension reduction of 1–10 g s⁻², we find that the Bond number is large and of the order of 10²–10³; surface diffusivity for most surfactants is of the order of $D^* = 10^{-10}–10^{-9} \text{ m}^2 \text{ s}^{-1}$ (Chang & Franses 1995), which leads to a very large surface Péclet number of the order of 10⁷–10⁹; furthermore, Suminagashi patterns can be transferred to paper within a minute, during which any boundary thickens under molecular (surface) diffusivity by less than 0.3 mm.

### E.1. Spreading steps

Each spreading step starts from initial conditions (2.7) with $F(x_0, y_0) = C_q(x_1, y_1, 0.5, 1 – \delta)$, and by imposing $L_1 = 13, L_2 = 11$, so that this is now a single deposit of exogenous
surfactant spreading into endogenous surfactant. For the first spreading step, we set $\delta = 0.01$, and solve (2.31) and (2.32). We store the final location of the boundary of this initial deposit, which we call $C_{1,(1)}$, where the first index labels each deposit edge, and the second index refers to how many total deposits have been released. Next, we set up a new problem by replacing $\delta$ by $\bar{\Gamma}$ from the previous problem, and we shift $(x_3, y_3)$ to a displaced position (using an algorithm described below as (E4)), and solve again (2.31) and (2.32), finding the edge of the new deposit, and the new locations of the edges of previously released deposits, repeating the above process $J$ times. This process is summarised below.

## E.2. Summary of the algorithm for the spreading steps

For $j = 1$ to $J$: solve (2.31) and (2.32) subject to initial conditions (2.7) with $\mathcal{F}(x_0, y_0) = C_q(x_1, y_1, 0.5, 1 - \delta)$, using the solution method outlined in § 2.3.2 with $L_1 = 13, L_2 = 11, \bar{\Gamma}_2 = \bar{\Gamma}_3 = \delta, r_1 = 0.5, r_2 = 0, r_3 = 0$, with $(x_1, y_1) = (x_{1,j}, y_{1,j})$ and $\delta = \delta_j$, where this last quantity is found using

$$ \delta_j = \bar{\Gamma}_{j-1}, \quad \bar{\Gamma}_0 = 0.05, \qquad (E1) $$

and where at each step $\bar{\Gamma}_j$ is computed using (2.4). We find all curves at step $j$ by

$$ C_{i,(j+1)} = X_j(C_{i,(j)}) \quad \text{for } i = 1, 2, \dots, j, \qquad (E2) $$

where $X_j$ is the solution of (2.31) and (2.32) at step $j$, and one new evolving curve is introduced at each step $j$ by

$$ C_{j,(j)} = \{(x, y) \mid (x - x_{1,j})^2 + (y - y_{1,j})^2 = 0.5\}. \qquad (E3) $$

The centres of the new deposits are chosen such that

$$ (x_{1,j+1}, y_{1,j+1}) = (x_{c,j} - 0.75, y_{c,j} - 0.2), \quad \text{where } (x_{c,j}, y_{c,j}) = \max_x (C_{j,(j)}), \qquad (E4) $$

until $\max_x(C_{j,(j)}) > 12.9$, and then

$$ (x_{1,j+1}, y_{1,j+1}) = (x_{c,j} - 0.05, y_{c,j} + 0.75), \quad \text{where } (x_{c,j}, y_{c,j}) = \min_y (C_{j,(j)}), \qquad (E5) $$

with

$$ (x_{1,1}, y_{1,1}) = (3.5, 8.5). \qquad (E6) $$

We vary the locations of the centres of the new deposits, as we have observed Suminagashi artists to do this. The final picture is created by plotting $C_{i,(J+1)}$ for all $i = 1$ to $J$.

## E.3. Blowing steps

Between some of the spreading steps, we introduce blowing steps, which are computed by finding a divergence-free map $X_b$ by time-stepping

$$ \frac{d}{dt}X_b = \sum_{n_s=1}^{\infty} A_{n_s} \left( \begin{array}{l} \frac{n_s \pi}{11} \sin \left( \frac{n_s \pi x}{13} \right) \cos \left( \frac{n_s \pi y}{11} \right) \\ -\frac{n_s \pi}{13} \cos \left( \frac{n_s \pi x}{13} \right) \sin \left( \frac{n_s \pi y}{11} \right) \end{array} \right), \qquad (E7) $$

for some integer $n_s$, between $t = 0$ and some $t_s$, where $A_{n_s}$ are the amplitudes of the modes, and where $X_b = [x_0, y_0]$ at $t = 0$. The resulting map $X_b$ is divergence-free (such that concentrations do not change) and satisfies the required boundary conditions. The functions on the right-hand side of (E7) form a basis, such that any divergence-free map can be obtained by choosing a set of amplitudes $A_{n_s}$ and some choice of $t_s$. The new location of deposit edge $i$, after spreading step $j$, is given by $X_b(C_{i,(j)})$.
Lagrangian surfactant dynamics in confined domains

surfactant spreading into endogenous surfactant. For the first spreading step, we set $\delta = 0.01$, and solve (2.31) and (2.32). We store the final location of the boundary of this initial deposit, which we call $C_{1,(1)}$, where the first index labels each deposit edge, and the second index refers to how many total deposits have been released. Next, we set up a new problem by replacing $\delta$ by $\bar{\Gamma}$ from the previous problem, and we shift $(x_3, y_3)$ to a displaced position (using an algorithm described below as (E4)), and solve again (2.31) and (2.32), finding the edge of the new deposit, and the new locations of the edges of previously released deposits, repeating the above process $J$ times. This process is summarised below.

## E.2. Summary of the algorithm for the spreading steps

For $j = 1$ to $J$: solve (2.31) and (2.32) subject to initial conditions (2.7) with $\mathcal{F}(x_0, y_0) = C_q(x_1, y_1, 0.5, 1 - \delta)$, using the solution method outlined in § 2.3.2 with $L_1 = 13, L_2 = 11, \bar{\Gamma}_2 = \bar{\Gamma}_3 = \delta, r_1 = 0.5, r_2 = 0, r_3 = 0$, with $(x_1, y_1) = (x_{1,j}, y_{1,j})$ and $\delta = \delta_j$, where this last quantity is found using

$$ \delta_j = \bar{\Gamma}_{j-1}, \quad \bar{\Gamma}_0 = 0.05, \qquad (E1) $$

and where at each step $\bar{\Gamma}_j$ is computed using (2.4). We find all curves at step $j$ by

$$ C_{i,(j+1)} = X_j(C_{i,(j)}) \quad \text{for } i = 1, 2, \dots, j, \qquad (E2) $$

where $X_j$ is the solution of (2.31) and (2.32) at step $j$, and one new evolving curve is introduced at each step $j$ by

$$ C_{j,(j)} = \{(x, y) | (x - x_{1,j})^2 + (y - y_{1,j})^2 = 0.5\}. \qquad (E3) $$

The centres of the new deposits are chosen such that

$$ (x_{1,j+1}, y_{1,j+1}) = (x_{c,j} - 0.75, y_{c,j} - 0.2), \quad \text{where } (x_{c,j}, y_{c,j}) = \max_x (C_{j,(j)}), \qquad (E4) $$

until $\max_x(C_{j,(j)}) > 12.9$, and then

$$ (x_{1,j+1}, y_{1,j+1}) = (x_{c,j} - 0.05, y_{c,j} + 0.75), \quad \text{where } (x_{c,j}, y_{c,j}) = \min_y (C_{j,(j)}), \qquad (E5) $$

with

$$ (x_{1,1}, y_{1,1}) = (3.5, 8.5). \qquad (E6) $$

We vary the locations of the centres of the new deposits, as we have observed Suminagashi artists to do this. The final picture is created by plotting $C_{i,(J+1)}$ for all $i = 1$ to $J$.

## E.3. Blowing steps

Between some of the spreading steps, we introduce blowing steps, which are computed by finding a divergence-free map $X_b$ by time-stepping

$$ \frac{d}{dt}X_b = \sum_{n_s=1}^{\infty} A_{n_s} \left( \begin{array}{l} \frac{n_s \pi}{11} \sin \left( \frac{n_s \pi x}{13} \right) \cos \left( \frac{n_s \pi y}{11} \right) \\ -\frac{n_s \pi}{13} \cos \left( \frac{n_s \pi x}{13} \right) \sin \left( \frac{n_s \pi y}{11} \right) \end{array} \right), \qquad (E7) $$

for some integer $n_s$, between $t = 0$ and some $t_s$, where $A_{n_s}$ are the amplitudes of the modes, and where $X_b = [x_0, y_0]$ at $t = 0$. The resulting map $X_b$ is divergence-free (such that concentrations do not change) and satisfies the required boundary conditions. The functions on the right-hand side of (E7) form a basis, such that any divergence-free map can be obtained by choosing a set of amplitudes $A_{n_s}$ and some choice of $t_s$. The new location of deposit edge $i$, after spreading step $j$, is given by $X_b(C_{i,(j)})$.
<figure><img src="image_22.png" /></figure>

Figure 12. Graphs showing the normalised absolute error (2.34) for several different values of $\delta$. The initial conditions are (2.7) and (2.9), with $C_q$ replaced by $C_c$ defined in (F1) except for the first solution, which is taken from Eul[0.002] and MA[0.002] and has a quadratic initial profile. (a) The logarithm of the Euclidean distance between the predictions of final location, for each initial particle location from a 221 × 261 grid, represented as a data point. The whiskers show 1.5 times the interquartile range above and below the quartiles, and particles outside this range are considered outliers and plotted as a cloud of points in blue, with the 90th percentile plotted as an orange line. (b) The same data as in (a), plotted in the form of a cumulative distribution function (CDF), with $\delta$ increasing in the direction of the arrow. Horizontal dashed lines indicate the median, the 75th percentile and the 90th percentile, respectively.

### E.4. Sequential depositions: the Suminagashi patterns

To illustrate the utility of the Monge–Ampère approximation, a solution of the Suminagashi algorithm is presented in figure 10. Here, for (E7), we choose $A_{ns} = 0$ for every $n_s$ except $N_s = 4$, and we choose $A_4 = 1$. We run 20 spreading steps, with a blowing step imposed after every four spreading steps with $t_s = 0.15$. The creation of figure 10(f) is shown after four spreading steps, after four spreading steps and one blowing step, after 12 spreading steps, after 12 spreading steps and three blowing steps, and finally after 20 spreading steps before the final blowing step. The final result is given in a monochrome colour scheme in figure 10(e).

## Appendix F. Difference in the prediction between Eulerian particle tracking and Monge–Ampère for different values of $\delta$

Figure 12 shows a series of box and whisker plots of the normalised absolute error (2.34), measuring the Euclidean distance between predictions for final particle location between Monge–Ampère and Eulerian particle-tracking, for multiple values of $\delta$. Except for the solution with $\delta = 0.002$, which has quadratic initial conditions (2.7) and (2.9), the rest of the solutions have cosine-shaped exogenous deposit profiles such that

$$C_c(\mathbf{x}_0; \mathbf{x}_c, r, \Gamma_{0,c} - \delta) = \begin{cases} \left(\frac{\Gamma_{0,c} - \delta}{2}\right) \left(1 + \cos\left(\frac{\pi (|\mathbf{x}_0 - \mathbf{x}_c|^2)}{r^2}\right)\right), & |\mathbf{x}_0 - \mathbf{x}_c| \le r, \\ 0, & |\mathbf{x}_0 - \mathbf{x}_c| > r, \end{cases} \quad (\text{F1})$$

replaces $C_q(\mathbf{x}_0; \mathbf{x}_c, r, \Gamma_{0,c} - \delta)$ in (2.7) and (2.9). The shape of the exogenous deposit profiles does not have a significant effect on results as shown in figure S7 of § S5 of the supplementary material, but we change the profile shape here for the sake of variety.
<figure><img src="image_28.png" /></figure>

Figure 12. Graphs showing the normalised absolute error (2.34) for several different values of $\delta$. The initial conditions are (2.7) and (2.9), with $C_q$ replaced by $C_c$ defined in (F1) except for the first solution, which is taken from Eul[0.002] and MA[0.002] and has a quadratic initial profile. (a) The logarithm of the Euclidean distance between the predictions of final location, for each initial particle location from a 221 × 261 grid, represented as a data point. The whiskers show 1.5 times the interquartile range above and below the quartiles, and particles outside this range are considered outliers and plotted as a cloud of points in blue, with the 90th percentile plotted as an orange line. (b) The same data as in (a), plotted in the form of a cumulative distribution function (CDF), with $\delta$ increasing in the direction of the arrow. Horizontal dashed lines indicate the median, the 75th percentile and the 90th percentile, respectively.

### E.4. Sequential depositions: the Suminagashi patterns

To illustrate the utility of the Monge–Ampère approximation, a solution of the Suminagashi algorithm is presented in figure 10. Here, for (E7), we choose $A_{ns} = 0$ for every $n_s$ except $N_s = 4$, and we choose $A_4 = 1$. We run 20 spreading steps, with a blowing step imposed after every four spreading steps with $t_s = 0.15$. The creation of figure 10(f) is shown after four spreading steps, after four spreading steps and one blowing step, after 12 spreading steps, after 12 spreading steps and three blowing steps, and finally after 20 spreading steps before the final blowing step. The final result is given in a monochrome colour scheme in figure 10(e).

## Appendix F. Difference in the prediction between Eulerian particle tracking and Monge–Ampère for different values of $\delta$

Figure 12 shows a series of box and whisker plots of the normalised absolute error (2.34), measuring the Euclidean distance between predictions for final particle location between Monge–Ampère and Eulerian particle-tracking, for multiple values of $\delta$. Except for the solution with $\delta = 0.002$, which has quadratic initial conditions (2.7) and (2.9), the rest of the solutions have cosine-shaped exogenous deposit profiles such that

$$C_c(\mathbf{x}_0; \mathbf{x}_c, r, \Gamma_{0,c} - \delta) = \begin{cases} \left(\frac{\Gamma_{0,c} - \delta}{2}\right) \left(1 + \cos\left(\frac{\pi (|\mathbf{x}_0 - \mathbf{x}_c|^2)}{r^2}\right)\right), & |\mathbf{x}_0 - \mathbf{x}_c| \le r, \\ 0, & |\mathbf{x}_0 - \mathbf{x}_c| > r, \end{cases} \quad (F1)$$

replaces $C_q(\mathbf{x}_0; \mathbf{x}_c, r, \Gamma_{0,c} - \delta)$ in (2.7) and (2.9). The shape of the exogenous deposit profiles does not have a significant effect on results as shown in figure S7 of § S5 of the supplementary material, but we change the profile shape here for the sake of variety.
Lagrangian surfactant dynamics in confined domains

For moderate $\delta$, the error between different methods remains small, but it increases rapidly as $\delta$ tends to zero. However, for $\delta$ as small as 0.075, the normalised absolute error is below 0.01 for 75% of particles, and even for $\delta = 0.002$, the error is below 0.05 for 90% of particles.

REFERENCES

AVERY, M.E. & MEAD, J. 1959 Surface properties in relation to atelectasis and hyaline membrane disease. *Am. J. Dis. Child* **97**, 517–523.

BENAMOU, J.-D. & BRENIER, Y. 2000 A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem. *Numer. Math.* **84** (3), 375–393.

BENAMOU, J.-D., FROESE, B.D. & OBERMAN, A.M. 2012 A viscosity solution approach to the Monge–Ampère formulation of the optimal transportation problem. arXiv:1208.4873.

BENAMOU, J.-D., FROESE, B.D. & OBERMAN, A.M. 2014 Numerical solution of the optimal transportation problem using the Monge–Ampère equation. *J. Comput. Phys.* **260**, 107–126.

BORGAS, M.S. & GROTBERG, J.B. 1988 Monolayer flow on a thin film. *J. Fluid Mech.* **193**, 151–170.

CAFFARELLI, L.A. & MCCANN, R.J. 2010 Free boundaries in optimal transport and Monge–Ampère obstacle problems. *Ann. Maths* **171**, 673–730.

CARRILLO, J.A., MATTHES, D. & WOLFRAM, M.-T. 2021 Lagrangian schemes for Wasserstein gradient flows. In *Handbook of Numerical Analysis* (ed. A. Bonito & R.H. Nochetto), vol. 22, chap. 4, 271–311. Elsevier.

CHAMBERS, A. 1991 *Suminagashi: The Japanese Art of Marbling: A Practical Guide*. Thames & Hudson.

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

DENG, Y., ZHENG, X., BAI, Y., WANG, Q., ZHAO, J. & HUANG, J. 2018 Surfactant-controlled ink drying enables high-speed deposition of perovskite films for efficient photovoltaic modules. *Nat. Energy* **3** (7), 560–566.

ESPINOSA, F.F., SHAPIRO, A.H., FREDBERG, J.J. & KAMM, R.D. 1993 Spreading of exogenous surfactant in an airway. *J. Appl. Physiol.* **75** (5), 2028–2039.

FROESE, B.D. & OBERMAN, A.M. 2011 Fast finite difference solvers for singular solutions of the elliptic Monge–Ampère equation. *J. Comput. Phys.* **230** (3), 818–834.

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

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

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

HAITSMA, J.J., LACHMANN, U. & LACHMANN, B. 2001 Exogenous surfactant as a drug delivery agent. *Adv. Drug Deliv. Rev.* **47** (2–3), 197–207.

HALLIDAY, H.L. 2008 Surfactants: past, present and future. *J. Perinatol.* **28** (1), S47–S56.

HALPERN, D., JENSEN, O.E. & GROTBERG, J.B. 1998 A theoretical study of surfactant and liquid delivery into the lung. *J. Appl. Physiol.* **85** (1), 333–352.

HENKEL, C., SNOEIJER, J.H. & THIELE, U. 2021 Gradient-dynamics model for liquid drops on elastic substrates. *Soft Matt.* **17** (45), 10359–10375.

HIDALGO, A., CRUZ, A. & PÉREZ-GIL, J. 2015 Barrier or carrier? Pulmonary surfactant and drug delivery. *Eur. J. Pharm. Biopharm.* **95**, 117–127.

IASELLA, S., SHARMA, R., GAROFF, S. & TILTON, R.D. 2024 Interaction of impinging Marangoni fields. *J. Colloid Interface Sci.* **653**, 807–820.

ISHII, T. & MURO, J. 1989 The fractal aspects of Suminagashi and marbling. *Thin Solid Films* **178**, 109–113.

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

JENSEN, O.E., HALPERN, D. & GROTBERG, J.B. 1994 Transport of a passive solute by surfactant-driven flows. *Chem. Engng Sci.* **49** (8), 1107–1117.

JOBE, A.H. 1993 Pulmonary surfactant therapy. *New Engl. J. Med.* **328** (12), 861–868.

JORDAN, R., KINDERLEHRER, D. & OTTO, F. 1998 The variational formulation of the Fokker–Planck equation. *SIAM J. Math. Anal.* **29** (1), 1–17.

KANTOROVICH, L.V. 1942 On the translocation of masses. *Dokl. Akad. Nauk SSSR* **37** (7–8), 227–229.

KANTOROVICH, L.V. 2006 On the translocation of masses. *J. Math. Sci.* **133** (4), 1381–1382.
Lagrangian surfactant dynamics in confined domains

For moderate $\delta$, the error between different methods remains small, but it increases rapidly as $\delta$ tends to zero. However, for $\delta$ as small as 0.075, the normalised absolute error is below 0.01 for 75% of particles, and even for $\delta = 0.002$, the error is below 0.05 for 90% of particles.

REFERENCES

AVERY, M.E. & MEAD, J. 1959 Surface properties in relation to atelectasis and hyaline membrane disease. *Am. J. Dis. Child* **97**, 517–523.

BENAMOU, J.-D. & BRENIER, Y. 2000 A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem. *Numer. Math.* **84** (3), 375–393.

BENAMOU, J.-D., FROESE, B.D. & OBERMAN, A.M. 2012 A viscosity solution approach to the Monge–Ampère formulation of the optimal transportation problem. arXiv:1208.4873.

BENAMOU, J.-D., FROESE, B.D. & OBERMAN, A.M. 2014 Numerical solution of the optimal transportation problem using the Monge–Ampère equation. *J. Comput. Phys.* **260**, 107–126.

BORGAS, M.S. & GROTBERG, J.B. 1988 Monolayer flow on a thin film. *J. Fluid Mech.* **193**, 151–170.

CAFFARELLI, L.A. & MCCANN, R.J. 2010 Free boundaries in optimal transport and Monge–Ampère obstacle problems. *Ann. Maths* **171**, 673–730.

CARRILLO, J.A., MATTHES, D. & WOLFRAM, M.-T. 2021 Lagrangian schemes for Wasserstein gradient flows. In *Handbook of Numerical Analysis* (ed. A. Bonito & R.H. Nochetto), vol. 22, chap. 4, 271–311. Elsevier.

CHAMBERS, A. 1991 *Suminagashi: The Japanese Art of Marbling: A Practical Guide*. Thames & Hudson.

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

DENG, Y., ZHENG, X., BAI, Y., WANG, Q., ZHAO, J. & HUANG, J. 2018 Surfactant-controlled ink drying enables high-speed deposition of perovskite films for efficient photovoltaic modules. *Nat. Energy* **3** (7), 560–566.

ESPINOSA, F.F., SHAPIRO, A.H., FREDBERG, J.J. & KAMM, R.D. 1993 Spreading of exogenous surfactant in an airway. *J. Appl. Physiol.* **75** (5), 2028–2039.

FROESE, B.D. & OBERMAN, A.M. 2011 Fast finite difference solvers for singular solutions of the elliptic Monge–Ampère equation. *J. Comput. Phys.* **230** (3), 818–834.

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

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

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

HAITSMA, J.J., LACHMANN, U. & LACHMANN, B. 2001 Exogenous surfactant as a drug delivery agent. *Adv. Drug Deliv. Rev.* **47** (2–3), 197–207.

HALLIDAY, H.L. 2008 Surfactants: past, present and future. *J. Perinatol.* **28** (1), S47–S56.

HALPERN, D., JENSEN, O.E. & GROTBERG, J.B. 1998 A theoretical study of surfactant and liquid delivery into the lung. *J. Appl. Physiol.* **85** (1), 333–352.

HENKEL, C., SNOEIJER, J.H. & THIELE, U. 2021 Gradient-dynamics model for liquid drops on elastic substrates. *Soft Matt.* **17** (45), 10359–10375.

HIDALGO, A., CRUZ, A. & PÉREZ-GIL, J. 2015 Barrier or carrier? Pulmonary surfactant and drug delivery. *Eur. J. Pharm. Biopharm.* **95**, 117–127.

IASELLA, S., SHARMA, R., GAROFF, S. & TILTON, R.D. 2024 Interaction of impinging Marangoni fields. *J. Colloid Interface Sci.* **653**, 807–820.

ISHII, T. & MURO, J. 1989 The fractal aspects of Suminagashi and marbling. *Thin Solid Films* **178**, 109–113.

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

JENSEN, O.E., HALPERN, D. & GROTBERG, J.B. 1994 Transport of a passive solute by surfactant-driven flows. *Chem. Engng Sci.* **49** (8), 1107–1117.

JOBE, A.H. 1993 Pulmonary surfactant therapy. *New Engl. J. Med.* **328** (12), 861–868.

JORDAN, R., KINDERLEHRER, D. & OTTO, F. 1998 The variational formulation of the Fokker–Planck equation. *SIAM J. Math. Anal.* **29** (1), 1–17.

KANTOROVICH, L.V. 1942 On the translocation of masses. *Dokl. Akad. Nauk SSSR* **37** (7–8), 227–229.

KANTOROVICH, L.V. 2006 On the translocation of masses. *J. Math. Sci.* **133** (4), 1381–1382.
KOLOURI, S., PARK, S.R., THORPE, M., SLEPCEV, D. & ROHDE, G.K. 2017 Optimal mass transport: signal processing and machine-learning applications. *IEEE Sig. Proc. Mag.* **34** (4), 43–59.

MAHAN, B. 2011 Graphical artist https://www.beamahan.com. Pictures taken from flickr (https://www.flickr.com/photos/beamahan/5618566534/in/album-72157626370020523/, https://www.flickr.com/photos/beamahan/5932561347/in/album-72157626370020523/) under licence CC BY-NC-ND 2.0 DEED: https://creativecommons.org/licenses/by-nc-nd/2.0/. Figure 1(b) has been cropped, figure 10(f) has been cropped and rotated.

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

MCNAIR, R., JENSEN, O.E. & LANDEL, J.R. 2022 Surfactant spreading in a two-dimensional cavity and emergent contact-line singularities. *J. Fluid Mech.* **930**, A15.

MCNAIR, R., TEMPRANO-COLETO, F., PEAUDECERF, F.J., GIBOU, F., LUZZATTO-FEGIZ, P., JENSEN, O.E. & LANDEL, J.R. 2023 Exogenous-endogenous surfactant interaction yields heterogeneous spreading in complex branching networks. *Phys. Rev. Lett.* (submitted), arXiv:2311.08133.

MEIRMANOV, A.M., PUKHNACHEV, V.V. & SHMAREV, S.I. 1997 Evolution equations and lagrangian coordinates. In *De Gruyter Expositions in Mathematics* (ed. O.H. Kegel, V.P. Maslov, W.D. Neumann & R.O. Wells), vol. 24. de Gruyter.

MONGE, G. 1781 Mémoire sur la théorie des déblais et des remblais. *Mem. Math. Phys. Acad. Royale Sci.*, 666–704.

OTTO, F. 2001 The geometry of dissipative evolution equations: the porous medium equation. *Commun. Part. Diff. Equ.* **26**, 101–174.

ROCKAFELLAR, R.T. 1970 Convex analysis. In *Princeton Landmarks in Mathematics*, vol. 28. Princeton University Press.

RODRIGUEZ, R.J. 2003 Management of respiratory distress syndrome: an update. *Respir. Care* **48** (3), 279–287.

ROUWET, D. & IORIO, M. 2017 The sedimentation of Suminagashi-like floating clay patterns at El Chichón crater lake (Chiapas, Mexico). *Geol. Soc. Lond. Spec. Publ.* **437** (1), 153–161.

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

THESS, A., SPIRN, D. & JÜTTNER, B. 1997 A two-dimensional model for slow convection at infinite Marangoni number. *J. Fluid Mech.* **331**, 283–312.

THIELE, U., ARCHER, A.J. & PISMEN, L.M. 2016 Gradient dynamics models for liquid films with soluble surfactant. *Phys. Rev. Fluids* **1** (8), 083903.

THIELE, U., SNOEIJER, J.H., TRINSCHEK, S. & JOHN, K. 2018 Equilibrium contact angle and adsorption layer properties with surfactants. *Langmuir* **34** (24), 7210–7221.
KOLOURI, S., PARK, S.R., THORPE, M., SLEPCEV, D. & ROHDE, G.K. 2017 Optimal mass transport: signal processing and machine-learning applications. *IEEE Sig. Proc. Mag.* **34** (4), 43–59.

MAHAN, B. 2011 Graphical artist https://www.beamahan.com. Pictures taken from flickr (https://www.flickr.com/photos/beamahan/5618566534/in/album-72157626370020523/, https://www.flickr.com/photos/beamahan/5932561347/in/album-72157626370020523/) under licence CC BY-NC-ND 2.0 DEED: https://creativecommons.org/licenses/by-nc-nd/2.0/. Figure 1(b) has been cropped, figure 10(f) has been cropped and rotated.

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

MCNAIR, R., JENSEN, O.E. & LANDEL, J.R. 2022 Surfactant spreading in a two-dimensional cavity and emergent contact-line singularities. *J. Fluid Mech.* **930**, A15.

MCNAIR, R., TEMPRANO-COLETO, F., PEAUDECERF, F.J., GIBOU, F., LUZZATTO-FEGIZ, P., JENSEN, O.E. & LANDEL, J.R. 2023 Exogenous-endogenous surfactant interaction yields heterogeneous spreading in complex branching networks. *Phys. Rev. Lett.* (submitted), arXiv:2311.08133.

MEIRMANOV, A.M., PUKHNACHEV, V.V. & SHMAREV, S.I. 1997 Evolution equations and lagrangian coordinates. In *De Gruyter Expositions in Mathematics* (ed. O.H. Kegel, V.P. Maslov, W.D. Neumann & R.O. Wells), vol. 24. de Gruyter.

MONGE, G. 1781 Mémoire sur la théorie des déblais et des remblais. *Mem. Math. Phys. Acad. Royale Sci.*, 666–704.

OTTO, F. 2001 The geometry of dissipative evolution equations: the porous medium equation. *Commun. Part. Diff. Equ.* **26**, 101–174.

ROCKAFELLAR, R.T. 1970 Convex analysis. In *Princeton Landmarks in Mathematics*, vol. 28. Princeton University Press.

RODRIGUEZ, R.J. 2003 Management of respiratory distress syndrome: an update. *Respir. Care* **48** (3), 279–287.

ROUWET, D. & IORIO, M. 2017 The sedimentation of Suminagashi-like floating clay patterns at El Chichón crater lake (Chiapas, Mexico). *Geol. Soc. Lond. Spec. Publ.* **437** (1), 153–161.

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

THESS, A., SPIRN, D. & JÜTTNER, B. 1997 A two-dimensional model for slow convection at infinite Marangoni number. *J. Fluid Mech.* **331**, 283–312.

THIELE, U., ARCHER, A.J. & PISMEN, L.M. 2016 Gradient dynamics models for liquid films with soluble surfactant. *Phys. Rev. Fluids* **1** (8), 083903.

THIELE, U., SNOEIJER, J.H., TRINSCHEK, S. & JOHN, K. 2018 Equilibrium contact angle and adsorption layer properties with surfactants. *Langmuir* **34** (24), 7210–7221.