# Exogenous-endogenous surfactant interaction yields heterogeneous spreading in complex branching networks

Richard Mcnair, Fernando Temprano-Coleto, François Peaudecerf, Frédéric Gibou, Paolo Luzzatto-Fegiz, Oliver E Jensen, Julien R. Landel

## ► To cite this version:

Richard Mcnair, Fernando Temprano-Coleto, François Peaudecerf, Frédéric Gibou, Paolo Luzzatto-Fegiz, et al.. Exogenous-endogenous surfactant interaction yields heterogeneous spreading in complex branching networks. Physical Review Letters, 2025, 134 (3), pp.034001. 10.1103/Phys-RevLett.134.034001. hal-04833558v2

HAL Id: hal-04833558

https://hal.science/hal-04833558v2

Submitted on 24 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Exogenous-endogenous surfactant interaction yields heterogeneous spreading in complex branching networks

Richard Mcnair¹, Fernando Temprano-Coleto²,³, François J. Peaudecerf⁴,

Frédéric Gibou⁵, Paolo Luzzatto-Fegiz⁵, Oliver E. Jensen¹, Julien R. Landel⁶,¹,*

¹Department of Mathematics, University of Manchester, Manchester M13 9PL, UK

²Andlinger Center for Energy and the Environment, Princeton University, Princeton, New Jersey 08544, USA

³Department of Mechanical and Aerospace Engineering,
Princeton University, Princeton, New Jersey 08544, USA

⁴Univ Rennes, CNRS, IPR (Institut de Physique de Rennes) – UMR 6251, F - 35000 Rennes, France

⁵Department of Mechanical Engineering, University of California, Santa Barbara, California 93106, USA

⁶Université Claude Bernard Lyon 1, Laboratoire de Mecanique des Fluides et d'Acoustique (LMFA),
UMR5509, CNRS, Ecole Centrale de Lyon, INSA Lyon, 69622 Villeurbanne, France

(Dated: December 10, 2024)

Experiments have shown that surfactant introduced to a liquid-filled maze can find the solution path. We reveal how the maze-solving dynamics arise from interactions between the added surfactant and endogenous surfactant present at the liquid surface. We simulate the dynamics using a nonlinear model solved with a discrete mimetic scheme on a graph. Endogenous surfactant transforms local spreading into a non-local problem with an omniscient view of the maze geometry, key to the maze-solving dynamics. Our results offer insight into surfactant-driven transport in complex networks such as lung airways.

Recent experiments [1] showed that surfactant, made visi-
ble by dye, can ‘solve’ a liquid maze (Fig. 1(a), Supplemen-
tary Movie 1 [2]). The surfactant propagates along the solu-
tion path, with little advance into lateral branches. This phe-
nomenon seems counterintuitive, as no mechanism appears to
draw the surfactant preferentially to the maze exit: one might
expect Marangoni stresses to spread the surfactant throughout
the whole maze. Yet, the surfactant proceeds with apparent
awareness of the whole geometry. Although no quantitative
theory exists, a hypothesis was suggested based on the poten-
tial effect of pre-existing, ‘endogenous’ surfactant in the maze
liquid [1]. The spread of the ‘exogenous’ (added) surfactant
would be influenced by the endogenous surfactant, in a man-
ner sensitive to the maze geometry. Inspired by these obser-
vations, we investigate the impact of endogenous surfactant
on the spreading of exogenous surfactant on thin liquid layers
confined in asymmetric branching networks.

Whether from natural or contaminant sources, endogenous
surfactants are usually undetectable *a priori*, but have macro-
scopic influence on flow dynamics across numerous applica-
tions [3–8]. Natural endogenous surfactant in human lungs
can affect exogenous surfactant therapies [9–14] for vari-
ous diseases, including acute respiratory distress syndrome
(ARDS) [15]. Although efficacy of surfactant replacement
therapy is indisputable for neonatal respiratory distress syn-
drome [16–19], there remains debate for other pathologies
[20–23], including ARDS, which has a 40% mortality rate in
the US [24]. To understand mixed results from clinical tri-
als [20], experimental and theoretical models analysed physi-
cal, chemical and biological processes in surfactant transport
through branching networks simulating lungs [13, 25, 26].
In distal lung airways, where surfactant therapy is effec-
tive, Marangoni transport dominates gravitational transport
[11, 13, 27, 28]. Endogenous surfactants can affect therapeutic
surfactant transport, preventing it from reaching distal airways

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Photograph of the liquid maze experiment [1]: surfactant deposited at the inlet finds the solution path with minimal penetration into lateral branches, as shown by red dye advected by Marangoni forces. (b) Results from our model (2) showing the concentration (blue-scale) of the surfactant field in the same geometry; the red front is the exogenous-endogenous boundary.</figcaption></figure>

[9–12, 14, 29, 30].

Pulmonary surfactant models have used simplifying geo-
metric assumptions to reduce analytical, computational and
experimental difficulties. Single-airway models [14, 31]
may suppress crucial dynamics associated with branching
[12, 32]. More complex models used a Weibel morphometry
[33], which yields a dichotomous, symmetric and self-similar
branching network [9–12, 27, 34, 35]; however, human lungs
are asymmetric and highly irregular [36–38]. Recent mod-
els [32, 39], which neglected Marangoni spreading, showed
that transport by air flow and gravity in asymmetric lungs can
yield heterogeneity in surfactant distribution. The maze ex-
periments of [1] were not designed to model lung surfactant
transport; nevertheless, the two problems share key physical
features. Understanding surfactant dynamics in mazes can
give insight about heterogeneities affecting surfactant spread-
ing in large complex networks, such as human lungs.

To test the hypothesis that endogenous surfactant explains
the maze-solving behavior of exogenous surfactant [1], and
provide new physical insight, we model surfactant transport on a thin film confined in a 1D branching network, replicating the maze geometry from [1]. The experimental maze consists of a network of grooves in an acrylic substrate, with larger 20 x 20 mm branches at the inlet and outlet (Fig. 1(a)). The large outlet creates an asymmetry in the geometry. Branches were filled with a 50% milk-cream mixture, which provides a dense, viscous liquid contrasting with red tracer dye added at the inlet. The aqueous dye is soluble in the milk and slightly buoyant, and thus remains mostly at the surface. A surfactant drop (0.3% w/w soap-water solution) was deposited at the inlet at $t = 0$, starting a Marangoni flow that entrains the dye. Owing to the solubility of all components, the milk, dye and surfactant solution form a single liquid phase. Within 60 s, one of the red dye fronts, assumed to follow the exogenous surfactant fronts, reached the outlet, through the solution path of length 155 mm (Sec. S1 [40]).

We assume the dynamics are governed by a balance of viscous and Marangoni stresses, that solubility and diffusion of surfactant are negligible, and that inertia and interface deformation are negligible, as shown quantitatively in [40]. Endogenous and exogenous contributions are modeled as a single surfactant concentration field $\Gamma$. The axial and vertical coordinates are $x$ and $z$, $t$ is time, $u|_s$ the surface axial velocity, and $\gamma$ the surface tension. Variables are normalized by appropriate physical scales [40], such that e.g. $0 \le z \le 1$. The key equations expressing surfactant conservation and stress balance at the surface are therefore, in nondimensional variables,

$$ \Gamma_t + (u|_s \Gamma)_x = 0, \quad \text{and} \quad \gamma_x = u|_s, \quad (1a,b) $$

where subscripts in $x, z, t$ denote partial derivatives, and $u|_s$ is the viscous stress at the surface. To introduce key ideas, ignore initially the effect of lateral ($y$) confinement, such that to leading order the viscous flow in the interior is governed by $u_{zz} = p_x$. Integrating twice in $z$ and assuming that the volume flux along the channel is uniformly zero, such that $\int_0^1 u dz = 0$, one finds the surface velocity $u|_s = p_x/6$ and viscous stress $u|_s = 2p_x/3$. Eliminating $p_x$ gives $u|_s = 4u|_s$. The nondimensional surfactant concentration and surface tension are related by $\gamma_x = -\Gamma_x$ [40]; substituting this relation, together with $u|_s = 4u|_s$, into (1b) yields $u|_s = -\frac{1}{4}\Gamma_x$. Eliminating $u|_s$ in (1a) yields an equation that only involves $\Gamma$,

$$ \Gamma_t = \frac{1}{4}(\Gamma \Gamma_x)_x \equiv \frac{1}{8}(\Gamma^2)_{xx}. \quad (2) $$

An alternative derivation of (2), starting from lubrication theory, is possible [40]. If we include the effect of lateral no-slip boundaries at $y = \pm W/2$, the factor of 1/4 in (2) and $u|_s$ is replaced by a positive coefficient $f(W) < 1/4$, calculated a priori as a function of $W$. This factor is a nonessential time scaling, as we rescale $t$ by experimental completion time $t_{ref}$; however, this analysis shows that narrow channels significantly reduce the surface velocity, helping to keep the flow within the viscous limit (for $[1] W \approx 0.94$ yielding $f \approx 0.15$).

Model (2) is solved numerically in a branching network simulating the experimental maze topology, illustrated in

<figure><img src="image_3.png" /><figcaption>FIG. 2. Schematic showing how the transport model (2) is implemented in a maze mimicking the experiment, with the branch numbering scheme indicated (Fig. S1 [40]). The maze topology is modelled by a 1D branching network (blue). Along each edge, the model results are used to advect the exogenous surfactant fronts (red).</figcaption></figure>

Fig. 2. At junctions, we impose continuity of surfactant concentration and flux, whilst no flux is imposed at dead-ends (Sec. S2 [40]). The dye fronts are computed by Lagrangian integration of the surface velocity $u|_s$. To model the initial conditions, we assume that exogenous surfactant is deposited instantaneously and uniformly at the inlet at $t = 0$. The normalised initial concentration is set to $\Gamma = 1$ at the inlet, whilst the rest of the maze has initial endogenous surfactant concentration $\Gamma = \delta$ everywhere, where $0 < \delta < 1$ is an empirical parameter. The branches are numbered as illustrated in Fig. 2 (see [40] for further details). Branch $i$ of length $L_i$ (normalised by the solution path length) has longitudinal coordinate $0 \le x_i \le L_i$. Internal branches have lengths representative of the experiment. The greater width of the inlet ($i = 0$) and outlet ($i = 37$) branches is treated by constructing 1D equivalents. Since the outlet has small variations in concentration, its length $L_{37}$ is calculated by dividing its surface area by the width of an interior branch. The inlet length $L_0$ is calculated by matching the amount of initial exogenous surfactant.

We estimate the ratio of endogenous to exogenous surfactant mass from the equilibrium state of the experiment, as encoded by the ratio of white versus red surface areas (Sec. S2.4 [40]). To avoid numerical instabilities, care is taken to regularize the initial surfactant concentration distribution between the inlet and the first branch (Sec. S2.2 [40]). The time at which the front enters branch 1 is labelled as $t = 0$. As noted earlier, in the model, we re-scale time so that the completion time (when the front enters the outlet) matches the experimental completion time $t_{ref}$ up to a multiplying factor $\tau_{end}$ chosen between $0.5 \le \tau_{end} \le 1.5$.

To solve the system of equations of the form (2), coupled between each of the branches of the maze, we use a mimetic finite-difference numerical method, which balances fluxes at junctions to machine precision. Mimetic finite-difference
methods [41–43] use tools from graph theory and discrete calculus to obtain finite-difference operators which reproduce desired features of vector calculus operators. As shown in Fig. 2, we represent the maze as a graph with vertices at every junction and dead end. Extra vertices are placed along each edge to refine the discretisation. We construct finite-difference mimetic operators on this graph using the incidence matrix of the network A, the columns of which correspond to the vertices of the graph, and the rows of which correspond to the edges. Each edge is assigned an orientation directed away from the maze inlet. Components of A are zero except at vertices where an edge leaves (−1), and where an edge enters (+1). Metric information is encapsulated by the diagonal edge and vertex length matrices L and V respectively. Each component of L is the length of an edge, and each component of V is one half the sum of lengths of every edge incident to the corresponding vertex — an important nuance ensuring mass conservation (see Sec. S3.1 [40]). The matrix L⁻¹A is therefore a finite difference gradient operator acting on scalar functions defined on the graph’s vertices, and is second order accurate at edge midpoints. Likewise, −V⁻¹Aᵀ is a divergence operator acting on functions defined at edge midpoints, and is second order accurate at the vertices. The concentration of surfactant at the vertices of the graph form the components of vector Γ(t). The system of differential equations of the form (2) solved on the graph is found by combining gradient and divergence operators into a Laplacian, and is therefore

