2.1. Chaotic advection in fracture networks

To analyse advective mixing in fracture networks, we first consider the propensity for these systems to admit complex mixing dynamics and use this assessment to develop a parsimonious representation of advective mixing. We focus on fracture networks at negligible Reynolds number, $Re \ll 1$ , with large fracture aspect ratio $\ell \ll L$ . Under these conditions, inertia-driven mixing mechanisms, such as those that arise at fracture intersections (Lee & Kang Reference Lee and Kang2020; Yang et al. Reference Yang, Chen, Lee and Kang2024), are absent. However, complex mixing dynamics may also arise in inertialess flows in topologically complex domains.

In the Stokes regime, steady flow in fracture networks may generate chaotic advection in a manner analogous to porous networks or granular media (Lester et al. Reference Lester, Metcalfe and Trefry2013, Reference Lester, Heyman, Méheust and Le Borgne2025a ). Pore and fracture networks both involve merging and branching of streamlines at pore junctions and fracture intersections, respectively, giving rise to localised stagnation points $\boldsymbol{x}_p$ , and associated unstable and stable 2-D hyperbolic manifolds; material surfaces that respectively undergo exponential stretching and contraction which are shown in figures 1(c), 3, 4, 5 and 7. If these manifolds intersect transversely, fluid elements undergo exponential stretching and folding, leading to chaotic advection (Ottino Reference Ottino1989). Indeed, several studies (Mourzenko et al. Reference Mourzenko, Yousefian, Kolbah, Thovert and Adler2002; Johnson, Brown & Stockman Reference Johnson, Brown and Stockman2006; Zhao et al. Reference Zhao, Li and Jiang2014) have examined steady 3-D Stokes flow in wide and rough fracture intersections, and found that although chaotic advection was not explicitly diagnosed, velocity fields were observed that exhibited 3-D structure consistent with non-trivial streamline braiding and thus chaotic advection (Boyland, Aref & Stremler Reference Boyland, Aref and Stremler2000). Conversely, tangential intersections between 2-D manifolds prohibit chaotic advection as the exponential expansion and contraction of hyperbolic manifolds exactly cancel. As fracture networks comprising narrow fractures impose strongly constrained geometries that favour near-tangential manifold intersections, the prevalence of chaotic advection under steady Stokes flow is unclear.

The common assumption of isotropic Darcy flow in 3-D DFNs further constrains the Lagrangian kinematics. Within fractures, steady 3-D Stokes flow may be validly approximated as a Hele-Shaw flow with variable aperture $h(\boldsymbol{x})\sim \ell$ across the fracture. Hence, for $\ell \ll L$ , the flow in fractures can be validly approximated as an isotropic heterogeneous Darcy flow (Bear Reference Bear1972)

(2.1) \begin{equation} \boldsymbol{v}(\boldsymbol{x})=-\frac {k(\boldsymbol{x})}{\mu }\boldsymbol{\nabla }p(\boldsymbol{x}), \end{equation}

where $\mu$ is the fluid viscosity and the permeability $k(\boldsymbol{x})=h(\boldsymbol{x})^2/12$ . For 2-D and 3-D DFNs, it is assumed that isotropic Darcy flow extends throughout the entire fracture network. This renders the velocity field helicity-free (Sposito Reference Sposito2001; Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2021), where from (2.1), the helicity density (Moffatt Reference Moffatt1969), defined as $\mathcal{H}(\boldsymbol{x}) \equiv \boldsymbol{v}\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{v})$ , vanishes identically as

(2.2) \begin{equation} \mathcal{H}(\boldsymbol{x})\equiv \boldsymbol{v}(\boldsymbol{x})\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{v}(\boldsymbol{x}))=\frac {1}{\mu ^2}k(\boldsymbol{x})\boldsymbol{\nabla }p(\boldsymbol{x})\boldsymbol{\cdot }(\boldsymbol{\nabla }k(\boldsymbol{x})\times \boldsymbol{\nabla }p(\boldsymbol{x}))=0. \end{equation}

Although anisotropic 3-D Darcy flow relaxes this constraint (Lester et al. Reference Lester, Metcalfe, Trefry and Dentz2025b ), such formulations have not, to our knowledge, been implemented in fracture network models. For isotropic Darcy flow, the condition $\mathcal{H}(\boldsymbol{x})=0$ implies that the flow is globally integrable (i.e. non-chaotic) in simply connected domains in the sense of Arnol’d (Reference Arnol’d1965). However, this integrability may not persist in multiply connected domains such as fracture networks.

Beyond fracture networks, the question of whether chaotic advection can arise under steady isotropic 3-D Darcy flow in topologically complex (multiply connected) domains is relevant to many applications, including transport processes in groundwater systems, resin transfer moulding and packed-bed reactors. Although Camacho & Neto (Reference Camacho and Neto1981) showed that closed (recirculating) streamlines in multiply connected domains can break global integrability, such streamlines cannot arise in Darcy flow (Bear Reference Bear1972). Nevertheless, multiply connected domains typically contain numerous stagnation points according to the Poincaré–Hopf theorem (Poincaré Reference Poincaré1881, Reference Poincaré1882; Hopf Reference Hopf1927). At such points, the Frobenius theorem (Frobenius Reference Frobenius1877) – which links zero helicity density to integrability – breaks down, meaning the flow is only locally integrable. Consequently, both steady 3-D Stokes flow and isotropic Darcy flow in multiply connected domains may admit chaotic advection through fluid deformation at stagnation points. However, the zero helicity density nature of isotropic Darcy flow means that streamline braiding (another potential mechanism for chaotic advection) can only occur in Stokes flow.

The degree of chaotic advection in fracture networks due to stagnation points may be estimated using a simple model for the dimensionless Lyapunov exponent $\lambda _\infty \equiv \hat {\lambda }_\infty \tau _v$ developed for random porous networks (Lester, Metcalfe & Trefry Reference Lester, Metcalfe and Trefry2013), where $\hat \lambda _\infty$ is the dimensional Lyapunov exponent and $\tau _v$ is the mean advection time between fracture intersections. In Appendix A.1 we apply this model to fracture networks with aspect ratio $\ell /L$ and obtain the estimate $\lambda _\infty \approx \ell /(8L) \ll 1$ , which is substantially smaller than those reported for random porous networks ( $\lambda _\infty \approx 0.1178$ ) or granular media ( $\lambda _\infty \approx 0.17$ – $0.21$ ) (Heyman et al. Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020; Souzy et al. Reference Souzy, Lhuissier, Méheust, Le Borgne and Metzger2020; Heyman, Lester & Le Borgne Reference Heyman, Lester and Le Borgne2021). In Appendix A.2 we show that under Stokes flow, vorticity is generated at no-slip boundaries which can lead to helicity generation in the bulk. This suggests that the strongly constrained geometry of fracture networks renders chaotic advection negligible, although further investigation is required to establish this conclusion more firmly.

2.2. Model scope and assumptions

Accordingly, we assume chaotic advection is negligible in large aspect ratio ( $\ell \ll L$ ) fracture networks. For simplicity of exposition, we also assume that intersections involve at most two fractures (precluding higher-order junctions) and that pressure gradients along intersections are negligible. To account for dead-end fractures which arise from pressure gradients along intersections (and play an important role in advective transport (Kang et al. Reference Kang, Hyman, Han and Dentz2020; Yoon et al. Reference Yoon, Hyman, Han and Kang2023)), we partition each intersection into ‘inflow’ and ‘outflow’ segments relative to the dead-end fracture. This yields a piecewise-constant approximation of the pressure gradient along intersections (§ 2.4).

Under these assumptions, advective mixing is both deterministic and reversible, and is governed entirely by streamline organisation within the network. To quantify this process, we consider both streamline routing (§§ 2.3, 2.4) and the steady transport of a continuously injected, non-diffusive tracer. As tracer concentration is conserved along streamlines, mixing is controlled solely by their transverse organisation and is independent of the Lagrangian travel time $\tau$ . Consequently, shear-driven mechanisms (e.g. Taylor–Aris dispersion) and purely advective longitudinal dispersion (Dentz & Hyman Reference Dentz and Hyman2023) are neglected. The present framework can, however, be extended (through characterisation of the Lagrangian travel time statistics (Davy et al. Reference Davy, Le Goc, Darcel, Pinier, Selroos and Le Borgne2024)) to incorporate time-dependent processes such as solute diffusion and reaction (§ 5).

