# A unified theory of pore-scale chaotic advection

Daniel Robert Lester, Joris Heyman, Yves Méheust, Tanguy Le Borgne

► To cite this version:

Daniel Robert Lester, Joris Heyman, Yves Méheust, Tanguy Le Borgne. A unified theory of pore-scale chaotic advection. Journal of Fluid Mechanics, 2025, 1017, pp.A13. 10.1017/jfm.2025.10477. insu-05214991

HAL Id: insu-05214991

https://insu.hal.science/insu-05214991v1

Submitted on 19 Aug 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
JM

PAPERS

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

# A unified theory of pore-scale chaotic advection

Daniel Robert Lester<fnref n="1" />, Joris Heyman<fnref n="2" />, Yves Méheust<fnref n="2" />,<fnref n="3" /> and Tanguy Le Borgne<fnref n="2" />

**Corresponding author:** Daniel Robert Lester, daniel.lester@rmit.edu.au

(Received 15 April 2025; revised 14 June 2025; accepted 11 July 2025)

Recent studies reveal the central role of chaotic advection in controlling pore-scale processes including solute mixing and dispersion, chemical reactions, and biological activity. These dynamics have been observed in porous media (PM) with a continuous solid phase (such as porous networks) and PM comprising discrete elements (such as granular matter). However, a unified theory of chaotic advection across these continuous and discrete classes of PM is lacking. Key outstanding questions include: (i) topological unification of discrete and continuous PM; (ii) the impact of the non-smooth geometry of discrete PM; (iii) how exponential stretching arises at contact points in discrete PM; (iv) how fluid folding arises in continuous PM; (v) the impact of discontinuous mixing in continuous PM; and (vi) generalised models for the Lyapunov exponent in both PM classes. We address these questions via a unified theory of pore-scale chaotic advection. We show that fluid stretching and folding (SF) in discrete and continuous PM arise via the topological complexity of the medium. Mixing in continuous PM manifests as discontinuous mixing through a combination of SF and cutting and shuffling (CS) actions, but the rate of mixing is governed by SF only. Conversely, discrete PM involves SF motions only. These mechanisms are unified by showing that continuous PM is analogous to discrete PM with smooth, finite contacts. This unified theory provides insights into the pore-scale chaotic advection across a broad class of porous materials and points to design of novel porous architectures with tuneable mixing and transport properties.

**Key words:** chaotic advection, porous media

<footnote n="1">¹School of Engineering, RMIT University, Melbourne, Australia</footnote>

<footnote n="2">²Université de Rennes, CNRS, Géosciences Rennes, Rennes, France</footnote>

<footnote n="3">³Institut Universitaire de France (IUF), Paris, France</footnote>
# 1. Introduction

