2.1. Background

We first consider the kinematics of steady 3-D pore-scale flow before upscaling to Darcy flow. Steady 3-D pore-scale flow is described by the Stokes equations

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

subject to no-slip conditions $\hat {\boldsymbol{v}}|_{\partial \varOmega _\textit{fs}}=\boldsymbol{0}$ at the pore boundary $\partial \varOmega _\textit{fs}$ . Here, $\mu$ is the fluid viscosity, $\hat {\boldsymbol{v}}$ is the pore-scale velocity, $\hat {p}$ the pore-scale pressure field, $\varOmega _f$ , $\varOmega _s$ respectively are the fluid (pore) and solid (grain) domains. For steady 3-D flow, the topological complexity of the boundary $\partial \varOmega _\textit{fs}$ generates pore-scale chaotic advection (Lester et al. Reference Lester, Metcalfe and Trefry2013). Upscaling by suitable averaging of (2.1) (denoted by $\langle \boldsymbol{\cdot }\rangle _p$ ) yields the Darcy equation Bear Reference Bear1972)

(2.2) \begin{equation} \boldsymbol{v}(\boldsymbol{x})=-\boldsymbol{K}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{\nabla } \phi (\boldsymbol{x}),\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}=0,\,\,\,\,\boldsymbol{x}\in \varOmega , \end{equation}

where $\varOmega \equiv \langle \varOmega _f\cup \varOmega _s\rangle _p$ denotes the Darcy-scale porous medium, $\boldsymbol{v}\equiv \langle \hat {\boldsymbol{v}}\rangle _p$ is the (upscaled) Darcy-scale velocity, $\phi \equiv \langle \hat {p}\rangle _p$ the upscaled potential field and $\boldsymbol{K}(\boldsymbol{x})$ is the hydraulic conductivity tensor field that captures viscous drag due to no-slip conditions on $\partial \varOmega _\textit{fs}$ . As this pore boundary generates chaotic advection, omission of $\partial \varOmega _\textit{fs}$ in (2.2) via upscaling also omits these pore-scale kinematics. Hence, homogeneous Darcy flow (i.e. $\boldsymbol{K}(\boldsymbol{x})= \text{const.}$ ) is non-chaotic at the Darcy scale while the pore-scale flow is in general chaotic. At the Darcy scale, the effects of pore-scale chaotic advection must be incorporated via coarse-grained models of e.g. solute mixing (Lester et al. Reference Lester, Dentz and Le Borgne2016b ), dispersion (Puyguiraud, Gouze & Dentz Reference Puyguiraud, Gouze and Dentz2021) and reaction (Aquino et al. Reference Aquino, Le Borgne and Heyman2023). As shall be shown, transverse solute dispersion is especially relevant to understanding chaotic advection in Darcy flow. For solute dispersion, this pore-scale chaotic advection generates mechanical dispersion under pure advection, whereas for diffusive solutes, the interplay of molecular diffusion with pore-scale velocity fluctuations generates hydrodynamic dispersion, which is a function of the pore-scale Péclet number, $\textit{Pe}_p$ (Bear Reference Bear1972; Bijeljic & Blunt Reference Bijeljic and Blunt2006, Reference Bijeljic and Blunt2007; Liu et al. Reference Liu, Xiao, Aquino, Dentz and Wang2024), that compares the advection and diffusion time scales over a characteristic pore length. Hydrodynamic dispersion due to pore-scale velocity fluctuations is also referred to as local-scale dispersion. For many applications the pore-scale transverse dispersion coefficient $D_{T,p}$ is assumed to scale with the mean velocity as

(2.3) \begin{equation} D_{T,p}=D_m+\alpha _{T,p}||\boldsymbol{v}||^n, \end{equation}

where $D_m$ is molecular diffusivity, $\alpha _{T,p}$ is the pore-scale transverse dispersivity and the empirical index $n\in [1,2]$ is often taken as unity.

In this study we focus solely on heterogeneous conductivity fields as the majority of porous materials are heterogeneous at the Darcy scale (Bear Reference Bear1972; Cushman Reference Cushman2013). At scales significantly larger than the largest correlation length scale $\ell$ of the conductivity field $\boldsymbol{K}(\boldsymbol{x})$ , the impact of Darcy-scale velocity fluctuations on solute spreading are captured via the macro-dispersion concept (Gelhar & Axness Reference Gelhar and Axness1983; Dagan Reference Dagan1987).

On the continuum scale, solute transport is typically advection dominated, that is, advection dominates over local-scale dispersion on the correlation scale of hydraulic conductivity. This is measured by the Darcy-scale Péclet number, $Pe$ , which compares the advection and dispersion times over correlation scale $\ell$ . As shall be shown, the scaling of the transverse macro-dispersion coefficient $D_{T,m}$ with the pore-scale transverse dispersion coefficient $D_{T,p}$ provides important insights regarding the prevalence of chaotic advection in heterogeneous Darcy flow.

For most porous materials the hydraulic conductivity tensor $\boldsymbol{K}(\boldsymbol{x})$ is anisotropic due to anisotropy of the underlying pore-scale structure, which occurs in, e.g. stratified formations in geological media and structured engineered media (Bear Reference Bear1972). Under conventional upscaling approaches the conductivity tensor $\boldsymbol{K}(\boldsymbol{x})$ (a second rank tensor) is symmetric and positive definite (Bear Reference Bear1972), with six independent components in 3-D space which vary spatially in heterogeneous formations but are all assumed to be smooth, continuous and finite. In the absence of boundaries and fluid sources and sinks, the positive–definite nature of $\boldsymbol{K}(\boldsymbol{x})$ renders Darcy-scale flow stagnation free, i.e. $\boldsymbol{v}(\boldsymbol{x})\neq \boldsymbol{0} \, \forall \boldsymbol{x}\in \varOmega$ (Bear Reference Bear1972). As $\boldsymbol{K}(\boldsymbol{x})$ is symmetric, it may be locally diagonalised into principal directions with three independent components as $\boldsymbol{K}^\prime (\boldsymbol{x})=\boldsymbol{R}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{K}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{R}^{-1}(\boldsymbol{x})$ , where $\boldsymbol{R}(\boldsymbol{x})$ is a rotation tensor field and

(2.4) \begin{equation} \boldsymbol{K}^\prime (\boldsymbol{x})\equiv \left ( \begin{array}{ccc} K_{11}^\prime (\boldsymbol{x}) & 0 & 0\\[4pt] 0 & K_{22}^\prime (\boldsymbol{x}) & 0\\[4pt] 0 & 0 & K_{33}^\prime (\boldsymbol{x}) \end{array} \right )\!. \end{equation}

Even if $\boldsymbol{R}(\boldsymbol{x})$ is spatially invariant, this conductivity field is fundamentally anisotropic if one of the diagonal components $K_\textit{ii}^\prime (\boldsymbol{x})$ differs from the others.

Conversely, many studies (Gelhar & Axness Reference Gelhar and Axness1983; Neuman & Zhang Reference Neuman and Zhang1990; Chaudhuri & Sekhar Reference Chaudhuri and Sekhar2005; Delgado Reference Delgado2007; Janković et al. Reference Janković, Steward, Barnes and Dagan2009; Beaudoin & de Dreuzy Reference Beaudoin and de Dreuzy2013; Boon et al. Reference Boon, Bijeljic, Niu and Krevor2016; Dartois, Beaudoin & Huberson Reference Dartois, Beaudoin and Huberson2018) consider the hydraulic conductivity tensor to be isotropic as $\boldsymbol{K}(\boldsymbol{x})=k(\boldsymbol{x})\boldsymbol{I}$ , where $k(\boldsymbol{x})=K_\textit{ii}^\prime (\boldsymbol{x})$ for $i=1:3$ , and (2.2) simplifies to

(2.5) \begin{equation} \boldsymbol{v}(\boldsymbol{x})=-k(\boldsymbol{x})\boldsymbol{\nabla }\phi ,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}=0,\,\,\,\,\boldsymbol{x}\in \varOmega . \end{equation}

The kinematics of isotropic Darcy flow markedly differ from that of anisotropic Darcy flow. If $k(\boldsymbol{x})$ is smooth and $\boldsymbol{v}$ does not admit stagnation points, then isotropic Darcy flow only admits simple kinematics because the helicity density $\mathcal{H}(\boldsymbol{x})\equiv \boldsymbol{v}\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{v})$ is identically zero

(2.6) \begin{equation} \mathcal{H}(\boldsymbol{x})=k\boldsymbol{\nabla }\phi \boldsymbol{\cdot }(\boldsymbol{\nabla }\phi \times \boldsymbol{\nabla } k)=0, \end{equation}

regardless of the heterogeneity of $k(\boldsymbol{x})$ . For inviscid flows, the total helicity $H\equiv \int _\varOmega \mathcal{H}(\boldsymbol{x})\,{\rm d}\boldsymbol{x}$ over the flow domain $\varOmega$ is a topological invariant that characterises the topological complexity (knottedness) of vortex lines of the flow (Woltjer Reference Woltjer1958; Moreau Reference Moreau1961; Moffatt Reference Moffatt1969), but in general the zero helicity density condition $\mathcal{H}(\boldsymbol{x})=0$ precludes chaotic dynamics in steady 3-D flows (Arnol’d Reference Arnol’d1965; Moffatt & Tsinober Reference Moffatt and Tsinober1992). Hence steady 3-D isotropic Darcy flows are non-chaotic and integrable and so admit two invariants $\psi _1$ , $\psi _2$ that act as streamfunctions (Yoshida Reference Yoshida2009; Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2022), leading to the Euler velocity representation

(2.7) \begin{equation} \boldsymbol{v}(\boldsymbol{x})=\boldsymbol{\nabla }\psi _1(\boldsymbol{x})\times \boldsymbol{\nabla }\psi _2(\boldsymbol{x}). \end{equation}

Streamlines of steady isotropic Darcy flow are confined to the intersections of streamsurfaces given by level sets of $\psi _1(\boldsymbol{x})$ , $\psi _2(\boldsymbol{x})$ . Such confinement prohibits transverse macro-dispersion (henceforth simply termed transverse dispersion) of the streamlines of $\boldsymbol{v}$ in the absence of local-scale dispersion, regardless of medium heterogeneity (Lester et al. Reference Lester, Dentz, Singh and Bandopadhyay2023).

2.2. Implications for hydraulic conductivity modelling

These kinematic constraints are illustrated in figure 1(c), which shows typical streamlines for isotropic Darcy flow generated by a strongly heterogeneous ( $\sigma ^2_{\textit{ln} K}=4$ ) hydraulic conductivity field. Although the streamlines are highly tortuous due to heterogeneity of the medium, they do not exhibit asymptotic transverse dispersion. In § 5.1 the transverse dispersivity of this flow is numerically computed to be effectively zero.

These kinematic constraints persist even if the scalar field $k(\boldsymbol{x})$ is statistically anisotropic (i.e. has different correlation structures in different principal directions), as the corresponding Darcy flow is still locally isotropic and $\mathcal{H}(\boldsymbol{x})=0$ . We remark on the application of these results to the so-called ‘helicity paradox’ (Cirpka et al. Reference Cirpka, Chiogna, Rolle and Bellin2015) that occurs when a statistically anisotropic but locally isotropic conductivity field is upscaled from the Darcy scale to the block field scale, resulting in an anisotropic block-scale conductivity field (Bear Reference Bear1972). This spuriously adds degrees of freedom to the Lagrangian kinematics and permits block-scale chaotic advection where none should exist based on the fully resolved Darcy-scale flow. While beyond the scope of this paper, our results highlight the need for upscaling methods that obey the appropriate kinematic constraints.

Figure 1. (a) Isosurfaces of the typical normalised heterogeneous log-conductivity field $f(\boldsymbol{x})=\textrm{ln}\, K(\boldsymbol{x})/\sigma ^2_{\textit{ln} K}$ used to model isotropic $k(\boldsymbol{x})\boldsymbol{I}$ and anisotropic $\boldsymbol{K}(\boldsymbol{x})$ conductivity tensors and (b) associated potential field $\phi (\boldsymbol{x})$ for anisotropic Darcy flow driven by a uniform mean potential gradient. Note coloured surfaces are isopotential surfaces $\phi =$ const. Associated streamlines for heterogeneous Darcy flow with (c) isotropic conductivity field ( $\delta =0$ in (5.1)) and (d) anisotropic conductivity field ( $\delta =1$ in (5.1)) with log-conductivity variance $\sigma ^2_{\textit{ln} K}=4$ and parameters $N=4$ , $N_i=2$ in (5.2).