Resolving streamline dynamics at both the aperture scale $\ell$ and fracture scale $L$ , we show that advective mixing is governed by streamline routing through fractures and their intersections (§ 2.4), leading to discontinuous mixing via cutting-and-shuffling (CS) of fluid elements (Sturman Reference Sturman2012; Park et al. Reference Park, Umbanhowar, Ottino and Lueptow2016; Smith et al. Reference Smith, Rudman, Lester and Metcalfe2017), rather than the stretching-and-folding (SF) mechanisms characteristic of chaotic advection in porous media (§ 2.6). We also assess the contribution of fluctuating fluid deformation superimposed on this process.

This discontinuous mixing is naturally represented by a graph-based formulation (§ 3), the mixing graph ${\mathcal{G}_M}$ , which extends the conventional DFN intersection graph ${\mathcal{G}_I}$ (Hyman et al. Reference Hyman, Hagberg, Srinivasan, Mohd-Yusof and Viswanathan2017). Comparison with fully resolved DFN simulations in § 4 demonstrates that ${\mathcal{G}_M}$ accurately captures advective mixing, supporting both the discontinuous mixing mechanism and its parsimonious graph representation. By linking mixing dynamics directly to network topology, ${\mathcal{G}_M}$ enables quantitative analysis of advective mixing in large random networks using both numerical and analytical approaches. This also provides a foundation for ab initio modelling of time-dependent transport processes, including solute diffusion, chemical reactions and colloidal transport.

2.3. Streamfunction coordinate system

To resolve advective mixing in 3-D DFNs as a deterministic and reversible process, we quantify streamline routing via a streamline coordinate system, the existence of which is guaranteed if $\mathcal{H}(\boldsymbol{x})=0$ , ensuring the existence of a holonomic (coherent) (Schutz Reference Schutz1980) coordinate frame. As the 3-D DFN velocity field $\boldsymbol{v}(\boldsymbol{x})$ is globally integrable, it may be expressed in terms of the dual streamfunctions (Lester et al. Reference Lester, Dentz, Singh and Bandopadhyay2023)

(2.3) \begin{equation} \boldsymbol{v}(\boldsymbol{x})=Q\boldsymbol{\nabla }\hat \psi _1(\boldsymbol{x})\times \boldsymbol{\nabla }\hat \psi _2(\boldsymbol{x}), \end{equation}

where $Q$ is the total volumetric flow rate through the fracture network, and the streamfunctions ( $\hat \psi _1,\hat \psi _2$ ) are orthogonal ( $\boldsymbol{\nabla }\hat \psi _1\boldsymbol{\cdot }\boldsymbol{\nabla }\hat \psi _2=0$ ) and normalised ( $\hat \psi _1\in [0,1],\hat \psi _2\in [0,1]$ ). As the dual streamfunctions are non-unique (Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2021), we define without loss of generality that the $\hat \psi _1$ and $\hat \psi _2$ streamfunctions respectively vary along the fracture width ( ${\sim} L$ ) and aperture ( ${\sim} \ell$ ).

The dual streamfunctions are analogous to the single streamfunction $\psi (x_1,x_2)$ for a 2-D incompressible flow, $\boldsymbol{v}=\boldsymbol{\nabla }\psi \times \hat {\boldsymbol{e}}_3$ , which may be extruded into a 3-D flow in the $x_3$ direction as $\boldsymbol{v}=\boldsymbol{\nabla }\psi _1\times \boldsymbol{\nabla }\psi _2$ , with $\psi _1=\psi$ and $\psi _2=x_3$ . Hence, the dual streamfunction representation provides a natural 3-D extension of 2-D incompressible flow (Bear Reference Bear1972). In both cases, the flow is globally integrable as there exist $d$ -1 (where $d$ is Eulerian dimension of the flow) invariants of the flow given by the streamfunctions (as $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\hat \psi _1=\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\hat \psi _2=0$ ).

As the streamfunctions $\hat \psi _1$ , $\hat \psi _2$ are invariant, streamlines are confined to coherent 2-D streamsurfaces $\hat \psi _1=$ const., $\hat \psi _2=$ const. (henceforth termed $\hat \psi _1$ - and $\hat \psi _2$ -streamsurfaces) and the flow is termed epi-2D (Yoshida & Morrison Reference Yoshida and Morrison2017). Hence, 1-D streamlines arise at the intersection of $\hat \psi _1$ - and $\hat \psi _2$ -streamsurfaces, and so each 1-D streamline may be uniquely labelled as the streamfunction pair $(\hat \psi _1,\hat \psi _2)$ .

From (2.1)–(2.3), the fluid pressure gradient $\boldsymbol{\nabla }p(\boldsymbol{x})$ is also orthogonal to the streamfunctions ( $\boldsymbol{\nabla }p\boldsymbol{\cdot }\boldsymbol{\nabla }\hat {\psi }_i=0$ , $i=1,2$ ), and so provides a structured, orthogonal 3-D Lagrangian coordinate system $\boldsymbol \zeta \equiv (\zeta _1,\zeta _2,\zeta _3)\equiv (\hat \psi _1,\hat \psi _2,p)$ . The transform from Eulerian space $\boldsymbol{x}=(x_1,x_2,x_3)$ to $\boldsymbol \zeta$ completely quantifies advective mixing. We note that the coordinate system $\boldsymbol \zeta$ naturally resolves advective mixing in rough fractures and intersections of variable aperture fractures, as the 2-D boundaries of 3-D fractures correspond to coherent $\hat \psi _1$ - or $\hat \psi _2$ -streamsurfaces. Streamline coordinates will be used extensively to quantify advective mixing in 3-D DFNs.

2.4. Intersection and fracture mixing

We consider simple 3-D DFNs (figure 2) as a basis to develop methods that extend to arbitrarily large networks. Despite their simplicity, these systems exhibit two distinct mixing mechanisms controlled by fracture geometry and connectivity. In the top configuration, streamline routing within fractures produces large-scale mixing at the fracture length scale $L$ , whereas in the bottom configuration, routing at intersections induces fine-scale mixing at the aperture scale $\ell$ . These define fracture mixing (at length scale $L$ ) and intersection mixing (at length scale $\ell$ ).

We denote the 2-D cross-section immediately adjacent to an intersection as a face $f_n$ ( $n=1{:}N_f$ ). Under the assumption of negligible longitudinal pressure gradients along intersections, the velocity field $\boldsymbol{v}(\boldsymbol{x})$ is orthogonal to all faces $f_n$ . Hence, each face is classified as an inlet $f_n^i$ or outlet $f_n^o$ relative to its adjacent intersection; by convention, we define an inflow face (i.e. that with fluid flowing into the network) to be an outlet face $f_n^o$ and likewise an outflow face (i.e. that with fluid flowing out of the network) to be an inlet face $f_n^i$ . Streamlines thus traverse the network periodically encountering outlet faces, fractures, inlet faces and intersections in sequence.

Figure 3.

(a) Schematic of two intersecting fractures $F_1$ , $F_2$ of characteristic aperture $\ell$ and width $L$ forming an intersection $I_1$ . Streamlines are separated by a continuous line of hyperbolic stagnation points $\boldsymbol{x}_p$ at the fracture intersection, and the associated stable (blue) and unstable (red) hyperbolic manifolds route surrounding streamlines (black) to their respective outlet faces. (b) Schematic detailing intersection $I_1$ shown in panel (a) depicting the inlet faces $f^i_1$ , $f^i_2$ associated with fracture $F_1$ and outlet faces $f^o_1$ , $f^o_2$ associated with fracture $F_2$ . Also shown are the local streamfunctions $\psi _{1,n}\in [0,1]$ , $\psi _{2,n}\in [0,1]$ , and associated $\psi _{1,n}$ - and $\psi _{2,n}$ -streamsurfaces that are respectively oriented along the fracture width $L$ and fracture aperture $\ell$ of each face. Similar to panel (a), the hyperbolic point $\boldsymbol{x}_p$ and associated manifolds route streamlines into downstream outlet faces. Local streamfunctions $\psi _2$ govern routing of streamlines, which leads to cutting and shuffling of fluid elements leaving $f^i_1$ , $f^i_2$ and arriving at $f^o_1$ , $f^o_2$ .

Figure 4.

Schematic of the complete set of distinct streamline topologies in a cross-section (a–d) X-, (e–f) T- and (g) L-shaped intersection. Each cross-section corresponds to a $\hat \psi _1$ -streamsurface, hence, the streamline topology is the same as that for a steady 2-D flow. Stagnation points $\boldsymbol{x}_p$ and corresponding (blue) stable and (red) unstable hyperbolic manifolds (which act as separating streamlines) act to route fluid parcels through inlet $f_n^i$ and outlet $f_n^o$ faces of the intersection. Black arrows indicate flow direction.