Fluid flow in porous media plays host to a broad range of chemical, physical, biological and geological processes (Dentz et al. 2011; Rolle & Le Borgne 2019; Valocchi, Bolster & Werth 2019). These processes are chiefly controlled by the transport, mixing and dispersion of solutes, nutrients, colloids and microorganisms in the fluid phase. Therefore, detailed knowledge of the mixing dynamics of these flows is required to quantify, understand and predict these phenomena. For example, incomplete mixing of solute plumes has been recognised (Gramling, Harvey & Meigs 2002; de Anna et al. 2014; Berkowitz et al. 2016; Wright, Zadrazil & Markides 2017; Valocchi et al. 2019) to significantly impact the propagation of chemical reactions. Hence, reactive transport models based on assumptions of complete mixing can lead to large errors (Dentz et al. 2011). In recent years, a new understanding of mixing in porous media has emerged through the characterisation of the pore-scale kinematics that drive mixing (Lester, Metcalfe & Trefry 2013; Lester, Dentz & Le Borgne 2016b; Kree & Villermaux 2017; Turuban et al. 2018, 2019; Heyman et al. 2020; Souzy et al. 2020, 2021). These studies have established that chaotic mixing – where fluid particles undergo chaotic orbits and fluid elements are stretched exponentially in time – is inherent to steady flow in almost all 3-D porous media. These ubiquitous kinematics have profound consequences for fluid-borne phenomena (Aref et al. 2017), including sustenance of chemical gradients at the pore-scale (Heyman et al. 2020), accelerated diffusive mixing (Lester, Metcalfe & Rudman 2014a), transport (Turuban et al. 2019) and dispersion (Lester et al. 2014b) of diffusive solutes, augmented clustering and deposition of colloidal particles (Ouellette, O'Malley & Gollub 2008; Sapsis & Haller 2010), singular enhancement of autocatalytic and competitive reactions (Károlyi et al. 2000; Tel et al. 2005), enhanced chemical signalling (Stocker 2012), metabolic pathways (Károlyi et al. 2002; Neufeld & Hernandez-Garcia 2009), and augmented alignment of particles (John & Mezić 2007). A complete understanding of the mechanisms of chaotic mixing in porous media is critical for characterisation and prediction of these phenomena. Such understanding also facilitates development of novel methods to control and optimise these processes directly in engineered systems or via interventions in natural systems.

Porous media may be broadly classified into two distinct classes based on continuity of the pore-scale solid phase. The first group, continuous porous media, has a solid phase which comprises a continuous, smooth medium, such as the open porous networks shown in figure 1(a–e), which include vascular networks, biological materials, tissue scaffolds, ceramic and metallic foams, static and microfluidic mixers, and catalyst supports. Chaotic mixing arises in continuous porous media through the continual splitting and recombination of fluid elements at pore branches and mergers, leading to efficient mixing of fluids. The second group, classed as discrete porous media, involve a solid phase which consists of a jammed array of discrete particles that have small point-like contacts. Such granular matter includes gravels and sands, packed and jammed media shown in figure 1(f–j). Although contacts between discrete particles with finite elasticity are always finite-sized (due to non-zero confining pressure), their contacts are often very small and so are considered to be point-wise herein. In this study, we show that transition from infinitesimal to very small but finite contacts makes no practical difference to the mixing dynamics. Chaotic mixing in discrete porous media arises due to the continual distortion of fluid elements as they flow through highly tortuous paths between grains, leading to rapid deformation and mixing of a dye plume, as shown in figure 1(j).

This classification into discrete and continuous covers most types of porous and permeable media except for fractured media, which have constrained mixing dynamics due
<figure><img src="image_3.png" /><figcaption>Figure 1. (a–e) Various porous networks as examples of continuous porous media: (a) tofu microstructure (Huang et al. 2018); (b) gyroidal tissue scaffold (Melchels, Feijen & Grijpma 2009); (c) ceramic foam (https://filterceramic.com/alu-ceramic-foam-filter); (d) vascular network of the heart (Huang et al. 2009); (e) mixing of dyes in a 3-D micromixer (Therriault, White & Lewis 2003). (f–j) Various granular materials as examples of discrete porous media: (f) granular sandstone (El Bied, Sulem & Martineau 2002); (g) corn kernels; (h) packed corks; (i) glass beads; (j) mixing of a continuously injected dye plume through a random glass bead pack (Heyman et al. 2020). Fluid is refractive index-matched with the beads and only a few beads are shown (grey) at 40 % of their true diameter.</figcaption></figure>

to the pseudo-two-dimensional (2-D) nature of the fractures, and heterogeneous systems such as granular assemblies of porous particles, which may be considered as multi-scale combinations of continuous and discrete porous media. The fundamental differences in pore-scale architecture between continuous and discrete porous media is an important distinction as the fluid dynamical features at the fluid/solid interface drives chaotic mixing in both continuous (Lester et al. 2013, 2016b) and discrete (Turuban et al. 2018, 2019) porous media. The characteristics of chaotic mixing in these distinct classes are illustrated in figure 2, which shows the evolution of particle trajectories and dye plumes in these media as well as the invariant structures (skin friction field, hyperbolic manifolds, critical lines and points) that govern chaotic mixing, as shall be explained in § 2.1.

Recent experiments have directly (Kree & Villermaux 2017; Heyman et al. 2020) and indirectly (Souzy et al. 2020; Heyman, Lester & Le Borgne 2021) observed chaotic mixing
<figure><img src="image_4.png" /><figcaption>Figure 2. (a–c) Characteristics of chaotic mixing in discrete porous media: (a) numerically reconstructed trajectories of tracer particles, taken from PIV experiments within a glass bead pack (adapted from Souzy et al. 2020); (b) numerically computed skin friction field $u(x)$ over the surface of a sphere for steady three-dimensional (3-D) Stokes flow within a bead pack (other spheres not shown) with node $x_p^n$ (green) and saddle $x_p^s$ (black) points, and one-dimensional (1-D) stable $W_{1-D}^s$ and unstable $W_{1-D}^u$ manifolds (black lines). Inset: the same sphere with streamlines shown close to the surface, indicating separation of streamlines in the vicinity of the two-dimensional (2-D) unstable manifold $W_{2-D}^U$. Image courtesy of Régis Turuban, Scuola Internazionale Superiore di Studi Avanzati, Italy; (c) sequences of experimental 2-D dye trace images for steady flow in a random bead pack at planes normal to the mean flow at different distances $x$ downstream from the injection point, measured in terms of the bead diameter $d$. These images show that bead contacts systematically trigger stretching and folding of fluid elements, leading to the formation of sharp cusps in the dye filament. Numbers label fixed spheres while arrows depict directions of fluid stretching (adapted from Heyman et al. 2020). (d–g) Characteristics of chaotic mixing in continuous porous media: (d) numerical simulation of Stokes flow mixing of a diffusive scalar in an archetypal element of an open (continuous) porous network involving a connected pore branch and merger, illustrating the formation of striated material distributions due to fluid stretching and folding which arises at (e) the saddle-type stagnation point ($x_p^s$) in the skin friction field; (f) experimental images of dyed fluid distribution near the ‘pore merger’ in a macroscopic analogue of the pore branch and merger shown in panel (d); (g) dyed fluids at the inlet (top) and outlet (bottom) of the macroscopic pore merger. Cross-section of the dye distribution exiting the pore merger (not shown) agrees well with the outlet scalar distribution shown in panel (d) (adapted from Lester & Chryss 2019).</figcaption></figure>
in discrete porous media. Figure 2(a–c) shows that chaotic mixing in granular media is generated at contact points between grains, leading to repeated stretching and folding motions of fluid elements that rapidly generate highly elongated and ramified material distributions within the pore space. Similarly, pore-scale chaotic mixing has been observed in the macroscopic analogue of an open porous network (Lester & Chryss 2019), and is generated at saddle-type stagnation points on the pore boundary (figure 2d–g). Both sets of observations are important as chaotic mixing generates complex, highly striated material distributions that can quickly obscure their mechanistic origins. In both cases, saddle or contact points appear to generate fluid stretching and folding motions in the fluid bulk, the hallmarks of chaotic mixing. Despite these detailed observations, there exist important outstanding challenges regarding the nature and origins of chaotic mixing in continuous and discrete porous media.

(i) An outstanding challenge is to obtain a unified framework for the description of chaotic mixing across both discrete and continuous porous media.

(ii) The non-smooth geometry at contact points in discrete porous media invalidates the topological theory (Lester et al. 2013) of chaotic mixing in continuous porous media.

(iii) It is unknown how exponential fluid stretching (a characteristic of chaotic mixing) arises at contact points in discrete porous media (Turuban et al. 2019).

(iv) Although the mechanisms of fluid stretching in continuous porous media are well understood, an understanding of how fluid folding (Thiffeault 2004) arises is incomplete.

(v) Experimental and numerical studies show evidence of discontinuous mixing (involving fluid cutting and shuffling) in continuous porous media that is not understood.

(vi) Although *ab initio* models for the Lyapunov exponent have been developed for continuous porous media (Lester et al. 2013), no analogue exists for discrete porous media.

These challenges (i)–(vi) highlight that the current understanding of chaotic mixing in both continuous and discrete porous media is incomplete. In this study, we address these outstanding questions and develop a unified description of chaotic mixing in discrete and continuous porous media. We highlight the differences and similarities between the mixing mechanisms in both classes of porous media, and connect these theories to the observations and dynamical structures shown in figure 2. This unified description of chaotic mixing in porous media provides deep insights into the generation of chaos, and facilitates prediction and optimisation of mixing and transport across a wide range of porous materials.

Throughout this study, for both continuous and discrete porous media, the approach we shall use is overtly topological in that the analysis is based on the topology of porous materials that inevitably arises from pore branches and mergers in the case of continuous porous media, and grain contacts in the case of discrete porous media. We shall show that topological complexity is an inherent characteristic of both classes of media and that this complexity inevitably gives rise to chaotic mixing at the pore-scale. As this approach is topological, it also is quite general and applies to topologically equivalent porous materials (i.e. all those that possess pore branches/mergers or grain contacts), regardless of their specific geometry. Indeed, recent experiments (Heyman et al. 2021) over a range of porous media have demonstrated the ubiquity of pore-scale chaotic mixing. Although the specific geometry may induce further mixing and transport effects, we show that the simplest embodiments of continuous and discrete porous media are sufficient to generate chaotic mixing.
D.R. Lester, J. Heyman, Y. Méheust and T. Le Borgne

To simplify exposition, we therefore limit scope to highly idealised pore networks and granular media, and ignore factors such as surface roughness, distributions of particle pore shapes and sizes. Although porous materials with these features are considerably more complex, they still share the basic classification into discrete or continuous porous media, based on whether the solid phase is a continuum or comprises discrete elements with point-like contacts. While this distinction may appear trivial for e.g. rough porous rocks or assemblies of rough granular matter, this distinction is important from the perspective of understanding chaotic mixing. As the stagnation points of the flow which drive chaotic mixing in continuous porous media become degenerate contact points in discrete porous media, it is unclear how chaotic mixing arises in these media. It is also important to note that while transport and dispersion in both granular matter and open porous networks may both be represented via pore network models (Bijeljic, Muggeridge & Blunt 2004), these models are not sufficient to resolve the mechanisms of pore-scale fluid mixing, where the distinction between discrete and continuous porous materials is significant.

The remainder of this paper is organised as follows. Topological equivalence of both porous media classes is established and a unified description of chaotic mixing is developed in § 2, addressing challenge (i) stated previously. The mechanisms of fluid stretching in both porous media classes are considered in § 3, addressing challenges (ii) and (iii). The mechanisms of fluid folding in both porous media classes are considered in § 4, addressing challenge (iv). The origins and implications of discontinuous fluid mixing are investigated in § 5, addressing challenge (v). In § 6, predictive models for fluid stretching in both porous media classes are developed, addressing challenge (vi), and conclusions are given in § 7.

## 2. Topology of discrete and continuous porous media

### 2.1. Background

To address challenges (i)–(vi), in this section (and Appendices A, B) we first briefly review the theory (Lester et al. 2013, 2016b) of fluid stretching in continuous porous media. To distinguish between the dynamical systems notion of ergodic mixing of non-diffusive fluid particles and physical mixing of a diffusive solute, throughout, we use ‘fluid mixing’ or ‘chaotic mixing’ to describe the former process and ‘diffusive mixing’ for the latter. We denote $v(x)$ as the steady divergence-free three-dimensional (3-D) fluid velocity field in the pore space $\Omega$ of the continuous porous medium, which satisfies the no-slip condition $v=0$ on the entire fluid/solid boundary $\partial\Omega$ and is governed by steady Stokes flow driven by a mean pressure gradient

$$ \mu \nabla^2 v(x) - \nabla p = 0, \quad \forall x \in \Omega, \qquad (2.1) $$

with $\langle \nabla p \rangle = \text{const.}$ We also define $u(x) = \partial v / \partial x'_3 = n \cdot \nabla v$ (where $n$ is the outward normal vector) as the skin friction field on $\partial\Omega$, where $x'_3$ is the local spatial coordinate normal to $\partial\Omega$. Note that while the velocity field $v$ is zero on $\partial\Omega$ (due to the no-slip condition), in general, the skin friction field $u$ is non-zero.

MacKay (1994) shows that critical points $x_p$ of the skin friction field (where $u(x_p) = 0$) such as stagnation or re-attachment points play a central role in the generation of chaotic mixing. These critical points generate local exponential fluid stretching if they are non-degenerate, i.e. if the eigenvalues $\eta_i$, $i = 1:3$ of the skin friction gradient tensor $A = \partial u / \partial x$ at $x = x_p$ are all non-zero. Without loss of generality, we denote the eigenvalues tangent to $\partial\Omega$ as $\eta_1, \eta_2$ (with $\eta_1 \le \eta_2$), and $\eta_3$ is the interior eigenvalue whose eigenvector is normal to $\partial\Omega$. The divergence-free condition means that the eigenvalues $\eta_i$
<figure><img src="image_5.png" /><figcaption>Figure 3. Schematic of the structure of the skin friction field **u** surrounding type I–IV critical points (black dots, summarised in Appendix A) on a portion (bounded by the dotted lines) of the fluid boundary ∂Ω and the associated stable Ws and unstable Wu manifolds. The interior 2-D manifolds for type III, IV critical points are shown as light blue surfaces. Arrows indicate the eigenvectors of the skin friction gradient tensor and the double arrows on the streamlines reflect the sum η₁ + η₂ + 2η₃ = 0. Adapted from Lester, Dentz & Le Borgne (2016b).</figcaption></figure>

satisfy (MacKay 1994)

$$ \eta_1 + \eta_2 + 2\eta_3 = 0, \qquad (2.2) $$

as $\sum_{i=1}^3 \eta_i = \nabla \cdot \mathbf{u} = -\eta_3$. Therefore, non-degenerate critical points $\mathbf{x}_p$ consist of one stable and two unstable eigenvalues or vice versa, and form either reattachment ($\eta_3 < 0$) or separation ($\eta_3 > 0$) points, as shown in figure 3, corresponding respectively to either the collision or emanation of an interior streamline with the critical point. As explained in Appendix A, saddle type critical points with $\eta_1 < 0$, $\eta_2 > 0$ generate 2-D stable $W_{2-D}^s$ and unstable $W_{2-D}^u$ hyperbolic manifolds which are 2-D material surfaces that respectively exponentially contract or expand into the fluid bulk. Unless symmetry conditions are imposed (Haller & Mezic 1998), these invariant 2-D hyperbolic manifolds intersect transversely in the fluid domain, forming a heteroclinic tangle, the hallmark of chaotic dynamics in Hamiltonian systems (Ott 2002). In practice, a single transverse intersection of stable $W_{2-D}^s$ and unstable $W_{2-D}^u$ 2-D manifolds implies a chaotic tangle of infinitely many (Ottino 1989), leading to strong fluid stretching and folding, and chaotic mixing. MacKay (1994) shows that interior 2-D unstable manifolds form a skeleton of the flow that comprises these surfaces of locally minimal transverse flux for diffusive solutes, and so organise both fluid and solute transport and mixing. Conversely, if the interior hyperbolic manifolds are one-dimensional (1-D), their impact on fluid transport is minimal.

Thus, only saddle points give rise to interior 2-D hyperbolic manifolds and chaotic mixing. The prevalence of saddle points on ∂Ω is related to the topological complexity of the continuous porous medium, given by the average number density ρ(χ) of the Euler characteristic χ(Ω), which is strongly negative (Vogel 2002; Scholz et al. 2012), reflecting the topological complexity inherent to all porous materials. For closed bounded manifolds Ω, the Euler characteristic of Ω is related (Armstrong et al. 2019) to that of its boundary δΩ as χ(δΩ) = 2χ(Ω), and hence the Euler characteristic of the pore boundary is also strongly negative. The Poincaré–Hopf theorem connects the Euler characteristic χ(δΩ) of the pore boundary to the sum of indices γp of stagnation points $\mathbf{x}_p$ in δΩ as

$$ \chi(\delta\Omega) = 2(1 - g) = \sum \gamma_p(\mathbf{x}_p) = n_n - n_s, \qquad (2.3) $$

where g is the topological genus of the pore boundary, γp = +1 for node points and γp = −1 for saddle points, and so $n_n$ and $n_s$ respectively denote the number of node and saddle points, and $|ρ(χ)|$ provides a lower bound for the number density of saddle points which naturally arise at pore branches and mergers (figure 4). Although the specific
<figure><img src="image_6.png" /><figcaption>Figure 4. Schematic of a (a) pore branch and (b) merger in continuous porous media. Non-degenerate critical points $x_p$ (black dots) generate 2-D hyperbolic stable $W_{2-D}^3$ and unstable $W_{2-D}^u$ manifolds (grey) which are surfaces of locally minimum transverse flux. The angles $\Delta$, $\delta$ characterise the relative orientation of pore branch and merger elements in the pore network. The red lines pertain to §5 and depict evolution of a continuously injected material line (red). Segments AB and CD of this material line are separated by the critical line $\zeta$ in the pore branch, and are advected through different branches of these pores. Dotted red lines indicate connected material elements that are not resolved by the spatial maps $M, M^{-1}$ defined in (B3).</figcaption></figure>

geometry of pore branches and mergers may vary, the basic topology (shown in figure 4) that drives chaotic mixing is universal to all continuous porous media.

For open porous networks, this leads to simple predictive models (detailed in Appendix B) for the infinite-time Lyapunov exponent $\lambda_{\infty}$ as a function of the geometry of the pore network. In this study, we consider $\lambda_{\infty}$ as the volume average $\langle \cdot \rangle$ of the infinite-time limit $\hat{\lambda}_{\infty}(X)$ of the finite-time Lyapunov exponent $\lambda(X, t)$ as

$$ \lambda_{\infty} \equiv \langle \hat{\lambda}_{\infty}(X) \rangle = \left\langle \lim_{t \to \infty} \lambda(X, t) \right\rangle, \quad (2.4) $$

where $X$ denotes Lagrangian space. For ordered porous media, there can exist a mixture of distinct regions of regular (non-chaotic) flow (where $\hat{\lambda}(X) = 0$) and chaotic mixing (where $\hat{\lambda}(X) > 0$), in which case $\lambda_{\infty}$ represents a volume-weighted average of $\hat{\lambda}(X)$. Conversely, for random porous media, fluid particle trajectories are ergodic (regardless of whether the system is chaotic), hence the infinite-time limit $\hat{\lambda}_{\infty}(X)$ is invariant and so $\lambda_{\infty} = \lambda_{\infty}(X)$.

## 2.2. Topological complexity of discrete and continuous porous media

The fundamental elements of continuous and discrete media are respectively the solid grain, and the connected pore branch and merger shown in figures 5(a) and 5(c). Connections of pore bifurcations (figure 5c) form an extensive 3-D pore network similar to those shown in figure 1(a–e), whereas assemblies of grains (figure 5a) form granular media similar to that shown in figure 1(f–j).

Although pore networks can differ with respect to pore topology, and pore size and shape distributions, they all share the basic feature of many branching and merging pores that can be represented by the basic pore bifurcation shown in figure 5(c). As discussed in §2.1, the connection of many such bifurcations renders porous networks topologically
<figure><img src="image_7.png" /><figcaption>Figure 5. Topological equivalence of the 2-D pore boundary $\delta\Omega$ separating the fluid (pore) $\Omega$ and solid domains of the fundamental elements of (a) discrete and (c) continuous porous media. The normal vector $n$ indicates the normal vector pointing into the fluid domain (pore) $\Omega$ from $\delta\Omega$, and $\delta\delta\Omega$ (green lines) is the 1-D boundary of the pore boundary $\delta\Omega$. $\delta\Omega$ is coloured according to its local Gaussian curvature $K$. (a) Pore boundary $\delta\Omega$ of single spherical grain (semi-transparent) with four contact points (black) associated with contacting grains and uniform positive curvature ($K = +1$) in discrete porous media. (b) Pore boundary of the same grain as panel (a) but with the cusp-shaped contact points smoothed to form a smooth pore boundary $\delta\Omega$ with finite-sized connections between contacting grains, forming boundaries $\delta\delta\Omega$. (c) Pore boundary $\delta\Omega$ for a connected pore branch and merger associated with continuous porous media.</figcaption></figure>

complex in that the topological genus $g$ of the pore-boundary $\delta\Omega$ is large, guaranteeing the existence of saddle points on the pore boundary $\delta\Omega$ that generate chaotic mixing in the pore space $\Omega$. For porous networks, the topological genus $g$ and Euler characteristic $\chi(\delta\Omega)$ of the pore network boundary are given by the total number $n_b$ of pore bifurcations as

$$g = 1 - \frac{\chi(\delta\Omega)}{2} = 1 + \frac{n_b}{2}. \qquad (2.5)$$

Hence, a straight pore ($n_b = 0$) with periodic boundary conditions (to avoid trivial complications associated with macroscopic boundaries of the porous medium) is topologically equivalent to a torus and so has genus $g = 1$, while a periodically connected pore branch and merger (termed a pore element) ($n_b = 2$) forms an additional handle and so has genus $g = 2$. Thus, the pair of bifurcations shown in figure 5(c) also has genus $g = 2$ and Euler characteristic $\chi(\delta\Omega) = -2$. As shown in figure 5(c), this genus $g = 2$ corresponds to a pore coordination number (the number of connected pore throats) $C_p = 4$. The addition of more pairs of pore branches and mergers in any topological configuration to form a pore network then increases the genus $g$ as per (2.5). Hence, all continuous porous media are topologically complex in that they possess a large topological genus $g$ number density, which is critical to chaotic mixing.

Similarly, for discrete porous media, close-packing of many grains such as those shown in figure 5(a) generates extensive 3-D granular assemblies that are also topologically complex, which has also been shown (Turuban et al. 2019) to be critical in generating chaotic mixing. Although various granular materials may differ with respect to particle sizes, shapes and number of contacts, they all share the same basic topology in that they comprise discrete particles with point-like contacts, except for specialised cases where discrete grains contact over a finite-sized area. These cases are not considered here for simplicity of exposition, but topological equivalence between this class of granular matter and continuous porous media also applies via the same approach outlined as follows. The major difference between continuous and discrete media is that the pore boundary $\delta\Omega$ of the latter is non-smooth, as the granular matter involves cusp-like contacts between discrete grains, regardless of grain shape and orientation.
D.R. Lester, J. Heyman, Y. Méheust and T. Le Borgne

Under periodic boundary conditions, the topological genus *g* and Euler characteristic *χ*(*δΩ*) of the grain boundary *δΩ* of granular assemblies is given in terms of the total number of grain contacts *n*<sub>c</sub> as (Turuban *et al.* 2019)

$$g = 1 - \frac{\chi(\delta\Omega)}{2} = n_c. \qquad (2.6)$$

As such, subject to periodic boundary conditions, the grain with four contact points shown in figure 5(c) has $n_c = 2$ (as a contact point is shared between two spheres), genus $g = 2$ and Euler characteristic $χ(δΩ) = -2$. From (2.6), granular assemblies are also topologically complex in that they possess a large topological genus *g* per unit volume.

## 2.3. Topological equivalence of discrete and continuous porous media

Topological equivalence between these seemingly disparate classes of porous media may be established by smoothing the cusp-shaped contacts between grains as shown in figure 5(b). By smoothing the cusp-shaped contact between grains, the discrete porous material now has the characteristics of a continuous porous material in that the pore boundary $δΩ$ is smooth and continuous. This smoothed grain has the same topology ($g = 2$, $χ(δΩ) = -2$) as the discrete grain shown in figure 5(a), but it is topologically equivalent (i.e. it can be morphed solely by stretching and deforming) to the coupled pore branch and merger shown in figure 5(c), but with the distinct difference that the fluid domain (indicated by the outward normal **n**) is external to the pore boundary $δΩ$ shown in figure 5(a,b), whereas the fluid domain is internal to the pore branch and merger shown in figure 5(c). Hence, discrete porous media with smoothed contacts may be considered as the phase inverse of continuous porous media.

Despite this inverse relationship, the pore boundaries $δΩ$ in continuous and smoothed discrete porous media are topologically equivalent, as reflected by their topological genus *g*. Thus, the Poincaré-Hopf theorem (2.3) also applies to the skin friction field **u** on the boundary $δΩ$ of smoothed discrete porous media, and so ensures the existence of saddle points and 2-D hyperbolic manifolds regardless of the orientation of **n**. As the pore boundary $δΩ$ of smoothed discrete and continuous porous media are topologically equivalent, they both have the same total negative Gaussian curvature *K* (see Appendix C for details), hence the same basic mechanism drives chaotic mixing in these different porous media classes. The question as to what extent this connection persists for discrete porous media with non-smooth contacts shall be addressed in §§ 3 and 4. Although discrete and continuous porous media have equivalent topology, their geometry is markedly different. This is explored in more detail in Appendix C and leads to generation of different types of critical points than those which arise in continuous porous media. However, as shall be shown in § 3, these do not alter the generation of chaotic mixing.

In summary, when the contacts in discrete porous media are smoothed, chaotic mixing arises via the same fundamental mechanism as continuous porous media. Such smoothed discrete porous media is topologically equivalent to continuous porous media, albeit with a phase inverse which does not impact the basic mechanism for chaotic mixing. Note that although simple representations of continuous and discrete porous media have been used in this section, these results extend to all topologically equivalent media, such as those shown in figure 1. What remains unknown is whether this mechanism for chaotic mixing in continuous porous media extends to non-smooth, discrete porous media. Turuban *et al.* (2018, 2019) and Heyman *et al.* (2020) have respectively identified that fluid stretching and folding is generated at the contact points of such discrete porous media. In §§ 3 and 4 respectively, the detailed mechanisms that govern fluid stretching and folding will be uncovered.
<figure><img src="image_8.png" /></figure>

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

Figure 6. Hyperbolic manifolds, critical points and lines in (a) continuous and (b) discrete porous media. Intersection (dotted green line) of stable (blue surface) and unstable (red surface) 2-D hyperbolic manifolds in (a) pore branch and merger (grey volume) in an open porous network and (b) a finite-volume numerical simulation (with residual $10^{-16}$) of 3-D Stokes flow through a body-centred cubic (bcc) lattice of monodisperse spheres (modified from Turuban et al. (2019), note only a few spheres in the lattice are shown for clarity). The manifold intersection connects the saddle points (black dots) of the branch/merger in panel (a) and the contact points (black dots) between spheres in panel (b). In panel (b), the 2-D manifolds emerge from the skin friction field of the two spheres labelled 1 and 4, and the green points indicate node points. The open plane indicates the orientation of transverse cross-sections in figure 9, and the black cell is the BCC unit cell.

### 3. Fluid stretching in continuous and discrete porous media

#### 3.1. *Comparison of fluid stretching in discrete and continuous porous media*

A schematic of fluid deformation in continuous porous media is shown in figure 6(a). This figure shows a connected pore branch/merger with offset angle $\Delta = \pi/2$ and the associated 1-D and 2-D stable and unstable manifolds which emanate from the saddle points $x_p$ situated at the pore branch and merger junctions. The intersection of the 2-D stable and unstable manifolds forms a 1-D critical line that connects the saddle points. If many of these pore branch/mergers are connected in series in a manner that maintains transverse connections, then multiple 2-D stable and unstable manifolds emanate respectively from the down- and up-stream saddle points into the element shown in figure 6(a), producing many more manifold intersections, a heteroclinic tangle and chaotic mixing. Conversely, if the stable and unstable 2-D manifolds connect smoothly and tangentially, they cancel each other out and so do not produce persistent exponential stretching.

As shown in figure 6(b), a similar mechanism arises in a finite-volume simulation of Stokes flow over a body-centred cubic (bcc) lattice of smooth spheres (Turuban *et al.* 2019) (note only a few spheres are shown in this figure). The same basic invariant features arise in this flow with contact points playing a similar role to saddle points, and stable and unstable 2-D manifolds emanate from the sphere's skin friction fields and intersect along 1-D critical lines that connect contact points between spheres. Multiple manifold intersections also occur from 2-D hyperbolic manifolds emanating from other spheres up- and down-stream of those shown in figure 6(b), leading to a heteroclinic tangle and chaos.

Figure 6(b) also shows isolated saddle points which lie along the 1-D intersection of the 2-D manifolds with the sphere surfaces. An important difference to continuous porous media is that fluid deformation at contact points is limited by the local cusp-shaped geometry, potentially rendering these points degenerate (meaning at least one of eigenvalues $\eta_i$ of A is zero) with non-exponential stretching. Conversely, the smoothed grains shown in figure 5(b) admit non-degenerate saddle points with exponential fluid stretching. Hence, it is unclear how chaotic mixing is generated in discrete porous media.
D.R. Lester, J. Heyman, Y. Méheust and T. Le Borgne

<figure><img src="image_10.png" /><figcaption>Figure 7. Schematic of evolving stable $W_{2-D}^S$ (blue surfaces) and unstable $W_{2-D}^U$ (red surfaces) manifolds as they are advected over a finite-sized connection between two spheres (grey). Black arrows indicate fluid stretching directions in the bulk fluid and skin friction field. Black and green points respectively indicate saddle and node points. These manifolds become degenerate in the neighbourhood of the connection (indicated by transition to grey colour) and exchange stability as they pass over, generating non-affine folding of material lines (green).</figcaption></figure>

## 3.2. Fluid stretching at contact points

Figure 7 depicts the stable and unstable manifolds associated with a smoothed contact between two spheres (analogous to figure 5c) whose centres are oriented normal to the far-field flow direction. The point-wise contact between these spheres has been smoothed to form a smooth ‘neck’ of diameter $a$. Due to the symmetries of this flow, the velocity field $v(x)$ in the symmetry plane between these two spheres has no transverse component. As shown in figure 8(a), this symmetry plane contains an inclusion from the smoothed contact ‘neck’ that is a diameter $a$ disc with saddle points $x_p$ on the upstream and downstream sides that are connected to critical lines. Figure 8(b) shows that as this finite contact shrinks to an infinitesimal point-like contact ($a \to 0$), these two critical points coalesce into a single critical point. As shown in figure 8(b), this critical point appears to be degenerate as the associated steady 2-D flow in the symmetry plane between the spheres cannot generate exponential stretching.

Turuban et al. (2019) proposed that this critical point is a topological saddle (Bakker 1991; Brøns & Hartnack 1999), which is degenerate (i.e. the skin friction gradient tensor $A$ at $x_p$ has zero eigenvalues) and so does not generate exponential fluid stretching as the manifolds are no longer hyperbolic. As a topological saddle has Poincaré index $\gamma = -2$ (Brøns & Hartnack 1999), then from (2.3), the Euler characteristic in discrete porous media is related to the number $n_d$ of degenerate topological saddles, nodes and regular saddles as

$$ \chi(\Omega) = 2(1 - g) = n_n - n_s - 2n_d. \qquad (3.1) $$

As the topological genus per grain $g$ is equal to the number of contacts $n_c$, then under the assumption that all contact points are topological saddles $n_c = n_d$, the number of non-degenerate nodes $n_n$ and saddle $n_s$ points on each grain is independent of the Euler
<figure><img src="image_11.png" /><figcaption>Figure 8. Streamlines (thin) and critical lines (thick) for 2-D velocity field $v_{2-D}(x)$ in the symmetry plane between two spheres connected by (a) a smoothed contact of diameter $a$ and (b) and infinitesimal contact point. Critical points are denoted as either a saddle point $x_p^s$ or a degenerate topological saddle $x_p^d$. Due to the no-slip condition, in both cases, the skin friction field $u(x) = \partial v / \partial x'_3$ on the grain boundaries local to this contact shares the same topology as the velocity field $v_{2-D}(x)$.</figcaption></figure>

characteristic $\chi(\Omega)$; $n_s = n_n - 2$. Hence, it was proposed by Turuban *et al.* (2019) that the simplest possible skin friction field in discrete porous media involves two isolated node points ($n_n = 2, \gamma = +1$) and no saddle points ($n_s = 0$), and therefore no hyperbolic 2-D manifolds in the fluid bulk and no chaotic mixing. Hence, it was concluded that chaotic mixing in discrete porous media is not ubiquitous as it requires bifurcation of the skin friction field beyond its simplest possible state.

However, this conclusion is not supported by the saddle-like flow structures observed in figure 9, which arise at the intersection of the non-degenerate stable $W_{2-D}^s$ and unstable $W_{2-D}^u$ 2-D manifolds that emanate from contact points. These flow structures clearly indicate the hyperbolic nature of the manifolds as they pass through contact points and intersect with each other (also shown in figure 6b), whereas degenerate points are not associated with exponential (hyperbolic) stretching (Turuban *et al.* 2019).

However, deeper inspection reveals that contact points are only degenerate with respect to the 2-D flow $v_{2-D}(x)$ in the symmetry plane shown in figure 7, but still exhibit exponential stretching in the direction transverse to the symmetry plane between contacting grains. In § 2.1, it was shown that isolated node points (type III and type IV critical points in figure 3) cannot generate 2-D hyperbolic manifolds in the fluid interior. However, the interaction between a node and contact point can produce an interior 2-D hyperbolic manifold, as shall first be shown for finite-sized contacts.

### 3.3. Hyperbolic manifolds at finite-sized contacts

Figure 7 shows that a saddle point (green) arises on each side (upstream and downstream) of a smoothed connection, and a pair of node points (black) also arise on each side (upstream and downstream) of the connected spheres. The 1-D manifold along the sphere surfaces connecting the upstream nodes and saddles is unstable with respect to these node points and stable with respect to the saddle points. The combination of this stable 1-D manifold with the stable 1-D manifold in the fluid interior that connects to the upstream saddle point generates the 2-D stable manifold $W_{2-D}^s$ that is shown in figure 7 without material (green) lines. Similarly, a downstream 2-D unstable manifold $W_{2-D}^u$ (figure 7 without material (green) lines) is generated that connects the downstream saddle and node points. Both of these 2-D manifolds are oriented parallel to the line connecting the sphere centres and so are termed parallel 2-D manifolds. Conversely, the manifolds transverse to this connecting line (shown in figure 7 with material (green) lines) are termed transverse 2-D manifolds.

The local cusp-shaped geometry near contact points constrains any (un)stable 2-D manifolds that are generated (up)downstream to pass over this connection as transverse
<figure><img src="image_12.png" /><figcaption>Figure 9. Series of cross-sections (transverse to the mean flow) with increasing downstream distance of Stokes flow through a bcc lattice of spheres between the two bead pairs (1–2) and (3–4) shown in figure 6(b). Purple lines show vector field lines of the 2-D velocity field transverse to the mean flow direction, thick and thin black lines indicate material lines advected by the flow, and red and blue lines respectively represent 2-D unstable and stable manifolds of the flow. Intersection of the 2-D stable and unstable manifolds forms a 1-D curve that connects the contact points of the (1–2) and (3–4) bead pairs. Modified from Heyman *et al.* (2020).</figcaption></figure>

2-D manifolds, as per figure 6(b). The analogous situation for a porous network is shown in figure 6(a), where the 2-D (un)stable manifold is generated by a (up)downstream saddle point and connects with an (down)upstream saddle point. This constrained geometry forces the transverse 2-D stable and unstable manifolds shown in figure 7 with material (green) lines to connect tangentially over the sphere connection. The 1-D intersection (green dotted line in figure 7) of the transverse and parallel 2-D manifolds are responsible for the saddle-type flow structures shown in figure 9.

### 3.4. Hyperbolic manifolds at contact points

If the connection between the spheres in figure 7 shrinks to an infinitesimal contact point, the up- and downstream saddle points coalesce into a single saddle point which is degenerate with respect to the transverse manifolds, but non-degenerate with respect to the parallel manifolds. The interior (un)stable 1-D manifolds that connect from (down)upstream to the contact point, and the exterior (un)stable 1-D manifolds that connect the (down)upstream node points to the contact point both persist as this contact shrinks to a point. Hence, the parallel 2-D stable and unstable hyperbolic manifolds persist for a contact point, and still impart exponential stretching to the fluid bulk.

Hence, node points and saddle-like contact points between grains generate 2-D hyperbolic manifolds and chaotic mixing in discrete porous media. This basic mechanism is very similar to that of continuous porous media, where connected contact and node
points play the role of saddle points in continuous porous media. As there is no need for isolated saddle points (away from contacts), chaotic mixing is inherent to discrete porous media.

# 4. Fluid folding in discrete and continuous porous media

## 4.1. *Fluid folding in discrete porous media*

In this section, we describe the mechanisms that drive fluid folding in continuous and discrete porous media. Although Heyman et al. (2020) show that contact points in discrete porous media generate folding of fluid elements, the detailed mechanism is unclear. Figure 7 shows that when a 2-D unstable transverse manifold impacts a contact point, this manifold becomes locally degenerate at the contact point and transitions to a stable 2-D transverse manifold that emanates from the contact. Conversely, a 2-D stable parallel manifold terminates at the contact point and upstream node points, and a 2-D unstable parallel manifold emanates downstream from the contact point and downstream node points.

This exchange in stability can generate strong folding of fluid elements in the symmetry plane between the spheres, as indicated by the green material lines in figure 7. Although the schematic in figure 7 is fore-aft symmetric, this can be broken by additional grains in the assembly or contacts that are skew to the mean flow direction, meaning that folding of material elements also manifests in the plane transverse to the mean flow direction. As shown in the sequence of panels in figure 9, subsequent fluid stretching elongates these folds into tight cusps and highly filamentous structures such as those observed in figures 2(c) and 9.

These insights provide a link between the folding dynamics described here and by Heyman et al. (2020) and recent studies (Turuban et al. 2018, 2019) that describe the mechanisms leading to pore-scale fluid stretching in discrete porous media. This means that chaotic mixing proceeds in discrete porous media via iterated stretching and folding of fluid elements as they pass over contact points between grains.

## 4.2. *Fluid folding in continuous porous media*

The maps $S_b, S_m$ (B8) can be used to uncover the folding mechanism in continuous porous media. The dye trace images in figure 10 show evidence of cutting, shuffling, stretching and folding of material lines as they are advected through the pore-space. To unravel these motions, we consider the evolution of a continuously injected line source through the pore branch and merger shown in figure 4. In this figure, the line source with segments marked AB and CD is continuously injected along the line $x_r = 0$ before bifurcating into two separate segments at the critical point $x_p$. The intersection of the 2-D stable manifold $W_{2-D}^s$ with the boundary $\partial\Omega$ represents a repelling (i.e. $\nabla \cdot \mathbf{u} > 0$) critical line $\zeta$ of the skin friction field $\mathbf{u}$. Any injected material line that crosses the mid-plane $y_r = 0$ of the inlet pore corresponding to $W_{2-D}^s$ is stretched and folded as it is advected downstream over $\zeta$, and fluid particles on this critical line are held there indefinitely, as shown by the temporal map $T^*$ (B6), which has a divergent residence time for particles entering along $y_r = 0$. Folding of fluid elements over this critical line effectively splits the fluid element as it is advected into two separate pores downstream. Thus, the indefinite holdup of fluid particles at the critical line $\zeta$ manifests as the distinct segments AB, CD in each outlet of the pore branch, yet these are connected via a continuous material line extending from the critical line.
<figure><img src="image_13.png" /><figcaption>Pore geometry (ii):</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Pore geometry (iii):</figcaption></figure>

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

$n=1$

$n=2$

$n=3$

$n=4$

<figure><img src="image_16.png" /><figcaption>Pore geometry (iv):</figcaption></figure>

$n=1$

$n=2$

$n=5$

$n=10$

Figure 10. Evolution of fluid elements (black points) with pore number *n* over the circular inlet plane of a pore branch (figure 4a) for various open porous networks, ranging from (i–iii) ordered to (iv) random pore networks. Red lines depict the web of discontinuities associated with cutting and shuffling of fluid elements.

Upon exiting their respective pores, these segments are reoriented at angle $\Delta$ before entering the branch merger. The solid red lines in figure 4 denote parts of the material line that are resolved by the spatial map $M$ (B3), and the dashed red lines indicate material elements that are not resolved by this map yet remain connected in 3-D due to continuity. The bifurcation and folding process is then reversed in the pore merger (figure 4b), where fluid elements from separate pores are merged together over the attracting (i.e. $\nabla \cdot u < 0$) critical line $\zeta$ formed by intersection of the 2-D unstable manifold $W_{2-D}^u$ and the pore boundary $\partial\Omega$. Fluid elements that pass near this attracting critical line also experience diverging residence times, leading to further folding of material elements prior to the merger of elements arriving from different inlet pores.

When projected onto the (2-D) pore outlet planes, this strong folding of fluid elements over the critical lines $\zeta$ in the 3-D pore space manifests as an effective ‘cutting’ of material elements, as shown in figure 10. In the following subsection, we consider the nature of such discontinuous mixing and its interplay with 3-D fluid stretching and folding. In addition to this strong folding that manifests as discontinuous mixing, the spatial map $M$ also imparts weaker folding that is illustrated by its affine $M_a$ (B4) and non-affine $M_n$ (B5) components.
Fluid cutting is also encoded in the non-affine map $M_n$, such that fluid elements with $y_r > 0$ ($y_r < 0$) are mapped to the upper (lower) branch in figure 4(a). Weak folding of material elements in these 2-D outlets arises from the nonlinear component of the non-affine map $M_n$ and its inverse $M_n^{-1}$, which map fluid particles along the straight line $y_r = 0$ to the curved lower pore boundary. Although these maps only induce weak folding, subsequent fluid stretching amplify these folds into the sharp kinks observed in figure 10. However, figure 10 shows that packing of highly elongated material lines into the 2-D pore outlet is dominated by cutting of material elements rather than weak non-affine folding.

Hence, strong folding of fluid elements in continuous porous media arises via a similar mechanism to that of discrete porous media. In continuous media, fluid elements are held up at pore junctions, whereas for discrete media, fluid elements are held up at contact points, both leading to strong folding. However, for continuous porous media, this hold-up generates discontinuous mixing when projected onto a 2-D plane transverse to the mean flow direction. Such discontinuous mixing does not occur in discrete media due to the infinitesimal contacts. In addition, the non-affine nature of $M$ imparts weak folding that may be amplified into sharp kinks by subsequent fluid stretching. As weak folding plays a secondary role, in § 5, we focus on strong fluid folding which manifests as discontinuous mixing.

# 5. Discontinuous mixing in continuous porous media

## 5.1. *Origins of discontinuous mixing*

Figure 10 shows the evolution of material distributions in an open porous network under stretching and folding (SF) and cutting and shuffling (CS) actions for the following ordered and random pore network types:

(i) ordered network corresponding to 3-D baker's map: $\Delta = \pi/2$, $\delta = 0$, $\lambda_{\infty} = \ln 2 \approx 0.693$;

(ii) ordered network, chaotic mixing: $\Delta = \pi/2$, $\delta = \pi/6$, $\lambda_{\infty} = \ln[(5\sqrt{3} + \sqrt{11})/8] \approx 0.403$;

(iii) ordered network, non-chaotic mixing: $\Delta = \pi/4$, $\delta = \pi/4$, $\lambda_{\infty} = 0$

(iv) random network, chaotic mixing: $\Delta \sim U[0, \pi]$, $\delta \sim U[0, \pi]$, $\bar{\lambda}_{\infty} \approx 0.118$.

As expected, the ordered network (i) shown in figure 10 corresponding to the baker's map exhibits the fastest mixing, whereas the random network (iv) is significantly slower. The evolution of these material distributions shown in figure 10 are produced via application of the advective maps $M, M^{-1}$ described in Appendix B for the various values of $\delta, \Delta$ described previously.

The origins of the CS actions shown in figure 10 are illustrated in figure 4(a), where the line segments AB, CD are advected to the different outlets of the pore merger element in an apparently discontinuous manner. However, as indicated by the dashed red lines in figure 4(b), these segments are actually connected by a material line that is stretched and folded as it is advected over saddle point $x_p$.

In many applications, this 3-D distribution of fluid elements is less relevant than the 2-D distribution in a planar cross-section corresponding to the pore outlet planes in figure 4 and images in figure 10. Material lines and sheets that are continuous in 3-D can manifest as discontinuous elements in this 2-D frame, even within the same pore. Henceforth, we consider mixing of fluid elements in these 2-D planes as being 'discontinuous', although we understand the fluid undergoes smooth and continuous 3-D deformation, unlike e.g. granular matter or plastic materials that deform via slip planes.
D.R. Lester, J. Heyman, Y. Méheust and T. Le Borgne

Hence, chaotic mixing in continuous porous media is akin to mixed cutting and shuffling (CS) and stretching and folding (SF) actions observed during chaotic mixing of discontinuously deforming systems (Juarez et al. 2010; Christov, Lueptow & Ottino 2011; Sturman 2012; Smith et al. 2016, 2017a,b, 2019), such as granular tumblers (Christov et al. 2011; Smith et al. 2017b), and microfluidic micro-mixers based on iterated pore branches and mergers (Sudarsan & Ugaz 2006). The motions imparted by the spatial maps $S_b, S_m$ as fluid elements flow through each pore branch and merger shown in figure 4(a,b) can be summarised as ‘cut-stretch-fold-rotate’ in the pore branch, and ‘compress-fold-glue-rotate’ in the pore merger.

For the pore branch, ‘cut’ refers to splitting of fluid elements along the critical line $\zeta$, ‘stretch’ refers to the affine stretching (by $M_a$) by a factor of 2 in the $y_r$ direction, ‘fold’ refers to non-affine folding (by $M_n$) of elements to conform to the outlet pore boundaries, ‘rotate’ refers to reorientation (by $R(\Delta)$) of fluid elements due to orientation of the downstream pore mergers. In the pore merger, ‘compress’ refers to the compression (by $M_a^{-1}$) of fluid elements by a factor of 1/2 in the $y_r$ direction, ‘fold’ refers to non-affine folding (by $M_n^{-1}$) to conform with the required semi-circular disc, ‘glue’ refers to gluing of these semi-circular discs to form the disc-shaped pore outlet and ‘rotate’ refers to reorientation (by $R(\delta)^{-1}$). As fluid elements are in general not smoothly reconnected at the gluing step, these actions lead to cutting and shuffling of material elements, as shown in figure 10. Conversely, for discrete porous media, the infinitesimal nature of grain contacts means that these discontinuities have zero extent and so only purely continuous deformation arises.

## 5.2. Coupled continuous and discontinuous mixing

The introduction of discontinuous CS actions significantly alters the mixing process as constraints associated with fluid continuity are relaxed. One impact is that the Poincaré index is no longer conserved, leading to the formation of new coherent Lagrangian structures such as leaky KAM surfaces around elliptic tori, and pseudo-elliptic and pseudo-hyperbolic points (Smith et al. 2017a,b). The natural mathematical framework to analyse CS motions is via piecewise isometries (PWIs) (Sturman 2012), where fluid elements are removed (cutting) and placed in different locations and orientations (shuffling). For CS-only systems, these actions generate a web of discontinuities (figure 10), which comprises the forward iterates of the cutting line (corresponding to the line $y_r = 0$ in the map $M$ associated with bifurcation of fluid elements along the critical line $\zeta$ in the pore branch element), and complete mixing occurs if this web becomes space-filling with time. Under these conditions, pure CS actions leads to weak ergodic mixing, and associated mixing measures decay at an algebraic rate in time.

Conversely, for pure SF actions in continuous systems, fluid stretching arises from affine deformations, whereas folding is associated with non-affine deformation. Together, these stretching and folding actions can generate chaotic advection which is characterised by strong ergodic mixing and mixing measures that decay at an exponential rate. To understand how combined SF and CS actions impact the mixing dynamics of continuous porous media in both ordered and random pore networks, we examine the four network geometries (i–iv) described in Appendix B.1, whose mixing dynamics is shown in figure 10. In all cases, these SF and CS actions lead to complete mixing of fluid elements, albeit at significantly different rates.

## 5.3. Measures of continuous and discontinuous mixing

Conventional tools such as Lyapunov exponent, scalar variance decay and entropy measures are not suitable to characterise discontinuous mixing as they are based on **1017** A13-18
dissipative or continuous mixing. The mix-norm measure (Mathew, Mezić & Petzold 2005) was developed to quantify multiscale mixing across diverse applications, including non-dissipative and discontinuous mixing. The mix-norm $\Phi(c) \in [0, 1]$ captures the extent of mixing of a zero mean concentration field $c(\mathbf{x})$ across different normalised length scales $s \in [0, 1]$ (where $s = 1$ corresponds to the velocity correlation length scale $\ell$) as

$$ \Phi(c) = \sqrt{\int_{0}^{1} \phi(c, s)^2 \mu(ds)}, \qquad (5.1) $$

$$ \phi(c, s) = \sqrt{\int_{p \in \Omega} d^2(c, \mathbf{p}, s) \mu(d\mathbf{p})}, \qquad (5.2) $$

$$ d(c, \mathbf{p}, s) = \frac{1}{\mathrm{vol}[B(\mathbf{p}, s)]} \int_{x \in B(\mathbf{p}, s)} c(\mathbf{x}) \mu(d\mathbf{x}), \qquad (5.3) $$

where $\mu$ denotes the Lebesgue measure, $d(c, \mathbf{p}, s)$ is the average of $c(\mathbf{x})$ over the ball $B(\mathbf{p}, s) := \{\mathbf{x} : |\mathbf{x} - \mathbf{p}| < s/2\}$ centred at point $\mathbf{p}$ in the fluid domain $\Omega$, $\phi(c, s)$ is the $L^2$-norm of $d$ over $\Omega$ and $\Phi(c)$ is the $L^2$-norm of $\phi(c, s)$ over $s \in [0, 1]$. As such, the $L^2$-norm $\phi(c, s)$ characterises the variance of the scalar field $c(\mathbf{x})$ after averaging over length scale $s$ and the mix-norm $\Phi(c)$ accounts for mixing across all scales by quantifying the $L^2$ norm of this measure over all $s$.

Hence, $\Phi = 1$ for scalar distributions that are completely segregated at the largest length scale of the system and $\Phi(c) = 0$ for mixtures that are well mixed at all scales. Under weak ergodic mixing characteristic of CS systems, the mix-norm $\Phi(c_n)$ (where $c_n$ is the concentration field after $n$ of iterations of the mixing process) decays algebraically with $n$, whereas under strong mixing characteristic of SF systems, $\Phi(c_n)$ decays exponentially with $n$. For SF systems that can be represented by a linear map such as $M_a$, the mix-norm decays exponentially at a rate of half the Lyapunov exponent (Mathew et al. 2005) as

$$ \Phi(c_n) = \Phi(c_0) \exp \left( -\frac{\lambda_\infty n}{2} \right). \qquad (5.4) $$

Although many studies have considered mixing due to purely CS or SF actions, only a handful (Smith et al. 2016, 2017, 2019a,b) have considered the impact of coupled CS and SF motions. Smith et al. (2017a,b) show that the combination of CS and SF motions can either act to retard or enhance the rate of strong mixing given by SF actions alone, depending upon whether hyperbolic and pseudo-elliptic points form a complete set of building blocks for Lagrangian structures (Smith et al. 2017b), similar to hyperbolic and elliptic points in SF-only systems. In this case, CS can retard mixing and the Lyapunov exponent in (5.4) forms an upper bound for the decay rate of the mix-norm. For mixed CS and SF systems, a key question regarding the impact of discontinuous deformation upon SF mixing is the rate and extent to which the web of discontinuities becomes space-filling. To quantify this, we consider the growth rate of line elements in the pore geometries (i)–(iv) whose mixing dynamics are shown in figure 10.

Figure 11(a) shows that the mean growth of the normalised length $\rho_n \equiv \langle \delta l_n / \delta l_0 \rangle$ of an ensemble of infinitesimal line elements in the ordered and random pore networks is well approximated by the theoretical estimate (B1) of the Lyapunov exponent $\lambda_\infty$, hence, exponential fluid stretching is captured by linearisations of the spatial maps $S_b, S_m$. For the non-chaotic case ($\Delta = \pi/4$, $\delta = \pi/4$), the Lyapunov exponent is zero and the relative length $\rho_n$ grows linearly as $\rho_n \approx \rho_0 + n/2$. As all of these pore geometries completely mix as $t \to \infty$, then the growth of the total length $l_n$ of the web of discontinuities follows the same stretching process as the fluid phase.
<figure><img src="image_17.png" /></figure>

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

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

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

Figure 11. (a) Relative elongation $\rho$ of infinitesimal material lines in (a) ordered and (b) random pore networks where black, red, grey and blue lines and points respectively correspond to the (i) the baker's map, (ii) chaotic, (iii) non-chaotic and (iv) 100 realisations of random pore geometries. Points represent numerical results from maps $S_b, \tilde{S}_m$, and lines represent stretching rates based upon the Lyapunov exponent for ordered (B1) and random (B2) networks, except for the non-chaotic case where stretching evolves linearly as $\rho_n = \rho_0 + \alpha n$. Thin blue lines in panel (b) represent stretching of $10^2$ realisations of the random network and the thick blue line the ensemble average. (c) Growth of the total length $l_n$ of the web of discontinuities with number $n$ of pore branches and mergers with same colour code as for panels (a) and (b). Points represent numerical results from maps $S_b, S_m$, and lines represent analytic predictions based on stretching rates for chaotic (5.5) and non-chaotic (5.6) cases. (d) Evolution of mix-norm $\Phi(c_n)$ in pore networks with chaotic and (inset) non-chaotic mixing, where points represent numerical results from maps $S_b, S_m$ and (5.1)–(5.3), and lines represent the mix-norm estimate (5.4) based upon pure SF motions and the Lyapunov exponent given by (B1) and (B2).

As an additional cut is added to the web of discontinuities at each iteration of the map $\Lambda_n$ (B10), the evolution of $l_n$ is related to the growth of the length $L_n$ of finite material lines with initial length $L_0 = 2$ (corresponding to the length of the cut along the chord $y_r = 0$) as $l_n = \sum_{i=0}^n L_i$. As the exponential stretching rate in chaotic pore-scale flows is log-normally distributed (Lester et al. 2016b) with mean stretching rate $\lambda_\infty$ and variance $\sigma_\lambda^2$, for pore geometries that admit chaotic mixing, these finite line elements grow exponentially as $L_n \approx 2 \exp(\Lambda_\infty n)$, where $\Lambda_\infty = \lambda_\infty + \sigma_\lambda^2/2$ and $\sigma_\lambda^2$ is found (Lester et al. 2016 b) to be $\sigma_\lambda^2 \approx 0.1144$ for random networks, whereas for the ordered networks, this
variance is negligible, $\sigma_{\lambda}^2 \approx 0$. Thus, the length $l_n$ of the web of discontinuities for pore networks that exhibit chaotic mixing evolves with $n$ as

$$l_n \approx 2 \frac{\exp(\Lambda_{\infty}(n+1)) - 1}{\exp(\Lambda_{\infty}) - 1} \qquad (5.5)$$

Conversely, for the non-chaotic case $\Delta = \delta = \pi/4$, the length of line elements grows linearly as $L_n \approx 2+n$ and so the total length of the web of discontinuities grows algebraically as

$$l_n \approx \frac{(n+1)(n+4)}{2} \qquad (5.6)$$

Figure 11(c) shows that for all pore geometries (i)–(iv), the growth rate of the web of discontinuities is well approximated by these scalings, indicating that the stretching of line elements mediates cutting and shuffling.

### 5.4. Interactions between continuous and discontinuous mixing

To determine how CS actions of these maps impact the mixing dynamics, we compute the mix-norm $\Phi(n)$ for all four cases, and the results are summarised in figure 11(c). For the chaotic mixing cases, figure 11(c) indicates that the mix-norm decays exponentially at a rate that is half the Lyapunov exponent in (5.4), suggesting that the CS actions do not contribute significantly to rate of mixing, but rather these are dominated by the exponential SF dynamics. For the non-chaotic case, the mix-norm decays linearly at a rate that is also half the growth rate $\alpha$, also suggesting that mixing is governed by fluid stretching as stretching also mediates the CS process. Although it is tempting to conclude that SF motions dominate the mixing process, it is important to note that the estimate (5.4) of the mix-norm only accounts for the amount of stretching experienced by fluid elements, and does not consider whether packing of stretched material lines into a finite-sized pore is achieved by either folding or cutting of this material.

Figure 10 clearly shows that material lines experience much more cutting (due to strong 3-D folding) rather than weak folding (arising from the non-affine map $M_n$), hence, cutting of material lines is integral to this packing process. Therefore, although CS plays a critical role in mixing by facilitating packing of elongated material elements into the pore-space, of the SF motions, only fluid stretching is required to achieve chaotic mixing as weak fluid folding plays a secondary role. As such, the theoretical estimates for the Lypaunov exponents in both ordered (B1) and random (B2) networks (which have been validated against line stretching simulations using the spatial maps (B8)) provide an accurate means of predicting the mixing rate in continuous porous media, and represent useful quantitative tools for the control and optimisation of mixing in engineered porous media. These insights into these CS and SF actions provide a complete description of the mechanisms that govern mixing in continuous porous media.

## 6. Unified prediction of Lyapunov exponent

While the prediction of the Lyapunov exponent for continuous porous media is well established (see Appendix B.2), current models (Heyman et al. 2020) for discrete porous media rely on experimental or numerical observations. Here, we derive a novel general formulation for the Lyapunov exponent that encompasses continuous and discrete porous media.
<figure><img src="image_21.png" /><figcaption>Figure 12. Illustration of the mechanism motivating the stretching model for discrete porous media, adapted from figure 9. (a) Immediately after a contact, the material line (black) is oriented at an angle $\delta$ from the stable manifold (blue line) due to asymmetry of the bead pack. This leads to a cusp of initial size $\epsilon$ that elongates along the unstable manifold (red line). (b) When approaching the next contact, the cusp has elongated to a length approximately equal to the grain diameter $d$.</figcaption></figure>

## 6.1. Lyapunov model in discrete porous media

In the deformation model for continuous media presented in Appendix B.2, fluid stretching and compression transverse to longitudinal flow direction is decoupled over the pore branch and merger. In a branch, fluid stretching (by a factor 2) in the transverse direction parallel to the branching is fully compensated by deceleration in the longitudinal direction, while the other transverse component is neutral. The reverse situation occurs symmetrically in a merger.

In discrete porous media, fluid stretching and compression act in both transverse directions (figure 9). To formulate a Lyapunov model for discrete media, we thus assume that compression and stretching occur simultaneously and perpendicularly in the transverse plane along the stable and unstable manifolds, respectively. Fluid deformation from one contact point to the next in the 2-D plane transverse to the mean flow direction is approximated by the deformation gradient tensor

$$F = \begin{pmatrix} \alpha & 0 \\ 0 & 1/\alpha \end{pmatrix}, \qquad (6.1)$$

where $\alpha \ge 1$ quantifies stretching along the unstable manifold. Note that the deformation represented by $F$ is equivalent to the deformation $D(\Delta)$ over a couplet in the pore network model (B9) with $\Delta = \pi/2$ and $\alpha = 2$. After passing the contact, the un/stable manifolds are inverted (figure 9a.3), which is reflected by a reorientation of $F$ of angle $\pi/2$ (transition from figure 9a.2 to figure 9a.4). However, the presence of other beads generates an asymmetry in the orientation of the hyperbolic manifolds upstream and downstream of the contact point, inducing a perturbation of angle $\delta$ (figure 9a.4). Therefore, the downstream (un)stable manifold is oriented at an angle $\varphi = \pi/2 + \delta$ from the upstream (un)stable manifold. Since the material line at the contact point is aligned with the upstream unstable manifold, it becomes oriented at an angle $\delta$ from the compression direction downstream of the contact point as per figures 12(a) and 9(a.4).

The asymmetry that generates the perturbation angle $\delta$ also induces a small cusp of initial size $\epsilon$ (figure 12a) that is amplified via compression and stretching to a fold of the order of the grain diameter $d$ at the following contact (figures 12(b) and 9a.5). Hence, the deformation between the two successive contact points is given by $\alpha \approx d/\epsilon$. The net
deformation $D_{2n}$ arising from subsequent stretching events over $2n + 1$ contact points may then be represented as

$$D_{2n} = \prod_{j=1}^{n} \left( F \cdot [R(\varphi) \cdot F \cdot R(\varphi)^{\top}] \right), \qquad (6.2)$$

where $R(\varphi)$ is the rotation matrix. The eigenvalues $r_i$ of $D_{2n}$ are

$$r_i = \left\{ \xi \pm \sqrt{\xi^2 - 1} \right\}^{n/2}, \quad \xi = \cos^2 \delta + \frac{\sin^2 \delta}{2\alpha^2} (1 + \alpha^4). \qquad (6.3)$$

The perturbation angle $\delta$ can be approximated from the initial cusp size $\epsilon$ as $\tan \delta \approx \epsilon/d$, hence, for small perturbation angle $\delta$, $\delta \approx \tan \delta$ and $\delta \approx \epsilon/d = 1/\alpha$. Under these assumptions and in the limit $\delta \to 0$, the eigenvalues $r_i$ are

$$r_i \approx \left( \frac{3 \pm \sqrt{5}}{2} \right)^{n/2}. \qquad (6.4)$$

Hence, the effective stretching from one contact point to the next ($n = 1$) is given by the largest eigenvalue $r \approx 1.618$. If $X_c$ is the average length of the critical line connecting contact points (Heyman et al. 2020), then the elongation $\rho$ over longitudinal distance $X$ is

$$\rho = \exp \left( \frac{X}{X_c} \ln(r) \right), \qquad (6.5)$$

and so the dimensionless Lyapunov exponent is

$$\lambda_{\infty} = \frac{d}{X_c} \ln(r), \qquad (6.6)$$

where $d$ is the mean grain size. This model (6.6) can be tested against simulations of Stokes flow through bcc sphere lattices (Turuban et al. 2018, 2019). In this system, the Lagrangian kinematics and distance between contact points is varied by changing the mean flow orientation $\theta$ with respect to the lattices symmetries (Turuban et al. 2019). The latter was estimated by measuring the distance of critical lines between successive contacts, described by an analytical model using simple geometrical principles. The prediction of (6.6) is in relatively good agreement with the measured Lyapunov exponents $\lambda_{\infty}$ for all angles $\theta$ using either numerical estimates or the analytical expression for $X_c$ (figure 13). It captures the evolution of the absolute value of the Lyapunov as a function of $\theta$, while previous predictions for this system were only relative (Turuban et al. 2019). Note that Heyman et al. (2020), who measured experimentally the Lyapunov exponent $\lambda_{\infty}$ in random loose granular packings, linked it to a folding frequency assuming $r = \ln(2)$. The different value of $r$ derived here should therefore introduce a minor adjustment to the estimated folding frequency in random packings

## 6.2. General formulation for the Lyapunov exponent

The models of (B1) and (6.6) for the Lyapunov exponent in continuous and discrete porous media, respectively, can be described by the general expression

$$\lambda_{\infty} = f \ln(r), \qquad (6.7)$$
<figure><img src="image_22.png" /><figcaption>Figure 13. Comparison of predicted and computed (Turuban et al. 2019) Lyapunov exponent $\lambda_{\infty}$ in BCC sphere lattices with various flow orientations $\theta$ with respect to the lattice symmetries. The numerically computed of Lyapunov exponents are shown as black circles. The predictions of (6.6) using the numerically measured distance $X_c$ and the analytical approximation $X_c(\theta)$ are shown as red squares and red dashed lines respectively.</figcaption></figure>

with $f$ the dimensionless frequency of folding events and $r$ the average incremental stretch at each event,

$$r = \lim_{n \to \infty} \left\langle \frac{\rho_{n+1}}{\rho_n} \right\rangle. \qquad (6.8)$$

For continuous media, the frequency $f$ of folding events is unity, as folding occurs at every pore branching or merging, but the magnitude of stretching $r \in [1, 2]$ in (B1) varies with the pore network parameters $\delta, \Delta$ (Lester et al. 2013). Conversely, for discrete media, the stretching increment is fixed at $r = 1.618$, but the frequency of contacts varies with the average length of critical lines between contact points as $f = d/X_c$ as per (6.6) which can be controlled by e.g. altering the flow orientation in crystalline lattices (Turuban et al. 2018).

These distinct controls on the Lyapunov exponent $\lambda_{\infty}$ are illustrated in figure 14 over the phase space $\{f, r\}$. For continuous media, the Lyapunov exponent can vary from zero to the theoretical maximum $\ln 2 \approx 0.6931$ for continuous 3-D systems given by the baker's map. For discrete media $\lambda_{\infty}$ ranges from zero to approximately 0.218 based on an upper bound of $d/X_c \approx 0.45$. While the considered granular packings and pore networks shown in figure 14 cover limited regions of the phase space $\{f, r\}$, more complex porous media, such as rocks or hierarchical media, may have mixed continuous/discrete porous media attributes and thus exhibit mean Lyapunov exponents that lie in between these extremes.

## 7. Conclusions

Despite fundamental differences in pore-scale geometry and topology, we have shown that discrete and continuous porous media share the same basic mechanism for the generation of chaotic mixing, leading to a unified description of pore-scale mixing in continuous and discrete porous media. The topological complexity inherent to all porous media generates exponential fluid deformation at critical points in the fluid/solid boundary that propagates into the fluid interior. This stretching manifests as 2-D hyperbolic un/stable manifolds that form a transverse heteroclinic connection in all but highly symmetric cases, generating fluid stretching and folding and chaotic mixing.
<figure><img src="image_23.png" /><figcaption>Figure 14. Distribution of Lypapunov exponents $\lambda_\infty$ predicted from unified theory (6.7) as a function of dimensionless folding frequency $f$ and net incremental stretching $r$. The Lyapunov exponent for granular packings and pore networks are located on the black and red lines, respectively, and random loose granular packing and the random pore network are indicated by a grey and orange dots, respectively.</figcaption></figure>

In § 2, continuous and discrete porous media were shown to be topologically equivalent by considering continuous porous media as discrete porous media with finite-sized connections between grains. Although this topological equivalence involves inversion of the fluid and solid phases, this does not impact chaotic mixing because it is generated by the total negative Gaussian curvature $K$ of the fluid–solid boundary $\delta\Omega$ in both classes of media. These media do however differ in that $\delta\Omega$ in continuous media has mostly an overall negative curvature and so admits saddle points, whereas discrete media has mostly positive curvature with node points, except at cusp-shaped contacts which have divergent negative curvature. This difference has significant implications for generation of chaotic mixing as saddle points inherently generate chaotic mixing, whereas node and contact points alone do not.

In § 3, both classes of media were shown to exhibit similar fluid stretching dynamics, controlled by interior stable and unstable 2-D manifolds that intersect transversely along 1-D critical lines that connect with boundary saddle points (continuous media) or contact points (discrete media). It was also established, contrary to prior understanding (Turuban et al. 2019), that while contact points between grains in discrete media are degenerate (in that the local velocity gradient tensor has a zero eigenvalue), in concert with node points, they can still generate hyperbolic 2-D manifolds and chaotic mixing. This establishes that chaotic mixing is ubiquitous in discrete porous media as contact points are inherent to these materials.

The mechanisms that drive fluid folding in both discrete and continuous porous media were uncovered in § 4. For discrete porous media, hold-up of fluid at contact points generates strong folding during advection, which manifests as weak folding transverse to the mean flow direction due to symmetry breaking. This weak folding is then amplified by fluid stretching and compression, leading to the formation of highly striated material distributions. For continuous porous media, strong folding is generated by fluid hold-up at
critical lines in pore branches. Although fluid mixing is continuous in the 3-D pore space, when projected to downstream 2-D planes transverse to the mean flow, this hold-up and strong folding manifests as discontinuous mixing, which is explored in § 5. Weak folding also occurs in continuous media due to advection over pore branches and mergers, which is also amplified by stretching and compression to form striated fluid distributions.

The impacts of discontinuous mixing due to fluid hold-up at pore branches in continuous media were considered in § 5. Although similar hold-up occurs in discrete porous media at contact points, their infinitesimal size means that they do not generate finite-sized discontinuities. Fluid elements in continuous media undergo both continuous (SF) and discontinuous (CS) actions, which may be summarised as ‘cut-stretch-fold-rotate’ in pore branches and ‘compress-fold-glue-rotate’ in pore mergers, leading to generation of a web of discontinuities along which cutting and shuffling occur. For systems that undergo CS only, complete mixing occurs if this web becomes space filling. We show that growth of the web of discontinuities can be accurately predicted by fluid stretching models, and highly stretched fluid elements are predominantly packed into finite-sized pores by CS rather than folding. As CS is mediated by SF, evolution of the mix-norm measure (Mathew et al. 2005) in continuous media is accurately predicted by estimates (5.4) based solely on the Lyapunov exponent.

In § 6, a generalised model for the dimensionless Lyapunov exponent in discrete porous media was presented. This model (6.6) encodes the fluid stretching mechanisms uncovered in § 3, and agrees well with numerical and experimental observations. Consequently, a unified model for the dimensionless Lyapunov exponent in both discrete and continuous porous media was also developed, which comprises the product of the magnitude of stretching increments *r* and their dimensionless frequency *f*. For continuous media, the frequency is unity, but the stretching increment varies due to the orientation of the pore branches and mergers in the network. Conversely, for discrete media, the stretching increment over a contact is fixed, but the frequency varies with the length of critical lines that connect contact points. These specific models cover a limited range of the parameter space {*f, r*}, whereas more complex media such as heterogeneous materials may span this space.

As such, the Lyapunov exponent can be used to quantify mixing in all porous media (figure 14) and provide a basis for optimisation of mixing and transport in engineered porous materials. Thus, predictive theories (Lester et al. 2013; Heyman et al. 2020) for the Lyapunov exponent play an important role in understanding and quantifying mixing in porous media, and recently developed characterisation techniques (Heyman et al. 2021) facilitate measurement of the Lyapunov exponent for a broad range of porous materials. In conjunction with previous studies (Lester et al. 2013, 2016b; Kree & Villermaux 2017; Turuban et al. 2018, 2019; Heyman et al. 2020; Souzy et al. 2020), these insights provide a complete description of the mechanisms that govern mixing in porous media, and so facilitate the development of improved theories of a wide range of physical, chemical and biological fluid-borne phenomena in porous media.

This study points towards fertile research directions involving experimental or computational methods. Although several experimental studies (Kree & Villermaux 2017; Heyman et al. 2020; Souzy et al. 2020; 2021) have focused on chaotic advection in discrete porous media, only limited experimental studies Lester & Chryss (2019) have examined mixing in continuous porous media. Examination of discontinuous mixing in these networks is an open problem, including verification of the mechanisms posited in this study. The extension to examine chaotic mixing in both continuous and discrete non-ideal porous media is an open question. Aside from the study of Heyman et al. (2021), no other known studies have resolved chaotic mixing in non-ideal porous materials, hence, the impacts of features such as surface roughness, grain/pore size and shape distributions
upon chaotic mixing are an open research topic. The results in this study also point to the exploration of novel porous architectures with tuneable mixing and transport properties. With the advent of rapid prototyping, the parameters controlling the mixing efficiency, as summarised in figure 14, could be optimised to develop and experimentally test pore-scale architectures with optimised transport properties for a variety of applications.

**Acknowledgements.** The authors thank the anonymous reviewers whose constructive comments have benefitted this paper.

**Funding.** This work was supported by the European Research Council (ERC) under the grant 101042466 (CHORUS).

**Declaration of interests.** The authors report no conflict of interest.

**Author contributions.** D.R.L. performed the simulations regarding discontinuous mixing, and D.R.L, J.H. and T.L-B. derived the theory. All authors contributed equally to analysing results and reaching conclusions, and in writing the paper.

# Appendix A. Critical point classification

The four different types of non-degenerate critical points $x_p$ on the pore boundary $\partial\Omega$ depicted in figure 3 are:

(i) separation point $\eta_3 > 0$, attractor for $u(x)$, $\eta_1 < 0$, $\eta_2 < 0$;

(ii) reattachment point $\eta_3 < 0$, repeller for $u(x)$, $\eta_1 > 0$, $\eta_2 > 0$;

(iii) separation point $\eta_3 > 0$, saddle for $u(x)$, $\eta_1 < 0$, $\eta_2 > 0$;

(iv) reattachment point $\eta_3 < 0$, saddle for $u(x)$, $\eta_1 < 0$, $\eta_2 > 0$.

Types (i) and (ii) critical points are node points, whereas types (iii) and (iv) critical points form saddle points. If, subject to technical conditions, the critical points are strongly hyperbolic (Surana, Grunberg & Haller 2006), the eigenvector directions indicated in figure 3 impart local exponential stretching of the fluid which persists away from the critical points as a continuation of these eigenvectors in the form of 1-D or 2-D stable $W^s$ and unstable $W^u$ hyperbolic manifolds. Hence, the stable and unstable manifolds from types (iii) and (iv) saddle points represent material surfaces that are respectively exponentially contracting or expanding in the fluid bulk. In addition, there also exist two types of degenerate critical points that arise at contact points:

(i) degenerate separation point $\eta_3 > 0$, attractor for $u(x)$, $\eta_1 < 0$, $\eta_2 = 0$;

(ii) degenerate reattachment point $\eta_3 < 0$, repeller for $u(x)$, $\eta_1 > 0$, $\eta_2 = 0$.

Although these points are degenerate with respect to the $\eta_2$ eigenvalue, they still exhibit exponential stretching/contraction in the $\eta_1$ and $\eta_3$ directions, as discussed in § 3.2.

# Appendix B. Fluid advection and stretching in porous networks

## B.1. *Advection mapping in continuous porous media*

As indicated by the angle $\theta = \Delta - \delta$ in figure 4, the 2-D stable and unstable manifolds emanating from the saddle points may either intersect smoothly ($\theta = 0$) or transversely ($\theta \neq 0$), leading respectively to either zero net stretching or chaotic mixing. In Appendix B, we show how simple but accurate maps of advective mixing in the pore branch and merger can be generated, leading to accurate predictions of fluid stretching in pore networks, as quantified by the dimensionless (infinite-time) Lyapunov exponent $\lambda_\infty \equiv \hat{\lambda}_\infty \ell / \langle v_1 \rangle$, where
ℓ the characteristic length scale of the medium (such as the length of a pore branch/merger) and 〈v₁〉 is the mean longitudinal velocity. The expression for λ∞ in ordered porous media with fixed δ, Δ is then

$$ \lambda_{\infty}(\delta, \Delta) = \ln|\xi + \sqrt{\xi^2 - 1}|, \quad \xi = \frac{9}{8} \cos \delta - \frac{1}{8} \cos(2\Delta + \delta), \qquad (B1) $$

where λ∞ = 0 for |ξ| ≤ 1, and maximum deformation (λ∞ = ln 2) occurs for Δ = (n + 1/2)π, δ = nπ, n = 0, 1, ... Chaotic mixing also occurs in random porous networks with uniformly distributed reorientation angles (Δ, δ ~ U[0, 2π]), where the Lyapunov exponent is given by the ensemble average

$$ \bar{\lambda}_{\infty} \equiv \frac{1}{4\pi^2} \int_{0}^{2\pi} \int_{0}^{2\pi} \lambda_{\infty}(\delta, \Delta) \, d\delta \, d\Delta \approx 0.1178, \qquad (B2) $$

which agrees very well with direct numerical simulations (Lester et al. 2013). Although the mechanisms that generate such exponential stretching in continuous porous media are now well understood, there exist much that is unknown regarding porous media more broadly, as summarised in questions (i–vi) which are addressed in the following sections.

This basic mechanism of fluid stretching has been used to develop a simple yet remarkably accurate model of advection in pore networks. Numerical simulations (Lester et al. 2014b, 2016a) of Stokes flow through the pore branch shown in figure 4(a) have shown that particle advection from a pore branch inlet to either outlet is well-approximated by the dimensionless spatial map M = M_n ∘ M_a,

$$ \mathcal{M} : (x_r, y_r) \mapsto \begin{cases} (x_r, 2y_r - \sqrt{1 - x_r^2}) & \text{if } y_r > 0, \\ (x_r, 2y_r + \sqrt{1 - x_r^2}) & \text{if } y_r \le 0, \end{cases} \qquad (B3) $$

which comprises affine M_a and non-affine M_n parts

$$ M_a : (x_r, y_r) \mapsto (x_r, 2y_r), \qquad (B4) $$

$$ M_n : (x_r, y_r) \mapsto \begin{cases} (x_r, y_r - \sqrt{1 - x_r^2}) & \text{if } y_r > 0, \\ (x_r, y_r + \sqrt{1 - x_r^2}) & \text{if } y_r \le 0, \end{cases} \qquad (B5) $$

where (x_r, y_r) are the local transverse coordinates within a pore branch inlet or outlet, such that x_r^2 + y_r^2 = 0, 1 respectively correspond to the pore centre and boundary. Note these maps are not area preserving due to 3-D effects, but the composite is area preserving. The inverse spatial map M⁻¹ = M_a⁻¹ ∘ M_n⁻¹ accurately approximates advection of fluid particles over a pore merger and the temporal T* map also accurately approximates the advection time t of fluid particles over a pore branch as

$$ T^* = T(\mathcal{M}), \quad T : t \mapsto t + \frac{1}{1 - x_r^2 - y_r^2}, \qquad (B6) $$

whereas the inverse temporal map T(M⁻¹) approximates advection of fluid particles over a pore merger. These maps may be extended to pore branches and mergers with respective orientations φ_b, φ_m in the x_r - y_r plane via the reoriented spatial maps

$$ S_b(\phi_b) = R(\phi_b) \circ M \circ R^{-1}(\phi_b), \qquad (B7) $$

$$ S_m(\phi_m) = R(\phi_m) \circ M^{-1} \circ R^{-1}(\phi_m), \qquad (B8) $$

where R is the rotation operator about the pore centre. The reorientation angle Δ between the pore branch and merger in figure 4 is Δ ≡ φ_m - φ_b. Rapid advection of fluid particles
in ordered and random model porous networks is performed via the map S over a coupled pore branch and merger (termed a couplet)

$$S = S_b(\varphi_b) \circ S_m(\varphi_m) = R(\varphi_b) \circ D(\Delta) \circ R^{-1}(\varphi_b), \qquad (B9)$$

where $D(\Delta) = M \circ R(\Delta) \circ M^{-1} \circ R^{-1}(\Delta)$. For a series of n concatenated couplets, the composite map $\Lambda_n$ is then

$$\Lambda_n = R(\varphi_{b,n}) \circ \left( \prod_{j=1}^{n} L(\delta_j, \Delta_j) \right) \circ R^{-1}(\varphi_{b,n}), \qquad (B10)$$

where $L(\delta, \Delta) \equiv R(\delta) \circ D(\Delta)$ and $\delta$ is the angle between couplets: $\delta_j = \varphi_{b,j+1} - \varphi_{b,j}$ for $j = 1:n-1$ and $\delta_n = \varphi_{b,n} - \varphi_{b,n}$. Previous studies (Lester et al. 2014a,b) have considered both random pore networks with a uniform distribution of the reorientation angles $\delta_j \sim U[0, 2\pi]$, $\Delta_j \sim U[0, 2\pi]$ and ordered pore networks with deterministic sequences of angles $\delta_j, \Delta_j$. Examples of application of these maps to simulate mixing in both ordered and random pore networks is shown in figure 10.

## B.2. Lyapunov models in continuous porous media

These advective maps can also generate accurate predictions of the fluid stretching in continuous porous media, as quantified by the dimensionless (infinite-time) Lyapunov exponent $\lambda_\infty$. Linearisation of $M$ in (B10) yields (Lester et al. 2013) the expression (B1). Hence, ordered 3-D porous networks represent extreme cases with respect to fluid stretching and deformation. As shown in figure 10, while a large class ($|\xi| \le 1$) of ordered media do not exhibit chaotic advection ($\lambda_\infty = 0$), some ordered networks exhibit maximal stretching ($\lambda_\infty = \ln 2$) associated with the baker's map. For random porous networks, the Lyapunov exponent is given by (B2). For such random pore networks, persistent fluid stretching arises because material lines undergoing stretching rotate towards the stretching direction which amplifies stretching, whereas material lines undergoing contraction rotate away from the contracting direction, retarding contraction. The Lyapunov exponent $\bar{\lambda}_\infty \approx 0.1178$ is thus a manifestation of the asymmetry between these stretching and contraction processes.

## Appendix C. Impact of pore boundary geometry on critical points

The different geometry of smoothed discrete and continuous porous media is reflected by the distribution of local Gaussian curvature $K = \kappa_1\kappa_2$ (where $\kappa_1 = 1/r_2$, $\kappa_2 = 1/r_2$ are the principal (maximum and minimum) local curvatures with corresponding radii $r_1, r_2$) over the surfaces shown in figure 5. The Gaussian curvature $K$ of a plane and a cylinder is zero (as at least one of the principal curvatures is zero), but is positive for a sphere and negative for a hyperbolic surface such as a saddle. However, as these fundamental elements share the same topology, their total curvature is equivalent, as reflected by the Gauss-Bonnet theorem,

$$\int_{\delta\Omega} K \, dA + \int_{\delta\delta\Omega} k_g \, ds = 2\pi \chi(\delta\Omega) = -4\pi, \qquad (C1)$$

where $k_g$ is the geodesic curvature of the 1-D boundary $\delta\delta\Omega$. For the pore branch/merger (figure 5c) and the smoothed grain (figure 5b), the net contribution of $\delta\delta\Omega$ to (C1) is zero due to periodicity, but for the discrete grain (figure 5a), the boundary $\delta\delta\Omega$ terminates at the contact points $x_c$, hence, there is a contribution of angle $-2\pi$ to (C1) at each contact point (Grinfeld 2013) due to the cusp-shaped angle $\pi$ formed with each contacting
sphere. As the spherical discrete grain in figure 5(a) has uniform Gaussian curvature $K = 1/r^2$, the first integral in (C1) is $4\pi$, whereas the second integral is $-8\pi$. Thus, whilst the discrete grains have positive curvature, grain contacts act as infinitesimal regions of infinite negative curvature, giving total negative curvature of $-4\pi$. The smoothed grain shown in figure 5(c) shows the finite contact regions have finite negative curvature, which yields the total negative curvature of $-4\pi$. Conversely, the branch/merger shown in figure 5(a) has predominantly negative or zero curvature.

Although the total curvature of these elements is preserved as per (C1), the distribution of Gaussian curvature $K$ impacts the number and type of critical points as saddle (node) points tend to arise in negative (positive) curvature regions. Hence, saddle points arise in continuous porous media near pore branch/merger junctions, whereas smoothed grains admit node points on the grain surface and saddle points near contact points.

Thus, the simplest arrangement of critical points for the branch/merger shown in figure 5(a) involves two saddle points $\mathbf{x}_p^s$ near the pore junctions, satisfying (2.3) as $\sum \gamma_p(\mathbf{x}_p^s) = -n_s = \chi(\delta\Omega)$. Conversely, the simplest arrangement of critical points for the grain with smoothed connections shown in figure 5(c) involves two node points $\mathbf{x}_p^n$ on the grain surface and two saddle points $\mathbf{x}_p^s$ per contact ($n_s = 4$ once periodicity of contacts is accounted for), satisfying (2.3) as $\sum \gamma_p(\mathbf{x}_p^s) = n_c - n_s = \chi(\delta\Omega)$. Therefore, despite topological equivalence, the different geometry of smoothed discrete porous media can generate different types of critical points than continuous porous media.

REFERENCES

DE ANNA, P., DENTZ, M., TARTAKOVSKY, A. & LE BORGNE, T. 2014 The filamentary structure of mixing fronts and its control on reaction kinetics in porous media flows. *Geophys. Res. Lett.* **41**, 4586–4593.

AREF, H., et al. 2017 Frontiers of chaotic advection. *Rev. Mod. Phys.* **89** (2), 025007.

ARMSTRONG, RYAN T., MCCLURE, J.E., ROBINS, V., LIU, Z., ARNS, C.H., SCHLÜTER, S. & BERG, S. 2019 Porous media characterization using Minkowski functionals: theories, applications and future directions. *Trans. Porous Med.* **130** (1), 305–335.

BAKKER, P.G. 1991 *Bifurcations in Flow Patterns*. vol. 2. Springer.

BERKOWITZ, B., DROR, I., HANSEN, S.K. & SCHER, H. 2016 Measurements and models of reactive transport in geological media. *Rev. Geophys.* **54** (4), 930–986.

BIJELJIC, B., MUGGERIDGE, A.H. & BLUNT, M.J. 2004 Pore-scale modeling of longitudinal dispersion. *Water Resour. Res.* **40** (11), W11501.

BRØNS, M. & HARTNACK, J.N. 1999 Streamline topologies near simple degenerate critical points in two-dimensional flow away from boundaries. *Phys. Fluids* **11** (2), 314–324.

CHRISTOV, I.C., LUEPTOW, R.M. & OTTINO, J.M. 2011 Stretching and folding versus cutting and shuffling: an illustrated perspective on mixing and deformations of continua. *Am. J. Phys.* **79** (4), 359–367.

DENTZ, M., LE BORGNE, T., ENGLERT, A. & BIJELJIC, B. 2011 Mixing, spreading and reaction in heterogeneous media: a brief review. *J. Cont. Hydrol.* **120-121**, 1–17.

EL BIED, A., SULEM, J. & MARTINEAU, F. 2002 Microstructure of shear zones in Fontainebleau sandstone. *Intl J. Rock Mech. Min.* **39** (7), 917–932.

GRAMLING, C.M., HARVEY, C.F. & MEIGS, L.C. 2002 Reactive transport in porous media: a comparison of model prediction with laboratory visualization. *Environ. Sci. Technol.* **36** (11), 2508–2514.

GRINFELD, P. 2013 *Introduction to Tensor Analysis and the Calculus of Moving Surfaces*. Springer.

HALLER, G. & MEZIC, I. 1998 Reduction of three-dimensional, volume-preserving flows with symmetry. *Nonlinearity* **11** (2), 319.

HEYMAN, J., LESTER, D.R. & LE BORGNE, T. 2021 Scalar signatures of chaotic mixing in porous media. *Phys. Rev. Lett.* **126**, 034505.

HEYMAN, J., LESTER, D.R., TURUBAN, R., MÉHEUST, Y. & LE BORGNE, T. 2020 Stretching and folding sustain microscale chemical gradients in porous media. *Proc Natl Acad Sci USA* **117** (24), 13359–13365.

HUANG, J., HUANG, K., YOU, X., LIU, G., HOLLETT, G., KANG, Y., GU, Z. & WU, J. 2018 Evaluation of tofu as a potential tissue engineering scaffold. *J. Mater. Chem. B* **6**, 1328–1334.

HUANG, J.-H., KIM, J., AGRAWAL, N., SUDARSAN, A.P., MAXIM, J.E., JAYARAMAN, A., UGAZ, V.M. 2009 Rapid fabrication of bio-inspired 3-D microfluidic vascular networks. *Adv. Mater.* **21** (35), 3567–3571.
JOHN, T. & MEZIĆ, I. 2007 Maximizing mixing and alignment of orientable particles for reaction enhancement. *Phys. Fluids* **19** (12), 123602.

JUAREZ, G., LUEPTOW, R.M., OTTINO, J.M., STURMAN, R. & WIGGINS, S. 2010 Mixing by cutting and shuffling. *EPL (Europhys. Lett.)* **91** (2), 20003.

KÁROLYI, G., PÉNTEK, Á., SCHEURING, I., TÉL, T. & TOROCZKAI, Z. 2000 Chaotic flow: the physics of species coexistence. *Proc. Natl Acad. Sci.* **97** (25), 13661–13665.

KÁROLYI, G., SCHEURING, I. & CZÁRÁN, T. 2002 Metabolic network dynamics in open chaotic flow. *Chaos* **12** (2), 460–469.

KREE, M. & VILLERMAUX, E. 2017 Scalar mixtures in porous media. *Phys. Rev. Fluids* **2**, 104502.

LESTER, D., METCALFE, G. & RUDMAN, M. 2014a Control mechanisms for the global structure of scalar dispersion in chaotic flows. *Phys. Rev. E* **90** (2), 022908.

LESTER, D.R., TREFRY, M.G. & METCALFE, G. 2016a Chaotic advection at the pore scale: Mechanisms, upscaling and implications for macroscopic transport. *Adv. Water Resour.* **97**, 175–192.

LESTER, D.R. & CHRYSS, A. 2019 Topological mixing of yield stress materials. *Phys. Rev. Fluids* **4**, 064502.

LESTER, D.R., DENTZ, M. & LE BORGNE, T. 2016b Chaotic mixing in three-dimensional porous media. *J. Fluid Mech.* **803**, 144–174.

LESTER, D.R., METCALFE, G. & TREFRY, M.G. 2013 Is chaotic advection inherent to porous media flow? *Phys. Rev. Lett.* **111**, 174101.

LESTER, D.R., METCALFE, G. & TREFRY, M.G. 2014b Anomalous transport and chaotic advection in homogeneous porous media. *Phys. Rev. E* **90**, 063012.

MACKAY, R.S. 1994 Transport in 3-D volume-preserving flows. *J. Nonlinear Sci.* **4** (1), 329.

MATHEW, G., MEZIĆ, I. & PETZOLD, L. 2005 A multiscale measure for mixing. *Physica D: Nonlinear Phenom.* **211** (1), 23–46.

MELCHELS, F.P.W., FEIJEN, J. & GRIPMA, D.W. 2009 A poly(d,l-lactide) resin for the preparation of tissue engineering scaffolds by stereolithography. *Biomaterials* **30** (23), 3801–3809.

NEUFELD, Z. & HERNANDEZ-GARCIA, E. 2009 *Chemical and Biological Processes in Fluid Flows: A Dynamical Systems Approach*. Imperial College Press.

OTT, E. 2002 *Chaos in Dynamical Systems*. 2nd edn. Cambridge University Press.

OTTINO, J.M. 1989 *The Kinematics of Mixing: Stretching, Chaos, and Transport*. Cambridge University Press.

OUELLETTE, N.T., O'MALLEY, P.J.J. & GOLLUB, J.P. 2008 Transport of finite-sized particles in chaotic flow. *Phys. Rev. Lett.* **101** (17), 174504.

ROLLE, M. & LE BORGNE, T. 2019 Mixing and reactive fronts in the subsurface. *Rev. Mineral. Geochem.* **85** (1), 111–142.

SAPSIS, T. & HALLER, G. 2010 Clustering criterion for inertial particles in two-dimensional time-periodic and three-dimensional steady flows. *Chaos* **20** (1), 017515.

SCHOLZ, C., WIRNER, F., GÖTZ, J., RÜDE, U., SCHRÖDER-TURK, G.E., MECKE, K. & BECHINGER, C. 2012 Permeability of porous materials determined from the euler characteristic. *Phys. Rev. Lett.* **109**, 264504.

SMITH, L.D., RUDMAN, M., LESTER, D.R. & METCALFE, G. 2016 Mixing of discontinuously deforming media. *Chaos* **26** (2), 023113.

SMITH, L.D., RUDMAN, M., LESTER, D.R. & METCALFE, G. 2017a Impact of discontinuous deformation upon the rate of chaotic mixing. *Phys. Rev. E* **95**, 022213.

SMITH, L.D., UMBANHOWAR, P.B., LUEPTOW, R.M. & OTTINO, J.M. 2019 The geometry of cutting and shuffling: an outline of possibilities for piecewise isometries. *Phys. Rep.* **802**, 1–22.

SMITH, L.D., UMBANHOWAR, P.B., OTTINO, J.M. & LUEPTOW, R.M. 2017b Mixing and transport from combined stretching-and-folding and cutting-and-shuffling. *Phys. Rev. E* **96**, 042213.

SOUZY, M., LHUISSIER, H., MÉHEUST, Y., LE BORGNE, T. & METZGER, B. 2020 Velocity distributions, dispersion and stretching in three-dimensional porous media. *J. Fluid Mech.* **891**, A16.

STOCKER, R. 2012 Marine microbes see a sea of gradients. *Science* **338**, 6107.

STURMAN, R. 2012 The role of discontinuities in mixing. *Adv. Appl. Mech.* **45**, 51–90.

SUDARSAN, A.P. & UGAZ, V.M. 2006 Multivortex micromixing. *Proc. Natl Acad. Sci.* **103** (19), 7228–7233.

SURANA, A., GRUNBERG, O. & HALLER, J.A. 2006 Exact theory of three-dimensional flow separation. Part 1. Steady separation. *J. Fluid Mech.* **564**, 57–103.

TEL, T., DE MOURA, A., GREBOGI, C. & KÁROLYI, G. 2005 Chemical and biological activity in open flows: a dynamical system approach. *Phys. Rep.* **413**, 91–196.

THERRIAULT, D., WHITE, S.R. & LEWIS, J.A. 2003 Chaotic mixing in three-dimensional microvascular networks fabricated by direct-write assembly. *Nat. Mater.* **2** (4), 265–271.

THIFFEAULT, J.-L. 2004 Stretching and curvature of material lines in chaotic flows. *Physica D: Nonlinear Phenom.* **198** (3), 169–181.
TURUBAN, R., LESTER, D.R., HEYMAN, J., BORGNE, T.L. & MÉHEUST, Y. 2019 Chaotic mixing in crystalline granular media. *J. Fluid Mech.* **871**, 562–594.

TURUBAN, R., LESTER, D.R., LE BORGNE, T. & MÉHEUST, Y. 2018 Space-group symmetries generate chaotic fluid advection in crystalline granular media. *Phys. Rev. Lett.* **120**, 024501.

VALOCCHI, A.J., BOLSTER, D. & WERTH, C.J. 2019 Mixing-limited reactions in porous media. *Trans. Porous Med.* **130** (1), 157–182.

VOGEL, H.-J. 2002 *Topological Characterization of Porous Media*, pp. 75–92. Springer.

WRIGHT, S.F., ZADRAZIL, I. & MARKIDES, C.N. 2017 A review of solid-fluid selection options for optical-based measurements in single-phase liquid, two-phase liquid-liquid and multiphase solid-liquid flows. *Exp. Fluids* **58**, 108.