These results provide insights into the connection between transverse macro-dispersion in steady Darcy flow and anisotropy of the hydraulic conductivity tensor $\boldsymbol{K}(\boldsymbol{x})$ . Unfortunately, the prevalence (or otherwise) of transverse macro-dispersion in field studies is an open question. While recent analyses (Zech et al. Reference Zech, Attinger, Bellin, Cvetkovic, Dietrich, Fiori, Teutsch and Dagan2019) of field experiments of solute transport (Sudicky, Cherry & Frind Reference Sudicky, Cherry and Frind1983; Rajaram & Gelhar Reference Rajaram and Gelhar1991; Hess, Wolf & Celia Reference Hess, Wolf and Celia1992) over a range of geological media clearly establish that longitudinal macro-dispersivity is significantly greater than pore-scale dispersivity, indicating the presence of large-scale conductivity heterogeneities at the Darcy scale, similar data are inconclusive with respect to transverse dispersivity. Hence, the prevalence of transverse macro-dispersion in groundwater flows in heterogeneous aquifers is currently unknown.

Although several numerical studies (Janković et al. Reference Janković, Steward, Barnes and Dagan2009; Beaudoin & de Dreuzy Reference Beaudoin and de Dreuzy2013; Dartois et al. Reference Dartois, Beaudoin and Huberson2018) have found that isotropic Darcy flow can generate purely advective transverse dispersion, these observations arise in systems with either non-smooth conductivity fields or impermeable inclusions, or from numerical schemes that violate the kinematic constraints associated with (2.7). Similar to how integration schemes that are not symplectic violate Hamiltonian constraints (i.e. invariance of the Hamiltonian), numerical schemes that are not explicitly designed to do so also violate the constraints associated with (2.7). Lester et al. (Reference Lester, Dentz, Singh and Bandopadhyay2023) show that numerical schemes that implicitly enforce these constraints yield zero transverse dispersivity.

Conversely, anisotropic conductivity fields $\boldsymbol{K}(\boldsymbol{x})$ generate non-zero helicity density

(2.8) \begin{equation} \mathcal{H}(\boldsymbol{x})=(\boldsymbol{K}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi )\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{K}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi )\neq 0, \end{equation}

and non-zero total helicity $H\neq 0$ . Hence, the streamlines of anisotropic Darcy flow are no longer constrained to the streamsurfaces of $\psi _1(\boldsymbol{x})$ , $\psi _2(\boldsymbol{x})$ , but rather freely wander through the medium (Lester et al. Reference Lester, Dentz, Singh and Bandopadhyay2023), giving rise to persistent transverse dispersion. This behaviour is shown in figure 1(d), where the streamlines associated with strongly heterogeneous ( $\sigma ^2_{\textit{ln} K}=4$ ) anisotropic Darcy flow spread transversely throughout the flow domain as they are advected longitudinally. In § 5.1 the transverse dispersion coefficient $D_T$ of this flow is quantitatively shown to be non-zero.

Note that non-zero helicity density $\mathcal{H}(\boldsymbol{x})\neq 0$ alone is not sufficient to ensure non-zero transverse dispersivity. Yoshida (Reference Yoshida2009) shows that only 3-D velocity fields given by the complete Clebsch parameterisation

(2.9) \begin{equation} \boldsymbol{v}(\boldsymbol{x})=\sum _{i=1}^{N=2} \alpha _i(\boldsymbol{x})\boldsymbol{\nabla }\beta _i(\boldsymbol{x})+\boldsymbol{\nabla }\varphi (\boldsymbol{x}), \end{equation}

where $\alpha _i(\boldsymbol{x})$ , $\beta _i(\boldsymbol{x})$ , $\varphi (\boldsymbol{x})$ are smooth continuous scalar fields, have non-zero total helicity $H$ and so admit chaotic advection. Conversely, the classical Clebsch parameterisation

(2.10) \begin{equation} \boldsymbol{v}(\boldsymbol{x})=\alpha _1(\boldsymbol{x})\boldsymbol{\nabla }\beta _1(\boldsymbol{x})+\boldsymbol{\nabla }\varphi (\boldsymbol{x})\quad \boldsymbol{x}\in \varOmega , \end{equation}

has been found to be incomplete (Yoshida Reference Yoshida2009), meaning that some 3-D vector fields cannot be represented by (2.10). Yoshida & Morrison (Reference Yoshida and Morrison2017) classify flows conforming to (2.10) as epi-2-D flows which are topologically two-dimensional with non-zero helicity density $\mathcal{H}(\boldsymbol{x})=\boldsymbol{\nabla }\varphi _1\boldsymbol{\cdot }(\boldsymbol{\nabla }\alpha _1\times \boldsymbol{\nabla }\beta _1)\neq 0$ but are helicity free, with $H=0$ . Conversely, realistic heterogeneous conductivity models yield 3-D velocity fields of the form (2.9) that admit chaotic advection and transverse dispersion.

To summarise, if transverse macro-dispersion is non-zero in heterogeneous aquifers in the advection-dominated limit, that is, in the absence of local-scale dispersion, then hydraulic conductivity models of heterogeneous porous media must be anisotropic to be consistent. Although the prevalence of transverse macro-dispersion in field studies is inconclusive, we note that many aquifers have an anisotropic geological structure and many groundwater models utilise anisotropic hydraulic conductivity fields to predict solute transport in groundwater applications. As non-zero transverse macro-dispersion requires $H\neq 0$ , this raises the potential for chaotic advection, which is inevitable in nonlinear continuous systems with arbitrary coefficients and three degrees of freedom (dof = 3) (Speetjens et al. Reference Speetjens, Metcalfe and Rudman2021). In the following section we firmly establish that transverse macro-dispersion and chaotic advection are intimately linked in Darcy flow.

3.1. Background

Figure 2. (a) Braid diagram depicting stretching of material elements (red) around fluid particle streamlines (black) numbered $n=1$ to $n=3$ from left to right in a 3-dof flow. For a steady 3-D flow braiding evolves in the longitudinal $x_1$ coordinate, and for an unsteady 2-D flow braiding evolves with respect to time $t$ . The pathline crossing events are characterised by the respective braid generators $\sigma _1$ and $\sigma _2^{-1}$ , which act to stretch the material element due to the topology of the braiding motions. (b) Stretching of material elements (brown) due to braiding motions (adapted from Thiffeault Reference Thiffeault2022) of streamlines (red circles) corresponding to the braid diagram in (a) that evolves with the longitudinal direction $x_1$ or time $t$ .

The link between dispersion and chaotic mixing can be explored via consideration of streamline braiding. As 1-D streamlines are invariants of steady 3-D flow, braiding of streamlines stirs fluid elements in a complex manner, as illustrated in figure 2. This concept has been used to quantify stirring in another 3-dof system – unsteady 2-D flow – where non-trivial braiding of 1-D pathlines can also stretch and fold the fluid continuum. For unsteady 2-D flow, braid group theory (Handel & Thurston Reference Handel and Thurston1985; Boyland, Aref & Stremler Reference Boyland, Aref and Stremler2000; Moussafir Reference Moussafir2006) has been developed to quantify the topological entropy $h_{\textit{braid}}$ of the pathline braiding motions, closely related to the Lyapunov exponent $\lambda _\infty$ (Thiffeault Reference Thiffeault2005). In § 3.2 we show that steady 3-D Darcy flow is topologically equivalent to 2-D unsteady flow, hence this mathematical framework applies to these flows.

The topological equivalence between stagnation-free steady 3-D flow and unsteady 2-D flow has been established since the earliest studies of chaotic advection (Aref Reference Aref1984; Bajer Reference Bajer1994), with particular focus on steady 3-D duct flows (Lester et al. Reference Lester, Kuan and Metcalfe2018b ; Speetjens et al. Reference Speetjens, Metcalfe and Rudman2021). Steady heterogeneous 3-D Darcy flow driven by a unidirectional mean pressure gradient shares many similarities with stagnation-free steady 3-D duct flows in that both flows are topologically equivalent to unsteady 2-D flow and chaotic advection arises due to transverse perturbations away from a fully integrable state. However, a fundamental difference between these flows is that heterogeneous Darcy flow is typically random and aperiodic, meaning the tools and techniques of Hamiltonian chaos for temporally periodic 2-D systems (Poincaré sections, analysis of periodic points, resolution of hyperbolic manifolds and Kolmogorov–Arnold–Moser islands) (Ottino Reference Ottino1989) do not apply to these flows. While tools such as Lagrangian coherent structures (LCS) (Haller Reference Haller2015) are well suited to uncover the advection structure of aperiodic flows, in this study we focus on braid group theory (or topological fluid mechanics) (Boyland et al. Reference Boyland, Aref and Stremler2000) as this approach, as shall be shown, provides a vital quantitative link between chaotic advection and transverse dispersion that is central to this study.

Braid group theory is used to quantify stirring in unsteady 2-D flows (Boyland et al. Reference Boyland, Aref and Stremler2000) via a well-developed mathematical framework (Artin Reference Artin1947; Moussafir Reference Moussafir2006; Thiffeault Reference Thiffeault2022) that efficiently encodes fluid stirring via a symbolic braiding representation of the topology of pathlines in space–time. This encoding can be used to compute the topological complexity of their braiding motions, and the Neilson–Thurston classification theorem (Handel & Thurston Reference Handel and Thurston1985) characterises braids as periodic, reducible or pseudo-Anosov types, the latter of which is interpreted as an indicator of chaotic dynamics (Boyland et al. Reference Boyland, Aref and Stremler2000). Due to topological equivalence, this framework can also be applied to streamline braiding in steady 3-D Darcy flow.

Figure 3. (a) Propagation of 1-D streamlines (black) in a steady unidirectional 3-D flow with mean flow in the longitudinal ( $x_1$ ) direction. (b) Schematic of braiding of streamlines (black circles) labelled $n-1$ , $n$ , $n+1$ in the transverse $x_2-x_3$ plane via clockwise (black) $\sigma _n$ , anti-clockwise (blue) $\sigma _n^{-1}$ braid generators.

For steady unidirectional 3-D flow, streamline braiding is encoded via a sequence of braid generators $\sigma _n^{\pm 1}$ (Artin Reference Artin1947). Figure 3(a) shows propagation in the $x_1$ direction of a set of streamlines of a steady unidirectional 3-D flow that undergo random ‘crossing’ motions with respect to the transverse $x_3$ coordinate. As shown in figure 3(b), these crossing motions can be characterised in terms of the braid generators $\sigma _n$ and $\sigma _n^{-1}$ , which respectively characterise clockwise and counter-clockwise crossing of streamlines $n$ and $n+1$ in the transverse $(x_2,x_3)$ plane. Note streamlines are labelled ( $\ldots ,n-1, n,n+1,\ldots$ ) with respect to their relative position in the $x_2$ coordinate, so these labels are updated after each crossing event. An ordered sequence of braid generators forms the braid word $\boldsymbol{b}=\sigma _{n_1}^{\pm 1}\sigma _{n_2}^{\pm 1}\ldots$ , which is determined by the sequence of streamline crossings in the $x_1$ direction. Hence, $\boldsymbol{b}$ completely characterises the topological braiding motions of the streamlines. As shown in figure 2(b), these braiding motions act to stretch and fold fluid elements (brown), leading to exponential stretching if they are of pseudo-Anosov type.

The rate of stretching generated by a braid word $\boldsymbol{b}$ is given by its topological braid entropy $h_{\textit{braid}}$ . In general, the topological entropy $\hat {h}$ of a dynamical system measures the rate of loss of information of the system about its initial conditions. For fluid flows $\hat {h}$ is closely related to the largest infinite-time Lyapunov exponent $\hat {\lambda }_\infty$ and forms an upper bound for $\hat {\lambda }_\infty$ (Thiffeault Reference Thiffeault2010). In practice, the topological entropy of a flow can be interpreted as the asymptotic exponential growth rate of the length $l(t)$ of a material line (Newhouse & Pignataro Reference Newhouse and Pignataro1993) as

(3.1) \begin{equation} \hat {h}=\lim _{t\rightarrow \infty }\frac {1}{t}\textrm{ln}\, \frac {l(t)}{l(0)}. \end{equation}