As shown in figure 3, intersection mixing is governed by streamline routing between inlet and outlet faces. Owing to the epi-2-D nature of isotropic Darcy flow, streamlines lie on $\hat \psi _1$ -streamsurfaces and exhibit the topology of steady 2-D planar flow. Within each surface, routing is determined by separatrices connected to the stagnation point $\boldsymbol{x}_p$ . This topology is invariant across $\hat \psi _1$ -streamsurfaces, such that zero-dimensional (0-D) stagnation points and 1-D separatrices correspond to 1-D stagnation lines and 2-D separating $\hat \psi _1$ -streamsurfaces spanning the 3-D intersection. Intersections of two fractures may also form X-, T- or L-shaped cross-sections, with four, three or two faces, respectively (figure 4). These admit seven distinct planar streamline topologies, distinguished by inlet–outlet configuration and stagnation structure. As only $\hat \psi _2$ varies within each $\hat \psi _1$ -streamsurface, it fully characterises streamline routing during intersection mixing.

Similar to intersection mixing, fracture mixing is driven by the routing of streamlines within fractures to different faces $f_n$ . As streamlines are confined to $\hat \psi _2$ -streamsurfaces (figure 5 a), streamline topology within each fracture is the same as that of a 2-D planar flow. Within each $\hat \psi _2$ -streamsurface, streamlines are routed by separating streamlines connected to stagnation points and other critical points (defined in §§ 2.5). Similar to intersection mixing, these 0-D critical points and 1-D separating streamlines manifest as 1-D critical lines and 2-D separating $\hat \psi _1$ -streamsurfaces across the full 3-D fracture.

As shown in figure 5(a), even simple fractures (e.g. $F_2$ ) connecting a single source ( $I_1$ ) and sink ( $I_2$ ) contain multiple separatrices that partition flow between faces, yielding more complex topologies than for intersection mixing (figure 4). Figure 5(b) illustrates the treatment of dead-end fractures under the assumption of constant pressure along intersections. By splitting the intersection $I_1$ in figure 5(a) into ‘inflow’ and ‘outflow’ components (with respect to the dead-end fracture $F_1$ ), streamlines are routed from $I_1$ to $I_3$ via the same mechanism as in non-dead-end fractures. Figure 6(a) shows typical streamlines (computed using DFN.lab) for an X-shaped intersection. The flux-weighted random streamline routing in DFN.lab distributes red and blue streamlines across both outlets, whereas deterministic routing (figure 4 b) directs all streamlines of a given type into at least one outlet. This discrepancy highlights the limitations of random routing for representing advective – and thus also finite- $Pe$ – mixing. Figure 6(b) shows the corresponding routing in a dead-end fracture, where $I_3$ is split into $I_3$ and $I_4$ ; the resulting intersection graphs are shown in figure 6(c), with intersections as nodes and fractures as edges.

Figure 5.

(a) Schematic of a simple fracture network (at the fracture scale $L$ ) comprising fractures $F_n$ (where $F_1$ and $F_3$ respectively are inlet and outlet fractures) and intersections $I_m$ (green lines), separating streamlines (black lines) and critical points (black dots). Separating streamlines route fluid ‘parcels’ from ‘source’ intersection faces to different ‘sink’ intersection faces. Stagnation (critical) points are denoted $\boldsymbol{x}_p$ (with corresponding hyperbolic stable (blue) and unstable (red) manifolds), while terminal and interior critical points are unlabelled. All separating streamlines either originate or terminate at stagnation or terminal critical points. (b) Same schematic as panel (a), but with a dead-end fracture ( $F_1$ ) which is treated by splitting intersection $I_1$ in panel (a) into an ‘outflow’ intersection $I_1$ and an ‘inflow’ intersection $I_3$ , with streamlines routed in $F_1$ from $I_1$ to $I_3$ .

Figure 6.

(a) Random, flux-weighted routing of streamlines through an X-shaped intersection, computed using DFN.lab. Note that random routing directs both red and blue inlet streamlines into both outlet fractures, which cannot occur for deterministic streamlines. (b) Random, flux-weighted routing of streamlines through a dead-end fracture $F_3$ , computed using DFN.lab. Note that intersection $I_3$ in panel (a) can be split into two separate fractures $I_3$ , $I_4$ . (c) Intersection graphs ${\mathcal{G}_I}$ depicting (top) the X-shaped intersection shown in panel (a) and (bottom) the X-shaped intersection with dead-end shown in panel (b).

Despite the differences in complexity of streamline topology and the presence of dead-end fractures, fracture mixing is very similar to intersection mixing in that it is driven by streamline routing. Hence, critical points, lines and surfaces control both intersection and fracture mixing in 3-D DFNs.

2.5. Critical points, lines and surfaces

The inherent topological complexity of fracture networks gives rise to critical features – stagnation points and lines, separating streamlines and streamsurfaces (figures 3, 4, 5) – that control streamline routing and advective mixing. In 3-D DFNs, stagnation points $\boldsymbol{x}_p$ arise at junctions within fractures, which are connected to separating streamlines that route fluid parcels to/from different intersections within the fracture. In addition to stagnation points, terminal critical points shown (figure 5 a) arise at the termini of 1-D intersections and differentiate flows into or out of different faces $f_n$ on each side of an intersection. Hence, in fracture flows, there are additional separating streamlines that originate or terminate at terminal critical points, and these separating streamlines can either originate or terminate at stagnation points, terminal points or midway along faces $f_n$ , forming the interior critical points shown in figure 5(a). Together, these critical points and their associated separating streamlines act to route fluid elements from parts of outlet faces $f_n^o$ to parts of inlet faces $f_n^i$ , leading to discontinuous mixing as discussed in the following.

2.6. Discontinuous mixing in fractured media

Although stretching and folding (SF) of fluid elements around stagnation points in 3-D DFNs cannot generate significant chaotic advection, it does generate discontinuous mixing (Sturman Reference Sturman2012; Smith et al. Reference Smith, Umbanhowar, Lueptow and Ottino2019), involving cutting-and-shuffling (CS) of fluid elements. CS is akin to shuffling of a deck of infinitely divisible cards, where an initially continuous set of streamlines is ‘cut’ into a set of streamline ‘parcels’ that are then ‘shuffled’ and recombined in a different arrangement. In the language of ergodic theory (Arnold & Avez Reference Arnold and Avez1968), mixing via CS is weak mixing, as compared with strong mixing under SF; however, CS mixing can still lead to rapid and complete (ergodic) advective mixing (Sturman Reference Sturman2012; Smith et al. Reference Smith, Umbanhowar, Lueptow and Ottino2019).

Figure 7.

(a) Stretching and folding (SF) due to hold-up of a fluid element (shaded, shown at successive times) at the stagnation point $\boldsymbol{x}_p$ of a fracture intersection. As the ‘legs’ of the fluid element are advected downstream into separate outlets, the element is stretched exponentially in time at $\boldsymbol{x}_p$ . Further downstream, the legs of the element are effectively disconnected (‘cut’) when their intersection with an arbitrary planar cross-section (vertical dashed line) $P_o$ is considered. (b) Merger of two fluid elements (dark and light shading) as they approach a stagnation point $\boldsymbol{x}_p$ of a fracture intersection. Within the planar cross-section $P_o$ , these elements are essentially merged next to each other (as the width of the separating fluid decays exponentially at $\boldsymbol{x}_p$ ) effectively ‘shuffling’ fluid elements together. In unison, the stretching and folding (SF) motions shown in panels (a) and (b) manifest as cutting and shuffling (CS) actions when considered within the cross-section $P_o$ .

Figure 7(a) illustrates how stretching and folding (SF) near a stagnation point $\boldsymbol{x}_p$ manifests as cutting and shuffling (CS): part of a fluid element is held up at $\boldsymbol{x}_p$ , stretched and folded across the intersection, and then ‘cut’ downstream as its legs disconnect in the cross-section $P_0$ . Conversely, two fluid elements passing near $\boldsymbol{x}_p$ (figure 7 b) are effectively merged, ‘shuffling’ them together. Thus, the local SF dynamics project as CS in $P_0$ . The fine-scale interleaving in figure 2(c) exemplifies CS at intersections, while hold-up at fracture critical points produces large-scale CS (figure 2 a). A similar mechanism occurs in open porous networks, though these admit both SF and CS (Lester et al. Reference Lester, Heyman, Méheust and Le Borgne2025a ).