$$ \frac{d\Gamma}{dt} = -\frac{f(W)}{2}V^{-1}A^{T}L^{-1}A\Gamma^{2}, \quad (3) $$

with Γ² the component-wise square of Γ. This mimetic finite-difference formulation automatically implements no-flux boundary conditions for all bounding vertices of the graph (dead ends), and continuity of concentration and flux for all junctions within the graph. A semi-implicit scheme solves this system of equations in time (Sec. S3.2 [40]). Simultaneously, we solve the equation tracking the dye fronts, which in branch i at $x_i = \Lambda_i(t)$ is $\partial\Lambda_i/\partial t = u_i(\Lambda_i, t) = -f(W)\Gamma_x(\Lambda_i, t)$. When a front reaches a junction, new fronts are created in the downstream branches.

To reduce computational costs and provide additional insight, we consider the limit of small concentration gradients, which arises in applications where concentration differences between exogenous and endogenous surfactants are small. We decompose the concentration into a time- and space-averaged component and a fluctuating component: $\Gamma = \bar{\Gamma} + \hat{\Gamma}(x, t)$, with $|\hat{\Gamma}| \ll \bar{\Gamma}$. Model (2) simplifies to the linear diffusion equation $\hat{\Gamma}_t = f(W)\bar{\Gamma}\hat{\Gamma}_{xx}$. In the discrete mimetic representation, the Laplacian operator $f(W)\bar{\Gamma}V^{-1}A^T L^{-1}A$ has orthogonal eigenvectors $\phi_n$ under an inner product $\langle\phi_l, \phi_m\rangle = \phi_l^T V \phi_m$ (Sec. S4 [40]). The vector of vertex concentrations can be decomposed in a sum of eigenmodes, namely $\Gamma(t) = \sum_{n=1}^N A_n \phi_n e^{-\lambda_n t}$, with amplitudes $A_n = \langle\phi_n, \Gamma(t=0)\rangle$; $\lambda_n$ are the corresponding eigenvalues. In the linear regime, the concentration in the maze can be approximated without the need for time integration, by instead using only a finite number N of eigenmodes.

<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Front location $\Lambda_M(t)$ in time along the solution path, normalised by overall length $L_M$, from experiments (red curves, experimental error shaded pink), nonlinear model (black curves), and 3-mode linear model (dashed curves). Inset shows lateral branch 4. (b) First 3 eigenmodes of the linear model along the solution path, with eigenvalues $\lambda$ shown. Thick vertical black lines show the start and end points. Dashed lines show junctions.</figcaption></figure>

We test our models by comparing their predictions to the dye front locations in the experiments, as shown in Fig. 3(a), which plots time versus location along the global solution path; the inset shows data for a lateral branch. The nonlinear model (shown by the black continuous line) captures the spreading dynamics. Its least-square error relative to the experiment is minimized by setting the order-one coefficients $\delta = 0.15$ and $\tau_{end} = 0.24$ [40]. In the experiments, fronts in lateral branches are observed to eventually recede after initially advancing; this is most visible near the inlet (Supplementary Movie [2]); an example is shown in the inset of Fig. 3(a). This receding effect occurs as the concentration at the beginning of the lateral branch eventually decreases over time, as the exogenous surfactant spreads across the rest of the maze, yielding a reversal of the Marangoni stress. The nonlinear model captures the flow reversal qualitatively, generally overpredicting relaxation time (inset of Fig. 3(a)). The nonlinear model predictions are not sensitive to the exact value of $\delta$: any value in the range $0.02 \le \delta \le 0.21$ captures the dynamics in the experiment.

Since the shapes of the eigenmodes depend only on the network geometry, transport dynamics in branching networks can be studied from the dominant eigenmodes and eigenvalues, without having to solve the full nonlinear transport model. The maze-solving dynamics described above are already reproduced by a linear model truncated to just three modes, as
shown by the dashed lines in Fig. 3(a). The front location on the global solution path is captured quantitatively. The spatial structure of each mode is in Fig. 3(b). The zeroth mode represents the surfactant concentration at steady state, $t \to \infty$. The first mode, $\lambda_1 = 0.75$, has a decreasing profile along the solution path, corresponding to the main concentration gradient driving the front, key to the maze-solving behaviour. The second mode, $\lambda_2 = 2.52$, presents a wavy profile, which, in combination with the first mode, controls the receding dynamics observed in lateral branches. The exact temporal evolution requires computation of the mode amplitudes, which depend on initial conditions. If three modes were used, $\delta \approx 0.03$ minimized the least-square error. Using a larger number of modes, e.g. of order 100, recovers $\delta \approx 0.15$ found by the nonlinear model. The main challenge for the linear model is to capture the nonlinear distribution of concentration at $t = 0$, better approximated using a larger number of modes.

To assess the importance of 2D dynamics neglected in the model, we perform a 2D simulation of the surfactant transport through branches 7–14 (see Fig. 4). We use COMSOL® to solve the 2D analogue of the 1D model (2), namely $\Gamma_t = \nabla \cdot (\Gamma \nabla \Gamma)/4$, in the maze geometry. A good match to experiment is achieved by imposing $\Gamma = 1$ at the entrance of branch 7 and $\Gamma = 0.7$ at the exit of branch 14, with initial condition $\Gamma = 0.7$ everywhere; as before, care is taken to regularize the initial $\Gamma$ field [40]. We use the particle-tracing module to track the exogenous–endogenous surfactant front. Dominant dynamical features are reproduced by the 2D model (Fig. 4). In the experiments, the exogenous surfactant tends to spread along the inside of sharp corners (blue arrows), whilst it spreads through the middle of the path as it enters new branches from a junction (green arrows). The 2D simulation captures qualitatively these behaviours, which are due to the compression of surfactant concentration level curves at inside corners (Fig. 4(b)). The nonlinear model appears consistent with surfactant dynamics on 2D confined surfaces. However, the 2D simulation does not capture some of the gaps that appear between the maze walls and dye (orange arrows). These features could be due to aspects omitted from the model such as geometrical imperfections, surface deformation at menisci or chemical heterogeneities.

Our model reveals the importance of exogenous–endogenous surfactant interaction in the maze-solving behaviour. Although the model makes assumptions about the endogenous surfactant, our results show that the model is robust; the maze-solving dynamics continue to be reproduced as the ratio of endogenous to exogenous concentration is varied across one order of magnitude. To confirm that this conclusion remains valid for systems comprising other liquids, surfactants, and tracers, we also designed and performed additional experiments, as reported in [40]. These new experiments show that the same maze solving dynamics are observed, as long as laminar flow is maintained. To this end, the original experiments of [1] leveraged the fact that soap in milk induces mild surface tension changes, as compared to typical surfactants in water [44]. As an example, we show

<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Experiment, showing branches 7–14. (b) 2D nonlinear simulation. Tracers (red dots) are advected from their starting position (blue dots) along Lagrangian trajectories (black lines). Experimental features are reproduced qualitatively: the front spreads inside sharp corners (blue arrows) and centrally at straight junctions (green arrows), owing to the bending of the concentration contours.</figcaption></figure>

that the maze experiment can be performed with another mild surfactant, namely 2-propanol in glycerol, and using tracer particles (rather than dye) for visualization [40]. Through the invisible endogenous surfactant, the exogenous surfactant is ‘aware’ of the whole network topology. The large area of the outlet compared to that of the lateral branches ensures that large concentration gradients remain along the solution path, thereby driving the exogenous surfactant to the outlet. In lateral branches, compression of endogenous surfactant limits access by exogenous surfactant. The elliptic nature of the spatial differential operator in (2) means that the surfactant spreads through the maze with an *omniscient* view of the geometry, since together the exogenous and endogenous surfactants fill the whole maze at all times.

Although our model is not designed to simulate surfactant transport in lungs, it may provide insight and novel methodologies for studying surfactant therapies for lungs. Its essential ingredients, endogenous surfactant [9–14] and network asymmetry [36–38], are characteristic of lungs. Exogenous surfactant delivering drugs to the distal airways and alveoli may spread non-uniformly through the lung, reducing drug efficacy. The model could inspire lung-scale models to test the impact of lung network asymmetry on surfactant drug delivery. Our mimetic finite-difference implementation on a network is amenable to scale-up, as we found from preliminary tests on a network with $O(10^4)$ branches based on real lung scans. If the linear regime applies, the transport model decomposes into a series of modes computed from the network topology, thus helping understand surfactant transport without computationally expensive simulations. Our model frame-
work provides a novel approach compared to single-branch
models [9–12, 34, 35], which overlook complex lung topol-
ogy, and molecular dynamic simulations [45–48], which are
numerically expensive to run over a whole lung.

Our model reveals how the combination of asymmetry in
the maze network and exogenous-endogenous surfactant in-
teractions are key to understanding confined transport prob-
lems in complex branching networks, such as lungs.

J.R.L. and O.E.J. acknowledge financial support from EP-SRC grant EP/T030739/1. P.L.-F. acknowledges support from NSF CAREER grant 2048234. F.T.-C. acknowledges support from a Distinguished Postdoctoral Fellowship from the Andlinger Center from Energy and the Environment at Princeton University. For the purpose of Open Access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

* Corresponding author: julien.landel@univ-lyon1.fr

[1] F. Temprano-Coleto, F. J. Peaudecerf, J. R. Landel, F. Gibou, and P. Luzzatto-Fegiz, Phys. Rev. Fluids 3, 100507 (2018).

[2] Supplementary videos at [URL will be inserted by publisher] of our nonlinear simulation show how the surfactant field in the lateral branches compresses and then relaxes at later times, leading to the receding behaviour observed.

[3] F. J. Peaudecerf, J. R. Landel, R. E. Goldstein, and P. Luzzatto-Fegiz, Proc. Nat. Acad. Sci. USA 114, 7254 (2017).

[4] A. Hourlier-Fargette, A. Antkowiak, A. Chateauminois, and S. Neukirch, Soft Matter 13, 3484 (2017).

[5] H. Manikantan and T. M. Squires, J. Fluid Mech. 892, 10.1017/jfm.2020.170 (2020).

[6] G. Koleski, A. Vilquin, J.-C. Loudet, T. Bickel, and B. Pouligny, Phys. Fluids 32, 092108 (2020).

[7] A. Faugaret, Y. Duguet, Y. Fraigneau, and L. Martin Witkowski, J. Fluid Mech. 900, A42 (2020).

[8] F. Temprano-Coleto, S. M. Smith, F. J. Peaudecerf, J. R. Landel, F. Gibou, and P. Luzzatto-Fegiz, Proc. Natl. Acad. Sci. USA 120, 10.1073/pnas.2211092120 (2023).

[9] F. F. Espinosa, A. H. Shapiro, J. J. Fredberg, and R. D. Kamm, J. Appl. Physiol. 75, 2028 (1993).

[10] J. Grotberg, D. Halpern, and O. Jensen, J. Appl. Physiol. 78, 750 (1995).

[11] D. Halpern, O. E. Jensen, and J. B. Grotberg, J Appl Physiol 85, 333 (1998).

[12] Y. L. Zhang, O. K. Matar, and R. V. Craster, Med. Eng. Phys. 25, 115 (2003).

[13] A. Z. Stetten, S. V. Iasella, T. E. Corcoran, S. Garoff, T. M. Przybycien, and R. D. Tilton, Curr. Op. Colloid Interface Sci. 36, 58 (2018).

[14] C. García-Mouton, M. Echaide, L. A. Serrano, G. Orel-
lana, F. Salomone, F. Ricci, B. Pioselli, D. Amidani,
A. Cruz, and J. Pérez-Gil, Pharmaceutics 15, 10.3390/pharma-
ceutics15010256 (2023).

[15] A. Anzueto, R. P. Baughman, K. K. Guntupalli, J. G. Weg, H. P. Wiedemann, A. A. Raventós, F. Lemaire, W. Long, D. S. Zaccardelli, and E. N. Pattishall, New Eng. J. Med. 334, 1417 (1996).

[16] M. E. Avery and J. Mead, Am. J. Dis. Child. 97, 517 (1959).

[17] A. H. Jobe, New Eng. J. Med. 328, 861 (1993).

[18] R. J. Rodriguez, Resp. Care 48, 279 (2003).

[19] H. Halliday, J. Perinat. 28, S47 (2008).

[20] S. V. Baudouin, New England Journal of Medicine 351, 853 (2004).

[21] Y. Y. Zuo, R. A. W. Veldhuizen, A. W. Neumann, N. O. Petersen, and F. Possmayer, Biochimica et Biophysica Acta-Biomembranes 1778, 1947 (2008).

[22] M. J. Sankar, N. Gupta, K. Jain, R. Agarwal, and V. K. Paul, J. Perinatol. 36, S36–S48 (2016).