For ergodic flows such as purely hyperbolic steady 3-D flow (which has a single positive Lyapunov exponent and does not admit non-hyperbolic features such as elliptic or parabolic LCS (Haller Reference Haller2015), the upper bound given by the topological entropy is exact (Yomdin Reference Yomdin1987; Newhouse Reference Newhouse1988; Newhouse & Pignataro Reference Newhouse and Pignataro1993; Matsuoka & Hiraide Reference Matsuoka and Hiraide2015; Catalan Reference Catalan2019), i.e.

(3.2) \begin{equation} \hat {h}=\hat {\lambda }_\infty . \end{equation}

Conversely, many studies dating over half a century (Cocke Reference Cocke1969; Girimaji & Pope Reference Girimaji and Pope1990; Hinch Reference Hinch1999; Kalda Reference Kalda2000; Duplat, Innocenti & Villermaux Reference Duplat, Innocenti and Villermaux2010; Villermaux Reference Villermaux2019) propose that as fluid elements undergo stretching, the relative length $\rho \equiv \delta l(t)/\delta l(0)$ of infinitesimal line elements is distributed log normally with log mean $\hat {\lambda }_\infty t$ and log variance $\hat {\sigma }^2_\lambda t$ . In § 5.2 and Appendix C we confirm that an ab initio derivation of fluid stretching in steady 3-D flows yields the same result if the velocity distribution does not yield anomalous transport. For this scenario a reasonable conclusion (Hinch Reference Hinch1999; Duplat et al. Reference Duplat, Innocenti and Villermaux2010; Villermaux Reference Villermaux2019) is that the growth rate of a material line of length $l(t)=\int _0^{l(0)}\delta l(X,t) {\rm d}X$ (where $X$ is the Lagrangian coordinate along the line) is given by the ensemble average of $\rho$ , which, via (3.1), yields

(3.3) \begin{equation} \hat {h}=\hat {\lambda }_\infty +\frac {\sigma ^2_\lambda }{2}, \end{equation}

as the variance of $\textrm{ln}\, \rho$ contributes to the growth of ${ln} l(t)$ . However, the discrepancy between (3.3) and (3.2) has recently been explained by Lester & Dentz (Reference Lester and Dentz2025), who show that in the asymptotic limit, temporal sampling of the finite-time Lyapunov exponent $\hat {\lambda }(\boldsymbol{X},t)$ (which converges to $\hat {\lambda }_\infty$ as $\sigma _\lambda /\sqrt {t}$ ) dominates over the ensemble average represented by (3.3), yielding (3.2). In § 4 we show that heterogeneous Darcy flow with random stationary conductivity $\boldsymbol{K}(\boldsymbol{x})$ is ergodic and purely hyperbolic.

The dimensionless topological entropy $h\equiv \hat {h}\ell /\langle v_1\rangle$ (with velocity correlation length $\ell$ ) may be approximated by considering the topological braid entropy $h_{\textit{braid}}$ of the braid word $\boldsymbol{b}$ comprising a set of $N_b$ braiding motions of $N_p$ streamlines of the flow (Thiffeault Reference Thiffeault2010). The braid entropy $h_{\textit{braid}}$ characterises the asymptotic growth of the number of distinguishable orbits of the braid (Thiffeault Reference Thiffeault2022). For steady 3-D flow this may be interpreted as the number of distinct linkages that a material line makes if it is fully contracted until it forms a loop $\ell _E$ that tightly contacts these streamlines, where the metric $|\ell _E|$ measures the number of linkages of the loop $\ell _E$ . For example, the evolving brown material element shown in figure 2(b) can be approximated as a number of straight linkages that either connect to or touch the red stirring rods. Under the initial configuration, the element can be approximated as one linkage ( $|\ell _E|=1$ ) connecting the $n=2$ and $n=3$ rods, but after the braid generator $\sigma _1$ , there are two linkages ( $|\ell _E|=2$ ) connecting the $n=1$ to $n=2$ to $n=3$ rods, and finally after $\sigma _1\sigma _2^{-1}$ there are three linkages ( $|\ell _E|=3$ ) connecting rods $n=1$ to $n=2$ , $n=2$ to $n=3$ and then $n=3$ to $n=2$ . The topological entropy $h_{\textit{braid}}$ of a braid word $\boldsymbol{b}$ is then

(3.4) \begin{equation} h_{\textit{braid}}=\lim _{k\rightarrow \infty }\frac {1}{k}\log \frac {|\boldsymbol{b}^k \ell _E|}{|\ell _E|}, \end{equation}

where $\boldsymbol{b}^k\ell _E$ indicates the braid $\boldsymbol{b}$ has operated $k$ times on the loop $\ell _E$ . Although non-trivial, various methods (Moussafir Reference Moussafir2006; Hall & Yurttaş Reference Hall and Yurttaş2009; Thiffeault Reference Thiffeault2022) are available to rigorously compute $h_{\textit{braid}}$ from a given braid word $\boldsymbol{b}$ . As the braid entropy $h_{\textit{braid}}$ only considers the complexity of streamline braiding, it forms a lower bound for $h$ , however, for many systems it converges to $h$ with increasing number of streamlines $N_p$ and braiding motions $N_b$ as

(3.5) \begin{equation} \lim _{N_p\rightarrow \infty }\lim _{N_b\rightarrow \infty }\frac {h_{\textit{braid}}}{t}\rightarrow \hat {h}=\hat {\lambda }_\infty . \end{equation}

Hence, for heterogeneous Darcy flow both the topological entropy $h$ and dimensionless Lyapunov exponent $\lambda _\infty \equiv \hat {\lambda }_\infty \ell /\langle v_1\rangle$ may be accurately estimated from the braid entropy $h_{\textit{braid}}$ of a sufficiently resolved set of streamlines (large $N_p$ ) recorded over a sufficient observation time (large $N_b$ ).

3.2. Topology of steady heterogeneous 3-D Darcy flow

The topology of steady 3-D flows $\boldsymbol{v}(\boldsymbol{x})=[v_1(\boldsymbol{x}),v_2(\boldsymbol{x}),v_3(\boldsymbol{x})]$ with one unidirectional (axial) velocity component (e.g. $v_1(\boldsymbol{x})\gt 0$ ) can be equated (Bajer Reference Bajer1994) to that of unsteady 2-D flow via the rescaling $\boldsymbol{v}^\prime (\boldsymbol{x})=\boldsymbol{v}(\boldsymbol{x})/v_1(\boldsymbol{x})$ , $t'=x_1$ , such that the $x_1$ coordinate acts as an analogue for time $t$ , hence the advection equation transforms as

(3.6) \begin{equation} \frac {\rm d}{{\rm d}t}[x_1,x_2,x_3]^\top =[v_1,v_2,v_3]^\top \Rightarrow \frac {\rm d}{{\rm d}t'}[x_2,x_3]^\top =\left [\frac {v_2(t',x_2,x_3)}{v_1(t',x_2,x_3)},\frac {v_3(t',x_2,x_3)}{v_1(t',x_2,x_3)}\right ]^\top . \end{equation}

This transformation brings several advantages with respect to analysis of the advective mixing properties of the flow, as axially periodic 3-D flows such as duct flows (Speetjens et al. Reference Speetjens, Metcalfe and Rudman2021) can be couched as temporally periodic 2-D flows, hence the tools and techniques of Hamiltonian chaos (Aref Reference Aref1984; Ottino Reference Ottino1989) can be used to study mixing in these systems. Although the analogous 2-D flow in (3.6) is no longer divergence free, the dynamics of this flow is still Hamiltonian in the stroboscopic frame (Bajer Reference Bajer1994). For aperiodic 3-D duct flows, conversion to aperiodic transient 2-D flows via (3.6) facilitates the use of tools such as LCS (Haller Reference Haller2015) to uncover the mixing dynamics of these flows. Indeed, this transformation has been used to study a wide class of both periodic and aperiodic steady 3-D duct flows (Speetjens et al. Reference Speetjens, Metcalfe and Rudman2006).

For strongly heterogeneous 3-D Darcy flows the transformation (3.6) breaks down when $v_1(\boldsymbol{x})$ is zero (Bajer Reference Bajer1994) due to flow reversal in the vicinity of low permeability structures. However, these flows are still fundamentally unidirectional as the potential field $\phi (\boldsymbol{x})$ is strictly monotonic decreasing along streamlines $\boldsymbol{\nabla }\phi (\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{v}(\boldsymbol{x})=-\boldsymbol{\nabla }\phi (\boldsymbol{x}){\times}\boldsymbol{K}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{\nabla }\phi (\boldsymbol{x})\lt 0$ due to the positive definite nature of $\boldsymbol{K}(\boldsymbol{x})$ (Bear Reference Bear1972). This leads to the intrinsic coordinate basis $\boldsymbol \xi =(\chi _1,\chi _2,\phi )$ , where $\chi _1$ , $\chi _2$ are arbitrary coordinates which are both orthogonal to $\phi$ , and the intrinsic velocity representation $ \boldsymbol{v}_\xi (\boldsymbol \xi )=[v_{\chi _1}(\boldsymbol \xi ),v_{\chi _2}(\boldsymbol \xi ),v_{\phi }(\boldsymbol \xi )]$ with $v_{\phi }(\boldsymbol \xi )\gt 0$ . Thus, the intrinsic velocity field $\boldsymbol{v}_\xi (\boldsymbol{\xi })$ may be rescaled as

(3.7) \begin{equation} \boldsymbol{v}_\xi ^\prime (\boldsymbol \xi )=\left [\frac {v_{\chi _1}(\boldsymbol \xi )}{v_{\phi }(\boldsymbol \xi )},\frac {v_{\chi _2}(\boldsymbol \xi )}{v_{\phi }(\boldsymbol \xi )},1\right ]\!,\quad v_{\phi }(\boldsymbol \xi )\gt 0\,\forall \,\boldsymbol \xi \in \varOmega . \end{equation}

This velocity field is topologically equivalent to that of unsteady 2-D flow, and the isopotential surfaces of steady 3-D Darcy flow (figure 1 b) are analogous to isochrones of unsteady 2-D flow. This means that steady 3-D Darcy flow always has an embedded Hamiltonian structure (and associated Lagrangian kinematics) as explained in § 4.3 of Speetjens et al. (Reference Speetjens, Metcalfe and Rudman2021), and so maintains a fundamental connection with the generation of chaotic advection in stagnation-free steady 3-D flows in general.

Figure 4. (a) Schematic of streamline (or pathline) braiding in a unidirectional steady 3-D flow (or unsteady 2-D flow). The set of four thick coloured streamlines (or pathlines) braid with each other as they evolve with the mean flow direction $x_3$ (or time $t$ ). This topological braiding motion stirs the fluid continuum (grey planes) and the length of the black rectangular material boundary grows exponentially with the number of braiding motions. Adapted from Thiffeault & Finn (Reference Thiffeault and Finn2006). This schematic also depicts streamline braiding in steady 3-D anisotropic Darcy flow with intrinsic coordinates $\boldsymbol \xi =(\chi _1,\chi _2,\phi )$ , where the grey planes depict isopotential surfaces. (b) Absence of braiding in isotropic Darcy flow in intrinsic coordinates $\boldsymbol \xi =(\chi _1,\chi _2,\phi )$ . As the velocity field is everywhere orthogonal to level sets of $\phi$ denoted by grey planes, streamlines in this coordinate system (which are simply straight lines) do not move laterally or undergo braiding.

As such, streamlines in steady 3-D Darcy flow are topologically equivalent to pathlines of unsteady 2-D flow and the mathematical framework developed for pathline braiding can be directly applied to steady 3-D Darcy flow. As per figure 4(b), the representation (3.7) also shows directly that isotropic Darcy flow cannot generate streamline braiding, where $(\chi _1,\chi _2)=(\psi _1,\psi _2)$ and (3.7) simplifies to $\boldsymbol{v}_\xi ^\prime (\boldsymbol \xi )=(0,0,1)$ because the velocity field is everywhere orthogonal to the level sets of $\phi$ , i.e. $\boldsymbol{v}\times \boldsymbol{\nabla }\phi =\boldsymbol{0}$ . These streamlines do not move relative to each other in the isopotential surfaces $\phi = \text{const.}$ , let alone undergo non-trivial braiding motions. In Cartesian coordinates (figure 1 c), the helicity-free condition constrains non-braiding yet tortuous streamlines to streamsurfaces $\psi _1(\boldsymbol{x})$ = const., $\psi _2(\boldsymbol{x})$ = const. Conversely, anisotropic Darcy flow admits non-zero transverse velocity components $v_{\chi _1}(\boldsymbol \xi )$ , $v_{\chi _2}(\boldsymbol \xi )$ and so streamlines can undergo the braiding motions as shown in figures 4(a) and 1(d). Although transverse motion of streamlines in the $(\chi _1,\chi _2)$ directions does not necessarily generate non-trivial braiding, below we show that non-trivial braiding is inherent to steady random 3-D flows.

3.3. Linking dispersion and chaos in streamline braiding flows

Lester, Metcalfe & Trefry (Reference Lester, Metcalfe and Trefry2024) recently examined the link between chaotic advection and transverse dispersion in randomly braiding flows via a simple 1-D streamline model. As per figure 3(a), this model consists of a set of $N_p$ streamlines in $\mathbb{R}^3$ space that repeat periodically in the transverse $x_2$ direction and propagate due to mean flow in the longitudinal $x_1$ direction. These streamlines have the same $x_3$ -coordinate ( $x_3=0$ ) and uniform spacing $\varDelta x_2=\ell$ in the $x_2$ -direction, corresponding to the velocity correlation length scale $\ell$ . At integer multiples of longitudinal distance $\varDelta _L/N_p$ , one pair ( $n$ , $n+1$ ) of neighbouring streamlines randomly exchange position in the $x_2-x_3$ plane via clockwise or counter-clockwise rotations, as labelled by the respective braid generators $\sigma _n$ , $\sigma _n^{-1}$ . Figure 5(a) shows the braid diagram for a typical braid word $\boldsymbol{b}$ consisting of $N_p=20$ streamlines and $N_b=20$ braid generators. For an array of $N_p$ streamlines, the distance $\varDelta _L/N_p$ corresponds to a braid frequency $\hat {f}$ per streamline of

(3.8) \begin{equation} \hat {f}\equiv \frac {\langle v_1\rangle }{\varDelta _L}. \end{equation}

The streamline braiding motions are defined by a braid word $\boldsymbol{b}$ comprising $N_b$ braid generators $\boldsymbol{b}=\sigma _{n_1}^{\pm 1}\sigma _{n_2}^{\pm 1}\ldots \sigma _{n_{N_b}}^{\pm 1}$ which is constructed by randomly choosing each generator from the set $\{\sigma _1,\sigma _1^{-1},\ldots ,\sigma _{N_p},\sigma _{N_p}^{-1}\}$ . Note that, due to periodicity, the $n=N_p$ pathline can braid with the $n=1$ pathline via the $\sigma _{N_p}$ and $\sigma _{N_p}^{-1}$ braid generators as an annular braid (Finn & Thiffeault Reference Finn and Thiffeault2007). Hence, streamlines undertake an unbounded random walk along the $x_2$ coordinate as they propagate longitudinally (figure 5 b, inset), while deforming the interstitial fluid as per figure 2. The topological braid entropy $h_{\textit{braid}}$ of these motions is efficiently computed via the braidlab package (Thiffeault & Budišić Reference Thiffeault and Budišić2013–Reference Thiffeault and Budišić2021). Figure 5(b) shows linear growth of topological braid entropy $h_{\textit{braid}}$ and transverse variance $\sigma _{x_2}^2$ of this system with $N_b$ .

Figure 5. (a) Braid diagram for 1-D streamline depicting evolution of $x_2$ coordinate of $N_p=20$ streamlines over $N_b=20$ random braid actions in the longitudinal $x_1$ direction, leading to non-trivial braiding and dispersion of streamlines. (b) Linear growth of topological braid entropy $h_{\textit{braid}}$ (black points) with $N_b$ , in agreement with (3.13) (green dashed line). Growth of transverse variance $\sigma _{x_2}^2(N_b)$ (blue line) with $N_b$ , in agreement with (3.11) (red line). Inset: Brownian motion of pathlines with increasing $N_b$ .

Lester et al. (Reference Lester, Metcalfe and Trefry2024) show that the topological braid entropy $h_{\textit{braid}}$ and transverse variance $\sigma _{x_2}^2$ of this system grow linearly with $N_b$ , and in the limit of large $N_b$ and $N_p$ , the scaled mean topological braid entropy $\langle h_{\textit{braid}}\rangle$ for $10^3$ realisations of this simple 1-D streamline model converges to the random braid entropy $\langle \lambda _\sigma \rangle$ as

(3.9) \begin{equation} \lim _{N_p\rightarrow \infty }\lim _{N_p\rightarrow \infty }\frac {N_p}{N_b}\langle h_{\textit{braid}}\rangle \rightarrow \frac {\varDelta _L}{\langle v_1\rangle }\hat {h}\equiv \langle \lambda _\sigma \rangle = 0.8529. \end{equation}

The random braid entropy $\langle \lambda _\sigma \rangle$ is a fundamental quantity in that it characterises the topological entropy of random braids in limit of large $N_p$ and $N_b$ . From (3.9), the random braid entropy is related to the dimensionless Lyapunov exponent $\lambda _\infty$ and topological entropy $h$ as

(3.10) \begin{equation} \lambda _\infty = h = \frac {\ell }{\varDelta _L}\langle \lambda _\sigma \rangle = f \langle \lambda _\sigma \rangle , \end{equation}

where $f\equiv \hat {f}\ell /\langle v_1\rangle$ is the dimensionless braiding frequency that characterises the frequency of braiding events for a single streamline. As $\langle \lambda _\sigma \rangle$ is a universal constant, $f$ is necessary and sufficient to completely characterise streamline braiding in random braiding flows. For streamlines with mean spacing $\ell$ , $f$ characterises the average longitudinal distance $\varDelta _L$ between braiding events in steady 3-D flow. Together, $\langle \lambda _\sigma \rangle$ and $f$ quantify the topological entropy $h$ of the 1-D streamline model shown in figure 3(a).

These quantities can also be connected to the advective transverse dispersion coefficient $D_T$ driving the growth of $\sigma _{x_2}^2$ shown in figure 5(b). As two of the $N_p$ streamlines make random jumps in the $x_2$ direction of $\pm \ell$ with each braiding event, the spatial variance $\sigma _{x_2}^2$ of streamlines grows linearly with $N_b$ as

(3.11) \begin{equation} \sigma _{x_2}^2(t)=\sigma _{x_2}^2(0)+2\ell ^2\frac {N_b}{N_p}\equiv \sigma _{x_2}^2(0)+2D_T t, \end{equation}

where the variance growth per braid event is $2\ell ^2/N_p$ as each event moves two of the $N_p$ streamlines by $\pm \ell$ . As $N_b=N_p\langle v_1\rangle t/\varDelta _L$ , then

(3.12) \begin{equation} D_T=\frac {\ell ^2\langle v_1\rangle }{\varDelta _L}, \end{equation}

which, as shown in figure 5(b), agrees well with model observations. Defining the advective transverse Péclet number as $\textit{Pe}_T\equiv D_T/(\langle v_1\rangle \ell )$ then yields $\textit{Pe}_T=1/f$ . Hence, for the 1-D streamline model the dimensionless Lyapunov exponent $\lambda _\infty$ and topological entropy $h$ are linearly related to the transverse dispersion coefficient $D_T$ as

(3.13) \begin{equation} \lambda _\infty =h=\frac {\langle \lambda _\sigma \rangle }{\textit{Pe}_T}=f \langle \lambda _\sigma \rangle . \end{equation}

This linear relationship arises as the 1-D streamline model generates transverse dispersion and chaotic advection in equal measure.

Two different 2-D extensions of this 1-D streamline model were also considered by Lester et al. (Reference Lester, Metcalfe and Trefry2024); one involving a 2-D streamline array analogous to the 1-D array in figure 5(a), and another involving 3-D streamlines that undertake independent random walks constructed by making a jump of fixed magnitude by random orientation in the $x_2-x_3$ plane at regular intervals in the $x_1$ coordinate. For both of these models the appropriately scaled topological braid entropy $h_{\textit{braid}}$ was also found to converge to the random braid entropy $\langle \lambda _\sigma \rangle$ , leading to a simple relationship between the Lyapunov exponent $\lambda _\infty$ of the flow and the transverse Péclet number $\textit{Pe}_T$ as

(3.14) \begin{equation} h=\lambda _\infty =f^{1/d}\langle \lambda _\sigma \rangle =\left (\frac {d}{\textit{Pe}_T}\right )^{1/d}\langle \lambda _\sigma \rangle ,\quad d=1,2, \end{equation}

which generalises (3.13) to $d=1$ or $d=2$ dimensional streamline arrays. The persistence of this relationship across these diverse models suggests that they all belong to the same universality class (Ódor Reference Ódor2004) associated with streamline braiding. This establishes that chaotic advection and transverse dispersion are intimately linked in heterogeneous Darcy flow because they are both driven by non-trivial streamline braiding. Equation (3.14) also points to development of methods to estimate $\lambda _\infty$ from transverse dispersivity data for specific systems; however, further research is required to extend this link to non-ideal systems.

3.4. Fluid deformation in heterogeneous Darcy flow

In addition to streamline braiding, it is also instructive to consider ab initio evolution of fluid deformation in heterogeneous Darcy flow from a kinematic perspective. This can be quantified by considering deformation as a random process. Le Borgne et al. (Reference Le Borgne, Dentz and Carrera2008a , Reference Le Borgne, Dentz and Carrerab ) established that for steady flow in random media, the velocity magnitude $v$ decorrelates with distance $s$ along streamlines and is described by a spatial Markov process with respect to spatial correlation length $\ell$ . These studies were applied to (multi-log Gaussian) hydraulic conductivity fields with a well-defined minimum correlation length scale, similar to those considered in § 5. It has also been shown (Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2022) that the velocity gradient in steady Darcy flow is also spatially Markovian, hence fluid deformation is characterised by sampling the dimensionless velocity gradient tensor $\boldsymbol \epsilon (t;\boldsymbol{X})\equiv \boldsymbol{\nabla }\boldsymbol{v}(\boldsymbol{x}(t;\boldsymbol{X}))^\top \ell /\langle v_1\rangle$ equidistantly with distance $s$ along streamlines. Rotation into the Protean coordinate frame $\boldsymbol{x}^\prime$ (Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018a ) renders the velocity gradient tensor $\boldsymbol \epsilon ^\prime (t;\boldsymbol{X})$ upper triangular (see Appendix B for details) and the fluid deformation gradient tensor $\boldsymbol{F}^\prime (t;\boldsymbol{X})$ which evolves as

(3.15) \begin{equation} \frac {{\rm d}\boldsymbol{F}^\prime (t;\boldsymbol{X})}{{\rm d}t}=\boldsymbol \epsilon ^\prime (t;\boldsymbol{X})\boldsymbol{\cdot }\boldsymbol{F}^\prime (t;\boldsymbol{X}),\quad \boldsymbol{F}^\prime (0;\boldsymbol{X})=\boldsymbol{1}, \end{equation}

is likewise. The diagonal elements $F_\textit{ii}^\prime$ represent principal stretches and the off-diagonal components $F_{\textit{ij}}^\prime$ (non-zero for $j\gt i$ ) represent shear deformations. Hence, the hyperbolic stretches that govern the Lyapunov exponents are given by the diagonal components $\epsilon _\textit{ii}^\prime$ , and are not conflated with shear deformations given by the off-diagonal components $\epsilon _\textit{ii}^\prime$ . From (3.15), the diagonal elements of $\boldsymbol{F}^\prime (t;\boldsymbol{X})$ grow exponentially as

(3.16) \begin{align} F_\textit{ii}^\prime (t;\boldsymbol{X}) =\exp \left [\int _0^t {\rm d}t^\prime \epsilon _\textit{ii}^\prime (t^\prime ;\boldsymbol{X})\right ]\!, \quad i=1,2,3. \end{align}

For $i=1$ , this yields

(3.17) \begin{equation} F_{11}^\prime (t;\boldsymbol{X}) =\exp \left [\int _0^t {\rm d}t^\prime \frac {\partial v}{\partial s}\right ] =\exp \left [\int _{v(0;\boldsymbol{X})}^{v(t;\boldsymbol{X})} {\rm d}v^\prime \frac {1}{v}\right ]=\frac {v(t;\boldsymbol{X})}{v(0;\boldsymbol{X})}, \end{equation}

where $s$ is the distance along the streamline labelled by $\boldsymbol{X}$ and $v(t;\boldsymbol{X})={\rm d}s/{\rm d}t$ . From (3.17), fluid stretching in the streamwise direction simply fluctuates with the local velocity as $F_{11}^\prime (t)=v(t)/v(0)$ , hence $\langle \epsilon ^\prime _{11}(t)\rangle =0$ . As $\sum _{i}\epsilon ^\prime _\textit{ii}=0$ for divergence-free flow flows, the Lyapunov exponent is then

(3.18) \begin{equation} \lambda _\infty =\langle \epsilon ^\prime _{22}\rangle =-\langle \epsilon ^\prime _{33}\rangle , \end{equation}

where $\langle \boldsymbol{\cdot }\rangle$ denotes spatial averaging along streamlines. Hence, $F_{22}^\prime \sim \exp (\lambda _\infty t)$ and $F_{33}^\prime \sim \exp (-\lambda _\infty t)$ , and so chaotic advection in steady 3-D flow is characterised by the single Lyapunov exponent $\lambda _\infty$ .

In Appendix C we show that (3.16) leads to an ab initio continuous time random walk (CTRW) for the relative stretch $\rho \equiv \delta l(\boldsymbol{X},t)/\delta l(\boldsymbol{X},0)$ of infinitesimal fluid line elements along streamlines as

(3.19) \begin{align} \begin{split} \textrm{ln}\, \rho _{n+1} &= \textrm{ln}\, \rho _n + \frac {\ell \,\epsilon _n}{v_n},\\ t_{n+1} &= t_n + \tau _n = t_n + \frac {\ell }{v_n},\\ s_{n+1} &=s_n+\ell , \end{split} \end{align}

where $v_n$ and $\tau _n$ respectively are the streamline velocity magnitude and advection time between $s_n$ and $s_{n+1}$ and $\epsilon _n$ is the relevant Protean velocity gradient component $\epsilon ^\prime _{22}$ . The distribution $\psi (\tau )$ of waiting times $\tau$ is related to the Lagrangian velocity distribution $p_s(v)$ as

(3.20) \begin{equation} \psi (\tau )=\frac {\ell p_s(\ell /\tau )}{\tau ^2}, \end{equation}

and $p_s(v)=v p_e(v)/\langle v \rangle _e$ (Dentz et al. Reference Dentz, Lester, Le Borgne and de Barros2016), where $p_e(v)$ is the Eulerian velocity distribution and $\langle \boldsymbol{\cdot }\rangle _e$ represents an Eulerian average. For heterogeneous Darcy flow the Lagrangian velocity distribution may be heavy-tailed in the small velocity limit (Berkowitz et al. Reference Berkowitz, Cortis, Dentz and Scher2006), i.e. $p_s(v)\propto v^{\beta -1}$ , $\beta \gt 1$ for $v\ll \langle v\rangle$ depending on the distribution of hydraulic conductivity (Hakoun, Comolli & Dentz Reference Hakoun, Comolli and Dentz2019). Thus, the waiting time distribution scales as (Dentz et al. Reference Dentz, Lester, Le Borgne and de Barros2016)

(3.21) \begin{equation} \psi (\tau )\approx \frac {\psi _0}{|\varGamma (-\beta )|} t^{-1-\beta }\quad \text{for}\,\,\tau /\langle \tau \rangle \gg 1, \end{equation}

where $\varGamma$ is the Gamma function and the $q$ th-order moments $\langle \tau ^q\rangle$ of $\psi (\tau )$ are finite for $q\lt \lfloor \beta \rfloor$ . Hence, normal transport arises if $\beta \gt 2$ , as the mean and variance of $\psi (\tau )$ are bounded, but anomalous transport arises for $1\lt \beta \lt 2$ as the variance of $\psi (\tau )$ is unbounded. For $\beta \gt 1$ , the advection time $t_n$ may be related to the number of increments $n$ as

(3.22) \begin{equation} t_n=\sum _{i=1}^n\tau _i\approx n\langle \tau \rangle = n\ell \left \langle \frac {1}{v_n}\right \rangle =\frac {n\ell }{\langle v\rangle _e}=n\tau _c, \end{equation}

where $\tau _c\equiv \ell /\langle v\rangle _e$ is the mean transition time. Similarly, the mean of ${ln} \rho$ also grows linearly in time as $\langle {ln} \rho \rangle =\langle \epsilon _n\rangle t=\lambda _\infty t$ . The rate of growth of the variance $\sigma ^2_{{ln} \rho }$ serves as an important input parameter for lamellar mixing theories (Le Borgne, Dentz & Villermaux Reference Le Borgne, Dentz and Villermaux2015; Villermaux Reference Villermaux2019) that predict evolution of the concentration probability distribution function (PDF) for advective–diffusive solute transport from fluid stretching behaviour. Rebenshtok et al. (Reference Rebenshtok, Denisov, Hänggi and Barkai2014) and Dentz et al. (Reference Dentz, Le Borgne, Lester and de Barros2015) show that in the asymptotic limit the coupled CTRW (3.19) generates $q$ th-order absolute central moments $\mu _q$ for $q\gt \beta$ which scale as

(3.23) \begin{equation} \mu _q(t)\equiv \lim _{t\rightarrow \infty }\langle |\textrm{ln}\, \rho -\lambda _\infty t |^q\rangle \sim \begin{cases} t^{1+q-\beta }\quad &\text{for}\,\,1\lt \beta \lt 2,\\ t^{q/2}\quad &\text{for}\,\,\beta \gt 2, \end{cases} \end{equation}

hence anomalous transport ( $\beta \lt 2$ ) in heterogeneous 3-D Darcy flow can also generate anomalous stretching dynamics. In Appendix C we show for the case $\beta \gt 2$ , $p_{{ln} \rho }(\textrm{ln}\, \rho )$ converges to a normal distribution with mean and variance

(3.24) \begin{align} \lim _{t\rightarrow \infty }\langle \textrm{ln}\, \rho \rangle =\lambda _\infty t \quad \lim _{t\rightarrow \infty }\sigma ^2_{{ln} \rho }=\sigma ^2_\lambda t, \end{align}

where $\sigma ^2_\lambda \equiv \sigma ^2_\epsilon \langle \tau ^2\rangle /\tau _c$ and $\sigma _\epsilon ^2$ is the variance of $p_{\epsilon }(\epsilon )$ . For the anomalous stretching case ( $1\lt \beta \lt 2$ ), the variance $\sigma ^2_{{ln} \rho }$ growth ranges from normal to ballistic as (Dentz et al. Reference Dentz, Le Borgne, Lester and de Barros2015)

(3.25) \begin{align} \lim _{t\rightarrow \infty }\sigma ^2_{{ln} \rho }=\frac {\psi _0\,\beta \,\varGamma (\beta +1)}{6-5\beta +\beta ^2}t^{3-\beta }. \end{align}

Under such anomalous stretching the equality (3.2) between $h$ and $\lambda _\infty$ persists as temporal averaging along streamlines still dominates (Lester & Dentz Reference Lester and Dentz2025), but convergence of (3.1) slows as $\beta \downarrow 1$ and physical processes such as mixing and dispersion may be impacted in this limit. In § 5, we find that the flows considered herein all exhibit normal transport. However, anomalous stretching dynamics form an important kinematic regime that warrants future investigation.

4.1. Chaotic advection in periodic Darcy flow

We first examine steady 3-D Darcy flows with a smooth, finite, heterogeneous and anisotropic hydraulic conductivity field $\boldsymbol{K}(\boldsymbol{x})$ that is $P$ -periodic in the direction of the mean potential gradient $\boldsymbol{g}\equiv -\langle \boldsymbol{\nabla }\phi \rangle /||\langle \boldsymbol{\nabla }\phi \rangle ||$ , i.e.

(4.1) \begin{equation} \boldsymbol{K}(\boldsymbol{x})=\boldsymbol{K}\left (\boldsymbol{x}+nP\, \boldsymbol{g}\right ),\quad n=1,2,\ldots .\end{equation}

The transform (3.6) ensures this flow has the same dynamics as a divergence-free temporally periodic flow (Bajer Reference Bajer1994; Speetjens et al. Reference Speetjens, Metcalfe and Rudman2021), in that the advection equations can be cast as a time-periodic Hamiltonian system (often referred to as a $1({1}/{2})$ degree of freedom Hamiltonian system). As such, the established tools and techniques of Hamiltonian chaos (Ottino Reference Ottino1989; Katok & Hasselblatt Reference Katok and Hasselblatt1995) can be used to understand the mechanisms leading to chaotic mixing in these flows, including construction of a stroboscopic map known as a Poincaré section of the flow, which is generated by recording the position of streamlines in the flow in the plane normal to $\boldsymbol{g}$ at integer multiples of $P$ . From the Brouwer fixed point theorem (Brouwer Reference Brouwer1911), there must exist period- $k$ points $\boldsymbol{x}_k$ (with $k=1,2,\ldots$ ) in the Poincaré section, such that streamlines at $\boldsymbol{x}_k$ are advected downstream to $\boldsymbol{x}_p+kP\,\boldsymbol{g}$ . These dynamics and features are exactly the same as those that arise in steady 3-D duct flows (Bajer Reference Bajer1994; Speetjens et al. Reference Speetjens, Metcalfe and Rudman2021).

Each periodic point $\boldsymbol{x}_k$ belongs to a corresponding periodic orbit (streamline). As per Hamiltonian dynamical systems theory (Katok & Hasselblatt Reference Katok and Hasselblatt1995), a periodic point may be classified as either a hyperbolic saddle point $\boldsymbol{x}_H$ if the local fluid deformation over the $k$ -period involves fluid stretching and contraction that is exponential in time, or an elliptic point $\boldsymbol{x}_E$ if the local fluid deformation involves rotation only. The number of hyperbolic ( $n_H$ ) and elliptic ( $n_E$ ) points is governed by the Poincaré–Hopf theorem, which may be expressed as

(4.2) \begin{equation} n_E-n_H=2(1-g), \end{equation}

where $g$ is the topological genus of the Poincaré section (which may be 0, 1 or 2 depending upon its periodicity). Elliptic points are associated with non-chaotic (non-mixing) regions of the flow known as Kolmogorov–Arnold–Moser (KAM) islands, where the local Lyapunov exponent is zero. For hyperbolic points, the directions of the exponentially stretching and contracting are respectively associated with hyperbolic unstable and stable manifolds emanating from $\boldsymbol{x}_k$ . If these manifolds intersect transversely (which arises in all but highly symmetric systems (Haller & Poje Reference Haller and Poje1998)), then a heteroclinic tangle results, leading to chaotic advection (Ottino Reference Ottino1989; Katok & Hasselblatt Reference Katok and Hasselblatt1995) in these regions, with positive Lyapunov exponent.

The transition from globally non-chaotic to globally chaotic dynamics in Hamiltonian systems is well established, with three predominant ‘routes to chaos’ identified (period doubling, quasi-periodicity and intermittency) as system parameters are altered. In the case of steady 3-D Darcy flows, these parameters are the degree of heterogeneity ( $\sigma ^2_{\textit{ln} K}$ ) or anisotropy ( $\delta$ in § 5) of the conductivity field $\boldsymbol{K}$ as either homogeneous or isotropic fields do not generate chaos. Under the period-doubling route, the increase of such control parameters leads to bifurcation of elliptic points into two elliptic and one hyperbolic point (thus maintaining (4.2)), and the stable and unstable manifolds associated with this hyperbolic point connect smoothly. However, with a further increase of the control parameters, these manifolds then intersect transversely, leading to a heteroclinic tangle, the hallmark of chaotic dynamics in Hamiltonian systems. This transition is associated with the formation of a chaotic ‘sea’ surrounding the two KAM ‘islands’ associated with the elliptic points. This process continues, with further bifurcation of elliptic points and growing chaotic regions associated with hyperbolic points until the domain is essentially filled with a chaotic sea and the KAM islands surrounding elliptic points are vanishingly small. In § 5 we demonstrate this route to chaos with increasing conductivity heterogeneity, including identification of periodic points and associated structures (KAM islands, hyperbolic manifolds).

4.2. Chaotic advection in random Darcy flow

For random Darcy flow, this picture is altered somewhat as there do not exist persistent features such as hyperbolic and elliptic periodic points, KAM islands and hyperbolic manifolds to elucidate the mixing dynamics of the flow. Instead, finite-time analogues of these features – elliptic, hyperbolic and parabolic LCS (Haller Reference Haller2015) – serve to organise fluid motion and uncover the Lagrangian dynamics of the flow over finite time scales. One major difference for random flows is that all particle trajectories are now ergodic (Liu, Muzzio & Peskin Reference Liu, Muzzio and Peskin1994; Poje, Haller & Mezić Reference Poje, Haller and Mezić1999; Kang et al. Reference Kang, Singh, Kwon and Anderson2008), meaning that statistics gathered over long times are equivalent to ensemble statistics as trajectories sample the entire phase space of the flow. This means that e.g. KAM islands cannot persist and all trajectories are chaotic, rather than a topologically distinct (in the Lagrangian frame) set of mixing and non-mixing regions. Hence, while elliptic and hyperbolic regions of the flow can be identified over finite periods, in the infinite time limit random flows are globally chaotic, and the finite-time Lyapunov exponent (FTLE) of each 3-D streamline converges to the unique global Lyapunov exponent of the flow (due to ergodicity). As the flow control parameter is increased from zero to finite values, the flow transitions form integrable (regular) dynamics to fully chaotic, and the Lyapunov exponent steadily increases from zero to finite values. A similar picture emerges if the flow is analysed via the braiding analogue of the FTLE, finite-time braiding exponent (FTBE) (Budišić & Thiffeault Reference Budišić and Thiffeault2015). Similar to the FTLE, for periodic flows the FTBE is zero inside KAM islands as the streamlines only undergo trivial braiding, and the FTBE is positive in chaotic regions due to non-trivial braiding with positive topological entropy. Conversely, for random chaotic Darcy flows the FTBE is positive everywhere at long times as all streamlines undergo non-trivial braiding due to ergodicity. As such, we make the conjecture that, similar to other aperiodic Hamiltonian systems (Liu et al. Reference Liu, Muzzio and Peskin1994; Poje et al. Reference Poje, Haller and Mezić1999), if an aperiodic Darcy flow exhibits chaotic trajectories, it must be globally chaotic due to ergodicity of streamlines. This conjecture has important implications for interpreting the results in § 5.

4.3. Chaotic advection from non-finite or non-smooth hydraulic conductivity

Although beyond the scope of this study, it is useful to briefly discuss the generation of chaos in steady Darcy flows with stagnation points or non-smooth hydraulic conductivity. For steady 3-D Darcy flow with non-finite conductivity fields (associated with impermeable inclusions) or boundary conditions or flow sources and sinks that generate stagnation points, chaotic advection can also arise in the absence of anisotropy of the hydraulic conductivity field. For these flows, as $v_\phi \downarrow 0$ at stagnation points $\boldsymbol{x}_p$ , the rescaling (3.7) breaks down (Bajer Reference Bajer1994), leading to chaotic advection via a similar mechanism to that of pore-scale flow (Lester et al. Reference Lester, Metcalfe and Trefry2013). Here, exponential stretching of fluid elements local to saddle-type stagnation points (Surana, Grunberg & Haller Reference Surana, Grunberg and Haller2006) (an analogue of hyperbolic periodic points) form stable and unstable hyperbolic manifolds in the flow. Under symmetry breaking, these hyperbolic manifolds form a heteroclinic tangle, the hallmark of chaos in continuous dynamical systems (Ottino Reference Ottino1989; Katok & Hasselblatt Reference Katok and Hasselblatt1995). For steady Darcy flows with isotropic but non-smooth heterogeneous conductivity fields $\boldsymbol{K}$ , the velocity gradient, vorticity and helicity are undefined over discontinuities in $\boldsymbol{\nabla }\boldsymbol{K}$ , leading to ‘leakage’ of streamlines from coherent streamsurfaces and the potential for chaos (Lester et al. Reference Lester, Dentz, Singh and Bandopadhyay2023). This mechanism is not well understood and it is currently unclear whether such discontinuities in $\boldsymbol{\nabla }\boldsymbol{K}$ are physically realistic, or a consequence of the breakdown of the Darcy approximation at small scales.

5.1. Onset of chaotic advection with anisotropy

The results in § 3 uncover a deep link between chaotic advection and transverse dispersion in heterogeneous Darcy flow, and show how anomalous transport generates anomalous stretching dynamics. To examine the onset of chaotic advection with medium anisotropy, we first consider Darcy flow in the simplest possible conductivity field that admits non-zero helicity

(5.1) \begin{equation} \boldsymbol{K}(\boldsymbol{x})=k_0(\boldsymbol{x})\boldsymbol{I}+\delta [k_\delta (\boldsymbol{x})-k_0(\boldsymbol{x})]\hat {\boldsymbol{e}}_1\otimes \hat {\boldsymbol{e}}_1, \end{equation}

where $k_0(\boldsymbol{x})\neq k_\delta (\boldsymbol{x})$ and $\delta \in [0,1]$ quantifies anisotropy of the conductivity tensor. Although this flow has been previously considered (Lester et al. Reference Lester, Metcalfe and Trefry2024) in the context of quantifying the link between stirring and dispersion, here, we focus on the onset of chaotic dynamics with $\delta$ . Note that, while similar conductivity tensor fields such as $\boldsymbol{K}(\boldsymbol{x})=k_0(\boldsymbol{x})(\boldsymbol{I}+\delta \hat {\boldsymbol{e}}_1\otimes \hat {\boldsymbol{e}}_1)$ have non-zero helicity density $\mathcal{H}(\boldsymbol{x})\neq 0$ , the resulting flows are shown (Appendix A) to be epi-two-dimensional (Yoshida & Morrison Reference Yoshida and Morrison2017) and hence non-chaotic (Arnol’d Reference Arnol’d1965; Holm & Kimura Reference Holm and Kimura1991).

Darcy flow over the conductivity field (5.1) is driven by a unit potential gradient $\boldsymbol{\nabla }\bar {\phi }=\{-1,0,0\}$ in a triply periodic unit cube (3-torus $\mathbb{T}^3$ ) $\varOmega :\boldsymbol{x}\in [0,1]\times [0,1]\times [0,1]$ , and the scalar log-conductivity fields $f(\boldsymbol{x})\equiv {ln} k(\boldsymbol{x})$ for the independent fields $k_0(\boldsymbol{x})$ , $k_\delta (\boldsymbol{x})$ are given by

(5.2) \begin{align} f(\boldsymbol{x})&=\varSigma ^2_{\textit{ln} K}\sum _{n=1}^{N_i}\sum _{i,j,k}^{N} \frac {A_{n,\textit{ijk}}}{\sqrt {i^2+j^2+k^2}} \cos \big[2\pi i \big(x_1+\chi ^1_{n,\textit{ijk}}\big)\big]\nonumber\\ &\cos \big[2\pi j \big(x_2+\chi ^2_{n,\textit{ijk}}\big)\big]\cos \big(2\pi k \big(x_3+\chi ^3_{n,\textit{ijk}}\big)\big), \end{align}

with $N=4$ and so the velocity correlation length scale is $\ell =1/(2N)$ . Here, $N_i=2$ is the number of realisations in each mode, and $A_{n,\textit{ijk}}$ and $\chi ^m_{n,\textit{ijk}}$ with $m=1,2,3$ are uniformly distributed random variables in $[0,1]$ . The coefficient $\varSigma _{\textit{ln} K}$ is chosen such that the log variance of the conductivity field is $||f^2||_\varOmega =\sigma ^2_{\textit{ln} K}=4$ . A typical scalar field $f(\boldsymbol{x})$ with these parameters is shown in figure 1(a). Note fields with $N=1$ do not generate a chaotic dynamics due to symmetry of the velocity field.

To solve the divergence-free condition $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}=0$ over $\varOmega$ , $\phi$ is decomposed into a mean and fluctuation as

(5.3) \begin{equation} \phi (\boldsymbol{x})=\bar \phi (\boldsymbol{x})+\tilde \phi (\boldsymbol{x}), \end{equation}

where $\bar \phi =-x_1$ . From (2.2), $\tilde {\phi }(\boldsymbol{x})$ is governed by

(5.4) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{K}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {\phi }(\boldsymbol{x}))- \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{K}\boldsymbol{\cdot }\hat {\boldsymbol{e}}_1)=0. \end{equation}