Although persistent exponential stretching is absent in 3-D DFNs, fluctuating stretching arises from local velocity variations. CS without deformation corresponds to piecewise isometries (PWIs), while CS with stretching yields piecewise smooth (PWS) transforms (Kreczak, Sturman & Wilson Reference Kreczak, Sturman and Wilson2017). In DFNs, this stretching is sub-exponential and bounded, reflecting non-chaotic flow. Sections 3 and 4 show that mixing is exactly described by PWS transforms which are quantified via a graph-based representation of advective mixing that is developed in the following section.

Figure 8.

(a) Intersection of local streamfunctions ( $\psi _{1,n},\psi _{2,n}$ ) within the face $f_n$ of the fracture network (with dimensions ${\sim} L\,\times {\sim} \ell$ ) that define a streamline with local velocity $\boldsymbol{v}_n$ . (b) Distribution of the local $\psi _{1,n}$ (vertical lines) and $\psi _{2,n}$ (horizontal lines) streamsurfaces over the face $f_n$ of a rough fracture of variable aperture. The mapping from local spatial coordinates $(x\prime_{1,n},x\prime_{2,n})$ in panel (b) to local streamfunction coordinates $(\psi _{1,n},\psi _{2,n})$ in panel (c) generates a structured orthogonal coordinate system over face $f_n$ . Note this structured grid is also depicted throughout the intersection shown in figure 3(b).

3.1. Global and local streamfunctions

Under this framework, the invariant streamfunctions $(\hat \psi _1$ , $\hat \psi _2)\in [0,1]\times [0,1]$ in (2.3) are termed global streamfunctions as they uniquely label each streamline over the entire fracture network. In addition to these global streamfunctions, we also define local streamfunctions $(\psi _{1,n},\psi _{2,n})$ (shown in figure 8) over each face $f_n$ of the network (e.g. figures 2 b and 2 d). These local streamfunctions characterise the local distribution of streamlines over $f_n$ and so facilitate quantification of advective mixing.

The local streamfunctions are also normal $(\psi _{1,n},\psi _{2,n})\in [0,1]\times [0,1]$ ) and orthogonal ( $\boldsymbol{\nabla }\psi _{1,n}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{2,n}=0$ ), and are quantified in terms of the local velocity field $\boldsymbol{v}_n$ as

(3.1) \begin{equation} \boldsymbol{v}_n(\boldsymbol{x})\equiv \boldsymbol{v}(\boldsymbol{x})\,\,\forall \,\,\boldsymbol{x}\in S_n, \end{equation}

where $S_n$ is the 2-D surface that spans the face $f_n$ . Consistent with the dual streamfunction representation (2.3), the local streamfunctions $(\psi _{1,n},\psi _{2,n})$ are then defined as

(3.2) \begin{equation} \boldsymbol{v}_n(\boldsymbol{x})\equiv Q_n\boldsymbol{\nabla }\psi _{1,n}(\boldsymbol{x})\times \boldsymbol{\nabla }\psi _{2,n}(\boldsymbol{x}),\quad n=1:N_{I},\quad \psi _{i,n}\in [0,1],\quad i=1,2, \end{equation}

where $Q_n$ is the volumetric flow rate through face $f_n$ . From (3.2), the face $f_n$ is bound by four stream surfaces: $\psi _{1,n}=0$ , $\psi _{1,n}=1$ , $\psi _{2,n}=0$ , $\psi _{2,n}=1$ . As $\boldsymbol{\nabla }\psi _{1,n}$ and $\boldsymbol{\nabla }\psi _{2,n}$ are both parallel to $S_n$ , the fluid flow rate $Q_n$ through this face is then conserved (Bear Reference Bear1972)

(3.3) \begin{equation} Q_n=\iint _{S_n}\boldsymbol{v}_n(\boldsymbol{x})\boldsymbol{\cdot }\text{d}\boldsymbol{S}=Q_n\int _0^1\int _0^1\text{d}\psi _{1,n}\,\text{d}\psi _{2,n}=Q_n. \end{equation}

The local streamfunctions $(\psi _{1,n},\psi _{2,n})$ represent a local re-labelling of the global streamfunctions over face $f_n$ . For example, if the $\hat \psi _1$ -streamsurface $\hat \psi _1=a=$ const. intersects face $f_n$ , this streamsurface also satisfies $\psi _{1,n}=b=$ const., but in general $a\neq b$ . Hence, a streamline may be labelled by the unique pair of global streamfunctions $(\hat \psi _{1},\hat \psi _{2})$ that are invariant over the entire network, but this streamline also traverses the fracture network via a sequence of faces $f_{n_1}, f_{n_2},\ldots$ with local streamfunction values $(\psi _{1,n_1},\psi _{2,n_1}),(\psi _{1,n_2},\psi _{2,n_2}),\ldots$ . As such, the global streamfunctions $(\hat \psi _1,\hat \psi _2)$ are used to label a given streamline, while the local streamfunctions $(\psi _{1,n},\psi _{2,n})$ are used to characterise the local position of that streamline in face $f_n$ . Hence, advective mixing can be quantified by the evolution of $(f_n,\psi _{1,n},\psi _{2,n})$ for a given streamline.

As illustrated in figure 8, $\psi _{1,n}$ varies along the fracture width ( ${\sim} L$ ) of face $f_n$ and $\psi _{2,n}$ varies along the fracture aperture ( ${\sim} \ell$ ), with streamlines arising at the intersection of local $\psi _{1,n}$ - and $\psi _{2,n}$ -streamsurfaces. We also introduce the local Eulerian spatial coordinates $(x^\prime_{1,n},x^\prime_{2,n})\in [0,{\sim} L]\times [0,{\sim} \ell ]$ (figure 8 b) that span $f_n$ and their dimensionless counterparts $(x_{1,n},x_{2,n})\in [0,1]\times [0,1]$ . Figure 8(c) shows that the mapping from local spatial coordinates $(x^\prime_{1,n},x^\prime_{2,n})$ to local streamfunction coordinates $(\psi _{1,n},\psi _{2,n})$ allows rough fractures with variable aperture to be mapped onto a structured orthogonal grid. As shown in figure 3(b), this means that even if rough, 3-D intersections are spanned and bounded by a structured orthogonal grid of $\hat \psi _1$ - and $\hat \psi _2$ -streamsurfaces, where the $\hat \psi _1$ -streamsurfaces are coherent 2-D cross-sections of the intersection.

For fracture networks with only one inlet fracture (i.e. $N_{I_i}=1$ ), the outlet face (where fluid enters the network) is defined as $f_n=f^o_1$ with flow rate $Q_1^o$ . Hence, the total network flow rate is $Q=Q^o_1$ , and from (2.3), (3.2), we can also define the global and local streamfunctions to coincide over face $f_n$ as

(3.4) \begin{align} \psi _{1,n=1}\equiv \hat \psi _1 && \psi _{2,n=1}\equiv \hat \psi _2. \end{align}

Conversely, for fracture networks with multiple inlet fractures ( $N_{I_i}\gt 1$ ), the corresponding outlet faces can be defined as $f_n=f^o_n$ with $n=1:I_i$ , and the total network flow rate $Q$ is

(3.5) \begin{equation} Q=\sum _{n=1}^{N_{I_i}}Q^o_{n}, \end{equation}

where $Q_{n,1}^o$ is the total flowrate through face $f_n^o$ . From (2.3), (3.2), we can uniquely define the global streamfunctions by spanning the $\hat \psi _1$ streamfunction over the $N_{I_i}$ ordered inlets as

(3.6) \begin{align} \psi _{1,m}\equiv \frac {1}{Q^o_{m,1}}\left (Q\,\hat \psi _1-\sum _{n=1}^{m-1}Q^o_{n,1}\right ) && \psi _{2,m}\equiv \hat \psi _2,\quad m\in [1:N_{I_i}], \end{align}

where $m$ is the index of the outlet face $f^o_m$ corresponding to $\psi _{1,m}$ , $\hat \psi _1$ . Note the flow rates in (3.6) arise from normalisation of the streamfunctions in (2.3), (3.2). The local streamfunctions $(\psi _{1,n},\psi _{2,n})$ for the remaining faces (i.e. $n\gt N_{I_i}$ ) in the network are defined by (3.2). Together, the local and global streamfunctions can be used to quantify intersection and fracture mixing.

3.2. Intersection and fracture mixing maps