[23] R. A. W. Veldhuizen, Y. Y. Zuo, N. O. Petersen, J. F. Lewis, and F. Possmayer, Expert Review of Respiratory Medicine 15, 597 (2021).

[24] G. D. Rubenfeld, E. Caldwell, E. Peabody, J. Weaver, D. P. Martin, M. Neff, E. J. Stern, and L. D. Hudson, New England Journal of Medicine 353, 1685 (2005), pMID: 16236739.

[25] R. Levy, D. B. Hill, M. G. Forest, and J. B. Grotberg, Integrative Comparative Biol. 54, 985 (2014).

[26] B. Baer, L. M. P. Souza, A. S. Pimentel, and R. A. Veldhuizen, Biochem. Pharmacol. 164, 64 (2019).

[27] F. F. Espinosa and R. D. Kamm, J. Appl. Physiol. 86, 391
(1999).

[28] J. B. Grotberg, Ann. Rev. Biomed. Eng. **3**, 421 (2001).

[29] O. E. Jensen and J. B. Grotberg, J. Fluid Mech. **240**, 259 (1992).

[30] J. L. Bull, L. K. Nelson, J. Walsh, J. T., M. R. Glucksberg, S. Schürch, and J. B. Grotberg, J. Biomech. Eng. **121**, 89 (1999).

[31] B. Baer, E. J. Veldhuizen, F. Possmayer, C. Yamashita, and R. Veldhuizen, Discov. Med. **26**, 207 (2018).

[32] A. Kazemi, B. Louis, D. Isabey, G. F. Nieman, L. A. Gatto, J. Satalin, S. Baker, J. B. Grotberg, and M. Filoche, PLOS Comp. Biol. **15**, 1 (2019).

[33] E. R. Weibel and D. M. Gomez, Science **137**, 577 (1962).

[34] M. Filoche, C.-F. Tai, and J. B. Grotberg, Proc. Nat. Acad. Sci. **112**, 9287 (2015).

[35] A. Copploe, M. Vatani, C. J. W., and H. Tavana, Ann. Biomed. Eng. **47**, 1435 (2019).

[36] K. Horsfield and G. Cumming, J. Appl. Physiol. **24**, 373 (1968).

[37] M. Tawhai, P. Hunter, J. Tschirren, J. Reinhardt, G. McLennan, and E. Hoffman, J. Appl. Physiol. **97**, 2310 (2004).

[38] A. Schmidt, S. Zidowitz, A. Kriete, T. Denhard, S. Krass, and H. Peitgen, Computerized Med. Imaging Graphics **28**, 203 (2004).

[39] S. Elias-Kirma, A. Artzy-Schnirman, H. Sabatan, C. Dabush, D. Waisman, and J. Sznitman, Journal of biomechanics **122**, 110458 (2021).

[40] See Supplemental Material at [URL will be inserted by publisher] for further details about the models, their numerical simulation and data analysis, which includes Refs [49–56].

[41] A. N. Tikhonov and A. A. Samarskii, USSR Comp. Math. & Math. Phys. **1**, 5 (1962).

[42] A. van der Schaft, Systems & Control Lett. **101**, 21 (2017).

[43] J. M. Hyman and M. Shashkov, Computers & Math. with Appl. **33**, 81 (1997).

[44] Q. Tian, L. Lai, Z. Zhou, P. Mei, Q. Lu, Y. Wang, D. Xiang, and Y. Liu, J. Agric. Food Chem. **67**, 6336 (2019).

[45] C. Laing, S. Baoukina, and P. D. Tieleman, Phys. Chem. Chem. Phys. **11**, 1916 (2009).

[46] S. Baoukina and P. D. Tieleman, Biophys. J. **100**, 1678 (2011).

[47] L. M. P. Souza, J. B. Nascimento, A. L. Romeu, E. D. Estrada-López, and A. S. Pimentel, Coll. Surf. B: Biointerfaces **167**, 345 (2018).

[48] L. M. P. Souza, M. C. Lima, L. F. S. Bezerra, and A. S. Pimentel, Chem. Phys. **563**, 111704 (2022).

[49] W. D. Rutz, C. Whitnah, and G. Baetz, J. Dairy Sci. **38**, 1312 (1955).

[50] A. Bakshi and D. Smith, J. Dairy Sci. **67**, 1157 (1984).

[51] D. Kristensen, P. Jensen, F. Madsen, and K. Birdi, J. Dairy Sci. **80**, 2282 (1997).

[52] O. E. Jensen and D. Halpern, J. Fluid Mech. **372**, 273 (1998).

[53] R. Mcnair, O. E. Jensen, and J. R. Landel, J. Fluid Mech. **930**, A15 (2022).

[54] M. Brio, J.-G. Caputo, and H. Kravitz, Appl. Numer. Math. **172**, 99 (2022).

[55] C. H. Chang and E. I. Franses, Colloids Surf. A **100**, 1 (1995).

[56] D. Ershov, M.-S. Phan, J. W. Pylvänäinen, S. U. Rigaud, L. Le Blanc, A. Charles-Orszag, J. R. W. Conway, R. F. Laine, N. H. Roy, D. Bonazzi, G. Duménil, G. Jacquemet, and J.-Y. Tinevez, Nat. Methods **19**, 829 (2022).
SUPPLEMENTAL MATERIAL

Exogenous–endogenous surfactant interaction yields heterogeneous spreading in complex branching networks

Richard Mcnair,<fnref n="1" /> Fernando Temprano-Coleto,<fnref n="2" /> François J. Peaudecerf,<fnref n="3" /> Frédéric Gibou,<fnref n="4" /> Paolo Luzzatto-Fegiz,<fnref n="4" /> Oliver E. Jensen,<fnref n="1" /> Julien R. Landel<fnref n="5" />,<fnref n="1" />,<fnref n="1" />

(Dated: December 10, 2024)

S1 Maze geometry and experimental parameters

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

Figure S1: Left: Picture at $t^* = 0$ s of the experiment conducted by Temprano-Coleto et al. [1], showing a numbering scheme for the branches. Right: Graph of the mobility $f(W)$ found from (S4g), with the sum truncated at $n = 1000$, plotted against the branch aspect ratio $W$ in blue, with $f(0.94)$ circled in red. The sum (S4g) converges to 2 significant figures for $n < 100$ for the values of $W$ shown. We take $W = 0.94$ to be the aspect ratio of the branches of the maze.

We model the transport dynamics in the maze experiment conducted by Temprano-Coleto et al. [1], as presented in their Fig. 2 and with the corresponding video showing the large maze at [1] (at 1:47 min). A picture of the experiment at what we define to be the initial time, $t^* = 0$ s (with a star denoting dimensional quantities), is shown in Fig. S1 (right), along with a branch numbering scheme. The experimental conditions and protocol are detailed in [1]. In table SI we report estimated ranges of values for the material parameters involved in the experiment. For the length of the solution path, the range of values is estimated by measuring the shortest and longest direct line through the solution path. The error associated with the measurement of lengths from image analysis is less than 1 mm. The parameter values related to the milk properties, the density $\rho^*$, the dynamic viscosity $\mu^*$ and the nominal surface tension $\gamma_0^*$ (assuming that only natural endogenous surfactant are present at the milk surface), are approximations for the milk-cream mixture used in the experiment [1]. The lowest surface tension of the milk-cream mixture $\gamma_c^*$ when the surfactant concentration is maximum during the experiment is estimated based on a crude approximation for soapy water. The Marangoni force $S^* = \gamma_0^* - \gamma_c^*$, calculated from the values $\gamma_0^*$ and $\gamma_c^*$, gives a sensible range of values. The viscosity of the milk-cream mixture is approximately 30

<footnote n="1">¹Department of Mathematics, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="2">²Andlinger Center for Energy and the Environment, Princeton University, Princeton, New Jersey 08544, USA</footnote>

<footnote n="3">³Institut de Physique, Université de Rennes, F-35000 Rennes, France</footnote>

<footnote n="4">⁴Department of Mechanical Engineering, University of California, Santa Barbara, California 93106, USA</footnote>

<footnote n="5">⁵Universite Claude Bernard Lyon 1, Laboratoire de Mecanique des Fluides et d'Acoustique (LMFA), UMR5509, CNRS, Ecole Centrale de Lyon, INSA Lyon, 69622 Villeurbanne, France</footnote>

<footnote n="1">¹Corresponding author: julien.landel@univ-lyon1.fr</footnote>
<table><caption>Table SI: Table showing the relevant dimensional parameters of the maze experiment conducted by Temprano-Coleto et al. [1]. The range of values for the Marangoni force, below the double lines, has been calculated from other values in the table. The parameter values related to the milk properties are approximations for the milk-cream mixture used in the experiment [1].</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Units</th><th>Estimated range in experiments</th></tr></thead><tbody><tr><td>Length of solution path</td><td>L<sub>M</sub><sup>*</sup></td><td>10<sup>-1</sup> m</td><td>1.36 – 1.84</td></tr><tr><td>Initial height of liquid film</td><td>h<sub>0</sub><sup>*</sup></td><td>10<sup>-3</sup> m</td><td>3.0 – 3.5</td></tr><tr><td>Acceleration due to gravity</td><td>g<sup>*</sup></td><td>m s<sup>-2</sup></td><td>9.8</td></tr><tr><td>Liquid density</td><td>ρ<sup>*</sup></td><td>10<sup>3</sup> kg m<sup>-3</sup></td><td>0.98 – 1.0[2]</td></tr><tr><td>Dynamic viscosity</td><td>μ<sup>*</sup></td><td>10<sup>-2</sup> Pa s</td><td>3 – 4 [3]</td></tr><tr><td>Surface tension of milk-cream mixture</td><td>γ<sub>0</sub><sup>*</sup></td><td>10<sup>-2</sup> N m<sup>-1</sup></td><td>3 – 5[4]</td></tr><tr><td>Surface tension of surfactant-laden interface</td><td>γ<sub>c</sub><sup>*</sup></td><td>10<sup>-2</sup> N m<sup>-1</sup></td><td>2.0 – 2.5</td></tr><tr><td>Marangoni force</td><td>S* = γ<sub>0</sub><sup>*</sup> – γ<sub>c</sub><sup>*</sup></td><td>10<sup>-2</sup> N m<sup>-1</sup></td><td>0.5 – 3.0</td></tr></tbody></table>

<table><caption>Table SII: Table of relevant dimensionless groups for the maze experiments [1], calculated using the dimensional parameter values from table SI.</caption><thead><tr><th>Dimensionless group</th><th>Definition</th><th>Model assumption</th><th>Estimated range in experiments</th></tr></thead><tbody><tr><td>Aspect ratio of liquid film ϵ</td><td>h<sub>0</sub><sup>*</sup>/L<sub>M</sub><sup>*</sup></td><td>≪ 1</td><td>1.7 × 10<sup>-2</sup> – 2.4 × 10<sup>-2</sup></td></tr><tr><td>Inverse capillary number C</td><td>ϵ<sup>2</sup>γ<sub>c</sub><sup>*</sup>/S<sup>*</sup></td><td>≪ 1</td><td>9.7 × 10<sup>-4</sup> – 1.0 × 10<sup>-2</sup></td></tr><tr><td>Bond number G</td><td>ϵ<sup>2</sup>ρ<sup>*</sup>g<sup>*</sup>L<sub>M</sub><sup>*</sup><sup>2</sup>/S<sup>*</sup></td><td>≫ 1</td><td>10 – 40</td></tr><tr><td>Reduced Reynolds number ϵ<sup>2</sup>Re</td><td>ϵ<sup>2</sup>ρ<sup>*</sup>S<sup>*</sup>h<sub>0</sub><sup>*</sup>/μ<sup>*</sup><sup>2</sup></td><td>≪ 1</td><td>1.6 × 10<sup>-3</sup> – 1.6 × 10<sup>-2</sup></td></tr></tbody></table>

to 50 times the viscosity of water. The main effect is to have slowed down the spreading dynamics of the surfactant Marangoni driven flow compared to experiments with water as the liquid medium. The red dye (red food color McCormick®) used as tracer in the experiments of [1] is water based and fully miscible in the milk-cream mixture. The dye is slightly buoyant in the milk-cream mixture, such that it remains at the surface to closely follow the exogenous surfactant spreading in the maze. We can notice that the fronts of the red dye remain fairly sharp in the experiments, as seen on the videos and images. This is due to the absence of external stirring or perturbations in the laminar flow, and slow molecular diffusion processes. The typical diffusion length is estimated to be less than 1 mm during the time of the experiments, based on a conservative estimate of the surface diffusivity of 5 × 10<sup>-9</sup> m<sup>2</sup>s<sup>-1</sup>, which is of the order of the spatial resolution in the images.