Methods to solve (5.4) and perform streamline tracking are detailed in Appendix D, and typical potential fields and streamlines are shown in figure 1. For each value of $\delta$ , a realisation of $\boldsymbol{K}(\boldsymbol{x})$ is generated, (5.4) is solved and $10^3$ 3-D streamlines are computed for distance $10^4 \ell$ , along with the transverse dispersion coefficient $D_T$ , Protean velocity gradient tensor $\boldsymbol \epsilon ^\prime$ and Lyapunov exponent $\lambda _\infty =\langle \epsilon _{22}^\prime \rangle$ , as described in Appendix C.

Typical streamfunctions $\psi _1(\boldsymbol{x})$ , $\psi _2(\boldsymbol{x})$ for the helicity-free flow $\delta =0$ given by (2.7) (see Appendix D) for solution details) are shown in figure 6(a). Streamlines (green) are confined to the intersections of the level sets of $\psi _1$ and $\psi _2$ , and their global behaviour is shown in figure 1(c). Despite their significant tortuosity, these confined streamlines do not exhibit persistent dispersion or non-trivial braiding. Conversely, for $\delta \gt 0$ the associated streamlines (red) in figure 6(a) are unconfined and so exhibit streamline braiding and transverse dispersion, see also figure 1(d) for $\delta =1$ .