We first develop a method to route streamlines between faces in the network. As a streamline is advected from face $f_j$ (with local streamfunctions $\psi _{1,j},\psi _{2,j}$ ) in the network, the next face $f_k$ it encounters (and local streamfunctions $\psi _{1,k},\psi _{2,k}$ ) is determined by streamline routing. The relationship between these faces and local streamfunctions can be efficiently encoded via two classes of algebraic maps we term mixing maps. If face $f_j$ is an inlet face, the streamline then passes through an intersection and $(f_k,\psi _{1,k},\psi _{2,k})$ is determined by an intersection mixing map ${\mathcal{M}_I}$ . Conversely, if face $f_j$ is an outlet face, the streamline then passes through a fracture and $(f_k,\psi _{1,k},\psi _{2,k})$ is determined by a fracture mixing map ${\mathcal{M}_F}$ . As streamlines are routed through either intersections or fractures with fixed global streamfunctions $(\hat \psi _1,\hat \psi _2)$ , the relationship between the local streamfunctions at faces $f_j$ and $f_k$ is given by the mappings between the global and local streamfunctions at these faces.

For the case of intersection mixing, we show in Appendix B that if the spatial distribution of areal fluxes,

(3.7) \begin{equation} q_n(x_1)\equiv \int _0^1 ||\boldsymbol{v}_n|| \,\text{d}\psi _{2,n}=Q_n\boldsymbol{\nabla }\psi _{1,n}, \end{equation}

between the inlet face $f_j$ and outlet face $f_k$ of an intersection are proportional, i.e. $q_j(x_1)/q_k(x_1)=$ const., then the local $\psi _1$ streamfunctions are equivalent, i.e. $\psi _{1,k}=\psi _{1,j}$ . Hence, streamline routing only alters the $\psi _2$ coordinate. This condition renders intersection mixing a PWI transform, where fluid elements are cut and shuffled with respect to $\psi _2$ without stretching in the $\psi _1$ direction. As such, the local streamfunctions evolve according to

(3.8) \begin{align} (\psi _{2,k},f_k)={\mathcal{M}_I}(\psi _{2,j},f_j) && \psi _{1,k}=\psi _{1,j}. \end{align}

If the flux condition is violated, i.e. $q_j(x_1)/q_k(x_1)\neq$ const., then intersection mixing is a PWS transform, where fluid stretching along the $\psi _1$ direction can occur in conjunction with cutting and shuffling with respect to $\psi _2$ . In this case, the intersection mixing map ${\mathcal{M}_I}$ acts on both local streamfunctions as

(3.9) \begin{equation} (\psi _{1,k},\psi _{2,k},f_k)={\mathcal{M}_I}(\psi _{1,j},\psi _{2,j},f_j). \end{equation}

We refer to (3.8) and (3.9) respectively as the PWI and PWS representations of ${\mathcal{M}_I}$ , where the PWI map is significantly simpler in that it can be completely quantified solely in terms of the total fluxes $Q_n$ through faces, whereas the PWS map requires the areal fluxes $q_n(x_1)$ as a function of $x_1$ along each face of the intersection (see Appendix B for details). Throughout §§ 2 and 3, we use ${\mathcal{M}_I}$ arbitrarily to denote either the PWI or PWS maps, and in § 4, we will directly test each representation by comparison against numerical simulation.

During fracture mixing, streamlines flow from the outlet face $f_k$ to the inlet face $f_l$ . The large aspect ratio $\ell \ll L$ renders the variation between $\psi _{2,j}$ and $\psi _{2,k}$ negligible ( $\mathcal{O}(\ell /L)$ ), and so only the local $\psi _1$ streamfunction is significantly altered

(3.10) \begin{align} (\psi _{1,l},f_l)&={\mathcal{M}_F}(\psi _{1,k},f_k) && \psi _{2,l}=\psi _{2,k}. \end{align}

For a given streamline, the overall map $\mathcal{M}$ which quantifies the evolution of local streamfunctions through the fracture network is given by a sequential composition of ${\mathcal{M}_F}$ , ${\mathcal{M}_I}$ as

(3.11) \begin{equation} \mathcal{M}={\mathcal{M}_F}\circ {\mathcal{M}_I}\circ \ldots {\mathcal{M}_F}\circ {\mathcal{M}_I}, \end{equation}

reflecting that streamlines sequentially travel through fractures and intersections. Hence, the maps ${\mathcal{M}_I}$ , ${\mathcal{M}_F}$ encode properties of the fracture network that govern the streamline routing process. Note that ${\mathcal{M}_F}$ , ${\mathcal{M}_I}$ are local in that they are specific to each face $f_n$ . From (3.11), the routing of a given streamline is characterised by the sequence of local streamfunctions and faces $f_{n_j}$ for $j=1,\ldots ,N$ it traverses through the fracture network

(3.12) \begin{equation} \{(\psi _{1,f_{n_1}},\psi _{2,f_{n_1}},f_{n_1}),(\psi _{1,f_{n_2}},\psi _{2,f_{n_2}},f_{n_2}),\ldots ,(\psi _{1,f_{n_N}},\psi _{2,f_{n_N}},f_{n_N})\}, \end{equation}

where

(3.13) \begin{equation} (\psi _{1,f_{n_i}},\psi _{2,f_{n_i}},f_{n_i})={\mathcal{M}_F}\circ {\mathcal{M}_I}\circ \ldots {\mathcal{M}_F}\circ {\mathcal{M}_I} (\psi _{1,f_{n_1}},\psi _{2,f_{n_1}},f_{n_1}). \end{equation}

This characterisation of local streamfunctions $(\psi _{1,n},\psi _{2,n})$ and faces $f_n$ completely quantifies advective mixing within the fracture network. Algebraic expressions for ${\mathcal{M}_I}$ , ${\mathcal{M}_F}$ based on respective flux balances within intersections and fractures are derived in Appendix B. In the following subsection, they are embedded in an efficient graph representation of mixing.

3.3. Intersection graph ${\mathcal{G}_I}$

Previous studies have used graph representations to quantify flow and transport in DFNs, via both fracture graphs ${\mathcal{G}_F}$ (where vertices represent fractures and edges represent intersections) and intersection graphs (vice versa) ${\mathcal{G}_I}$ . To accommodate the mixing mechanisms identified in § 2, we extend the intersection graph ${\mathcal{G}_I}$ to form the mixing graph ${\mathcal{G}_M}$ . The intersection graph ${\mathcal{G}_I}$ is defined by the set of $N_I$ vertices given by the intersections $I_m$ and the set of $N_e$ graph edges $e_n$ (with weights $Q_n$ ) given by flows between intersections within each fracture; hence, ${\mathcal{G}_I}$ is a weighted acyclic directed graph (WADG). Here, we use modified nomenclature where the $n$ th inlet (outlet) face of intersection $I_m$ is denoted $f_{m,n}^i$ , $n=1:N_{I_i}$ with flow rate $Q_{m,n}^i$ ( $f_{m,n}^o$ , $n=1:N_{I_o}$ with flow rate $Q_{m,n}^o$ ), where $N_m^i\in [1,3]$ ( $N_m^o\in [1,3]$ ) is the total number of inlet (outlet) faces of intersection $I_m$ . Using this nomenclature, all interior intersections $I_m$ satisfy flow conservation

(3.14) \begin{equation} Q_m=\sum _{n=1}^{N_m^i}Q^i_{m,n}=\sum _{n=1}^{N_m^o}Q^o_{m,n}\quad \text{for}\quad m=N_{I_i}+1:N_I-N_{I_o}. \end{equation}

Note that each inlet into the overall fracture network is considered to be an intersection $I_m$ with one outlet face ( $N_m^o=1$ ) and each outlet is considered to be an intersection $I_m$ with only one inlet face ( $N_m^i=1$ ). From (3.5), the total flow rate through the network is

(3.15) \begin{equation} Q= \sum _{m=1}^{N_{I_i}}Q^o_{m,1}=\sum _{m=N_I-N_{I_o}}^{N_I}Q^i_{m,1}. \end{equation}

These properties mean that the intersection graph ${\mathcal{G}_I}$ may also be classified as a zero-capacity flow network (Ahuja, Magnanti & Orlin Reference Ahuja, Magnanti and Orlin1993). Figure 9(a) shows the intersection graph ${\mathcal{G}_I}$ for the simple DFN shown in figure 1(b), which comprises $N_F=5$ fractures and $N_I=8$ intersections, which comprise $N_{I,i}=2$ inlet intersections, $N_{I_b}=4$ interior intersections and $N_{I_o}=2$ outlet intersections.

Figure 9.

