2.1. Background

To address challenges (i)–(vi), in this section (and Appendices A, B) we first briefly review the theory (Lester et al. Reference Lester, Metcalfe and Trefry2013, Reference Lester, Dentz and Le Borgne2016b ) 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 $\boldsymbol{v}(\boldsymbol{x})$ as the steady divergence-free three-dimensional (3-D) fluid velocity field in the pore space $\varOmega$ of the continuous porous medium, which satisfies the no-slip condition $\boldsymbol{v}=\boldsymbol{0}$ on the entire fluid/solid boundary $\partial \varOmega$ and is governed by steady Stokes flow driven by a mean pressure gradient

(2.1) \begin{equation} \mu {\nabla }^2\boldsymbol{v}(\boldsymbol{x})-{\boldsymbol{\nabla }} p=0,\quad \forall \boldsymbol{x}\in \varOmega , \end{equation}

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

MacKay (Reference MacKay1994) shows that critical points $\boldsymbol{x}_p$ of the skin friction field (where $\boldsymbol{u}(\boldsymbol{x}_p)=\boldsymbol{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 $\boldsymbol{A}\equiv \partial \boldsymbol{u}/\partial \boldsymbol{x}$ at $\boldsymbol{x}=\boldsymbol{x}_p$ are all non-zero. Without loss of generality, we denote the eigenvalues tangent to $\partial \varOmega$ as $\eta _1$ , $\eta _2$ (with $\eta _1\leqslant \eta _2$ ), and $\eta _3$ is the interior eigenvalue whose eigenvector is normal to $\partial \varOmega$ . The divergence-free condition means that the eigenvalues $\eta _i$ satisfy (MacKay Reference MacKay1994)

(2.2) \begin{equation} \eta _1+\eta _2+2\eta _3=0, \end{equation}

Figure 3. Schematic of the structure of the skin friction field $\boldsymbol{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 $\partial \varOmega$ and the associated stable $\mathcal{W}^s$ and unstable $\mathcal{W}^u$ 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 $\eta _1+\eta _2+2\eta _3=0$ . Adapted from Lester, Dentz & Le Borgne (Reference Lester, Dentz and Le Borgne2016b ).

as $\sum _{i=1}^3\eta _i={\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}=-\eta _3$ . Therefore, non-degenerate critical points $\boldsymbol{x}_p$ consist of one stable and two unstable eigenvalues or vice versa, and form either reattachment ( $\eta _3\lt 0$ ) or separation ( $\eta _3\gt 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\lt 0$ , $\eta _2\gt 0$ generate 2-D stable $\mathcal{W}^s_{2\text{-D}}$ and unstable $\mathcal{W}^u_{2\text{-D}}$ 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 Reference Haller and Mezic1998), 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 Reference Ott2002). In practice, a single transverse intersection of stable $\mathcal{W}^s_{2\text {-D}}$ and unstable $\mathcal{W}^u_{2\text {-D}}$ 2-D manifolds implies a chaotic tangle of infinitely many (Ottino Reference Ottino1989), leading to strong fluid stretching and folding, and chaotic mixing. MacKay (Reference MacKay1994) 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.

Figure 4. Schematic of a (a) pore branch and (b) merger in continuous porous media. Non-degenerate critical points $\boldsymbol{x}_p$ (black dots) generate 2-D hyperbolic stable $\mathcal{W}^s_{2\text{-D}}$ and unstable $\mathcal{W}^u_{2\text{-D}}$ manifolds (grey) which are surfaces of locally minimum transverse flux. The angles $\varDelta$ , $\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 $\mathcal{M}$ , $\mathcal{M}^{-1}$ defined in (B3).

Thus, only saddle points give rise to interior 2-D hyperbolic manifolds and chaotic mixing. The prevalence of saddle points on $\partial \varOmega$ is related to the topological complexity of the continuous porous medium, given by the average number density $\rho (\chi )$ of the Euler characteristic $\chi (\varOmega )$ , which is strongly negative (Vogel Reference Vogel2002; Scholz et al. Reference Scholz, Wirner, Götz, Rüde, Schröder-Turk, Mecke and Bechinger2012), reflecting the topological complexity inherent to all porous materials. For closed bounded manifolds $\varOmega$ , the Euler characteristic of $\varOmega$ is related (Armstrong et al. Reference Armstrong, McClure, Robins, Liu, Arns, Schlüter and Berg2019) to that of its boundary $\delta \varOmega$ as $\chi (\delta \varOmega )=2\chi (\varOmega )$ , and hence the Euler characteristic of the pore boundary is also strongly negative. The Poincaré–Hopf theorem connects the Euler characteristic $\chi (\delta \varOmega )$ of the pore boundary to the sum of indices $\gamma _p$ of stagnation points $\boldsymbol{x}_p$ in $\delta \varOmega$ as

(2.3) \begin{equation} \chi (\delta \varOmega )=2(1-g)=\sum \gamma _p(\boldsymbol{x}_p)=n_n-n_s, \end{equation}

where $g$ is the topological genus of the pore boundary, $\gamma _p=+1$ for node points and $\gamma _p=-1$ for saddle points, and so $n_n$ and $n_s$ respectively denote the number of node and saddle points, and $|\rho (\chi )|$ provides a lower bound for the number density of saddle points which naturally arise at pore branches and mergers (figure 4). Although the specific 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 \boldsymbol{\cdot }\rangle$ of the infinite-time limit $\hat {\lambda }_\infty (\boldsymbol{X})$ of the finite-time Lyapunov exponent $\lambda (\boldsymbol{X},t)$ as

(2.4) \begin{equation} \lambda _\infty \equiv \langle \hat {\lambda }_\infty (\boldsymbol{X})\rangle =\left \langle \lim _{t\rightarrow \infty }\lambda (\boldsymbol{X},t)\right \rangle , \end{equation}

where $\boldsymbol{X}$ denotes Lagrangian space. For ordered porous media, there can exist a mixture of distinct regions of regular (non-chaotic) flow (where $\hat {\lambda }(\boldsymbol{X})=0$ ) and chaotic mixing (where $\hat {\lambda }(\boldsymbol{X})\gt 0$ ), in which case $\lambda _\infty$ represents a volume-weighted average of $\hat {\lambda }(\boldsymbol{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 (\boldsymbol{X})$ is invariant and so $\lambda _\infty =\lambda _\infty (\boldsymbol{X})$ .

Figure 5. Topological equivalence of the 2-D pore boundary $\delta \varOmega$ separating the fluid (pore) $\varOmega$ and solid domains of the fundamental elements of (a) discrete and (c) continuous porous media. The normal vector $\boldsymbol{n}$ indicates the normal vector pointing into the fluid domain (pore) $\varOmega$ from $\delta \varOmega$ , and $\delta \delta \varOmega$ (green lines) is the 1-D boundary of the pore boundary $\delta \varOmega$ . $\delta \varOmega$ is coloured according to its local Gaussian curvature $K$ . (a) Pore boundary $\delta \varOmega$ 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 \varOmega$ with finite-sized connections between contacting grains, forming boundaries $\delta \delta \varOmega$ . (c) Pore boundary $\delta \varOmega$ for a connected pore branch and merger associated with continuous porous media.

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 5 c) form an extensive 3-D pore network similar to those shown in figure 1(a–e), whereas assemblies of grains (figure 5 a) 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 complex in that the topological genus $g$ of the pore-boundary $\delta \varOmega$ is large, guaranteeing the existence of saddle points on the pore boundary $\delta \varOmega$ that generate chaotic mixing in the pore space $\varOmega$ . For porous networks, the topological genus $g$ and Euler characteristic $\chi (\delta \varOmega )$ of the pore network boundary are given by the total number $n_b$ of pore bifurcations as

(2.5) \begin{equation} g=1-\frac {\chi (\delta \varOmega )}{2}=1+\frac {n_b}{2}. \end{equation}

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 \varOmega )=-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. Reference Turuban, Lester, Heyman, Borgne and Méheust2019) 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 \varOmega$ of the latter is non-smooth, as the granular matter involves cusp-like contacts between discrete grains, regardless of grain shape and orientation.

Under periodic boundary conditions, the topological genus $g$ and Euler characteristic $\chi (\delta \varOmega )$ of the grain boundary $\delta \varOmega$ of granular assemblies is given in terms of the total number of grain contacts $n_c$ as (Turuban et al. Reference Turuban, Lester, Heyman, Borgne and Méheust2019)

(2.6) \begin{equation} g=1-\frac {\chi (\delta \varOmega )}{2}=n_c. \end{equation}

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 $\chi (\delta \varOmega )=-2$ . From (2.6), granular assemblies are also topologically complex in that they posses 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 $\delta \varOmega$ is smooth and continuous. This smoothed grain has the same topology ( $g=2$ , $\chi (\delta \varOmega )=-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 $\boldsymbol{n}$ ) is external to the pore boundary $\delta \varOmega$ 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 $\delta \varOmega$ 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 $\boldsymbol{u}$ on the boundary $\delta \varOmega$ of smoothed discrete porous media, and so ensures the existence of saddle points and 2-D hyperbolic manifolds regardless of the orientation of $\boldsymbol{n}$ . As the pore boundary $\delta \varOmega$ 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. (Reference Turuban, Lester, Le Borgne and Méheust2018, Reference Turuban, Lester, Heyman, Borgne and Méheust2019) and Heyman et al. (Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020) 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.

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 $\varDelta =\pi /2$ and the associated 1-D and 2-D stable and unstable manifolds which emanate from the saddle points $\boldsymbol{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.

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. (Reference Turuban, Lester, Heyman, Borgne and Méheust2019), 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.

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. Reference Turuban, Lester, Heyman, Borgne and Méheust2019) (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 $\boldsymbol{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.

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 5 c) 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 $\boldsymbol{v}(\boldsymbol{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 $\boldsymbol{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\rightarrow 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.

Figure 7. Schematic of evolving stable $\mathcal{W}^S_{2\text{-D}}$ (blue surfaces) and unstable $\mathcal{W}^U_{2\text{-D}}$ (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).

Figure 8. Streamlines (thin) and critical lines (thick) for 2-D velocity field $\boldsymbol{v}_{2\text{-D}}(\boldsymbol{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 $\boldsymbol{x}_p^s$ or a degenerate topological saddle $\boldsymbol{x}_p^d$ . Due to the no-slip condition, in both cases, the skin friction field $\boldsymbol{u}(\boldsymbol{x})\equiv \partial \boldsymbol{v}/\partial x_3^\prime$ on the grain boundaries local to this contact shares the same topology as the velocity field $\boldsymbol{v}_{2\text{-D}}(\boldsymbol{x})$ .

Turuban et al. (Reference Turuban, Lester, Heyman, Borgne and Méheust2019) proposed that this critical point is a topological saddle (Bakker Reference Bakker1991; Brøns & Hartnack Reference Brøns and Hartnack1999), which is degenerate (i.e. the skin friction gradient tensor $\boldsymbol{A}$ at $\boldsymbol{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 Reference Brøns and Hartnack1999), 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

(3.1) \begin{equation} \chi (\varOmega )=2(1-g)=n_n-n_s-2n_d. \end{equation}

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 characteristic $\chi (\varOmega )$ ; $n_s=n_n-2$ . Hence, it was proposed by Turuban et al. (Reference Turuban, Lester, Heyman, Borgne and Méheust2019) 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 $\mathcal{W}^s_{2\text{-D}}$ and unstable $\mathcal{W}^u_{2\text{-D}}$ 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 6 b), whereas degenerate points are not associated with exponential (hyperbolic) stretching (Turuban et al. Reference Turuban, Lester, Heyman, Borgne and Méheust2019).

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. (Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020).

However, deeper inspection reveals that contact points are only degenerate with respect to the 2-D flow $\boldsymbol{v}_{2\text{-D}}(\boldsymbol{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 $\mathcal{W}_{2\text{-D}}^S$ that is shown in figure 7 without material (green) lines. Similarly, a downstream 2-D unstable manifold $\mathcal{W}_{2\text{-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 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.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. (Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020) 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. (Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020) and recent studies (Turuban et al. Reference Turuban, Lester, Le Borgne and Méheust2018, Reference Turuban, Lester, Heyman, Borgne and Méheust2019) 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 $\mathcal{S}_b$ , $\mathcal{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 $\boldsymbol{x}_p$ . The intersection of the 2-D stable manifold $\mathcal{W}_{2\text{-D}}^s$ with the boundary $\partial \varOmega$ represents a repelling (i.e. ${\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}\gt 0$ ) critical line $\zeta$ of the skin friction field $\boldsymbol{u}$ . Any injected material line that crosses the mid-plane $y_r=0$ of the inlet pore corresponding to $\mathcal{W}_{2\text{-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 $\mathcal{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 10. Evolution of fluid elements (black points) with pore number $n$ over the circular inlet plane of a pore branch (figure 4 a) 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 $\varDelta$ 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 $\mathcal{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 4 b), where fluid elements from separate pores are merged together over the attracting (i.e. ${\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}\lt 0$ ) critical line $\zeta$ formed by intersection of the 2-D unstable manifold $\mathcal{W}_{2\text{-D}}^u$ and the pore boundary $\partial \varOmega$ . 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 $\mathcal{M}$ also imparts weaker folding that is illustrated by its affine $\mathcal{M}_a$ (B4) and non-affine $\mathcal{M}_n$ (B5) components.

Fluid cutting is also encoded in the non-affine map $\mathcal{M}_n$ , such that fluid elements with $y_r\gt 0$ ( $y_r\lt 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 $\mathcal{M}_n$ and its inverse $\mathcal{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 $\mathcal{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.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:

1. (i) ordered network corresponding to 3-D baker’s map: $\varDelta =\pi /2$ , $\delta =0$ , $\lambda _\infty =\ln 2\approx 0.693$ ;

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

3. (iii) ordered network, non-chaotic mixing: $\varDelta =\pi /4$ , $\delta =\pi /4$ , $\lambda _\infty =0$

4. (iv) random network, chaotic mixing: $\varDelta \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 $\mathcal{M}$ , $\mathcal{M}^{-1}$ described in Appendix B for the various values of $\delta$ , $\varDelta$ 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 $\boldsymbol{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.

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. Reference Juarez, Lueptow, Ottino, Sturman and Wiggins2010; Christov, Lueptow & Ottino Reference Christov, Lueptow and Ottino2011; Sturman Reference Sturman2012; Smith et al. Reference Smith, Rudman, Lester and Metcalfe2016, Reference Smith, Umbanhowar, Ottino and Lueptow2017 Reference Smith, Rudman, Lester and Metcalfea , Reference Smith, Umbanhowar, Ottino and Lueptowb , Reference Smith, Umbanhowar, Lueptow and Ottino2019), such as granular tumblers (Christov et al. Reference Christov, Lueptow and Ottino2011; Smith et al. Reference Smith, Umbanhowar, Ottino and Lueptow2017b ), and microfluidic micro-mixers based on iterated pore branches and mergers (Sudarsan & Ugaz Reference Sudarsan and Ugaz2006). The motions imparted by the spatial maps $\mathcal{S}_b$ , $\mathcal{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 $\mathcal{M}_a$ ) by a factor of 2 in the $y_r$ direction, ‘fold’ refers to non-affine folding (by $\mathcal{M}_n$ ) of elements to conform to the outlet pore boundaries, ‘rotate’ refers to reorientation (by $R(\varDelta )$ ) of fluid elements due to orientation of the downstream pore mergers. In the pore merger, ‘compress’ refers to the compression (by $\mathcal{M}_a^{-1}$ ) of fluid elements by a factor of 1/2 in the $y_r$ direction, ‘fold’ refers to non-affine folding (by $\mathcal{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. Reference Smith, Umbanhowar, Ottino and Lueptow2017 Reference Smith, Rudman, Lester and Metcalfea , Reference Smith, Umbanhowar, Ottino and Lueptowb ). The natural mathematical framework to analyse CS motions is via piecewise isometries (PWIs) (Sturman Reference Sturman2012), 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 $\mathcal{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 dissipative or continuous mixing. The mix-norm measure (Mathew, Mezić & Petzold Reference Mathew, Mezić and Petzold2005) was developed to quantify multiscale mixing across diverse applications, including non-dissipative and discontinuous mixing. The mix-norm $\varPhi (c)\in [0,1]$ captures the extent of mixing of a zero mean concentration field $c(\boldsymbol{x})$ across different normalised length scales $s\in [0,1]$ (where $s=1$ corresponds to the velocity correlation length scale $\ell$ ) as

(5.1) \begin{align} \varPhi (c)&=\sqrt {\int _0^1\phi (c,s)^2 \mu ({\textrm d}s)}, \end{align}

(5.2) \begin{align} \phi (c,s)&=\sqrt {\int _{p\in \varOmega }d^2(c,\boldsymbol{p},s) \mu ({\textrm d}\boldsymbol{p})}, \end{align}

(5.3) \begin{align} d(c,\boldsymbol{p},s)&=\frac {1}{\text{vol}[B(\boldsymbol{p},s)]}\int _{x\in B(\boldsymbol{p},s)}c(\boldsymbol{x})\mu ({\textrm d}\boldsymbol{x}), \end{align}

where $\mu$ denotes the Lebesgue measure, $d(c,\boldsymbol{p},s)$ is the average of $c(\boldsymbol{x})$ over the ball $B(\boldsymbol{p},s):=\{\boldsymbol{x}:|\boldsymbol{x}-\boldsymbol{p}|\lt s/2\}$ centred at point $\boldsymbol{p}$ in the fluid domain $\varOmega$ , $\phi (c,s)$ is the $L^2$ -norm of $d$ over $\varOmega$ and $\varPhi (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(\boldsymbol{x})$ after averaging over length scale $s$ and the mix-norm $\varPhi (c)$ accounts for mixing across all scales by quantifying the $L^2$ norm of this measure over all $s$ .

Hence, $\varPhi =1$ for scalar distributions that are completely segregated at the largest length scale of the system and $\varPhi (c)=0$ for mixtures that are well mixed at all scales. Under weak ergodic mixing characteristic of CS systems, the mix-norm $\varPhi (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, $\varPhi (c_n)$ decays exponentially with $n$ . For SF systems that can be represented by a linear map such as $\mathcal{M}_a$ , the mix-norm decays exponentially at a rate of half the Lyapunov exponent (Mathew et al. Reference Mathew, Mezić and Petzold2005) as

(5.4) \begin{equation} \varPhi (c_n)=\varPhi (c_0)\exp \left (-\frac {\lambda _\infty n}{2}\right ). \end{equation}

Although many studies have considered mixing due to purely CS or SF actions, only a handful (Smith et al. Reference Smith, Rudman, Lester and Metcalfe2016, Reference Smith, Rudman, Lester and Metcalfe2017, Reference Smith, Umbanhowar, Lueptow and Ottino2019 Reference Smith, Rudman, Lester and Metcalfea , Reference Smith, Umbanhowar, Ottino and Lueptowb ) have considered the impact of coupled CS and SF motions. Smith et al. (Reference Smith, Rudman, Lester and Metcalfe2017a ,Reference Smith, Umbanhowar, Ottino and Lueptow 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. Reference Smith, Umbanhowar, Ottino and Lueptow2017b ), 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) 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 $\mathcal{S}_b$ , $\mathcal{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 $\mathcal{S}_b$ , $\mathcal{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 $\varPhi (c_n)$ in pore networks with chaotic and (inset) non-chaotic mixing, where points represent numerical results from maps $\mathcal{S}_b$ , $\mathcal{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).

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 $\mathcal{S}_b$ , $\mathcal{S}_m$ . For the non-chaotic case ( $\varDelta =\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\rightarrow \infty$ , then the growth of the total length $l_n$ of the web of discontinuities follows the same stretching process as the fluid phase.

As an additional cut is added to the web of discontinuities at each iteration of the map $\varLambda _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. Reference Lester, Dentz and Le Borgne2016b ) 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 (\varLambda _\infty n )$ , where $\varLambda _\infty =\lambda _\infty +\sigma _\lambda ^2/2$ and $\sigma _\lambda ^2$ is found (Lester et al. Reference Lester, Dentz and Le Borgne2016 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

(5.5) \begin{equation} l_n\approx 2\frac {\exp \left (\varLambda _\infty (n+1)\right )-1}{\exp \left (\varLambda _\infty \right )-1}. \end{equation}

Conversely, for the non-chaotic case $\varDelta =\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

(5.6) \begin{equation} l_n\approx \frac {(n+1)(n+4)}{2}. \end{equation}

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 $\varPhi (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 $\mathcal{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.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

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

where $\alpha \geqslant 1$ quantifies stretching along the unstable manifold. Note that the deformation represented by $F$ is equivalent to the deformation $\mathcal{D}(\varDelta )$ over a couplet in the pore network model (B9) with $\varDelta =\pi /2$ and $\alpha =2$ . After passing the contact, the un/stable manifolds are inverted (figure 9 a.3), which is reflected by a reorientation of $F$ of angle $\pi /2$ (transition from figure 9 a.2 to figure 9 a.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 9 a.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).

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$ .

The asymmetry that generates the perturbation angle $\delta$ also induces a small cusp of initial size $\epsilon$ (figure 12 a) 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 9 a.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

(6.2) \begin{equation} D_{2n}=\prod _{j=1}^n\left(F\boldsymbol{\cdot }\big[R(\varphi )\boldsymbol{\cdot }F\boldsymbol{\cdot }R(\varphi )^\top \big]\right), \end{equation}

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

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

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 \rightarrow 0$ , the eigenvalues $r_i$ are

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

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. Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020), then the elongation $\rho$ over longitudinal distance $X$ is

(6.5) \begin{equation} \rho =\exp \left ( \frac {X}{X_c} \ln (r) \right ), \end{equation}

and so the dimensionless Lyapunov exponent is

(6.6) \begin{equation} \lambda _{\infty }=\frac {d}{X_c}\ln (r), \end{equation}

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. Reference Turuban, Lester, Le Borgne and Méheust2018, Reference Turuban, Lester, Heyman, Borgne and Méheust2019). 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. Reference Turuban, Lester, Heyman, Borgne and Méheust2019). 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. Reference Turuban, Lester, Heyman, Borgne and Méheust2019). Note that Heyman et al. (Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020), 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

Figure 13. Comparison of predicted and computed (Turuban et al. Reference Turuban, Lester, Heyman, Borgne and Méheust2019) 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.

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

(6.7) \begin{equation} \lambda _{\infty }= f \ln (r), \end{equation}

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

(6.8) \begin{equation} r= \lim _{n\to \infty } \left \langle \frac {\rho _{n+1}}{\rho _{n}} \right \rangle . \end{equation}

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$ , $\varDelta$ (Lester et al. Reference Lester, Metcalfe and Trefry2013). 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. Reference Turuban, Lester, Le Borgne and Méheust2018).

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.

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.