Figure 6(b) shows that the mean helicity magnitude $\langle |\mathcal{H}|\rangle$ increases from $\mathcal{H}(\boldsymbol{x})=0$ everywhere for $\delta =0$ to a plateau $\langle |\mathcal{H}|\rangle \approx 4.74$ at $\delta =0.9$ , then suddenly increases to $\langle |\mathcal{H}|\rangle \approx 7.16$ at $\delta =1$ . This sharp increase as $\delta \uparrow 1$ is attributed to the loss of correlation between the $K_{11}$ and $K_{22}=K_{33}$ fields. The linear growth of $\langle |\mathcal{H}|\rangle$ for $\delta \ll 1$ is explained by decomposing the potential field $\phi$ as $\phi (\boldsymbol{x})=\phi _0(\boldsymbol{x})+\delta \,\phi _\delta (\boldsymbol{x})$ . To leading order $\boldsymbol{v}(\boldsymbol{x})$ is then

(5.5) \begin{align} \boldsymbol{v}(\boldsymbol{x})=\boldsymbol{v}_0(\boldsymbol{x})+\delta \boldsymbol{v}_\delta (\boldsymbol{x}),=-k_0\boldsymbol{\nabla }\phi _0-\delta [k_0\boldsymbol{\nabla }\phi _\delta -k_\delta \hat {\boldsymbol{e}}_1(\hat {\boldsymbol{e}}_1\boldsymbol{\cdot }\boldsymbol{\nabla }\phi _0)]+\mathcal{O}(\delta ^2),\end{align}