(a) Intersection graph ${\mathcal{G}_I}$ of a simple DFN depicting intersections $I_m$ (white vertices) and connecting flows as directed edges $e_n$ (red arrows). (b) Expansion of intersections $I_m$ in panel (a) into intersection clusters $C_m$ comprising inlet faces $f^i_{m,n}$ (grey vertices) and outlet faces $f^o_{m,n}$ (white vertices), and connecting flows as directed edges (blue arrows). (c) Mixing graph ${\mathcal{G}_M}$ formed by replacing the $N_I$ intersection vertices $I_m$ of ${\mathcal{G}_I}$ shown in panel (a) with the clusters $C_m$ shown in panel (b), and connecting the inlet faces $f^i_{m,n}$ (grey vertices) and outlet faces $f^o_{m,n}$ (white vertices) via streamline connections within fractures. Red arrows depict connection of faces via intersections and streamline routing via the intersection map ${\mathcal{M}_I}$ and blue arrows depict connection of faces via fractures and streamline routing via the intersection map ${\mathcal{M}_F}$ .

3.4. Mixing graph ${\mathcal{G}_M}$

As the intersection graph ${\mathcal{G}_I}$ does not resolve the intersection faces $f_n$ , it must be extended to form the mixing graph ${\mathcal{G}_M}$ . Each intersection $I_m$ is accompanied by $N^i_m$ inlet faces $f^i_{m,n}$ ordered in the direction of increasing $\psi _{2,n}$ as $n=1:N^i_m$ , and $N^o_m$ outlet faces $f^o_{m,n}$ ordered in the same manner as $n=1:N^o_m$ . Hence, an intersection cluster $C_m$ can be formed for each intersection $I_m$ (shown in figure 9 b) that comprises the inlet faces $f^i_{m,n}$ and outlet faces $f^o_{m,n}$ connected by $I_m$ . Then, ${\mathcal{G}_M}$ is formed by replacing each of the intersection vertices $I_m$ in ${\mathcal{G}_I}$ with its intersection cluster $C_m$ . Note that the inlet ( $I_1$ , $I_2$ ) and outlet intersections ( $I_7$ , $I_8$ ) in figure 9(b) only contain a single outlet face $f^o_{m,n}$ or inlet face $f^i_{m,n}$ , and so their intersection clusters $C_m$ comprise a single vertex. The directed graph edges between inlet $f^i_{m,n}$ and outlet $f^o_{m,n}$ faces (blue arrows in figure 9) are termed intersection edges as they are associated with flow through intersections, and the directed edges between outlet $f^o_{m,n}$ and inlet $f^i_{m,n}$ faces (red arrows in figure 9) are termed fracture edges and are associated with flow through fractures. Similar to ${\mathcal{G}_I}$ , ${\mathcal{G}_M}$ is both a WADG and a flow network.

As shown in figure 9(c), streamlines are advected through the graph vertices via the periodic sequence (outlet face, inlet face, outlet face, …). The intersection maps ${\mathcal{M}_I}$ and fracture maps ${\mathcal{M}_F}$ respectively apply at the intersection and fracture edges of the graph, forming a dynamic process known as a sequential dynamical system (SDS) (Mortveit & Reidys Reference Mortveit and Reidys2007), a class of discrete dynamical system embedded in graph structures such as ${\mathcal{G}_M}$ . Here, the topology of the SDS (fracture network) is encoded by ${\mathcal{G}_M}$ , and the dynamics (streamline routing) are encoded by the maps ${\mathcal{M}_{I}}$ , ${\mathcal{M}_F}$ respectively acting at intersection edges and fracture edges. The evolution of a scalar concentration field under the action of ${\mathcal{G}_M}$ is shown in figure 10, which highlights the CS actions that control mixing.

3.5. Quantification of ${\mathcal{G}_M}$ and computational overhead

The information required to construct the mixing graph ${\mathcal{G}_M}$ under the PWI or PWS models is summarised as follows. First, the fractures and intersections of the network must be identified as well as the intersection graph ${\mathcal{G}_I}$ (figure 9 a), including identification of dead-end fractures and splitting intersections into ‘inflow only’ or ‘outflow only’ regions with respect to the dead-end fracture (see § 2.4 for details). For each intersection, the relevant faces $f_n$ must be identified along with their volumetric flow rates $Q_n$ and the relative ordering of these faces around the intersection as detailed in Appendix B. This information is sufficient to quantify the intersection mixing map ${\mathcal{M}_I}$ under the PWI model. For the PWS model, the areal fluxes $q_n(x_1)$ along the faces are also required to quantify ${\mathcal{M}_I}$ . For the fracture mixing map ${\mathcal{M}_F}$ , the separating streamlines in each fracture must be identified, along with the streamline topology (figure 1 c) that routes streamlines from face to face within the fracture and the associated fluxes $Q_{j,n}$ , as detailed in Appendix B. This topological information allows construction of the clusters ${\mathcal{C}_m}$ associated with each intersection (figure 9 b), facilitating construction of the mixing graph ${\mathcal{G}_M}$ (figure 9 b) and quantification of advective mixing via the maps ${\mathcal{M}_I}$ , ${\mathcal{M}_F}$ .

Hence, the information required to quantify the PWI (PWS) model can be readily determined from DFN models by identifying all faces $f_n$ and computing the flow rates $Q_n$ ( $q_n(x_1)$ ), and computing the separating streamlines within the fractures and the associated fluxes $Q_{j,n}$ between these streamlines. As shown in figure 1(c), these separating streamlines can be identified and constructed by integrating tracer streamlines forward and backward in time from stagnation points and the end points of intersections. Hence, the mixing graph ${\mathcal{G}_M}$ can be constructed from minimal information from the DFN, and resolves the complex advective mixing dynamics with minimal computational overhead. This information (distribution of fluxes, streamline topology, etc.) can also be characterised in a statistical sense, facilitating development of stochastic models of advective mixing in large random DFNs.

Streamline routing fracture networks is computed efficiently via ${\mathcal{G}_M}$ . For a given streamline, we denote $n_F$ , $n_I=n_F+1$ respectively as the number of fractures and intersections routed from network inlet to outlet. From the operations in Appendix B, the number of floating point operations (FLOPS) for the fracture map ${\mathcal{M}_F}$ is $n_o+n_i+2$ . For the intersection map ${\mathcal{M}_I}$ , the PWI model also incurs $n_o+n_i+2$ FLOPS, whereas the PWS model incurs $n_o+n_i+4$ function valuations ( $q_n(x_1)$ , $\varPsi _{1,n}(x_1)$ , $X_{1,n}(\psi _1)$ ) of $p^2$ FLOPS each if they involve order $p$ spline functions. Thus this advection of $N$ streamlines requires $N(\langle n_F\rangle +p^2\langle n_I\rangle )(\langle n_o\rangle +\langle n_i\rangle +1)$ FLOPS, where $\langle n\rangle$ denotes the average of $n$ and $p=1$ for the PWI model.

3.6. Scalar transport and mixing measures

To quantify advective mixing in 3-D DFNs, we consider mixing of a non-diffusive scalar with a steady concentration distribution field $c$ that is defined in terms of the inlet concentration distribution (i.e. over the network inlets) as $c_0(\hat \psi _1, \hat \psi _2)$ . For simplicity of exposition, we consider the inlet distribution $c_0(\hat \psi _1, \hat \psi _2)=H(1/2-\hat \psi _1)$ (where $H$ is the Heaviside step function), but note that the complete mixing properties of the fracture network depends upon the mixing characteristics across all possible $c_0$ . The local solute concentration distribution $c_n(\psi _{1,n},\psi _{2,n})$ at each inlet face $f_n$ , $n=1:N_{I_i}$ of the fracture network may be given in terms of the local streamfunctions as

(3.16) \begin{equation} c_n(\psi _{1,n},\psi _{2,n})\equiv c_0 \left[\hat \psi _1(\psi _{1,n},n),\hat \psi _2(\psi _{2,n},n)\right],\quad n=1:N_{I_i}. \end{equation}

As the scalar concentration is invariant along streamlines, then for all faces,

(3.17) \begin{equation} c_n(\psi _{1,n},\psi _{2,n})=c_0\left [ \hat \psi _1(\psi _{1,n},\psi _{2,n},n), \hat \psi _2(\psi _{1,n},\psi _{2,n},n)\right ],\quad n=1:N_I. \end{equation}

Figure 10.