In table SII we report estimated values for non-dimensional groups. The model we use to simulate the experiment is based on a one-dimensional (1D) lubrication theory approximation to Stokes flow (see section S2). Inertial forces are dominated by viscous forces when the reduced Reynolds number ϵ<sup>2</sup>Re is small, where ϵ ≪ 1 is the aspect ratio of the main flow path, which is consistent with the values taken from table SI. The equation we use for the transport dynamics assumes that gravity dominates over surface tension and surface tension gradients to maintain a flat surface. Any curvature effects of the film surface associated for instance with the contact angle at the side contact lines are assumed negligible. Accordingly, the Bond number G is found to be large and the inverse capillary number C small, which is consistent with these assumptions.

The dimensionless length L<sub>i</sub> of branch i, for i = 0, 1, 2, . . . 37, is shown in table SIII, where we define one unit as the length of the solution path L<sub>M</sub><sup>*</sup> = 155 mm. These lengths were calculated by counting pixels from the video of the experiment [1], with all dimensional lengths measured with less than 1 mm accuracy. The maze video was taken directly overhead, minimizing optical distortion.
<table><caption>Table SIII: Table of the branch lengths of the maze, where one unit represents the length of the solution path ($L_M^* = 155$ mm). The branch are numbered based on the scheme showed in Fig. S1 (left). The length of the inlet branch $L_0$ is defined in section S2.2 as a function of $\delta$, the ratio of endogenous to exogenous initial surfactant concentration.</caption><tr><td>Branch No. (i)</td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td><td>11</td><td>12</td></tr><tr><td>L<sub>i</sub>/10<sup>-2</sup></td><td>L<sub>0</sub>/10<sup>-2</sup></td><td>9.6</td><td>19</td><td>6.0</td><td>14</td><td>6.5</td><td>6.4</td><td>10</td><td>16</td><td>3.4</td><td>6.3</td><td>9.6</td><td>3.4</td></tr><tr><td>Branch No. (i)</td><td>13</td><td>14</td><td>15</td><td>16</td><td>17</td><td>18</td><td>19</td><td>20</td><td>21</td><td>22</td><td>23</td><td>24</td><td>25</td></tr><tr><td>L<sub>i</sub>/10<sup>-2</sup></td><td>6.5</td><td>13</td><td>3.2</td><td>3.2</td><td>9.9</td><td>6.7</td><td>6.6</td><td>3.2</td><td>6.5</td><td>13</td><td>6.5</td><td>13</td><td>6.5</td></tr><tr><td>Branch No. (i)</td><td>26</td><td>27</td><td>28</td><td>29</td><td>30</td><td>31</td><td>32</td><td>33</td><td>34</td><td>35</td><td>36</td><td>37</td><td></td></tr><tr><td>L<sub>i</sub>/10<sup>-2</sup></td><td>6.5</td><td>3.3</td><td>3.2</td><td>6.3</td><td>6.9</td><td>3.1</td><td>3.2</td><td>9.5</td><td>9.9</td><td>3.3</td><td>7.9</td><td>71</td><td></td></tr></table>

## S2 The surfactant transport equation

The spreading of insoluble surfactant at the free surface of a thin two-dimensional (2D) layer of incompressible viscous fluid under the action of surface-tension gradients can be described by a nonlinear diffusion equation, as we briefly explain. We use 2D Cartesian coordinates with horizontal coordinate $x^*$. A liquid of viscosity $\mu^*$ is confined between a horizontal solid wall and a free surface, at height $h^*(x^*, t^*)$ on which a surfactant with concentration $\Gamma^*(x^*, t^*)$ is present at time $t^*$. The surface tension $\gamma^*$ of the free surface is assumed to diminish linearly as a function of $\Gamma^*$ for $\delta\Gamma_c^* \le \Gamma^* \le \Gamma_c^*$ from a maximum value $\gamma_0^*$ to a minimum value $\gamma_c^*$, with $\delta$ the ratio of endogenous to exogenous initial surfactant concentration such that $\delta\Gamma_c^*$ is the endogenous surfactant concentration initially. We define $S^* = \gamma_0^* - \gamma_c^*$. We consider the evolution of $\Gamma^*$ and $h^*$ over horizontal distances $L_M^*$ that are much larger than the characteristic film height $h_0^*$, and define $\epsilon = h_0^*/L_M^*$. Adopting scales appropriate to lubrication theory, we define $x = x^*/L_M^*$, and $t = \epsilon S^*t^*/(L_M^*\mu^*)$. We define the dimensionless film thickness and surfactant concentration as $h(x,t) = h^*/h_0^*$ and $\Gamma(x,t) = \Gamma^*/\Gamma_c^*$, respectively, and the surface tension as $\gamma = (\gamma^* - \gamma_c^*)/S^* = 1 - \Gamma$. Then, when the reduced Reynolds number is sufficiently small ($\epsilon^2\rho^*S^*h_0^*/\mu^*^2 \ll 1$), the flow is governed by the coupled evolution equations [5]

$$h_t - \frac{1}{2} (h^2 \Gamma_x)_x - \frac{1}{3} G (h^3 h_x)_x + \frac{1}{3} C (h^3 h_{xxx})_x = 0, \quad (S1a)$$

$$\Gamma_t - (\Gamma h \Gamma_x)_x - \frac{1}{2} G (\Gamma h^2 h_x)_x + \frac{1}{2} C (\Gamma h^2 h_{xxx})_x = 0, \quad (S1b)$$

where $C = \epsilon^2\gamma_c^*/S^*$ is the inverse capillary number and $G = \epsilon^2\rho^*g^*L_M^2/S^*$ the Bond number, which are dimensionless parameters representing the strength of surface tension and gravity relative to surface tension gradients, respectively. Surface diffusion of surfactant and surface rheological stresses are neglected. Spreading of surfactant from some initial condition can generate deflections of the free surface, although $h$ never exceeds 2 and cannot reach zero in the absence of disjoining pressure, with gravity suppressing deflections at large times [5]. The limit $G \gg 1$ is consistent with the experiment (Table SII). Expanding $h$ and $\Gamma$ in powers of $1/G$ about the base state $h = 1$ and imposing no-flux conditions at the ends of the domain, (S1), to leading order in $\epsilon$, reduces to the nonlinear diffusion equation [6]

$$\Gamma_t = \frac{1}{4} (\Gamma \Gamma_x)_x, \quad (S2)$$

which can be rearranged as $\Gamma_t = (\Gamma^2)_{xx}/8$. This captures the essential features of self-induced spreading of an insoluble surfactant on a thin film via Marangoni forces.

### S2.1 Surfactant transport along a channel

The factor of $1/4$ in (S2) is appropriate for a thin film, but it takes no account of the lateral no-slip boundary conditions that constrain spreading along a channel in a maze having rectangular cross section with finite width-to-height aspect ratio $W$. We refine the approximation to determine the mobility $f(W)$ of the surfactant as follows, where $f \to 1/4$ in the thin-film limit in which $W \to \infty$. We again assume (i) no out-of-plane deflection of the air-liquid interface and no externally imposed pressure gradient, so that the volume flux along the channel is uniformly zero, and (ii) the channel is long compared to its width and height so that rapid lateral Marangoni spreading quickly eliminates lateral surfactant concentration
gradients [7]. Under (i) and (ii), the axial velocity component $u^*(x^*, y^*, z^*, t^*)$ is driven by a surfactant gradient $-\Gamma_{x^*}^*(x^*, t^*)$ moving the bulk liquid forward and a counter pressure gradient $-p_{x^*}^*(x^*, t)$ driving a return flow in the bulk (arising via incompressibility from the no-interfacial-deflection condition and the fact that the channels are closed at all ends), where $y^*$ and $z^*$ are coordinates spanning the cross-section of the channel in the transverse and vertical directions, respectively. Adopting scalings introduced above, the leading-order axial velocity component satisfies $\nabla_{\perp}^2 u \equiv u_{yy} + u_{zz} = p_x$ across the cross-section of the channel, with $u = 0$ on the bottom and lateral walls of the channel ($z = 0, y = \pm W/2$) and the shear stress condition $u_z = -\Gamma_x$ on the interface at $z = 1$. Linearity allows the field to be decomposed as $u(x, y, z, t) = -p_x(x, t)u^b(y, z) - \Gamma_x(x, t)u^m(y, z)$ where $\nabla_{\perp}^2 u^b = -1$ and $u_z^b = 0$ on $z = 1$, and $\nabla_{\perp}^2 u^m = 0$ and $u_z^m = 1$ on $z = 1$. Averaging $u$ over the rectangular cross-section defines volume fluxes $Q^b(W)$ and $Q^m(W)$, and condition (i) then requires $0 = -Q^b p_x - Q^m \Gamma_x$, which determines $p_x$ in terms of $\Gamma_x$. Finally, the surfactant field is advected by the transversely-averaged surface velocity $\bar{u}_s \equiv W - 1 \int_{-W/2}^{W/2} u(x, y, 1, t) dy$, which can be expressed as $\bar{u}_s = -p_x \bar{u}_s^b - \Gamma_x \bar{u}_s^m$. Eliminating $p_x$ recovers the nonlinear diffusion equation

$$ \Gamma_t = f(W)(\Gamma\Gamma_x)_x, \quad \text{where} \quad f(W) = \bar{u}_s^m - (\bar{u}_s^b Q^m / Q^b). \qquad (S3) $$

In the limit of large $W$, $u^b = z - \frac{1}{2}z^2$, $u^m = z$, $\bar{u}_s^b = \frac{1}{2}$, $\bar{u}_s^m = 1$, $Q^b = 1/3$ and $Q^m = 1/2$ giving $f = \frac{1}{4}$.

For finite $W$, we use separation of variables to obtain the following expressions for $u^b$, $u^m$, $\bar{u}_s^b$, $\bar{u}_s^m$, $Q^b$, $Q^m$ and $f$:

$$ u^b(y, z) = -\frac{1}{2}(z^2 - 2z) - \sum_{n=1}^{\infty} \frac{2 \sin \left[ \left(n - \frac{1}{2}\right) \pi z \right] \left[ \sinh \left[ \left(n - \frac{1}{2}\right) \pi (y + W/2) \right] - \sinh \left[ \left(n - \frac{1}{2}\right) \pi (y - W/2) \right] \right]}{\sinh \left[ \left(n - \frac{1}{2}\right) \pi W \right] \left(n - \frac{1}{2}\right)^3 \pi^3}; \qquad (S4a) $$

$$ u^m(y, z) = - \sum_{n=1}^{\infty} \frac{W(-1)^n \cos\left(\frac{2(n-\frac{1}{2})\pi}{W}y\right) \sinh\left(\frac{2(n-\frac{1}{2})\pi z}{W}\right)}{(n-\frac{1}{2})^2 \pi^2 \cosh\left(\frac{2(n-\frac{1}{2})\pi}{W}\right)}; \qquad (S4b) $$

$$ \bar{u}_s^b = \frac{1}{2} + \sum_{n=1}^{\infty} \frac{2(-1)^n [2 \cosh [(n - \frac{1}{2})\pi W] - 2]}{\sinh [(n - \frac{1}{2})\pi W] (n - \frac{1}{2})^4 \pi^4 W}; \qquad (S4c) $$

$$ \bar{u}_s^m = \sum_{n=1}^{\infty} \frac{W \tanh \left[ 2(n - \frac{1}{2})\pi/W \right]}{(n - \frac{1}{2})^3 \pi^3}; \qquad (S4d) $$

$$ Q^b = 1/3 - \sum_{n=1}^{\infty} \frac{4(\cosh[(n - \frac{1}{2})\pi] - 1)}{\sinh[(n - \frac{1}{2})\pi W](n - \frac{1}{2})^5 \pi^5 W}; \qquad (S4e) $$

$$ Q^m = \sum_{n=1}^{\infty} \frac{W^2(\cosh[2(n - \frac{1}{2})\pi/W] - 1)}{2(n - \frac{1}{2})^4 \pi^4 \cosh[2(n - \frac{1}{2})\pi/W]}; \qquad (S4f) $$

$$ f(W) = \sum_{n=1}^{\infty} \frac{W \tanh \left( \frac{2(n-\frac{1}{2})\pi}{W} \right)}{(n-\frac{1}{2})^3 \pi^3} - \frac{\sum_{n=1}^{\infty} \frac{W^2(\cosh[2(n-\frac{1}{2})\pi/W]-1)}{2(n-\frac{1}{2})^4 \pi^4 \cosh[2(n-\frac{1}{2})\pi/W]} \left( \frac{1}{2} + \sum_{n=1}^{\infty} \frac{2(-1)^n[2\cosh[(n-\frac{1}{2})\pi W]-2]}{\sinh[(n-\frac{1}{2})\pi W](n-\frac{1}{2})^4 \pi^4 W} \right)}{1/3 - \sum_{n=1}^{\infty} \frac{4(\cosh[(n-\frac{1}{2})\pi W]-1)}{\sinh[(n-\frac{1}{2})\pi W](n-\frac{1}{2})^5 \pi^5 W}}. \qquad (S4g) $$