where $\boldsymbol{v}_0(\boldsymbol{x})$ is helicity free and the helical perturbation $\boldsymbol{v}_\delta (\boldsymbol{x})$ is generated by the difference $k_0-k_\delta$ , as is shown by inserting (5.5) into $\mathcal{H}(\boldsymbol{x})\equiv \boldsymbol{v}\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{v})$ and truncating to order $\delta$ to give

(5.6) \begin{align} \mathcal{H}(\boldsymbol{x})&=\delta [k_0\boldsymbol{\nabla }\phi \boldsymbol{\cdot }\big(\hat {\boldsymbol{e}}_1\times \boldsymbol{\nabla }(k_\delta -k_0)(\boldsymbol{\nabla }\phi \boldsymbol{\cdot }\hat {\boldsymbol{e}}_{1})\big)\nonumber\\ &\quad +(k_\delta -k_0)(\boldsymbol{\nabla }\phi \boldsymbol{\cdot }\hat {\boldsymbol{e}}_{1}\otimes \hat {\boldsymbol{e}}_1\boldsymbol{\cdot }(\boldsymbol{\nabla } k_0\times \boldsymbol{\nabla }\phi ))]+\mathcal{O}(\delta ^2). \end{align}

Figure 6. (a) Perturbation of $\delta =0.1$ streamlines (red) from $\delta =0$ (zero helicity) streamlines (green) and associated $\psi _1$ streamsurface (blue) for the conductivity field given in (5.1). Similar perturbation of $\delta =0.1$ streamlines away from the $\psi _2$ streamsurfaces (not shown) also occurs. (b) Growth of mean absolute helicity $\langle |\mathcal{H}|\rangle$ with $\delta$ , inset shows $|\mathcal{H}|$ fields for $\delta =0.9, 1$ (adapted from Lester et al. Reference Lester, Metcalfe and Trefry2024). (c) Growth of transverse dispersion coefficient $D_T/\langle v_1\rangle \ell$ with $\delta$ from simulations (red points) and fitted exponential (5.7) (red curve). (d) Growth of Lyapunov exponent $\lambda _\infty$ with perturbation parameter $\delta$ from simulations (black points) and fitted exponential (5.8) (blue curve). Also shown (red dotted curve) is the dimensionless Lyapunov exponent $\lambda _\infty$ predicted from fitted exponential in (b) and (3.14). (c), (d) Adapted from (Lester et al. Reference Lester, Metcalfe and Trefry2024).

Figure 6(c) shows $D_T$ increases exponentially with $\delta$ as

(5.7) \begin{equation} \frac {D_T}{\langle v_1\rangle \ell }=\frac {1}{\textit{Pe}_T}\approx 4.343\times 10^{-6}\big(e^{10.214\,\delta }-1\big), \end{equation}