Schematic of the evolution of the 2-D concentration profile $c_n(\psi _{1,n},\psi _{2,n})$ (coloured rectangles) at each face $f_n$ of the fracture network shown in figure 1(b), represented in terms of the mixing graph ${\mathcal{G}_M}$ shown in figure 9(c). As shown, the inlet concentration field at the fracture inlets $f^o_{3,1}$ and $f^o_{4,1}$ comprises red and blue coloured segments that only vary in the $\psi _1$ direction. Dashed green horizontal (vertical) lines represent splitting or ‘cutting’ due to intersection (fracture) mixing, and solid green horizontal (vertical) lines represent merging or ‘shuffling’ due to intersection (fracture) mixing.

Figure 10 shows evolution of $c_n$ across the intersection faces $f_n$ for the simple fracture network shown in figure 1(b) (under the PWI model for ${\mathcal{M}_I}$ ). This figure illustrates how advective mixing acts to cut and shuffle the local concentration field $c_n$ in both the $\psi _{1,n}$ and $\psi _{2,n}$ directions via the respective fracture ( ${\mathcal{M}_F}$ ) and intersection ( ${\mathcal{M}_I}$ ) mixing mechanisms.

To examine mixing over the entire fracture network, it is useful to construct the outlet concentration distribution $c_{\textit{out}}$ as the aggregated concentration distribution over the network outlets. To do so, we aggregate all inlet faces $f^i_{m,1}$ of the outlet intersections $I_m$ with $m=N_I - N_{I_o} + 1:N_I$ by defining a global outlet streamfunction $\hat \psi _1^{\textit{out}} \in [0,1]$ that concatenates the local $\psi _1$ streamfunctions $\psi _{1,m}$ weighted by their respective fluxes $Q_m$ as

(3.18) \begin{equation} \hat \psi _1^{\textit{out}} \equiv \frac {1}{Q} \left ( \sum _{j=N_I - N_{I_o} + 1}^{m-1} Q_j + Q_m \psi _{1,m} \right ),\qquad m=N_I - N_{I_o} + 1: N_I. \end{equation}

As these outlets are concatenated along the $\psi _1$ coordinate, the global $\psi _2$ streamfunction at the outlet is simply $\hat \psi _2^{\textit{out}}=\psi _{2,m}$ . Hence, the outlet global streamfunctions can be expressed in terms of the global streamfunctions as

(3.19) \begin{align} \hat \psi _1^{\textit{out}}= \hat \psi _1^{\textit{out}}\left( \hat \psi _1, \hat \psi _2\right), && \hat \psi _2^{\textit{out}}= \hat \psi _2^{\textit{out}}\left( \hat \psi _1, \hat \psi _2\right). \end{align}

The outlet concentration profile $c_{\textit{out}}$ is then

(3.20) \begin{align} c_{\textit{out}}\left(\hat \psi _1,\hat \psi _2\right)\equiv c_0\left[\hat \psi ^{\textit{out}}_1\left(\hat \psi _1,\hat \psi _2\right),\hat \psi ^{\textit{out}}_2\left(\hat \psi _1,\hat \psi _2\right)\right], \end{align}

which provides a convenient basis for quantification of mixing over the entire fracture network. Figure 11 shows how the inlet concentration $c_0(\hat \psi _1,\hat \psi _2)$ is mapped to the outlet concentration $c_0(\hat \psi _1^{\textit{out}},\hat \psi _2^{\textit{out}})$ for the concentration field shown in figure 10.

Figure 11.

Mapping of the inlet concentration $c_0(\hat \psi _1,\hat \psi _2)$ shown in figure 10 to the outlet concentration $c_{\textit{out}}(\hat \psi _1,\hat \psi _2)\equiv c_0(\hat \psi ^{\textit{out}}_1,\hat \psi ^{\textit{out}}_2)$ . The resultant concentration profile is a direct result of the CS actions mediated by the separating streamsurfaces (indicated by green dashed lines) in the DFN.

In Appendix C, we show that under purely advective (non-dissipative) mixing, typical mixing measures such as concentration variance do not alter between $c_0$ and $c_{\textit{out}}$ . Instead we use a multi-scale mixing measure, the (squared) mix-norm $\varPhi ^2_{2\textit{-D}}$ (Mathew, Mezić & Petzold Reference Mathew, Mezić and Petzold2005), which robustly measures the degree of mixing by integrating over the scale-dependent concentration variance $\phi ^2_{2-D}(c_{\textit{out}},s)$ as

(3.21) \begin{equation} \varPhi ^2_{2\textit{-D}}\equiv \frac {1}{L}\int _0^L \phi ^2_{2\textit{-D}}(c_{\textit{out}},s)\,\text{d}s, \end{equation}

see Appendix C for details.

The 2-D concentration field $c_{\textit{out}}(\hat \psi _1,\hat \psi _2)$ and mixing measures $\phi ^2_{2\textit{-D}}$ , $\varPhi ^2_{2\textit{-D}}$ provide a complete quantification of advective mixing over the fracture network. In Appendix C, we also show that a coarse-grained 1-D ( $\hat \psi _1$ -only) representation also provides a highly accurate representation of advective mixing for $\ell \ll L$ . In some instances, this coarse-grained representation is preferable to simplify visualisation of mixing and connect with observable concentration fields. Following Appendix C, a simple coarse-graining approach is employed which consists of averaging $c_0$ and $c_{\textit{out}}$ with respect to $\hat \psi _2$ to yield

(3.22) \begin{equation} c(\hat \psi _1)\equiv \int _0^1 c\left(\hat \psi _{1},\hat \psi _{2}\right)\,\text{d}\hat \psi _2,\quad c=c_0, c_{\textit{out}}, \end{equation}

where $c(\hat \psi _1)$ is differentiated from $c(\hat \psi _1,\hat \psi _2)$ by its argument. Similarly, the global streamfunctions can be coarse-grained as $\hat \psi _2^{\textit{out}})\mapsto \hat \psi _1^{\textit{out}}$ . Note that this averaging yields a simplified representation of advective mixing rather than a physical mixing process. The associated 1-D scale-dependent variance $\phi ^2_{1\textit{-D}}(s)$ and mix-norm $\varPhi ^2_{1\textit{-D}}$ are defined analogously as their respective 2-D counterparts in (C5), (C6), and are related to the 2-D measures as

(3.23) \begin{align} \phi ^2_{1\textit{-D}}(s)= \begin{cases} 0\quad &0\lt s\leqslant \ell ,\\ \phi ^2_{2\textit{-D}}(s)\quad &\ell \lt s\leqslant L, \end{cases} && \varPhi ^2_{1\textit{-D}}=\varPhi ^2_{2\textit{-D},L}. \end{align}

Hence, the 1-D representation provides an accurate representation of mixing dynamics and the associated 1-D mixing measures $\phi ^2_{1\textit{-D}}$ , $\varPhi ^2_{1\textit{-D}}$ provide accurate estimates of their 2-D counterparts $\phi ^2_{2\textit{-D}}$ , $\varPhi ^2_{2\textit{-D}}$ .

Many 2-D DFN models perform similar averaging on physical grounds, based on the assumption of complete mixing at fracture intersections. The validity of this assumption is quantified by the two Péclet numbers associated with the $\psi _1$ and $\psi _2$ coordinates, defined as

(3.24) \begin{align} Pe_\|\equiv \frac {v_c L}{D_m}, && Pe_\bot \equiv \frac {v_c\ell }{D_m}, \end{align}

where $v_c$ is the characteristic velocity scale and $D_m$ is the diffusion coefficient. Under such models, solute mixing is considered to be advection-dominated in the $\psi _1$ direction ( $Pe_\|\gg 1$ ) and diffusion-dominated in the $\psi _2$ direction ( $Pe_\bot \ll 1$ ), justifying the ‘well-mixed’ assumption commonly evoked at fracture intersections (Hyman & Jiménez-Martínez Reference Hyman and Jiménez-Martínez2018). However, such models cannot be readily generalised beyond these limits. Note that the averaging process (3.22) is not equivalent to such physical diffusion as it is applied as a post-processing step to the outlet concentration distribution $c_{\textit{out}}$ rather than a continual diffusion process throughout the fracture network.

Figure 12.

(a,c) Isometric and (b,d) top views of test (a, b) DFN A and (c, d) DFN B.

4.1. Test DFNs and numerical method