Fig. S1 (left) shows that for $W = 0.94$, which is the approximate aspect ratio of the channel in the experiment, $f \approx 0.15$, which indicates a 40% reduction approximately in the surfactant mobility in comparison to the wide-channel limit.

In addition to gravitational suppression of interfacial deflections, as used to derive (S2) from (S1), surface tension can also have a significant role in suppressing deflections in a channel flow via pinning of contact lines to the channel edges. Any interfacial deflection will have two components of curvature: a short length comparable to the channel width; and a longer axial scale induced by competition between viscous and capillary forces, as represented in (S2a). The former can be expected to dominate over the
latter, as it acts over shorter lengthscales, and it will act alongside gravity in maintaining a near-uniform liquid depth.

Transport of a single surfactant species in a shallow channel (S2) can be expressed dimensionally as $\Gamma_{t*}^* + (u^*\Gamma^*)_{x*} = 0$, where $u^* = -\frac{1}{4}(h_0^*A^*/\mu^*)\Gamma_{x*}^*$. Here $A^* = S^*/\Gamma_c^*$ is the activity of the surfactant (the magnitude of the slope of the assumed-linear relationship between surface tension and surfactant concentration). The equivalent model for a mixture of two surfactant species with concentrations $\Gamma_1^*$ and $\Gamma_2^*$ is $\Gamma_{i,t}^* + (u^*\Gamma_i^*)_{x*} = 0$ for $i = 1, 2$ with $u^* = -\frac{1}{4}(h_0^*/\mu^*)(A_1^*\Gamma_{1,x*}^* + A_2^*\Gamma_{2,x*}^*)$, again assuming a linear equation of state. Defining $\Gamma^* = (A_1^*\Gamma_1^* + A_2^*\Gamma_2^*)/(A_1^* + A_2^*)$, $A^* = A_1^* + A_2^*$ and summing the evolution equations weighted respectively by $A_1^*$ and $A_2^*$ recovers exactly the evolution equation for a single species $\Gamma^*$, equivalent to (S2). While for simplicity we assume that endogenous and exogenous surfactant have identical properties in the present analysis, the model accommodates two non-diffusing species with different activities.

## S2.2 Surfactant transport in a maze

We now extend the 1D model to show that exogenous surfactant added to pre-existing endogenous surfactant can simulate the maze-solving behavior seen in the experiments [1]. This behavior consists of the exogenous surfactant spreading from the end of branch 0 (see numbering scheme in Fig. S1, left) to the start of branch 37, while also not spreading to the end of the lateral branches.

To model the maze-solving behavior of the out-of-equilibrium surfactant field, we solve the simplified time-dependent surfactant transport equation (S3) in a network of connected 1D domains (1DDs). Equation (S3) solves the transport of both the exogenous and endogenous surfactants, as they are assumed to form a single field. Each 1DD represents one of the 36 branches of the maze, where the length of the ith 1DD, $L_i$, is the length of the ith branch. These lengths are given in table III. Additionally, we simulate the square branches 0 and 37 at the inlet and outlet of the maze, respectively, with two more 1DDs of longer lengths. The length of the outlet 1DD, $L_{37}$, is calculated to simulate the surface area of branch 37 in the experiment. Counting the number of pixels $P_{37}^*$ constituting branch 37 in the image in Fig. S1 (left) and, for example $P_1^*$, the number of pixels that constitute branch 1, we define $L_{37} = L_1P_{37}^*/P_1^*$. The length of the inlet 1DD simulating branch 0, $L_0$, is estimated independently to represent the mass of exogenous surfactant $M_{ex}^*$ introduced in this branch at the start of the experiment, as discussed below.

We use subscripts to number the 1DDs, so that for example $\Gamma_0(x_0, t)$ denotes the concentration in the inlet 1DD, for $0 \le x_0 \le L_0$ along the branch, and at time $t$. Each 1DD receives a coordinate direction such that, for example, along 1DD $i$ the coordinate direction is $x_i$, and this coordinate increases in the direction away from the junction nearest to the inlet 1DD. We solve the governing equation (S3) for the surfactant transport equation throughout the 1D network along with the following initial, continuity and boundary conditions, which simulate the conditions in the maze experiment,

$$ \Gamma_0(x_0, 0) = 1 \quad \text{for all } 0 \le x_0 \le L_0, \qquad (S5) $$

$$ \Gamma_1(x_1, 0) = \begin{cases} \frac{1}{2}(1 - \delta) \cos(\pi x_1 / x_f) + \frac{1}{2}(1 + \delta) & 0 \le x_1 \le x_f, \\ \delta & x_f < x_1 \le L_1, \end{cases} \qquad (S6) $$

$$ \Gamma_i(x_i, 0) = \delta \quad \text{for all } 0 \le x_i \le L_i, \quad \text{for } i = 2, 3, \dots, 37. \qquad (S7) $$

$$ \frac{\partial}{\partial x_0} \Gamma_0(0, t) = 0, \quad \frac{\partial}{\partial x_{37}} \Gamma_{37}(L_{37}, t) = 0, \qquad (S8) $$

$$ \frac{\partial}{\partial x_i} \Gamma_i(L_i, t) = 0 \qquad (S9) $$

for all 1DDs $i$ not connecting to a junction at $L_i$,

$$ \Gamma_j(0, t) = \Gamma_i(L_i, t) \qquad (S10) $$

for every 1DD $i$ terminating at a junction, and for all 1DDs $j \in J(i)$ where $J(i)$ is the set of all branches such that $x_j = 0$ locates the same junction as $x_i = L_i$. Finally,

$$ \frac{\partial}{\partial x_i} \Gamma_i(L_i, t) = \sum_{j \in J(i)} \frac{\partial}{\partial x_j} \Gamma_j(0, t), \qquad (S11) $$
for the same set of 1DDs as (S10). The initial, continuity and boundary conditions stated above, (S5)–(S11), closely model the experimental conditions described in section S1. The initial condition in the 1DD representing the inlet branch, $i = 0$, is that the surfactant concentration is uniformly equal to 1 (S5). All the other 1DDs have a uniform surfactant concentration of smaller value $\delta < 1$, which represents the initial ratio of concentrations of endogenous surfactant to exogenous surfactant (S7). The only exception is in the 1DD which connects the inlet 1DD to the rest of the maze, where a short section of length $x_f = 0.032$ smoothly connects the two uniform surfactant concentrations with half a cosine wave (S6). The dynamics are driven by the gradient of the surfactant concentration, so we impose the initial condition to be $C^1$ continuous so that the initial velocity is defined everywhere. The smooth profile in branch 1 has been added to avoid numerical issues at early times. Its influence on the overall transport dynamics is negligible. The boundary condition at unconnected ends of 1DDs in the network is no flux of surfactant, such as at the start of the inlet and end of the outlet branches (S8), and at the end of the peripheral 1DDs (S9). At junctions where multiple 1DDs are connected, the conditions are continuity of surfactant concentration (S10), and continuity of flux of surfactant (S11).

## S2.3 Boundary between the exogenous and endogenous surfactants

As shown in [8], in the absence of surface diffusion, the front of exogenous surfactant, added to a surface with a pre-existing endogenous surfactant concentration, moves like a material element transported at the surface velocity generated through the surfactant-induced Marangoni stress. To simulate the experiment, the evolution of the front location $x_j = \tilde{\Lambda}_j(t)$ (the tilde denotes a simulated solution to avoid confusion with the experimental data $\Lambda_j$) of exogenous surfactant in 1DD $j$, given (S3), satisfies

$$ \frac{d\tilde{\Lambda}_j(t)}{dt} = -f(W) \frac{\partial}{\partial x_j} \Gamma(x_j = \tilde{\Lambda}_j(t), t), \qquad (S12) $$

for all 1DDs $j$ where a front exists. When $\Lambda_j(t) = L_j$, the front ceases to exist in 1DD $j$, and appears in all 1DDs $i \in J(j)$, where $J(j)$ is the set of 1DDs originating at the head junction of 1DD $j$. The initial conditions are that initially a front only exists in 1DD 1, and

$$ \tilde{\Lambda}_1(0) = x_f, \qquad (S13) $$

and we impose that $\tilde{\Lambda}_i(t_b) = 0$ whenever the front in 1DD $i$ appears at time $t = t_b$. The equation (S12) is solved simultaneously with (S3). The maze simulation is completed when $\tilde{\Lambda}_{36} = L_{36}$.

## S2.4 Key non-dimensional parameters in the model

The model has three unknown parameters to be determined empirically. The first is $\delta$, the ratio of initial endogenous surfactant concentration to the reference exogenous surfactant concentration. The second is $\tau_{end}$, which is the ratio of simulation completion time to experiment completion time $t_{ref}$. We use $\tau_{end}$ to improve the comparison to the experiment of the dynamic behavior of exogenous surfactant along the solution path (see section S6 for further details), but $\tau_{end}$ does not affect the model dynamics significantly. The third a priori unknown parameter is $L_0$, the length of the inlet 1DD. As detailed below, the length $L_0$ is used to simulate the mass of exogenous surfactant $M_{ex}^*$ in the inlet branch 0.

Examining the video of the experiment [1], we observe complex dynamics occurring in the inlet branch 0 at early times. The exogenous surfactant is not added instantaneously, nor is it added to the entire area initially occupied by the dye. Part of branch 0 is unoccupied by the dye, and remains unoccupied throughout the duration of the experiment. Moreover, our quantitative analysis of how the dye front moves at early times from branch 0 to branch 1 shows a very fast dynamics, with a short time scale characteristically different from the dynamics throughout the internal branches of the maze. We assume that this rapid dynamics is associated with the inlet branch being a large square, rather than a narrow channel, such that the dynamics is mostly 2D, rather than 1D as assumed in our model. To avoid adding unnecessary complexity to our model, we do not model this early time dynamics at the inlet branch as it has limited impact on the long-time dynamics throughout the maze. We find that the overall maze dynamics is well-captured by a 1D inlet branch, whose length is slightly different from the geometric
value. Nevertheless, we still assume that the inlet branch is occupied uniformly by exogenous surfactant at $t = 0$. Thus, the inlet branch length is chosen to represent the mass of exogenous surfactant deposited at $t = 0$, rather than its geometric size. The actual size of branch 0 is not crucial to explaining the key behavior of the experiment (unlike branch 37). The role of branch 0 is to act as a (finite) reservoir of exogenous surfactant.

Owing to conservation of mass, we can estimate the exogenous mass of surfactant $M_{ex}^*$, which is used to determine $L_0$, by examining the video of the experiment at late times. At late times, the surface tension throughout the maze approaches a constant value. Under the modelling assumption that the material parameters of the exogenous surfactant are the same as the endogenous surfactant, this uniform surface tension means that the concentrations of both species of surfactant also evolve to a constant value, which we call $\Gamma_u^*$. Furthermore, under the assumptions that the surfactant species do not mix and that the red dye closely follows the exogenous-endogenous interface, the endogenous surfactant occupies the area of the maze which is white, and the exogenous surfactant occupies the area of the maze which is red. Calling the dimensional masses of exogenous and endogenous surfactant $M_{ex}^*$ and $M_{en}^*$ respectively, we make the approximations that

$$M_{en}^* \approx \int_{A_w^{\infty*}} \Gamma_u^* dS^* = \Gamma_u^* A_w^{\infty*}, \quad \text{and} \quad M_{ex}^* \approx \int_{A_r^{\infty*}} \Gamma_u^* dS^* = \Gamma_u^* A_r^{\infty*}, \qquad (S14)$$

where $A_w^{\infty*}$ and $A_r^{\infty*}$ are the surface areas of the maze not occupied by the dye, and occupied by the dye respectively when the system reaches a steady state at the end of the video. The concentration $\Gamma_u^*$ cancels when calculating the ratio of masses, which we exploit to impose the condition that the ratio of dimensionless masses in our model to be the same as the ratio of masses in the experiment. The non-dimensional exogenous surfactant mass relates to the dimensional mass by $M_{ex}^* = M_{ex} \Gamma_c^* L_M^* W^*$, where $W^*$ is the dimensional width of the maze branches, with a similar identity holding for the endogenous surfactant mass $M_{en}$. Thus, $M_{ex}/M_{en} = A_r^{\infty*}/A_w^{\infty*}$ provides an identity for $M_{ex}$ in terms of $\delta$