which is consistent with (2.7) as $D_T=0$ for $\delta =0$ (Lester et al. Reference Lester, Dentz, Singh and Bandopadhyay2023), indicating finite transverse dispersion arises for weak perturbations away from heterogeneous isotropic Darcy flow. Figure 6(d) shows that $\lambda _\infty$ also increases exponentially with $\delta$ as

(5.8) \begin{equation} \lambda _\infty \approx 0.00381\big(e^{5.283\,\delta }-1\big), \end{equation}

and is also non-zero for small $\delta \gt 0$ , indicating that chaotic advection occurs for weak perturbations away from isotropic Darcy flow. For $\delta =0$ the streamfunction formulation (2.7) enforces zero Lyapunov exponent (Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2022). Figure 6(d) also shows that insertion of the fitted exponential (5.8) for $D_T$ into a similar relationship to (3.14)

(5.9) \begin{equation} \lambda _\infty =2\sqrt {\frac {D_T}{\langle v_1\rangle \ell }} ,\end{equation}

yields excellent agreement with the measured Lyapunov exponent in (5.7). The different proportionality between $\lambda _\infty$ and $D_T$ to (3.14) is expected as the domain $\varOmega$ is finite and so does not belong to a universality class (Ódor Reference Ódor2004). These results indicate that chaotic advection is inherent to anisotropic heterogeneous Darcy flow. If the conjecture raised in § 4.2 is true, this statement also holds for weakly anisotropic random heterogeneous conductivity fields.

5.2. Onset of chaotic advection with medium heterogeneity

Figure 7. (Top row) Poincaré sections and (bottom row) FTLE fields with superposed Poincaré sections (colours used to distinguish different streamlines) of steady 3-D anisotropic Darcy flow ranging from weakly ( $\sigma ^2_{\textit{ln} K}=2^{-4}$ ) to moderately ( $\sigma ^2_{\textit{ln} K}=1$ ) heterogeneous conductivity fields. Blue and red points respectively denote elliptic ( $\boldsymbol{x}_E$ ) and hyperbolic ( $\boldsymbol{x}_H$ ) points. Note for more heterogeneous conductivity fields $\sigma ^2_{\textit{ln} K}\geqslant 2$ (not shown), the KAM islands shrink to infinitesimal size, indicating essentially globally chaotic dynamics.

To examine the impact of medium heterogeneity upon anisotropic Darcy flow, we consider flow generated by the anisotropic diagonal conductivity tensor $\boldsymbol{K}(\boldsymbol{x})$ in (2.4), where the scalar conductivity fields $K_\textit{ii}$ are generated in the same manner as $k_0$ , $k_\delta$ but with $N_i=4$ . The heterogeneity of these fields $K_\textit{ii}$ is varied from weakly ( $\sigma ^2_{\textit{ln} K}\lt 1$ ) to strongly ( $\sigma ^2_{\textit{ln} K}\gt 1$ ) heterogeneous over the range of log variances $\sigma ^2_{\textit{ln} K}=$ $(2^{-10},2^{-8},\ldots , 2^{-2})$ $\cup$ $(1/2,1,2,3,4)$ . For each value of $\sigma ^2_{\textit{ln} K}$ , a realisation of $\boldsymbol{K}(\boldsymbol{x})$ is generated and (5.4) is solved, and streamlines, $D_T$ and $\lambda _\infty$ are computed in the same manner as described in § 5.1.

Figure 8. (a) Variation of normalised mean longitudinal velocity $\langle v_1\rangle /v_h$ (red points) with log variance $\sigma ^2_{\textit{ln} K}$ and linear fit (5.14), (grey dashed line) in anisotropic Darcy flow. (Inset) Variation of small velocity scaling index $\beta$ with $\sigma ^2_{\textit{ln} K}$ . (b) The PDFs of normalised velocity magnitude $p_v(v/\langle v\rangle )$ as a function of log variance $\sigma ^2_{\textit{ln} K}$ and fitted log-normal distribution (black lines). (c) Variation of dimensionless Lyapunov exponent $\lambda _\infty$ (red dots) with log variance $\sigma ^2_{\textit{ln} K}$ and nonlinear fit (5.19), (grey dashed line). (d) Variation of dimensionless stretching variance $\sigma ^2_\lambda$ with log variance $\sigma ^2_{\textit{ln} K}$ and nonlinear fit (5.16), (grey dashed line). (Inset) Variation of Protean velocity gradient variance $\sigma ^2_\epsilon$ with log variance $\sigma ^2_{\textit{ln} K}$ and linear fit (5.15), (grey dashed line).

Figure 7 shows both contour plots of the FTLE $\lambda (t,\boldsymbol{X})$ and Poincaré sections and low-order ( $k=1,2,3$ ) periodic points for the cases $\sigma ^2_{\textit{ln} K}=2^{-4},2^{-1},1$ , where the log variance $\sigma ^2_{\textit{ln} K}$ acts as a control parameter that governs perturbation from the non-chaotic (integrable) state corresponding to homogeneous anisotropic Darcy flow with $\sigma ^2_{\textit{ln} K}=0$ . The FTLE field is computed using the LCS Tool package (Onu, Huhn & Haller Reference Onu, Huhn and Haller2015 for the analogous 2-D unsteady flow via the transform (3.6)), where the FTLE field is computed over one period of the flow $t\in [0,1]$ as

(5.10) \begin{equation} \lambda (t;\boldsymbol{X})=\frac {1}{2t}{ln} \sigma _{\textit{max}}(t;\boldsymbol{X}), \end{equation}

for $\boldsymbol{X}\in x_2\times x_3 = [0,1]\times [0,1]$ where $\sigma _{\textit{max}}(t;\boldsymbol{X})$ is the largest eigenvalue of the Cauchy–Green tensor $\boldsymbol{C}(t;\boldsymbol{X})\equiv \boldsymbol{F}(t;\boldsymbol{X})^\top \boldsymbol{\cdot }\boldsymbol{F}(t;\boldsymbol{X})$ . The Poincaré sections indicate that even for very weak heterogeneity ( $\sigma ^2_{\textit{ln} K}=2^{-4}$ ), the Lagrangian topology is rich, consisting of a large number of elliptic and hyperbolic points, with a greater number density than that suggested by the correlation length $\ell$ . For such weak heterogeneity, the Lagrangian dynamics are close to integrable, with stable and unstable manifolds (not shown) associated with hyperbolic points intersecting almost tangentially, corresponding to weakly chaotic dynamics, as reflected by the small Lyapunov exponent reported in figure 8(c). The integrable nature of this weakly heterogeneous flow is also established by considering the non-zero diagonal components $K_\textit{ii}(\boldsymbol{x})$ of $\boldsymbol{K}(\boldsymbol{x})$ in the limit $\sigma ^2_{\textit{ln} K}\rightarrow 0$ , such that $K_\textit{ii}(\boldsymbol{x})=\exp (\sigma _{\textit{ln} K}f_\textit{ii}(\boldsymbol{x}))\approx 1+\sigma _{\textit{ln} K}f_\textit{ii}(\boldsymbol{x})$ , hence

(5.11) \begin{equation} \boldsymbol{K}(\boldsymbol{x})\approx \boldsymbol{I}+\sigma _{\textit{ln} K}\boldsymbol{f}(\boldsymbol{x}),\quad \phi (\boldsymbol{x})\approx x_1+\sigma _{\textit{ln} K}\tilde \phi (\boldsymbol{x}), \end{equation}

where $\boldsymbol{f}(\boldsymbol{x})\equiv \sum _{i=1}^3 f_\textit{ii}(\boldsymbol{x})\hat {\boldsymbol{e}}_i\otimes \hat {\boldsymbol{e}}_i$ , and the potential fluctuation $\tilde \phi$ and velocity field satisfy

(5.12) \begin{align} {\nabla }^2\tilde \phi +\boldsymbol{\nabla } f_{11}(\boldsymbol{x})\boldsymbol{\cdot }\hat {\boldsymbol{e}}_1=0, && \boldsymbol{v}(\boldsymbol{x})=-\hat {\boldsymbol{e}}_1-\sigma _{\textit{ln} K}\big[f_{11}(\boldsymbol{x})\hat {\boldsymbol{e}}_1+\boldsymbol{\nabla }\tilde \phi\big]+\mathcal{O}\big(\sigma ^2_{\textit{ln} K}\big). \end{align}

Hence, to leading order, the perturbative flow is solely driven by $f_{11}(\boldsymbol{x})$ and the helicity density scales with $\sigma ^2_{\textit{ln} K}$ as

(5.13) \begin{equation} \mathcal{H}(\boldsymbol{x})=\sigma ^2_{\textit{ln} K}\,\hat {\boldsymbol{e}}_1\big(\boldsymbol{\nabla }\tilde \phi \times \boldsymbol{\nabla } f_{11}\big)+\mathcal{O}\big(\sigma ^3_{\textit{ln} K}\big), \end{equation}

so this flow is integrable in the limit $\sigma ^2_{\textit{ln} K}\rightarrow 0$ , and the perturbation that scales with $\sigma ^2_{\textit{ln} K}$ in (5.12) is heterogeneous and anisotropic. Hence, this perturbation represents a non-integrable perturbation away from a fully integrable state, and so shares the same basic features as non-integrable perturbation of integrable Hamiltonian systems. With increasing medium heterogeneity, ( $\sigma ^2_{\textit{ln} K}=2^{-1}$ ), the chaotic dynamics strengthen, with clearly visible stochastic layers (chaotic regions) between distinct KAM islands apparent in figure 7(b), but no bifurcation of low-order periodic points are apparent. As $\sigma ^2_{\textit{ln} K}$ grows to 1, evidence of period-doubling bifurcations in low-order periodic points become apparent as KAM islands break up into smaller islands. Other surviving integrable regions shrink further, before transitioning to essentially globally chaotic dynamics at $\sigma ^2_{\textit{ln} K}=2$ (not shown). These dynamics are also reflected by the FTLE distributions, which show some weakly negative regions that correspond to KAM islands as the analogous 2-D flow (3.6) is not divergence free but it is Hamiltonian (Bajer Reference Bajer1994). For weak ( $\sigma ^2_{\textit{ln} K}=2^{-4}$ ) to moderate ( $\sigma ^2_{\textit{ln} K}=2^{-1}$ ) medium heterogeneity, the basic structure of the FTLE field is invariant up to a scaling factor, reflecting influence of the hydraulic conductivity field upon the fluid stretching dynamics. This behaviour, including the transition from an integrable state to global chaos under a non-integrable perturbation (in this case an anisotropic and heterogeneous perturbation) is entirely consistent with the embedded Hamiltonian structure of steady 3-D Darcy flow (Speetjens et al. Reference Speetjens, Metcalfe and Rudman2021) discussed in § 3.2.

Figure 8(a) shows that the scaled mean longitudinal velocity $\langle v_1\rangle /v_h$ (where $v_h$ is the velocity for a homogeneous medium with the same mean permeability) increases linearly with log variance $\sigma ^2_{\textit{ln} K}$ as expected for small and moderate values of $\sigma ^2_{{ln} k}$ (Renard & De Marsily Reference Renard and de Marsily1997)

(5.14) \begin{equation} \frac {\langle v_1\rangle }{v_h}= 1+\alpha _1\,\sigma ^2_{\textit{ln} K}, \end{equation}

with $\alpha _1\approx 0.29178$ and coefficient of determination $R^2=0.999$ . Figure 8(b) shows that the Eulerian velocity PDF $p_v(v)$ roughly follows a log-normal distribution for all but strongly heterogeneous flows, and converges toward a delta function in the homogeneous limit $\sigma ^2_{\textit{ln} K}\rightarrow 0$ . For $\sigma ^2_{\textit{ln} K}\geqslant 1$ , the velocity PDF scales as a power law in the small velocity limit as $p_v(v)\propto v^{\beta -1}$ , and the inset in figure 8(a) shows that the index $\beta \gt 2$ and approaches $\beta \downarrow 2$ in the strongly heterogeneous regime, leading to normal transport for all $\sigma ^2_{\textit{ln} K}\geqslant 1$ . The persistence of normal transport in the strongly heterogeneous regime is attributed to the low probability of all three $K_\textit{ii}$ in (2.4) being simultaneously small. We also note that different random models such as Gamma-distributed conductivity fields have a greater propensity to generate anomalous transport (Hakoun et al. Reference Hakoun, Comolli and Dentz2019).