The mixing measures $c_{\textit{out}}$ , $\varPhi ^2_{1\textit{-D}}$ , $\phi ^2_{1\textit{-D}}(s)$ used to validate the advective mixing theory and associated mixing graph ${\mathcal{G}_M}$ by comparison with numerical simulations performed via the 2-D DFN simulation code DFN.lab (Le Goc et al. Reference Le Goc, Pinier, Darcel, Lavoine, Doolaeghe, de Simone, De Dreuzy and Davy2019). We consider two different test DFNs and test ${\mathcal{G}_M}$ under both the PWI (3.8) and PWS (3.9) variants of the intersection map ${\mathcal{M}_I}$ . The first DFN considered (denoted DFN A) is shown in figures 12(a) and 12(b), and corresponds to the DFN shown in figure 1(b), which is less disordered than the second DFN considered (denoted DFN B) shown in figures 12(c) and 12(d). For both DFNs, flow in each fracture is modelled as a steady planar 2-D Darcy flow with uniform transmissivity. A prescribed pressure differential between the DFN inlets and outlets fixes the total volumetric flow at $Q=1$ . The inlet concentration is set to $c_0(\hat \psi _1,\hat \psi _2)=0$ in one inlet fracture and $c_0(\hat \psi _1,\hat \psi _2)=1$ in the other. The two outlets collect the discharge and all other boundaries are no-flux. DFN.lab is used to solve the flow and advect approximately $10^6$ streamlines that are seeded along the network inlets in a flux-weighted manner. Note that while a large number of streamlines are required to robustly test the accuracy of the mixing graph, only a small number of critical streamlines are required to quantify the mixing graph as described in § 3.5.

4.2. Construction of mixing graph ${\mathcal{G}_M}$

From these numerical results, the topology of the mixing graph ${\mathcal{G}_M}$ is then determined by the connectivity (via streamlines) between fractures and intersections, and the edge weights of ${\mathcal{G}_M}$ are determined by the flow rates $Q_n$ between faces. The intersection ${\mathcal{M}_I}$ and fracture ${\mathcal{M}_F}$ maps that operate respectively on intersection edges and fracture edges are determined from streamline routing. Note that for the PWI model (3.8), the intersection maps ${\mathcal{M}_I}$ are determined solely from streamline routing data and the local $\psi _1$ streamfunction is invariant over an intersection. Conversely, $\psi _1$ can vary under the PWS model (3.10), as is quantified by the flux distribution $q_n$ over face $f_n$ as detailed in Appendix B. Two-dimensional DFN models such as DFN.lab do not resolve the $\psi _2$ coordinate, but instead perform random flux-weighted streamline routing at intersections, which is physically consistent with well-mixed conditions in the $\psi _2$ coordinate and the 1-D mixing representation outlined in Appendix C. Conversely, the PWI and PWS models for ${\mathcal{G}_M}$ perform completely deterministic streamline routing based on the mixing maps ${\mathcal{M}_I}$ , ${\mathcal{M}_F}$ . We note that while the methods detailed in § 3 allow both $\psi _1$ and $\psi _2$ to be resolved, in this section we only compare 1-D mixing results between DFN.lab and the mixing graph ${\mathcal{G}_M}$ .

Development and application of the mixing graph ${\mathcal{G}_M}$ and the mixing maps ${\mathcal{M}_I}$ (PWI model) and ${\mathcal{M}_F}$ for DFN A results in the 2-D concentration fields shown in figure 10. This clearly illustrates the impact of intersection and fracture mixing mechanisms upon the local concentration field $c_n(\psi _{1,n},\psi _{2,n})$ at each face $f_n$ , and the corresponding 1-D mixing plot is shown in figure 18(d).

4.3. Comparison of advective mixing in test DFNs

Figure 13 compares the 1-D mixing plots (see details in Appendix D) for both test DFNs generated by DFN.lab and ${\mathcal{G}_M}$ under both the PWI and PWS models. Figure 13(a) and table 1 show that the PWI model fairly accurately captures the mixing dynamics generated by DFN.lab for DFN A in terms of both the outlet concentration distribution $c_{\textit{out}}(\hat \psi _1)$ and mapping $\hat \psi _1\mapsto \hat \psi _1^{\textit{out}}$ , although some minor differences are apparent. In these plots, CS of fluid elements is clearly shown by the cut and shuffled 1-D lines in the $\hat \psi _1-\hat \psi _1^{\textit{out}}$ plane, as well as the impact of intersection mixing which dilutes the 1-D outlet concentration field as $0\lt c_{\textit{out}}(\hat \psi _1)\lt 1$ . An inherent property of the PWI model is that it produces straight lines in the $\hat \psi _1-\hat \psi _1^{\textit{out}}$ plane as fluid parcels only undergo CS in the $\psi _1$ coordinate without deformation (stretching or compression); hence, the mapping $\hat \psi _1\mapsto \hat \psi _1^{\textit{out}}$ is piecewise linear. Conversely, the DFN.lab results show some evidence of streamline deformation, as indicated by the slight curvature of lines in the $\hat \psi _1-\hat \psi _1^{\textit{out}}$ plane.

Figure 13.

Comparison of the 1-D mixing plots for (a, b) DFN A and (c,d) DFN B generated by DFN.lab and ${\mathcal{G}_M}$ under the (a,c) PWI and (b,d) PWS models. Colour bars denote inlet $c_0(\hat \psi _1)$ and outlet $c_{\textit{out}}(\hat \psi _1)$ concentration distributions, and shaded bands indicate overlap regions of intersection mixing leading to $0\lt c_{\textit{out}}(\hat \psi _1)\lt 1$ .

Figure 13(b) and table 1 shows that for DFN A, the PWS model predicts the outlet concentration $c_{\textit{out}}(\hat \psi _1)$ and streamfunctions $\hat \psi _1^{\textit{out}}$ almost exactly due to resolution of fluid deformation during advective mixing. The difference between the PWI and PWS models is more pronounced for DFN B, consistent with the more disordered nature of this DFN. Figure 13(c) shows a significant level of discrepancy between the PWI model and the DFN.lab results, whereas the PWS model (figure 13 d) accurately reproduces the mixing dynamics generated by DFN.lab. As many fracture networks typically exhibit a much larger level of disorder than DFN B, the impact of fluid deformation upon advective mixing is likely to be more pronounced than that indicated in figures 13(c) and 13(d).

Table 1.

Comparison of 1-D mix-norm values $\varPhi ^2_{1\textit{-D}}$ for the test DFNs between the DFN.lab results and the PWI and PWS models for ${\mathcal{G}_M}$ .

Figure 14 shows that that while the PWS model results very closely match the DFN.lab results for both test DFNs, the PWI model results slightly deviate at small to intermediate scales for DFN A and in a more pronounced fashion for DFN B. This is consistent with the $\varPhi ^2_{1\textit{-D}}$ values in table 1, which indicates that the PWI model for DFN A and B incurs respective mix-norm errors of 1.4 % and 10.5 %, whereas the PWS model is highly accurate. The excellent level of agreement between DFN.lab and the PWS model confirms that advective mixing in fracture networks follows a PWS discontinuous mixing process.

Figure 14.

Comparison of multi-scale variance $\phi ^2_{1\textit{-D}}(s)$ for (a) DFN A and (b) DFN B between DFN.lab results (blue dots) and the PWS (orange crosses) and PWI (grey crosses) models for ${\mathcal{G}_M}$ .

4.4. Advective mixing in iterated fracture networks

To highlight the utility of the graph-based approach, we also consider iterated mixing in DFN B, where the outlet concentration distribution $c_{\textit{out}}(\hat \psi _1)$ for one iteration of the network forms the inlet concentration distribution $c_0(\hat \psi _1)$ for the following iteration. The 1-D mixing plots in figure 15 show that under both the PWI and PWS models, this DFN exhibits rapid and complete mixing with increasing iteration number $n$ , where the outlet concentration $c_{\textit{out}}\rightarrow 1/2$ as $n\rightarrow \infty$ . Similarly, figure 16 shows that the 1-D mix-norm $\varPhi ^2_{1\textit{-D}}(n)$ for both models converges towards zero with $n$ as a power-law ( $\varPhi ^2_{1\textit{-D}}(n)\sim n^{-7}$ ), reflecting weak ergodic mixing (algebraic decay) of the PWI and PWS transforms. Figure 16 also shows that $\varPhi ^2_{1\textit{-D}}(n)$ is quite similar between the PWI and the PWS models, suggesting that although the PWI model does not resolve the detailed mixing dynamics, this has minimal impact on quantifying mixing in larger fracture networks.

Figure 15.

Evolution of 1-D mixing plots for iterated mixing in DFN B under the (a–d) PWS and (e–h) PWI mixing graph representations after $n=2$ , $n=4$ , $n=8$ and $n=16$ iterations.).

Figure 16.

Decay of 1-D mix-norm $\varPhi ^2_{1\textit{-D}}(n)$ as a function of iteration number $n$ for the orthogonal DFN under the PWS and PWI models.