$$M_{ex} = M_{en} \frac{A_r^{\infty*}}{A_w^{\infty*}} = \delta \left[ \sum_{i=1}^{N} L_i - x_f \right] \frac{A_r^{\infty*}}{A_w^{\infty*}}, \qquad (S15)$$

where the quantity $[\sum_{i=1}^N L_i - x_f]$ is the total length of all the 1DDs in the model not initially occupied by the exogenous surfactant. From the initial conditions, the exogenous surfactant mass satisfies

$$M_{ex} = \int_{0}^{L_0} \Gamma_0(x_0, 0)dx_0 + \int_{0}^{x_f} \Gamma_1(x_1, 0)dx_1 \\ = \int_{0}^{L_0} 1 dx_0 + \int_{0}^{x_f} \left[ \frac{1}{2}(1-\delta)\cos(\pi x_1/x_f) + \frac{1}{2}(1+\delta) \right] dx_1 = L_0 + \frac{1+\delta}{2}x_f. \quad (S16)$$

Using (S15) and (S16), the dimensionless length of the inlet 1DD ($i = 0$) can be estimated as

$$L_0 = \delta \left[ \sum_{i=1}^{N} L_i - x_f \right] \frac{A_r^{\infty*}}{A_w^{\infty*}} - \frac{1+\delta}{2} x_f, \qquad (S17)$$

where $A_r^{\infty*}$ and $A_w^{\infty*}$, are measured from the video, the dimensionless lengths $L_i$, $i = 1, \dots, N$ are taken from table SIII. We note that the precise value of $x_f$ (chosen as 0.032) has negligible impact on the maze dynamics. The remaining dimensionless parameters $\delta$ and $\tau_{end}$ are determined through comparison with the experiment (see section S6).

# S3 Numerical methods

## S3.1 Spatial discretisation

We model the transport dynamics in the maze by the time-dependent equation (S3) along a connected network of 1DDs, subject to the initial, continuity and boundary conditions (S5)–(S11) presented in
section S2.2. The endogenous-exogenous surfactant interface, assumed to be the front of the red dye, is calculated through (S12) under the initial condition (S13). We solve a discretised approximation of the continuum formulation of the problem on a network having the same topology as the maze, exploiting the methodology of mimetic finite differences. We represent the maze as a directed graph. Every junction and every unconnected end of a 1DD is represented by a vertex, and within each 1DD we add vertices to refine the discretisation with oriented edges between them. The number of vertices $N_i$ within the ith 1DD is chosen such that these vertices are equally spaced at a distance $\Delta x_i = L_i/(N_i + 1)$, where $\Delta x_i$ is as close as possible to a chosen global discretisation value $\Delta x$ such that it gives an integer $N_i$. This is accomplished by

$$ \Delta x_i = \frac{L_i}{\text{round}(L_i/\Delta x)}. \qquad (S18) $$

As a convention, we choose the orientation of every edge to point away from the vertex nearest to the inlet. Following the orientation of each edge, its vertex nearest to the inlet is designated as tail vertex, and its other vertex is designated as head vertex. The topology of the directed graph is encapsulated by its signed incidence matrix which is defined by

$$ A_{ij} = \begin{cases} +1 & \text{Edge } i \sim \text{Vertex } j \mid \text{Edge } i \text{ points towards Vertex } j, \\ -1 & \text{Edge } i \sim \text{Vertex } j \mid \text{Edge } i \text{ points away from Vertex } j, \\ 0 & \text{Edge } i \not\sim \text{Vertex } j, \end{cases} \qquad (S19) $$

with $i$ and $j$ two integers such that $1 \le i \le m$, $1 \le j \le n$ for a graph with $m$ edges and $n$ vertices, and where $\sim$ and $\not\sim$ mean ‘connected to’ and ‘not connected to’, respectively. Each row of $A$ corresponds to an edge of the graph and has exactly two non-zero entries: +1 corresponding to the head vertex of the edge, and -1 corresponding to the tail vertex of the edge. We designate $A_u$ the unsigned incidence matrix, which is the defined via $A_{u(ij)} = |A_{ij}|$. The labelling of edges and vertices is arbitrary. A subset of adjacent vertices is represented by a vector (or chain) $\mathbf{D}$ of length $n$ whose components are 1 for every vertex in the subset and zero otherwise. Non-zero components of the chain $\mathbf{AD}$ identify the oriented edges bounding the domain represented by $\mathbf{D}$.

To define physically meaningful operators on the graph, we include metric information describing the physical system represented by the graph. We define the edge length metric tensor as the diagonal matrix of size $m \times m$

$$ \mathbf{L} = \text{diag}(\Delta x_0, \dots, \Delta x_0, \dots, \Delta x_1, \dots, \Delta x_1, \dots \Delta x_2, \dots), \qquad (S20) $$

where $\Delta x_0, \Delta x_1$, etc. are the lengths of the edges of the graph, each obtained from (S18), ordered the same as the rows of $A$. We define the associated vertex length metric tensor as

$$ \mathbf{V} = \text{diag} \left( \frac{1}{2} \mathbf{A}_u^T \mathbf{L} \mathbf{1}_e \right), \qquad (S21) $$

where $\mathbf{1}_e$ is the $m$-vector (or chain) $(1, 1, \dots, 1)$ identifying all the edges of the graph. The quantity $V_{jj}$ is half the sum of edge lengths of all edges connected to vertex $j$. We define $\mathbf{\Gamma}_v$ as the $n$-vector (co-chain) of surfactant concentrations defined on the $n$ vertices of the graph, and $\mathbf{\Gamma}_v^2$ as the component-wise square of this vector. The incidence matrix can be used as an operator which can act on $\mathbf{\Gamma}_v^2$. For instance, $A\mathbf{\Gamma}_v^2$ is a differencing operation on $\mathbf{\Gamma}_v^2$ that returns a co-chain of variables defined on the edges of the graph, where the value on each edge is the difference of the square of the concentration on the bounding vertices of the edge. Thus, using (S20) we can define a gradient operator on $\mathbf{\Gamma}_v^2$ as $L^{-1}A\mathbf{\Gamma}_v^2$. This is equivalent to a second-order-accurate central-difference approximation of the gradient of $\mathbf{\Gamma}_v^2$ defined at the midpoint of each edge. We can therefore approximate the continuous flux $-f(W)(\mathbf{\Gamma}_v^2)_x/2$ in the rearranged version of (S3) by the discrete expression

$$ \mathbf{q}_e = -\frac{f(W)}{2} L^{-1} A \mathbf{\Gamma}_v^2. \qquad (S22) $$

The flux $\mathbf{q}_e$ is defined on the edges of the graph. The vector $\mathbf{V}\mathbf{\Gamma}_v$ can be interpreted as the mass of surfactant associated with each vertex. Then $\langle \mathbf{D}, \mathbf{\Gamma}_v \rangle = \mathbf{D}^T \mathbf{V}\mathbf{\Gamma}_v$ can be interpreted as the mass of surfactant in the simply connected domain represented by the chain $\mathbf{D}$. Mass conservation at each individual vertex (in the absence of sources and sinks) is expressed as

$$ \frac{d}{dt} \mathbf{V}\mathbf{\Gamma}_v = \mathbf{A}^T \mathbf{q}_e. \qquad (S23) $$
Here, $A^T q_e$ returns a co-chain of variables defined on the vertices of the graph giving the net amount of flux entering each vertex. Indeed, each row of $A^T$ corresponds to a vertex of the graph, where the components in the row are either +1 or -1 for each edge (corresponding to its column) pointing towards or pointing away from that vertex, respectively. Equivalently, from (S23), $-V^{-1}A^T$ is the analogue of the divergence operator applied to all the vertices. Therefore, for a fixed metric V in time, the discrete analogue of the surfactant transport equation (S3) in a maze of a given topology is

$$ \frac{d\Gamma_v}{dt} = V^{-1}A^T q_e = -\frac{f(W)}{2}V^{-1}A^T L^{-1}A\Gamma_v^2. \qquad (S24) $$

Left-multiplying (S24) by $D^T V$ gives

$$ \frac{d}{dt}\langle D, \Gamma_v \rangle = (AD)^T q_e. \qquad (S25) $$

We recognise $(AD)^T q_e$ as the net flux entering the boundary of the domain represented by $D$. This provides a discrete analogue of the divergence theorem on the graph. When implemented numerically, this formulation conserves mass to machine precision.

## S3.2 Temporal discretisation

We solve (S24) using a semi-implicit time-stepping finite-difference scheme. At the $(k+1)$th time step, we approximate $\Gamma_v^2$ as $Q\Gamma_v^{(k+1)}$ where $Q = \text{diag}(\Gamma_v^{(k)})$. Therefore, (S24) can be approximated by

$$ \frac{\Gamma_v^{(k+1)} - \Gamma_v^{(k)}}{\Delta t} = -\frac{f(W)}{2}V^{-1}A^T L^{-1}AQ\Gamma_v^{(k+1)}. \qquad (S26) $$

The solution at the $(k+1)$th time step is then obtained by solving the linear system

$$ \left( I + \frac{f(W)}{2} \Delta t V^{-1} A^T L^{-1} A Q \right) \Gamma_v^{(k+1)} = \Gamma_v^{(k)}, \qquad (S27) $$

where $I$ is the $n \times n$ identity matrix. As discussed, the boundary conditions (S9) to (S11) are implemented automatically within the formulation (S27). The initial conditions are imposed by projecting (S5)-(S7) onto the vertices of the graph, whereby we obtain the initial concentration vector $\Gamma_v^{(0)}$. This is done using the same order for the components of $\Gamma_v^{(0)}$ as we have ordered the columns of $A$ (this ordering is arbitrary but must be consistent). The numerical scheme is found to converge with second-order spatial accuracy following $\Delta x^2$, with $\Delta x \to 0$. This is consistent with Brio et al. [9], who found that the rows of (S24) are a second-order accurate approximation of the Laplacian at each vertex of a graph.

Equation (S12) is solved simultaneously with (S27) to keep track of the fronts of exogenous surfactant. We define $\Lambda^{(k)}$ as the vector of locations of exogenous surfactant fronts. This vector has the same number of components as the number of edges in the graph. To impose the initial conditions (S13), every component of $\Lambda^{(0)}$ is zero except the $\lfloor x_f/\Delta x_1 \rfloor$ components which represent the first $\lfloor x_f/\Delta x_1 \rfloor$ edges in 1DD 1, where $\lfloor a \rfloor$ means the floor of $a$, returning the largest integer smaller than $a$. These components of $\Lambda^{(0)}$ are set to $\Delta x_1$. We set the $\lfloor x_f/\Delta x_1 \rfloor + 1$ component in 1DD 1 equal to $x_f - \Delta x_1 \lfloor x_f/\Delta x_1 \rfloor$.

We also equip our scheme with an indicator vector $I^{(k)}$ which has the same number of components as there are edges. Each component is either 1 if a front of the exogenous surfactant is present in that edge or 0 if not. We set the $\lfloor x_f/\Delta x_1 \rfloor + 1$ component of $I^{(0)}$ in 1DD 1 to 1, and every other component to 0. The numerical scheme which solves equation (S12), relating time-step $k$ to $k-1$, is

$$ \Lambda^{(k)} = \Lambda^{(k-1)} - f(W) I^{(k)} \circ L^{-1} A \Gamma_v^{(k-1)}, \qquad (S28) $$

where $\circ$ is the component-wise product of vectors.

The $j$th component of $I^{(k-1)} \circ \Lambda^{(k)}$, if it is non-zero, corresponds to an exogenous surfactant front within the $j$th edge of the graph. After time-step $k$ we check whether every non-zero component of $I^{(k-1)} \circ \Lambda^{(k)}$ is within some small tolerance $\varepsilon$ of $\Delta x_i$, where $i$ is the index of the 1DD where the edge corresponding to the component is located. If $|\Delta x_i - \Lambda_j^{(k)}| < \varepsilon$ for this component $j$, we 'hand over' this
surfactant front to the connecting edges by setting the component of $I^{(k)}$ for that edge to 0, and the forward connecting edges to 1. The other components of $I^{(k)}$ remain the same as $I^{(k-1)}$. To reduce error we also set the components corresponding to the forward connecting edges of $\Lambda^{(k)}$ to $\Delta x_i - \Lambda_j^{(k)}$. We also take into account receding effects by testing whether or not the fronts have receded to the start of any edge, in which case we pass the front location back to the previous edge following a similar algorithm. The location of the exogenous surfactant front within any 1DD at time-step $k$ is given by the sum of the components of $\Lambda^{(k)}$ corresponding to the edges representing that 1DD. If we call $J(i)$ the set of indices corresponding to the edges discretising 1DD $i$, and we assume a uniform time step $\Delta t$,