Figure 8(c) shows that $\lambda _\infty$ may converge to a constant value with increasing log variance $\sigma ^2_{\textit{ln} K}$ , and is well fitted by the simple nonlinear function (5.19) as described below. Although the numerical methods used do not allow flows in more heterogeneous media to be computed, asymptotic convergence of $\lambda _\infty$ with increasing $\sigma ^2_{\textit{ln} K}$ is argued on physical grounds. We consider the limit $\sigma ^2_{\textit{ln} K}\rightarrow \infty$ and partition $\boldsymbol{K}$ over $\varOmega$ into infinitely permeable $\varOmega _\infty$ ( $\text{tr}\,\boldsymbol{K}(\boldsymbol{x})\rightarrow \infty$ ) and impermeable $\varOmega _0$ ( $\text{tr}\,\boldsymbol{K}(\boldsymbol{x})\rightarrow 0$ ) subdomains. As $\varOmega$ is the 3-torus $\mathbb{T}^3$ with topological genus $g=3$ , the boundary $\partial \varOmega _{\infty 0}$ between the subdomains is also topologically complex and thus admits saddle-type stagnation points $\boldsymbol{x}_0$ on $\partial \varOmega _{\infty 0}$ as per the Poincaré–Hopf theorem. As discussed in § 4.3, these saddle points generate chaotic advection via the same mechanism as for pore-scale Stokes flow (Lester et al. Reference Lester, Metcalfe and Trefry2013). Due to the linearity of Darcy flow, the rescaled fluid velocity field $\boldsymbol{v}/\langle v_1\rangle$ and associated saddle points and Lyapunov exponent are all invariant in the limit $\sigma ^2_{\textit{ln} K}\rightarrow \infty$ .

As shown in the inset of figure 8(d), the velocity gradient variance $\sigma ^2_\epsilon$ (from data presented in § 5.3) grows linearly with conductivity log variance as

(5.15) \begin{equation} \sigma ^2_\epsilon \approx \alpha _2\sigma ^2_{\textit{ln} K}, \end{equation}

where $\alpha _2=1.1656$ with $R^2=0.999$ , and $\sigma ^2_\lambda$ grows exponentially with $\sigma ^2_{\textit{ln} K}$ as

(5.16) \begin{equation} \sigma ^2_\lambda \approx \alpha _3\big(e^{\alpha _4\sigma ^2_{\textit{ln} K}}-1\big), \end{equation}

where $\alpha _3=2.2429$ and $\alpha _4=0.8415$ with $R^2=0.999$ . These results show that, although the Lyapunov exponent $\lambda _\infty$ converges toward the upper bound ${ln} 2$ with increasing medium heterogeneity, the corresponding variance $\sigma ^2_\lambda$ increases due to growth of the second moment $\langle \tau ^2\rangle$ . For the case of anomalous transport ( $1\lt \beta \lt 2$ ), this variance grows super-linearly as $\sigma ^2_{{ln} \rho }\sim t^{3-\beta }$ Lester et al. (Reference Lester, Dentz, Le Borgne and de Barros2018a ).

In addition to direct computation of $\lambda _\infty$ , the topological braid entropy $h_{\textit{braid}}$ of the streamlines is computed directly via the E-tec routine (Roberts et al. Reference Roberts, Sindi, Smith and Mitchell2019), which is limited to moderately heterogeneous ( $0.1\leqslant \sigma ^2_{\textit{ln} K}\leqslant 2$ ) systems due to lack of convergence for $\sigma ^2_{\textit{ln} K}\lt 0.1$ and flow reversal for $\sigma ^2_{\textit{ln} K}\gt 1$ . In principle, streamlines for the strongly heterogeneous cases could be placed in the intrinsic coordinates $\boldsymbol \xi =(\chi _1,\chi _2,\phi )$ , however, determination of $\chi _1$ , $\chi _2$ is difficult as the non-zero helicity density $\mathcal{H}(\boldsymbol{x})\neq 0$ means that mutual Lie derivatives of the flow do not vanish and so there does not exist an intrinsic holonomic basis (Schutz Reference Schutz1980) (a coherent orthogonal coordinate system) for $(\chi _1,\chi _2)$ . Figure 9(a) shows that $\lambda _\infty$ and $h_{\textit{braid}}$ agree to within 5 %, verifying (3.2).

Figure 9. (a,b) Lyapunov exponents and (c,d) dispersion coefficients in (a,c) weakly ( $\sigma ^2_{\textit{ln} K}\leqslant 1$ ) and (b,d) strongly ( $\sigma ^2_{\textit{ln} K}\geqslant 1$ ) heterogeneous porous media with anisotropic conductivity given by (2.4). Red points in (a,b) indicate Lyapunov exponents computed from the Protean frame (Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018a ), black points indicate topological entropy computed using E-tec (Roberts et al. Reference Roberts, Sindi, Smith and Mitchell2019). Dashed lines in (a,c) and (b,d) respectively indicate power law (5.17), (5.21), and linear (5.18), (5.22), fits to the Lyapunov exponents and dispersion coefficients in the weakly and strongly heterogeneous regimes.

In the weakly heterogeneous regime ( $\sigma ^2_{\textit{ln} K}\leqslant 1$ ), $\lambda _\infty$ varies almost linearly with log variance as

(5.17) \begin{equation} \lambda _\infty \approx \alpha _6\big(\sigma ^2_{\textit{ln} K}\big)^{\alpha _5}, \end{equation}

with $\alpha _6=0.154$ , $\alpha _5=1.129$ and $R^2=0.995$ , suggesting that chaotic advection arises in weakly heterogeneous anisotropic conductivity fields. Figure 9(b) shows that in the strongly heterogeneous regime, $\lambda _\infty$ scales linearly with $\sigma ^2_{\textit{ln} K}$ as

(5.18) \begin{equation} \lambda _\infty \frac {\langle v_1\rangle }{v_h}\approx \alpha _7\,\sigma ^2_{\textit{ln} K}, \end{equation}

with $\alpha _7\approx 0.198$ and $R^2=0.999$ . Dividing by (5.14) yields

(5.19) \begin{equation} \lambda _\infty \approx \frac {\alpha _7 \sigma ^2_{\textit{ln} K}}{1+\alpha _1\sigma ^2_{\textit{ln} K}}, \end{equation}

which, as shown in figure 8(c), provides an excellent fit of $\lambda _\infty$ as a function of $\sigma ^2_{\textit{ln} K}$ across the entire regime. Note that (5.19) is compatible with (5.17) in the weakly heterogeneous limit $\sigma _{\textit{ln} K}^2\ll 1$ under the approximation $\alpha _5\approx 1$ as $\alpha _6\approx \alpha _7/(1+\alpha _1)=0.1538$ . Equation (5.19) also implies that $\lambda _\infty$ converges in the strongly heterogeneous limit

(5.20) \begin{equation} \lim _{\sigma ^2_{\textit{ln} K}\rightarrow \infty }\lambda _\infty \approx \frac {\alpha _7}{\alpha _1}= 0.677. \end{equation}

This value is close to the upper bound $h={\rm {ln}} 2\approx 0.693$ for steady 3-D flow (Dinh & Sibony Reference Dinh and Sibony2008), indicating strong chaotic mixing arises in strongly heterogeneous anisotropic Darcy flow.

Figure 9(c) shows that the transverse dispersion coefficient scales nonlinearly in the weakly heterogeneous regime as

(5.21) \begin{equation} \frac {D_T}{\langle v_1\rangle \ell }\sim \big(\sigma ^2_{\textit{ln} K}\big)^{2.01}, \end{equation}

and figure 9(d) shows that transverse dispersion coefficient scales linearly in the strongly heterogeneous regime as

(5.22) \begin{equation} \frac {D_T}{v_h\ell }\approx 0.0234\,\sigma ^2_{\textit{ln} K}, \end{equation}

with $R^2=1.000$ in both cases. In the weakly heterogeneous regime, the scalings of $D_T$ (5.21) and $\lambda _\infty$ (5.17) with $\sigma ^2_{\textit{ln} K}$ agree with the relationship (3.14) linking dispersion and chaotic advection in random flows. However, as expected this relationship does not persist in the strongly heterogeneous regime due to the onset of flow reversal. Similar to the findings in § 5.1, these results indicate that chaotic advection is inherent to heterogeneous anisotropic Darcy flow. If the conjecture raised in § 4.2 is true, this statement also holds for weakly heterogeneous random anisotropic conductivity fields.

5.3. Velocity gradient statistics

In addition to the Lyapunov exponent, the Protean velocity gradient $\boldsymbol \epsilon ^\prime$ provides important information regarding the deformation structure of heterogeneous Darcy flow. For all computed flows $\boldsymbol \epsilon ^\prime$ is sampled along $10^3$ streamlines at fixed spatial increment $\ell$ for a distance of $10^4\ell$ . Figures 10(a,c) and 10(b,d) respectively show the PDFs of the principal stretches $\epsilon _\textit{ii}^\prime$ and shears $\epsilon _{\textit{ij}}^\prime$ for strongly heterogeneous ( $\sigma ^2_{\textit{ln} K}=4$ ) Darcy flow for (a,b) fully anisotropic ( $\delta =1$ ) and (c,b) isotropic ( $\delta =0$ ) conductivity fields (5.1). For all flows computed the divergence error $|\sum _{i=1}^3\epsilon ^\prime _\textit{ii}|\lt 10^{-6}$ and the streamwise mean stretch $|\langle \epsilon ^\prime _{11}(t)\rangle |\lt 10^{-4}$ is also close to zero. For all flows the distributions of both $\epsilon ^\prime _\textit{ii}$ and $\epsilon ^\prime _{\textit{ij}}$ for $j\gt i$ , $i=1:3$ are broad and exhibit a sharp cutoff due to the finite nature of $\varOmega$ . In all cases, the standard deviation is large $\sigma _{\epsilon ^\prime _\textit{ii}}\gg |\langle \epsilon ^\prime _\textit{ii}\rangle |$ , however, the large number of independent observations ( $10^7$ ) reduces the standard error to $\sigma _{\epsilon _\textit{ii}}/\sqrt {n}\sim 10^{-3}$ , yielding accurate estimates of the principal stretches.

Figure 10(a,c) shows that the PDF of $\epsilon ^\prime _{11}$ is significantly broader for the isotropic case. This is attributed to the fact that in anisotropic media the local permeability varies with the local velocity orientation, leading to a narrower velocity distribution as the flow has a greater flexibility to minimise dissipation. The off-diagonal shears $\epsilon _{\textit{ij}}^\prime$ in figure 10(b,d) are all similarly distributed except the longitudinal shears $\epsilon _{12}^\prime$ , $\epsilon _{13}^\prime$ of the anisotropic flow are more broadly distributed as these streamlines are not confined to streamsurfaces and hence have higher curvature. The Protean transverse shear $\epsilon _{23}^\prime$ in the isotropic case is zero as this is related to the helicity density as $\mathcal{H}(\boldsymbol{x})=v^\prime _i\,\varepsilon _{ijk}\,\epsilon ^\prime _{jk}=v\,\epsilon _{23}^\prime$ (Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018a ).

Figure 10. The PDFs of (a,c) diagonal $\epsilon _\textit{ii}$ and (b,d) off-diagonal velocity gradient components $\epsilon _{\textit{ij}}$ for (a,b) heterogeneous anisotropic Darcy flow (5.1) with $\sigma _{\textit{ln} K}^2=4$ , $\delta =1$ and (c,d) isotropic Darcy flow (5.1) with $\sigma _{\textit{ln} K}^2=4$ , $\delta =0$ .

For isotropic flow, the Lyapunov exponent is effectively zero ( $|\langle \epsilon ^\prime _\textit{ii}\rangle |\lt 10^{-4}$ ), whereas for anisotropic flow it is slightly larger than the theoretical upper bound $\lambda _\infty ={ln} 2$ , $(\langle \epsilon ^\prime _{11}\rangle ,\langle \epsilon ^\prime _{22}\rangle ,\langle \epsilon ^\prime _{33}\rangle )$ = (0.0001, 0.7148, −0.7149). The broad nature of the velocity gradient PDFs leads to a relatively large stretching variance, $\sigma ^2_\lambda /\lambda _\infty \sim 10^2$ as shown in figure 8(d). This magnitude is consistent with the finding (discussed in § 3.1) that the ensemble average (3.3) $h=\lambda _\infty +\sigma ^2/2$ is not correct, as this would not yield the agreement between $h_{\textit{braid}}$ and $\lambda _\infty$ observed in figure 9(a) The broad distribution of all components of $\boldsymbol \epsilon ^\prime$ has significant impacts upon the range of fluid processes hosted in these flows, as shall be explored further in § 6.