$$ \tilde{\Lambda}_i(k\Delta t) = \sum_{j \in J(i)} \Lambda^{(k)}(j). \qquad (S29) $$

# S4 Linearisation and modal decomposition of the transport model

## S4.1 Linear modes

We can approximate the nonlinear diffusion equation (S3) in the limit of small gradients of surfactant concentration as a linear diffusion equation. If we assume that the surfactant concentration has the form $\Gamma = \bar{\Gamma} + \hat{\Gamma}(x, t)$ where $\hat{\Gamma} \ll \bar{\Gamma}$, the nonlinear diffusion equation projected onto the discrete representation of the maze (S24) becomes, to leading order,

$$ \frac{d\hat{\Gamma}}{dt} = -\kappa V^{-1} A^T L^{-1} A \hat{\Gamma}, \qquad (S30) $$

where $\kappa = f(W)\bar{\Gamma}$, and we identify the Laplacian operator

$$ \kappa V^{-1} A^T L^{-1} A. \qquad (S31) $$

If we define $U = V^{1/2}\hat{\Gamma}$ then (S30) becomes

$$ U_t = -\kappa V^{-1/2} A^T L^{-1} A V^{-1/2} U. \qquad (S32) $$

The operator $\kappa V^{-1/2} A^T L^{-1} A V^{-1/2}$ is symmetric, and therefore has orthogonal eigenvectors, so that $U_i^T U_j = \delta_{ij}$ for normalised eigenvectors $U_i$ and $U_j$, with distinct eigenvalues $\lambda_i$ and $\lambda_j$ (as the maze is not symmetric, all eigenvalues are distinct), and where $\delta_{ij}$ is the Kronecker delta. Thus, we can express $\hat{\Gamma}$ in terms of orthonormal eigenvectors $\phi_i = V^{-1/2} U_i$ of (S31) that satisfy

$$ \phi_i^T V \phi_j = \langle \phi_i, \phi_j \rangle = \delta_{ij}. \qquad (S33) $$

This defines the natural inner product for the Laplacian operator (S31), which is the operator governing the linearised dynamics of the transport in the maze. The Laplacian (S31) is positive semi-definite, with a single zero eigenvalue corresponding to a constant eigenvector representing steady state concentrations as $t \to \infty$.

We can use a truncated sum of the first $N$ eigenvectors to construct an approximate solution for the vector of surfactant concentrations defined on the vertices of the graph,

$$ \Gamma(t) = \sum_{i=0}^{N} A_i \phi_i e^{-\lambda_i t}, \qquad (S34) $$

where the coefficients $A_i$ are the amplitudes of the eigenvectors $\phi_i$, and $\lambda_i$ are the corresponding eigenvalues of (S31). Utilising the orthogonality of these eigenvectors with respect to this inner product we can project an initial concentration profile $\Gamma(0)$ onto this basis. The amplitudes for the $i$th mode will be given by

$$ A_i = \langle \phi_i, \Gamma(0) \rangle. \qquad (S35) $$

We find the modes look similar to cosine waves through the solution path of the maze. The gradients of the cosine waves have discontinuities at the locations of junctions, which is to be expected due to the internal boundary condition (S11). The first 20 eigenvalues are plotted in Fig. S2(a).
<figure><img src="image_7.png" /></figure>

Figure S2: (a) The first 20 eigenvalues of the Laplacian (S31) (circles). In magenta is the approximation $\pi^2n^2/L_W$, where $L_W = 3.36$ is the nondimensional length of the solution path. In blue is the approximation $\pi^2n^2/L_T$, where $L_T = 5.07$ is the sum of the lengths of all the branches in the maze. The first expression approximates the dominant eigenvalues well, but the higher order modes are much better approximated by the second expression. (b) Results of the parameter optimisation. The blue line shows $E_{sp}(\delta)/\min E_{sp}(\delta)$ calculated from (S41), the normalised square error between the front location in the experiment and simulation along the solution path. The red line shows $E_{pb}(\delta)/\min E_{pb}(\delta)$ calculated from (S42), the normalised square error between experiment and simulation for the location of the fronts in the lateral branches. Both error calculations give minima which are close to each other. The vertical dashed line represents the minimum of the sum of the two normalised data sets, which is the optimal value we use for $\delta$ at 0.15.

## S4.2 Approximation of transport with the linear modes

We can capture qualitatively the main behaviours observed in the experiment (maze solving, receding in dead-end sections) with only three eigenmodes. We are however unable to capture the initial conditions even approximately with only three modes. To replicate the key behaviours we therefore need to adjust $\delta$ and $L_0$. Using a fitting parameter $\Phi$, we replace $L_0$ with $\Phi L_0$. We perform a parameter sweep over $1 \le \Phi \le 40$ and $0.01 \le \delta \le 0.1$ to find a minimum error between experiment and simulation using the same error quantification as in section S6. We find a minimum of the combined data sets at $\delta = 0.03$ and $\Phi = 25$.

## S5 Experimental data

Experimental data to compare with the model predictions were obtained through image analysis of the video of the experiment [1]. The video of the maze experiment was converted into a time series of $N + 1$ JPEG images, with each image taken at regular times $T_i$, where $i = 0, 1, 2, \dots, N$. The time $T_0$ is the start time and the corresponding image is given in Fig. S1 (left), and $T_N$ is the completion time $t_{ref}$ (chosen when the dye front just enters the outlet branch 37). We imported these images into MATLAB as matrices of dimension $h \times v \times 3$ where $h$ and $v$ are the number of rows and columns of pixels constituting the image, and for each pixel, the three numbers in the third index take integer values from 0 to 256, corresponding to the 8-bit light intensity value of the red, green and blue components at this pixel. We normalise these data by dividing by 256, such that ideal white is created by the 3-tuple [1, 1, 1], and ideal red by [1, 0, 0].

In each image of the time series, we average the pixel values in the transverse direction across the branches, for every branch in the maze, in order to obtain a 1D network of pixels representing the maze. For instance, let us say a horizontal (east-west orientation) section of a branch in a given image has $P_w$ pixels in the east-west direction and $P_l$ pixels in length in the north-south direction, then its image is a
$P_w \times P_l \times 3$ matrix P. We transform this matrix into a $P_l \times 3$ matrix p defined by

$$ \mathbf{p}(j, k) = \frac{1}{P_w} \sum_{i=1}^{P_w} \mathbf{P}(i, j, k). \qquad (S36) $$

Likewise, for a vertical section of a branch, its image is a $P_l \times P_w \times 3$ matrix, which we transform by

$$ \mathbf{p}(i, k) = \frac{1}{P_w} \sum_{j=1}^{P_w} \mathbf{P}(i, j, k). \qquad (S37) $$

At corners we neglect the transverse averaging, and instead trace a line of pixels diagonally from the middle of one branch section incident to the corner, indexing this line with the first index of a matrix $c(i, k)$, where $k$ are once again the 3 RGB components. Junctions are taken to belong wholly to the branch nearest the inlet. We then concatenate the matrices of the form (S36) and (S37), and corner sections. For example, the $i$th branch of a maze at time $T_j$ consists of first a horizontal section, which we turn into a matrix of form (S36), which we call $\mathbf{p}_{h1}$, followed by a corner matrix $\mathbf{c}_1$, followed by a vertical section of the form (S37), which we call $\mathbf{p}_{v1}$, followed by a corner $\mathbf{c}_2$, and so on; the matrix representing the branch is then

$$ \mathbf{B}_i^{T_j} = [\mathbf{p}_{h1} \ \mathbf{c}_1 \ \mathbf{p}_{v1} \ \mathbf{c}_2 \ \mathbf{p}_{h2} \ \mathbf{c}_3 \ \mathbf{p}_{v2} \ \dots]^\top. \qquad (S38) $$

The image of all 36 branches in the maze are thus transformed into matrices of the form (S38), where the first index represents the pixel location along the branch with its index increasing in the direction away from the junction nearest to the inlet, and the second index gives the transverse averaged pixel intensities at that location. This operation was repeated for each image in the time series from $T_0$ to $T_N$.

For each image in the time series, we compute upper and lower bounds for the locations of the fronts of exogenous surfactant in each branch, as a measure of the spread of the front. We define the upper bound $\Lambda^u$ as the pixel with the smallest index which has both green and blue components above 180, and the lower bound $\Lambda^l$ as the pixel with the minimum index where both green and blue were above 120. The difference in these values represent approximately a quartile of colour intensity. For larger values than 180, or smaller values than 120, we find a smooth approximation of the front location is unobtainable due to noise. Indeed, in both the video images, and transverse averaged 1D network, the pixel components creating the white of the milk and the red of the dye are noisy, such that any given pixel usually has RGB components in the range $[0.9, 0.9, 0.9] \le I_{milk} \le [1, 1, 1]$ and $[0.9, 0, 0] \le I_{red} \le [1, 0.1, 0.1]$. If the number of rows in matrix $\mathbf{B}_i$ is $P_B$, then these locations are

$$ \Lambda_i^u(T_j) = \max_{D:(B_i^{T_j}(D,2)>180)\cup(B_i^{T_j}(D,3)>180)} \frac{W}{P_B} L_i, \quad \text{and} \quad \Lambda_i^l(T_j) = \max_{D:(B_i^{T_j}(D,2)>120)\cup(B_i^{T_j}(D,3)>120)} \frac{W}{P_B} L_i. \quad (S39) $$

The difference $|\Lambda_i^l(T_j) - \Lambda_i^u(T_j)|$ gives us an experimental error for the exogenous front location. We define $\Lambda_i(T_j)$ as the mean of the upper and lower bounds

$$ \Lambda_i(T_j) = \frac{1}{2} (\Lambda_i^u(T_j) + \Lambda_i^l(T_j)). \qquad (S40) $$

The initial location of the unique exogenous front as used in (S7), (S13), and (S17) is given by $\Lambda_1(0) = x_f = 0.032$.

# S6 Comparison of the model predictions with the experiment

We compare the predictions for the time evolution of the fronts of the exogenous surfactant $\tilde{\Lambda}_i$, obtained from solving the MFD formulation (S27) and (S28), with the experimental data $\Lambda_i(T_j)$. We analyze two distinct behaviours. The first is the front location of exogenous surfactant along the solution path, and the second is the late time behavior of the front locations in the lateral branches. Hence, we define two separate error quantification procedures for the two behaviours and perform a parameter sweep over $\delta$ to find the values which minimise both errors.
## S6.1 Minimisation of the error along the solution path

We define $E_{sp}(\delta)$ as the difference between simulation and experiment for the front location along the solution path, defined by branches 1, 3, 5, 7, 9, 12, 14, 18, 22, 30, 32, 34 and 36 (Fig. S1, left). We choose a set of $M$ equally spaced points along the solution path, and for each point $\chi_i$ find the time taken $\tilde{\tau}_i$ for the exogenous surfactant front to arrive in a simulation run with a given $\delta$. Then, we use the data set obtained from the experiment of exogenous front locations to find an analogous data set from the experiment. The data set found in S5 maps points equally spaced in time $T_j$ to the spatial location of the exogenous front in branch $i$, $\Lambda_i(T_j)$. Using MATLAB®'s 'interp1' function we interpolate between these data points to find the data set mapping the same set of spatial points $\chi_i$ along the solution path to the time of arrival at those points $\tau_i$.

We obtain two normalised, comparable data sets from simulation and experiment $\tilde{\tau}_i/\tilde{\tau}_M$ and $\tau_i/\tau_M$, respectively, for $i = 1, 2, \dots, M$. We impose that experiment and simulation start at the same time, but we allow the finish time to vary by scaling the simulation data set by a parameter $0.5 \le \tau_{end} \le 1.5$. This secondary parameter provides a slightly better fit, but has negligible impact on the fundamental dynamics underpinning the model. To obtain a positive definite error we compute

$$E_{sp}(\delta) = \min_{0.5 \le \tau_{end} \le 1.5} \sum_{i=1}^{M} \left( \frac{\tau_i}{\tau_M} - \tau_{end} \frac{\tilde{\tau}_i(\delta)}{\tilde{\tau}_M(\delta)} \right)^2. \qquad (S41)$$

We compute $E_{sp}(\delta)$ for a range of values $0 < \delta < 1$. The normalised error $E_{sp}(\delta)$ is shown in Fig. S2(b) (blue). The minimum of $E_{sp}(\delta)$ is found at $\delta \approx 0.075$.

## S6.2 Minimisation of the error in the lateral branches

We quantify the error $E_{pb}(\delta)$ between experiment and simulation for the front location in the lateral branches. We consider branches, $I = \{2, 4, 6, 8, 10, 11, 13, 15, 19, 23, 24, 31, 33\}$ (Fig. S1, left). We neglect branch 35 where the experimental front does not penetrate. We record the location of the front inside each of the lateral branches in the set $I$ at the set of time points $1 \le j \le M$ described in subsection S6.1, in both simulation $\tilde{\Lambda}_i(\tilde{\tau}_j)$ and experiment $\Lambda_i(\tau_j)$. We compute the error as

$$E_{pb}(\delta) = \sum_{i \in I} \sum_{j=1}^{M} \frac{(\Lambda_i(\tau_j) - \tilde{\Lambda}_i(\tilde{\tau}_j))^2}{L_i^2}. \qquad (S42)$$

The normalised error $E_{pb}/\min E_{pb}$ is plotted in Fig. S2(b) (red) for a range of $\delta$. The minimum is at $\delta \approx 0.195$.

## S6.3 Optimal value for the endogenous concentration

The errors $E_{sp}(\delta)$ and $E_{pb}(\delta)$, defined in (S41) and (S42) respectively, give well-defined minima with different optimal values for $\delta$. We added these normalised data sets together, and found a combined minimum at $\delta \approx 0.15$ with $\tau_{end} \approx 1.24$, which are the values used in the simulation results shown in this study. We note that all values $0.02 \le \delta \le 0.175$ would give a reasonable fit between the model predictions and the experiments along the solution path, and all values $0.125 \le \delta \le 0.21$ would give a reasonable fit of the dynamics in the lateral branches.

## S7 Dynamics in lateral branches

In addition to the results presented in Fig. 3(a) (inset) in the main paper, we show the front dynamics in the first 9 lateral branches in the experiment and in the simulation using the optimal value of $\delta$ calculated in section S6, Fig. S3. We can observe the same receding effect of the front locations in the simulation and in the experiment. This is explained physically by the finite mass of surfactant available from the inlet branch (in contrast to a fixed concentration at the inlet) spreading through the solution path causing non-monotonic behavior of the surfactant concentration closer to the inlet. The concentration increases
<figure><img src="image_8.png" /><figcaption>Figure S3: Evolution of the front in time in the first 9 lateral branches for both the optimal prediction from the nonlinear model (black curve, $\delta = 0.15$) and the experimental data (red curve) (S40). The red shaded area shows the experimental error (S39) for the front location. We observe receding behavior in both simulation and experiment as the non-monotonic behavior of surfactant concentration in the early branches of the maze cause a Marangoni flow in the opposite direction in these branches. The receding in the experiment seems to occur over shorter timescales than we achieve in the simulation. We note that the 3-mode linear model (see section S4) has not been plotted in this figure, as we find that for most of these lateral branches the linear model produces no or little flow. This explains why, overall, the full model yields a more quantitative fit to the experimental data throughout all the branches of the maze, based on the global fitting criterion explained in section S6.</figcaption></figure>

substantially initially in the first lateral branches of the solution path, but then decreases at later times. This late time behavior causes a receding Marangoni flow in the opposite direction in branches 2, 4, 6, 8, 10, 11, 13, and 15.

## S8 Two-dimensional dynamics

We use COMSOL to simulate the spreading of surfactant in some of the individual branches, and use the particle-tracing feature to track the dynamics of the front of exogenous surfactant. We do this for the section which consists of part of branches 7 and 14, and branches 8, 9, 10, 11, 12 and 13 (Fig. S1, left). In this section, we solve the 2D equation

$$ \frac{\partial \Gamma}{\partial t} = f(W) \nabla \cdot (\Gamma \nabla \Gamma), \qquad (S43) $$

implementing no flux boundary conditions on all boundaries except at the inlet ($\Gamma = 1$) and outlet $\Gamma = \delta$. The initial conditions are $\Gamma = (1 - \delta) \cos(\pi x/2)/2 + (1 + \delta)/2$ for $0 \le x_7 \le 2$, and $\Gamma = \delta$ everywhere else.
<figure><img src="image_9.png" /><figcaption>Figure S4: Sequence of images showing an experiment to solve a liquid maze with surfactant using different experimental conditions than in the original experiments of [1] (see §S9.1 and experiment #1 in table SIV). (a) A pipette is used to introduce 20 µl of isopropyl alcohol (undyed). The maze was previously filled with glycerol, and seeded with pliolite particles (nominal diameter around 100 µm), which float at the surface. (b)-(d) show flow evolution. Tracks for each particle, starting from $t = 0$, are shown in color. (Length scale: the sides of the large entrance square inlet and outlet measure approximately 20 mm.)</figcaption></figure>

## S9 Additional experiments

Our main conclusion is that the maze solving behavior is primarily due to the interaction between the exogenous and endogenous surfactants. We sought to demonstrate that the maze-solving behavior is robust, and does not rely on the specific liquids and surfactant used in the original experiments of [1], which are also shown in our Letter. The original experiments of [1] used a mixture of milk and cream as the liquid medium, and the exogenous surfactant consisted of an aqueous soap solution. The flow was visualized using an aqueous red dye solution.

Here we show that a similar maze-solving behavior can be observed as long as the following requirements are met:

* The exogenous surfactant solution, as well as any passive tracers, must be less dense than the background liquid. This ensures that the surfactant and tracer move along the surface;

* The background liquid must be of sufficiently high viscosity to keep the flow essentially laminar; and
<figure><img src="image_10.png" /><figcaption>Figure S5: Sequence of images showing another experiment to solve a liquid maze with surfactant using different experimental conditions than in the original experiments of [1] (see §S9.2 and experiment #2 in table SIV). (a) A pipette is used to introduce a mixture of 50% isopropyl alcohol (2-propanol) and 50% dye (by mass). The background liquid is a solution of 10% 2-propanol and 90% glycerol. The amount added is 1 ml. (b)-(d) show flow evolution. (Length scale: the sides of the large entrance square inlet and outlet measure approximately 20 mm.)</figcaption></figure>

* The exogenous surfactant must be mild enough (as measured by the slope of the surface tension change with concentration), such that the associated Marangoni flow does not immediately lead to complete coverage of the interface.

The experiments of [1] leveraged the fact that soap in milk or cream behaves as a much milder surfactant than soap in water (see for example [10] for data and discussion on the interaction between milk casein and Sodium Dodecyl Sulfate from soap). Because we now wish to use a background fluid other than milk, we choose isopropyl alcohol (i.e. 2-propanol), which acts as a mild surfactant in water, as measured by the ratio between adsorption and desorption coefficients [11]. To maximize viscosity and density in the background liquid, we used glycerol mixtures. We report here two representative experiments, which are summarized in table SIV and are shown in Figs S4 and S5.

## S9.1 Glycerol background and 2-propanol added, visualized with particles

The first experiment is shown in Fig. S4, and aimed to also remove possible complications arising from the presence of dye. We therefore used glycerol (99% purity) as the maze liquid, and 2-propanol (99% concentration) as the exogenous surfactant. The flow was visualized by plastic particles (Pliolite, nominal
diameter approximately 100 µm) that were dispersed on the surface before the endogenous surfactant was introduced. The particles had nominal density 1.04·10⁻³ kg/m³, such that they floated over the glycerol, whose density was around 1.25·10⁻³ kg/m³. A pipette was used to introduce 20 µl of 2-propanol at *t* = 0.

In the experiments of [1], the bottom of the channels was transparent and a light source was placed below the maze, such that the white milk stood out against the black maze. In our new experiment shown in Fig. S4, the particles are white, and we cannot use the same lighting technique of [1], as the particles would be invisible against the background lighting. Instead, we manufactured a matt black maze, and illuminated it from above, such that the particles stood out against the background, as shown in Fig. S4.

In order to calculate particle tracks, images were analyzed using the Trackmate plugin in FIJI [12]. Selected snapshots are shown in panels (a) through (d) of Fig. S4. The tracks show how the surface flow in the main branch is directed primarily to the maze exit; this is especially evident in the final image, namely Fig. S4(d). The tracks also visualize motion in the side branches and in the reservoir. For example, we can see that particles in the destination reservoir (on the left of each panel in Fig. S4) began moving well before any particles from the main branch approached the end of the maze, as shown in Fig. S4(c,d); this observation is consistent with the long-range interactions associated with our transport model.

## S9.2 Background and added liquids both comprising 2-propanol and glycerol

To further illustrate the robustness of the maze-solving behavior arising from exogenous/endogenous surfactant interaction, we performed another experiment where the background and added liquids both contained glycerin and the same surfactant (2-propanol), but in different proportions. A small amount of red dye was introduced in the added liquid as a passive tracer. Specifically, the maze liquid was a solution (by mass) of 10% 2-propanol and 90% glycerol, whereas the added liquid consisted of 50% 2-propanol and 50% glycerol.

The experiment is depicted in Fig. S5. The dyed fluid reached the outlet in approximately one minute, yielding a pattern similar to that shown in Fig. 1 of our Letter, although with fuzzier edges, as shown in Fig. S5(d). This result further supports our proposed model for the surfactant-solving mechanism.

There are additional interesting observations that can be made from Fig. S5. At late stages, a thin ridge of dye is visible ahead of each dye front. We hypothesize that this is due to a secondary effect, which partially analogous to the ‘tears of wine’ phenomenon, as well as to the ‘coffee ring’ effect. Namely, we believe the effect may be driven by alcohol evaporation near the exogenous surfactant front, resulting in surface tension that is locally slightly higher than that of the following liquid. The locally higher tension draws dyed liquid to the front, where dye accumulates as the alcohol evaporates.

## S9.3 Summary of new experiments

In summary, for the different systems tested, the same underlying maze solving behavior is still observed. This provides further experimental evidence to support the main conclusion of this study. We expect that the maze solving behavior is possible with a wide range of liquids, surfactants and tracers, as long as they meet the conditions listed above – namely that the surfactant must be mild, that the exogenous surfactant and tracer must float above the maze liquid, and that the flow remains laminar (e.g. the viscosity must be high enough). Therefore, we anticipate that analogous interactions between exogenous and endogenous surfactants in confined geometries can be important in other systems where surfactant concentrations are out-of-equilibrium.

## References

[1] F. Temprano-Coleto, F. J. Peaudecerf, J. R. Landel, F. Gibou, and P. Luzzatto-Fegiz, Phys. Rev. Fluids 3, 100507 (2018).

[2] W. D. Rutz, C. Whitnah, and G. Baetz, J. Dairy Sci. 38, 1312 (1955).

[3] A. Bakshi and D. Smith, J. Dairy Sci. 67, 1157 (1984).
<table><caption>Table SIV: Table describing the experimental conditions of the qualitative experiments conducted to test the main conclusion of our modelling study: that the maze solving dynamics reported in [1] is primarily due to the interaction between exogenous and endogenous surfactants.</caption><thead><tr><th>Experiment #</th><th>Maze medium</th><th>Exogenous surfactant</th><th>Passive tracers</th></tr></thead><tbody><tr><td>1</td><td>99% glycerol</td><td>2-propanol</td><td>Pliolite 100 µm</td></tr><tr><td>2</td><td>90% glycerol + 10% 2-propanol</td><td>50% 2-propanol + 50% glycerol</td><td>Dye</td></tr></tbody></table>

[4] D. Kristensen, P. Jensen, F. Madsen, and K. Birdi, *J. Dairy Sci.* **80**, 2282 (1997).

[5] O. E. Jensen and J. B. Grotberg, *J. Fluid Mech.* **240**, 259 (1992).

[6] O. E. Jensen and D. Halpern, *J. Fluid Mech.* **372**, 273 (1998).

[7] R. Mcnair, O. E. Jensen, and J. R. Landel, *J. Fluid Mech.* **930**, A15 (2022).

[8] J. Grotberg, D. Halpern, and O. Jensen, *J. Appl. Physiol.* **78**, 750 (1995).

[9] M. Brio, J.-G. Caputo, and H. Kravitz, *Appl. Numer. Math.* **172**, 99 (2022).

[10] Q. Tian, L. Lai, Z. Zhou, P. Mei, Q. Lu, Y. Wang, D. Xiang, and Y. Liu, *J. Agric. Food Chem.* **67**, 6336 (2019).

[11] C.-H. Chang and E. I. Franses, *Colloids Surf. A Physicochem. Eng. Asp.* **100**, 1 (1995).

[12] D. Ershov, M.-S. Phan, J. W. Pylvänäinen, S. U. Rigaud, L. Le Blanc, A. Charles-Orszag, J. R. W. Conway, R. F. Laine, N. H. Roy, D. Bonazzi, G. Duménil, G. Jacquemet, and J.-Y. Tinevez, *Nat. Methods* **19**, 829 (2022).