2 Pore-scale chaotic advection in 3D porous media

In this section we summarize the mechanisms that generate pore-scale chaotic advection as described in Lester et al. (Reference Lester, Metcalfe and Trefry2013), and further explain the role of critical points of the skin friction field. Topological complexity is a defining feature of all porous media, from granular and packed media to fractured and open networks; these materials are all typified by a highly connected pore space. Such topological complexity is characterized by the Euler characteristic $\unicode[STIX]{x1D712}$ (related to the topological genus $g$ as $\unicode[STIX]{x1D712}=2(1-g)$ ) which measures the connectivity of the pore space (Vogel Reference Vogel, Mecke and Stoyan2002) as

where $N$ is the number of pores, $C$ is the number of redundant connections and $H$ is the number of completely enclosed cavities. For porous media it is more meaningful to consider the average number density of these quantities, where from computer tomography studies (Vogel Reference Vogel, Mecke and Stoyan2002) it is found that typically $N$ is large while $C$ and $H$ are small. Hence, the number density of the Euler characteristic $\unicode[STIX]{x1D712}$ is typically large and negative (Vogel Reference Vogel, Mecke and Stoyan2002; Scholz et al. Reference Scholz, Wirner, Götz, Rüde, Schröder-Turk, Mecke and Bechinger2012), reflecting topologically complex media.

When a continuous fluid is advected through such media, a large number of critical points (zeros) $\boldsymbol{x}_{p}$ arise in the skin friction field $\boldsymbol{u}(\boldsymbol{x})$ at the boundary $\unicode[STIX]{x2202}\mathscr{D}$ of the fluid domain $\mathscr{D}$ . The skin friction field $\boldsymbol{u}(\boldsymbol{x})$ is defined as

where $\boldsymbol{v}(x_{1},x_{2},x_{3})$ is the velocity field, $x_{3}$ is the coordinate normal to the fluid boundary $\unicode[STIX]{x2202}\mathscr{D}$ and the $x_{1}$ , $x_{2}$ coordinates are tangent to this boundary (de Winkel & Bakker Reference de Winkel and Bakker1988; Surana, Grunberg & Haller Reference Surana, Grunberg and Haller2006; Chong et al. Reference Chong, Monty, Chin and Marusic2012). The critical points $\boldsymbol{x}_{p}$ of the skin friction field $\boldsymbol{u}$ (where $\boldsymbol{u}(\boldsymbol{x}_{p})=\mathbf{0}$ ) represent the boundary counterparts to the stagnation points $\boldsymbol{x}_{0}$ of the interior flow field $\boldsymbol{v}$ (where $\boldsymbol{v}(\boldsymbol{x}_{0})=\mathbf{0}$ ). It should be noted that the majority of stagnation points in pressure-driven Stokes flow with immobile boundaries are exterior points, as mechanisms to generate flow reversal (e.g. Moffatt vortices (Moffatt Reference Moffatt1964)) are rare.

The Poincaré–Hopf theorem states that the sum of the indices $\unicode[STIX]{x1D6FE}_{p}$ of critical points $\boldsymbol{x}_{p}$ is related to the topological genus $g$ and Euler characteristic $\unicode[STIX]{x1D712}$ as

where the index $\unicode[STIX]{x1D6FE}_{p}=-1$ for saddle points, $\unicode[STIX]{x1D6FE}_{p}=1$ for node points and $\unicode[STIX]{x1D6FE}_{p}=0$ for centres of the skin friction field. Hence, $|\unicode[STIX]{x1D712}|$ represents a lower bound for the number density of saddle points under steady 3D Stokes flow, which, as shall be shown, impart significant fluid stretching to the fluid interior. A relevant question is what the impact of this stretching is with respect to mixing in porous media.

Critical points play an important role with respect to the Lagrangian dynamics of 3D steady flows, as it is at these points that the formal analogy between transport in steady 3D volume-preserving flows and one-degree-of-freedom Hamiltonian systems breaks down (Bajer & Moffatt Reference Bajer and Moffatt1990; Bajer Reference Bajer1994), and such points are widely (Mezić & Wiggins Reference Mezić and Wiggins1994; Wiggins Reference Wiggins2010) implicated in the creation of Lagrangian chaos. Associated with these critical points are stable $\mathscr{W}^{s}$ and unstable $\mathscr{W}^{u}$ manifolds which project into the fluid interior, and these may be either one-dimensional (1D) or 2D, as shown in figure 1. MacKay (Reference MacKay1994, Reference MacKay, Chandre, Leoncini and Zaslavsky2008) proposes that the manifolds that are 2D represent surfaces of minimal transverse flux, and so the set of these co-dimension 1 manifolds forms the ‘skeleton’ of the flow which organizes transport and mixing throughout the fluid domain. While 1D manifolds also play a role, they are less important as they do not separate the fluid domain.

Figure 1.

Schematic of type I–IV critical points of the skin friction field $\boldsymbol{u}$ on the boundary $\unicode[STIX]{x2202}\mathscr{D}$ and associated stable and unstable manifolds $\mathscr{W}^{s}$ , $\mathscr{W}^{u}$ . The double arrows reflect the eigenvalue sum $\unicode[STIX]{x1D702}_{1}+\unicode[STIX]{x1D702}_{2}+2\unicode[STIX]{x1D702}_{3}=0$ in (2.5).

The dimensionality of these interior manifolds $\mathscr{W}^{s}$ , $\mathscr{W}^{u}$ is dictated by the type of critical point (saddle, node, centre) and hence the index $\unicode[STIX]{x1D6FE}_{p}$ . Fluid deformation local to critical points $\boldsymbol{x}_{p}$ is quantified by the linearized skin friction tensor $\unicode[STIX]{x1D63C}=\unicode[STIX]{x2202}\boldsymbol{u}/\unicode[STIX]{x2202}\boldsymbol{x}|_{\boldsymbol{x}_{p}}$ which is the boundary counterpart of the strain-rate tensor $\unicode[STIX]{x1D641}=\unicode[STIX]{x2202}\boldsymbol{v}/\unicode[STIX]{x2202}\boldsymbol{x}|_{\boldsymbol{x}_{0}}$ which quantifies deformation local to the stagnation points $\boldsymbol{x}_{0}$ . The components $a_{ij}$ of the skin friction tensor $\unicode[STIX]{x1D63C}$ are given by the expansion

and de Winkel & Bakker (Reference de Winkel and Bakker1988) and Chong et al. (Reference Chong, Monty, Chin and Marusic2012) have derived the normal form of $\unicode[STIX]{x1D63C}$ around a critical point $\boldsymbol{x}_{p}$ in an incompressible flow as

For incompressible flow the eigenvalues $\unicode[STIX]{x1D708}_{i}$ of $\unicode[STIX]{x1D641}$ satisfy $\sum _{i=1}^{3}\unicode[STIX]{x1D708}_{i}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0$ , whereas the eigenvalues $\unicode[STIX]{x1D702}_{i}$ of $\unicode[STIX]{x1D63C}$ satisfy $\sum _{i=1}^{3}\unicode[STIX]{x1D702}_{i}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=-\unicode[STIX]{x1D702}_{3}$ , where $\unicode[STIX]{x1D702}_{1}$ , $\unicode[STIX]{x1D702}_{2}$ are the eigenvalues in the skin friction boundary (with $\unicode[STIX]{x1D702}_{1}\leqslant \unicode[STIX]{x1D702}_{2}$ ), and $\unicode[STIX]{x1D702}_{3}$ is the interior eigenvalue. Hence, critical points consist of one stable and two unstable eigenvalues, or vice versa, and $\unicode[STIX]{x1D702}_{3}>0$ for separation points and $\unicode[STIX]{x1D702}_{3}<0$ for reattachment points. As (2.5) indicates that the boundary $\unicode[STIX]{x2202}\mathscr{D}$ is invariant, there exist four critical point types as shown in figure 1:

1. (1) type I separation point $\unicode[STIX]{x1D702}_{3}>0$ , attractor for $\boldsymbol{u}$ with $\unicode[STIX]{x1D702}_{1},\unicode[STIX]{x1D702}_{2}<0$ ;

2. (2) type II reattachment point $\unicode[STIX]{x1D702}_{3}<0$ , repeller for $\boldsymbol{u}$ with $\unicode[STIX]{x1D702}_{1},\unicode[STIX]{x1D702}_{2}>0$ ;

3. (3) type III separation point $\unicode[STIX]{x1D702}_{3}>0$ , saddle for $\boldsymbol{u}$ with $\unicode[STIX]{x1D702}_{1}<0$ , $\unicode[STIX]{x1D702}_{2}>0$ ;

4. (4) type IV reattachment point $\unicode[STIX]{x1D702}_{3}<0$ , saddle for $\boldsymbol{u}$ with $\unicode[STIX]{x1D702}_{1}<0$ , $\unicode[STIX]{x1D702}_{2}>0$ .

Types I, II are node points, while types III, IV are saddle points, and so the associated interior manifolds are 2D, forming barriers to transport in the fluid interior. Figure 2 illustrates these features for pore branch (a) and merge (b) elements (the union of which is termed a couplet) in an open porous network. Here, the saddle-type reattachment (type IV) and separation point (type III) points arise due to the basic topology of the pore branch and merger, and these generate the respective interior 2D manifolds $\mathscr{W}_{2D}^{u}$ , $\mathscr{W}_{2D}^{s}$ . While the geometry of pore branches and mergers may vary significantly, the basic topology shown in figure 2 is universal to almost all pore branches and mergers.

Figure 2.

Schematic of pore branch (a) and merge (b) elements, where the black dots show the critical separation and attachment points respectively and the grey surfaces denote the associated 2D stable and unstable manifolds representing 2D surfaces of locally minimum transverse flux, the ‘skeleton’ of the flow. The angle $\unicode[STIX]{x1D703}$ denotes the relative orientation of the pore branch and merger elements.

The interaction of the 2D interior manifolds ( $\mathscr{W}_{2D}^{s},\mathscr{W}_{2D}^{u}$ ) controls fluid transport and mixing. If these 2D manifolds intersect transversely a heteroclinic tangle and chaotic dynamics result, whereas a smooth heteroclinic connection between them yields a 2D streamsurface that acts as a separatrix between saddle points. For the pore junction shown in figure 2, significant fluid stretching occurs transverse to the 2D interior manifolds, and folding of material elements arises from advection over the critical points $\boldsymbol{x}_{p}$ . These stretching and folding actions are the constituent motions of the Smale horseshoe map, a hallmark of chaotic dynamics in continuous systems. If a connected pore branch and merger are oriented symmetrically (as reflected by zero offset angle $\unicode[STIX]{x1D703}$ in figure 2), then this symmetry results in a smooth heteroclinic connection between the 2D manifolds. Conversely, when these pore elements are oriented with non-zero offset angle $\unicode[STIX]{x1D703}\neq 0$ , this symmetry is broken, resulting in transverse heteroclinic intersection of the 2D manifolds and hence chaotic streamline dynamics. The disordered nature of random porous media ensures that $\unicode[STIX]{x1D703}\neq 0$ and these symmetries are broken, leading to chaotic advection under steady flow conditions. While porous media may exhibit other features such as surface roughness or pore tortuosity which may also generate chaotic advection (Jones, Thomas & Aref Reference Jones, Thomas and Aref1989; Ottino & Wiggins Reference Ottino and Wiggins2004), a large density of saddle-type critical points is generic to all complex media; hence, the mechanisms described above represent a lower bound for chaotic dynamics inherent to porous media flow.

3 Three-dimensional open porous network model

To study the impact of such topology-induced chaotic advection upon transport and mixing in porous media, we consider a model 3D open porous network which consists of a random network of connected pore junctions and mergers shown in figure 2. We consider a 3D random network of these pore branches and mergers which is the simplest representation of a topologically complex open porous network that is homogeneous at the macroscale. Due to the spatial symmetries and reversibility of Stokes flow in these elements, we shall use a series of linear spatial transformations (reflection, scaling, reorientation, inversion) to map fluid transport and deformation in Stokes flow in the elements in figure 2 to construct the entire network. In this way we can develop a highly efficient model of fluid transport and deformation over a random 3D network.

Figure 3.

(a) Schematic of pore branches (black) and mergers (white) between mapping planes, with branch pores (black) and merge pores (white) shown. (b) (Middle) evolution of $10^{5}$ non-diffusive particles in the pore network from continuous injection in a single pore; (upper left) typical transverse distribution of non-diffusive coloured fluid particles with longitudinal pore number $n$ ; (lower right) typical distribution of residence times (pink $=$ maximum, orange $=$ minimum) with longitudinal pore number  $n$ .

To compose a 3D network over the semi-infinite domain $\unicode[STIX]{x1D6FA}:\boldsymbol{x}\in \mathbb{R}^{3},z>0$ , we use the pore branch and merge elements shown in figure 2 to connect a series of so-called ‘mapping planes’ oriented parallel to the $(x,y)$ plane, distributed along the $z$ axis at integer multiples of the pore element length $\ell$ . Each mapping plane consists of an infinite number of reflections in the $x$ , $y$ directions of a periodic square $(x,y)\in [0,L]\times [0,L]$ which contains $3J$ randomly located non-overlapping pores, as shown in figure 3. The location of the $j$ th pore in the $n$ th mapping plane is labelled $\boldsymbol{r}_{n,j}=(x_{n,j},y_{n,j},n\ell )$ . Of these $3J$ pores, $J$ are randomly labelled as ‘branch’ pores (i.e. pores that branch next in the positive $z$ direction), and the remaining $2J$ pores are labelled as ‘merge pores’ (i.e. those that are about to merge). A branch pore at plane $n$ is connected to two merge pores at plane $n+1$ by a pore branch element located between these planes, and conversely two branch pores at plane $n$ are connected to a single merge pore at plane $n$ by a pore merge element as per figure 2.

To accommodate connections between the merge pore locations $\boldsymbol{r}_{m_{1}},\boldsymbol{r}_{m_{2}}$ and the branch pore location $\boldsymbol{r}_{b_{1}}$ at adjacent mapping planes, the pore branch and merge elements are rotated and stretched such that the $\hat{\boldsymbol{e}}_{y}$ , $\hat{\boldsymbol{e}}_{z}$ vectors are aligned along $\boldsymbol{r}_{m_{1}}-\boldsymbol{r}_{m_{2}}$ , $\boldsymbol{r}_{b_{1}}-(\boldsymbol{r}_{m_{1}}+\boldsymbol{r}_{m_{2}})/2$ respectively, and the element height and length are rescaled to $\Vert \boldsymbol{r}_{m_{1}}-\boldsymbol{r}_{m_{2}}\Vert$ , $\Vert \boldsymbol{r}_{b_{1}}-(\boldsymbol{r}_{m_{1}}+\boldsymbol{r}_{m_{2}})/2\Vert$ respectively. The orientation angle $\unicode[STIX]{x1D711}$ in the $(x,y)$ plane of a pore branch/merge element is then

and the relative angle $\unicode[STIX]{x1D703}$ between subsequent elements is defined as $\unicode[STIX]{x1D703}_{n}=\unicode[STIX]{x1D711}_{n+1}-\unicode[STIX]{x1D711}_{n}$ . Hence, the set of non-overlapping pore locations $\boldsymbol{r}_{n,j}$ for $n=0:\infty$ , $j=1:3J$ completely defines a realization of the 3D open porous network, and the ensemble of all realizations of the porous network is ergodic and stationary. The resultant random orientation angle $\unicode[STIX]{x1D711}$ is uniformly distributed and is uncorrelated both across and within mapping planes; hence, transport within the 3D open porous network may be described by a Markov process in space.

To model transport within the porous network, we solve the fluid velocity field $\boldsymbol{v}$ for Stokes flow $\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+\unicode[STIX]{x1D735}p=0$ within the pore branch element $\unicode[STIX]{x1D6FA}$ driven by a differential pressure $p|_{z=0}-p|_{z=L}=\unicode[STIX]{x0394}p$ , subject to the no-slip boundary condition $\boldsymbol{v}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}=0$ . The flow field is calculated numerically to order $10^{-16}$ root mean square accuracy using the finite-volume computational fluid dynamics (CFD) package ANSYS-CFX with double precision calculation. The advection dynamics for passive tracer particles within the branch element $\unicode[STIX]{x1D6FA}$ under Stokes flow can be represented by the spatial $\mathscr{M}^{\ast }$ and temporal $\mathscr{T}^{\ast }$ maps

where $\{x_{r,n},y_{r,n}\}$ are the $x,y$ particle positions relative to the inlet of $\unicode[STIX]{x1D6FA}$ (where $x_{r,n}^{2}+y_{r,n}^{2}=0,1$ correspond to the pore centre and boundary respectively), and $\{x_{r,n+1},y_{r,n+1}\}$ are the positions relative to the outlet of $\unicode[STIX]{x1D6FA}$ . It should be noted that $x_{r,n+1},y_{r,n+1}$ do not reflect the specific outlet pore the particle travels to, as a direct consequence of the reflection symmetry of the flow about the plane $y_{r}=0$ . Similarly, the temporal map $\mathscr{T}^{\ast }$ describes the residence time from the inlet to outlet planes. From the CFD results, the exact spatial $\mathscr{M}^{\ast }$ and temporal $\mathscr{T}^{\ast }$ maps are remarkably well approximated (within error $\unicode[STIX]{x1D716}\sim 10^{-3}$ ) by composites of the following simple analytic maps: where for a pore branch $\mathscr{M}^{\ast }\approx \mathscr{M}$ , $\mathscr{T}^{\ast }\approx \mathscr{T}\circ \mathscr{M}$ , and due to the reflection symmetry of Stokes flow $\mathscr{M}^{\ast }\approx \mathscr{M}^{-1}$ , $\mathscr{T}^{\ast }\approx \mathscr{T}\circ \mathscr{M}^{-1}$ for a pore merger. These approximate maps greatly simplify the description of the advection dynamics and preserve the essential features of advective transport, namely stretching in a pore branch (and contraction in a pore merger) of fluid elements by a factor of 2 in the $y$ direction, and the Poiseuille-like distribution of $\mathscr{T}$ (scaled such that the minimum advection time is unity) due to the no-slip boundary conditions. The composite $\mathscr{T}^{\ast }\approx \mathscr{T}\circ \mathscr{M}$ accurately captures the residence time distribution because the advection map $\mathscr{M}$ accurately captures the streamline geometry, and the no-slip boundary condition within  $\unicode[STIX]{x1D6FA}$ dominates temporal advection. Hence, a simple Poiseuille distribution concatenated with $\mathscr{M}$ accurately approximates the temporal map $\mathscr{T}^{\ast }$ . While the advective maps $\mathscr{M}$ , $\mathscr{M}^{-1}$ are not area-preserving, they are implicitly volume-preserving when deformation in the $z$ direction is considered.

While the symmetry of the pore branch and merger results in zero fluid stretching (as reflected by the concatenation $\mathscr{M}^{-1}\circ \mathscr{M}=I$ , where $I$ is the identity operator), these elements at arbitrary orientations in the $xy$ plane can break this symmetry and lead to persistent fluid stretching. The reoriented advective maps $\mathscr{S}(\unicode[STIX]{x1D711})$ , $\mathscr{S}^{-1}(\unicode[STIX]{x1D711})$ define transport in pore branches and mergers at arbitrary angles,

where $R(\unicode[STIX]{x1D711})$ is a rotation operator and $\mathscr{S}^{-1}(\unicode[STIX]{x1D711}_{n+1})\circ \mathscr{S}(\unicode[STIX]{x1D711}_{n})\neq I$ for $\unicode[STIX]{x1D703}_{n}=\unicode[STIX]{x1D711}_{n+1}-\unicode[STIX]{x1D711}_{n}\neq 0$ . Coupled with the residence time distribution $\mathscr{T}$ , the composite maps $\mathscr{S}$ , $\mathscr{S}^{-1}$ quantify advective transport relative to a local co-moving reference frame through the 3D porous network, and so the position and residence time of fluid particles can be rapidly propagated via these approximate maps, as shown in figure 3. Hence, the sequence of alignment angles $\unicode[STIX]{x1D711}_{1},\unicode[STIX]{x1D711}_{2},\ldots$ for a particular trajectory can impart either regular or chaotic advection at the pore scale (Lester et al. Reference Lester, Metcalfe and Trefry2013), and the associated infinite-time Lyapunov exponent for ordered porous media spans $\unicode[STIX]{x1D706}_{\infty }\in [0,\ln 2]$ . It should be noted that $\unicode[STIX]{x1D706}_{\infty }$ is the Lyapunov exponent per pore branch, related to the Lyapunov exponent per unit time $\unicode[STIX]{x1D706}_{\infty }^{\ast }$ as $\unicode[STIX]{x1D706}_{\infty }^{\ast }=\unicode[STIX]{x1D706}_{\infty }/\unicode[STIX]{x0394}t_{a}$ , where $\unicode[STIX]{x0394}t_{a}$ is the mean residence time per pore branch. Conversely, for random media $\unicode[STIX]{x1D703}_{n}$ is uniformly distributed over $[-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]$ and so always generates chaotic advection.

This basic pore network model has been used (Lester et al. Reference Lester, Metcalfe and Trefry2013, Reference Lester, Metcalfe and Trefry2014) to study the impact of pore-scale chaotic advection upon dispersion in porous media, i.e. the spatial spreading of transport particles. Figure 3(b) illustrates the evolution of $10^{5}$ non-diffusive particles propagated by the composite maps $\mathscr{S}$ , $\mathscr{S}^{-1}$ through a realization of the random network. The lower right subfigure in figure 3(b) illustrates the creation of fine-scale structure via chaotic advection within the pore scale (from an initially segregated distribution of half blue/half red points across all pores), and the upper left subfigure shows the evolution of the residence time distribution evolved via the composite temporal maps $\mathscr{S}\circ \mathscr{T}$ , $\mathscr{S}^{-1}\circ \mathscr{T}$ . It should be noted that there is no visual evidence of folding of fluid elements in the $xy$ plane, as folding occurs around critical points $\boldsymbol{x}_{p}$ as fluid elements are advected downstream in the $z$ direction.

4 Fluid advection and deformation in 3D porous networks

5 Stretching CTRW

To quantify pore-scale mixing in terms of the stretching and advection dynamics considered above, we propose to represent the stochastic stretching process as a random walk in the space of time and elongation. As the orientation angles $\unicode[STIX]{x1D711}_{s}$ , $\unicode[STIX]{x1D711}_{c}$ for the pore branch and merge elements are random, the elongation rates are uncorrelated across successive pores. Hence, the successive material line elongations follow a Markov process defined by (4.9). Similarly, because of the random reorientation, the travel times across successive pores are also weakly correlated (Lester et al. Reference Lester, Metcalfe and Trefry2013). Hence, both the series of elongations and the travel times follow Markov processes in space. This implies that the deformation of a fluid material element can be modelled as a CTRW, whereby fluid elements undergo a series of punctuated stretching and compression events occurring in random times as they propagate through the pore network. While most CTRW models consider random spatial increments in random times to model spatial transport processes, we develop a stretching CTRW model that considers random elongation events occurring in random times to represent the stochastic stretching process.

The approach we shall employ to model fluid mixing is to consider both the deformation and the advection of an infinitesimal fluid element in the absence of diffusion as random walk processes conditioned upon the transport dynamics derived in the previous section. We then account for molecular diffusion by utilizing the deformation history in the CTRW to solve a local ADE, and the ensemble of such elements informs the evolution of the concentration PDF in the pore space. It should be noted that as the diffusion is fully coupled to the deformation we account explicitly for stretching enhanced diffusion.

The two-step stretching process (4.9) over a couplet provides a natural discretization for a CTRW of fluid deformation, and the temporal map $\mathscr{T}$ also quantifies the residence time. For this random walk model we do not need to resolve the particle location in terms of pore-space coordinates $(x_{r},y_{r})$ , but rather only consider the distributions of the travel time $\unicode[STIX]{x0394}t$ and stretch $\unicode[STIX]{x1D70C}_{s},\unicode[STIX]{x1D70C}_{c}$ derived from the relationships in § 4. The resultant two-step CTRW over a couplet is then

where the advection time between pore elements $\unicode[STIX]{x0394}t_{n}$ is distributed as $\unicode[STIX]{x1D713}(\unicode[STIX]{x0394}t)$ (4.18), and the random alignment angles $\unicode[STIX]{x1D711}_{2n}$ , $\unicode[STIX]{x1D711}_{2n+1}$ are uniformly distributed over $[-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]$ . Thus, the stretching CTRW quantifies the elongation of the length $\unicode[STIX]{x1D6FF}l_{n}$ of a material element at subsequent downstream positions $z_{n}$ , and the deformation time through the temporal random walk $t_{n}$ . For simplicity, we approximate the deformation rates during a stretching or compression transition as being constant. While the multiplicative two-step nature of the CTRW (5.1) does not permit analytic solution, analytic solutions are possible for an equivalent additive one-step CTRW over a couplet in terms of $s=\ln \unicode[STIX]{x1D6FF}l$ ,

It should be noted that the fluid stretching for this one-step CTRW is not uniform over a couplet, but rather is oscillatory due to stretching and compression over a pore branch or merger respectively. While such oscillations do not impact the overall stretching rate $\unicode[STIX]{x1D706}_{\infty }$ , such oscillatory stretching may generate a significant diffusive signature over the couplet. We shall return to this issue in § 6 where we develop a pore-scale mixing model based upon this one-step CTRW (5.2).

Figure 5.

(a) Schematic of the lamella formed by the 2D material sheet $\mathscr{H}_{2D}$ arising from continuous injection of a 1D material line in the inlet plane (b). The lamella $\mathscr{H}_{2D}$ (solid blue line) forms the backbone of the diffusive concentration field $c(\boldsymbol{x})$ (blurred blue region), shown here intersecting a pore outlet in the $xy$ plane, injected along the same material line in the inlet plane, resulting in a Gaussian concentration distribution transverse to $\mathscr{H}_{2D}$ . Here, $(\unicode[STIX]{x1D702},\unicode[STIX]{x1D701})$ are the material coordinates normal and tangential to $\mathscr{H}_{2D}$ in the $xy$ plane.

6 Scalar mixing and fluid stretching CTRW

To model pore-scale solute mixing we employ a 1D lamellar mixing model (Ranz Reference Ranz1979; Duplat & Villermaux Reference Duplat and Villermaux2008; Le Borgne et al. Reference Le Borgne, Dentz and Villermaux2015), where fluid stretching is quantified by the CTRW model developed in § 5. Due to exponential fluid stretching, an evolving concentration field conforms to a highly striated lamellar structure, where these lamellae align with the dominant stretching direction(s) of the flow and are simultaneously stretched longitudinally and compressed transversely.

The continuous injection of a 1D material line $\mathscr{L}_{1D}$ of length $l_{0}$ across a single pore of the inlet plane ( $z=n=0$ ) of the 3D network forms a continuous 2D material sheet $\mathscr{H}_{2D}$ shown in figure 5(a) which propagates throughout the pore space of the network and branches into multiple pores. We denote this 2D material sheet $\mathscr{H}_{2D}$ as a single lamella, and define the material coordinates $\{\unicode[STIX]{x1D702},\unicode[STIX]{x1D701}\}$ in the $xy$ plane, where $\unicode[STIX]{x1D701}$ is the coordinate along the lamella and $\unicode[STIX]{x1D702}$ is the transverse coordinate. We consider the steady concentration field $c(\boldsymbol{x})=c(\unicode[STIX]{x1D702},\unicode[STIX]{x1D701},z)$ which arises from the continuous injection of a concentration strip along the 1D material line $\mathscr{L}_{1D}$ with a Gaussian transverse profile

with maximum initial concentration $c_{m}(0)=c_{0}/\sqrt{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}_{0}^{2}}$ and variance $\unicode[STIX]{x1D70E}_{0}^{2}$ . As shown in figure 5(b), the lamella $\mathscr{H}_{2D}$ then forms the ridge of the 3D solute concentration field $c(\boldsymbol{x})$ throughout the pore network, and the concentration transverse to $\mathscr{H}_{2D}$ follows a Gaussian profile. For simplicity of exposition we consider a uniform solute injection applied at the inlet plane $z=0$ across the entire network. This means that we assume the same ‘initial’ concentration profile (6.1) across the lamellae in all pores located at the inlet plane at $z=0$ . Thus, the composition of the concentration field does not depend on the position transverse to the mean flow direction and only depends on the longitudinal distance $z$ . We shall extend these results in § 7 to the case of a heterogeneous solute injection across the inlet plane.

The macroscopically homogeneous initial condition described above results in an ensemble of 2D lamellae $\mathscr{H}_{2D}$ which are homogeneously distributed throughout the pore network. Hence, lamella segments from different inlet pores may coexist within a single pore, and the spacing between lamella segments decays in the longitudinal direction as the exponential stretching rate; however, these segments never overlap due to their invariant nature. Conversely, the associated 3D concentration field $c(\boldsymbol{x})$ has finite width and so will eventually overlap and coalesce. While a pore may contain segments of different lamellae, this redistribution of lamellar segments between pores does not impact the concentration distribution at a statistical level due to the homogeneous inlet condition described above. Hence, the statistical concentration distribution is the same across all pores for fixed $n$ and we do not need to resolve the redistribution of lamellar segments across pores. This homogeneous inlet condition greatly simplifies the exposition of the method, but in principle may be extended to inhomogeneous inlet conditions by explicitly tracking this redistribution, facilitating solution of both the evolving pore-scale and the macroscopic concentration distributions throughout the network.

The steady 3D concentration field $c(\boldsymbol{x})$ in the network is governed by the steady ADE

where $D_{m}$ is the molecular diffusivity and $\boldsymbol{n}$ is the outward normal from the boundary $\unicode[STIX]{x2202}\mathscr{D}$ . Under the assumption that the spatial concentration gradients are small in the $z$ direction, we ignore longitudinal diffusion and pose the steady 3D ADE (6.2) as an equivalent unsteady 2D ADE in the $xy$ plane,

where $v_{z}$ is the $z$ component of $\boldsymbol{v}(\boldsymbol{x})$ and $\bot$ denotes the $xy$ components. The $z$ coordinate is also parameterized in terms of the Lagrangian travel time $t=t(z)$ along a fluid particle trajectory as

The lamellar mixing model is based upon posing (6.3) in terms of the material coordinates $\{\unicode[STIX]{x1D702},\unicode[STIX]{x1D701}\}$ in the $xy$ plane; hence, $c(x,y,t)\mapsto c(\unicode[STIX]{x1D702},\unicode[STIX]{x1D701},t)$ . Under exponential fluid stretching, concentration gradients along the lamella ( $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D701}}$ ) decay exponentially while the gradients normal to the lamella ( $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}$ ) are maintained; thus, level sets of the concentration field identify with the lamella. This alignment occurs regardless of the inlet condition, and so the lamellar approach is quite generic. Under these conditions the ADE (6.3) simplifies to the 1D lamellar transport (Ranz Reference Ranz1979; Duplat & Villermaux Reference Duplat and Villermaux2008; Villermaux Reference Villermaux2012; Le Borgne et al. Reference Le Borgne, Dentz and Villermaux2013)

where the stretching rate $\unicode[STIX]{x1D6FE}(t,\unicode[STIX]{x1D701})\equiv \unicode[STIX]{x2202}_{\unicode[STIX]{x1D701}}v_{\unicode[STIX]{x1D701}}=-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}v_{\unicode[STIX]{x1D702}}$ arises from a first-order expansion of the $\unicode[STIX]{x1D702}$ component of the velocity field, and the stretching rate $\unicode[STIX]{x1D6FE}(t,\unicode[STIX]{x1D701})$ can be expressed in terms of the lamella elongation in the $xy$ plane $l(t,\unicode[STIX]{x1D701})$ as

It should be noted that the 1D nature of (6.5) renders it only valid up to the coalescence of lamellae in the pore space. While methods (Duplat & Villermaux Reference Duplat and Villermaux2008; Villermaux Reference Villermaux2012; Le Borgne et al. Reference Le Borgne, Dentz and Villermaux2013) exist to extend beyond coalescence and are readily applicable to this problem, they are beyond the scope of the present study.

Due to negligible transport in the $\unicode[STIX]{x1D701}$ direction, equation (6.5) can be posed with respect to a single fluid element trajectory within the lamella, and describes the evolution of the associated transverse concentration $c(\unicode[STIX]{x1D702},t)$ profile as

This particle-based description of fluid mixing is compatible with the CTRW stretching model. As such, we do not seek to directly solve the evolution of an entire lamella, but rather solve the evolving concentration profile $c(\unicode[STIX]{x1D702},t_{n})$ over a representative ensemble of points, which recovers the same statistics of the ensemble of all lamellae in the pore space due to ergodicity.

If we consider fluid advection and stretching as stochastic processes, then for given $z_{n}$ or $n$ the advection time $t_{n}=t(z_{n})$ is distributed randomly according to the waiting time distribution (4.18). Hence, $c(\unicode[STIX]{x1D702},t_{n})$ is a random variable through its dependence on both the random advection time $t_{n}$ and the random deformation $\unicode[STIX]{x1D6FE}_{n}=\unicode[STIX]{x1D6FE}(t_{n})$ . Under the assumption of independent lamellae, the PDF of concentration at a cross-section at $z_{n}$ is obtained by (i) sampling the concentration values across the lamella and (ii) sampling between the lamellae. This interpretation allows the lamellar mixing (6.7) to describe the evolution of the transverse concentration profile with the advection time and stretching history of any part of any lamella in the entire 3D porous network.

We non-dimensionalize (6.7) by introducing the variables $c^{\prime }=c/c_{m}(0)$ , $t^{\prime }=t/\unicode[STIX]{x0394}t_{a}$ , $\unicode[STIX]{x1D702}^{\prime }=\unicode[STIX]{x1D702}/(\sqrt{2}\unicode[STIX]{x1D70E}_{0})$ , with the characteristic advection time $\unicode[STIX]{x0394}t_{a}=\ell /v_{0}$ , where $\ell$ is the length of a pore and $v_{0}$ is the mean pore velocity. We also define the Péclet number as $Pe=\unicode[STIX]{x1D70E}_{0}^{2}v_{0}/(2D_{m}\ell )$ , and henceforth drop the primes for simplicity of notation. Under these scalings (6.7) is non-dimensionalized, with $D_{m}$ now replaced by $1/Pe$ , and the initial condition (6.1) is now $c(\unicode[STIX]{x1D702},0)=\exp (-\unicode[STIX]{x1D702}^{2})$ . To solve (6.7), following Ranz (Reference Ranz1979) we define the reduced coordinate $\unicode[STIX]{x1D702}_{0}$ and operational time $\unicode[STIX]{x1D70F}(t)$ respectively as

where the solution $c(\unicode[STIX]{x1D702},t)=g[\unicode[STIX]{x1D702}_{0}(t),\unicode[STIX]{x1D70F}(t)]$ is equivalent to a transformation into the characteristic system of (6.7). Making this substitution, equation (6.7) simplifies to the diffusion equation

with initial condition $g(\unicode[STIX]{x1D702}_{0},0)=\exp (-\unicode[STIX]{x1D702}_{0}^{2})$ . Hence, the concentration profile $c(\unicode[STIX]{x1D702},t)=g[\unicode[STIX]{x1D702}_{0}(t),\unicode[STIX]{x1D70F}(t)]$ across the strip is given by the Gaussian distribution

and so the concentration profile is determined by the total stretch $\unicode[STIX]{x1D70C}(t)$ which acts to narrow lamellae, and the stretching history $\unicode[STIX]{x1D70F}(t)$ which broadens and dilutes the lamellae. As $\unicode[STIX]{x1D70F}(t)$ depends upon the stretching history, stretching oscillations in the two-step process (5.1) over a couplet are captured by $\unicode[STIX]{x1D70F}(t)$ , and so must be accounted for if significant. From (6.10) the maximum concentration $c_{m}(t)=c(0,t)$ is now

Equations (6.8)–(6.11) show that the pore-space mixing explicitly depends upon $Pe$ and the stretching history $\unicode[STIX]{x1D70C}(t)$ which is used to encode the operational time $\unicode[STIX]{x1D70F}(t)$ . The fluid deformation CTRW model developed in § 4 yields expression (4.16) for the PDF of $\unicode[STIX]{x1D70C}(t_{n})$ , and so this must be extended to generate analytic expressions for the PDF of $\unicode[STIX]{x1D70F}(t_{n})$ , the concentration PDF and associated measures.

Equations (6.8)–(6.11) also form the basis of the diffusive strip method (Meunier & Villermaux Reference Meunier and Villermaux2010) which facilitates efficient numerical solution of the 1D lamellar ADE (6.7) across a broad range of $Pe$ . To test the accuracy of the analytic solutions for the concentration PDF described above, we employ this method in conjunction with the advective $\mathscr{M}$ and temporal $\mathscr{T}$ maps to simulate deformation of a 2D material filament within the 3D pore space. The diffusion is calculated as a post-processing step for various $Pe$ via (6.10) given computation of $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D70F}$ along the strip from the stretching history. This method is capable of accurately capturing the stretching, diffusion and coalescence of the scalar field up to $n\sim 100$ pores, beyond which the exponential growth of the lamellar structure is too large for feasible computation. Figure 6 shows a typical distribution of the operational time $\ln \unicode[STIX]{x1D70F}$ over the cross-section of a pore at $n=40$ using this method, as well as the concentration distribution $c(\boldsymbol{x},t_{n})$ calculated at $Pe=10^{8}$ .

Figure 6.

(a) Typical distribution of the logarithm of the operational time $\ln \unicode[STIX]{x1D70F}$ (black, $\ln \unicode[STIX]{x1D70F}=50$ ; white, $\ln \unicode[STIX]{x1D70F}=5$ ) representing the backbone of the lamellar concentration field, with detail (dashed box) shown in (b). The associated scalar field $c(\boldsymbol{x})$ for $Pe=10^{8}$ shown in (c) is calculated via the diffusive strip method (Meunier & Villermaux Reference Meunier and Villermaux2010).

In order to obtain a closed-form expression for the PDF of $\unicode[STIX]{x1D70F}(t_{n})$ from the CTRW model, we approximate the evolution of $\unicode[STIX]{x1D70C}(t)$ between $t=0$ and $t=t_{n}$ as a constant exponential stretching process,

in terms of the average stretching rate $\unicode[STIX]{x1D6E4}_{n}=\ln (\unicode[STIX]{x1D70C}_{n})/t_{n}$ after $n$ steps. From (6.12), the operational time $\unicode[STIX]{x1D70F}(t_{n})$ is then approximated as

This expression can be further simplified by noticing that the PDF (4.16) of $\ln (\unicode[STIX]{x1D70C}_{n})$ is Gaussian, with mean $n\unicode[STIX]{x1D6EC}_{\infty }/2$ and variance $n\unicode[STIX]{x1D70E}^{2}/2$ . Thus, we may approximate the denominator in (6.13) by its mean value $\ln (\unicode[STIX]{x1D70C}_{n})\approx n\unicode[STIX]{x1D6EC}_{\infty }/2$ ,

the accuracy of which is shown in figure 7(a). In essence, equation (6.14) solves the pore-scale fluid mixing, involving exponential growth of $\unicode[STIX]{x1D70F}_{n}$ with $n$ . Oscillations associated with the two-step stretching process (5.1) appear to play a minor role as $\unicode[STIX]{x1D70F}_{n}$ is well approximated by the one-step stretching process (5.2) characterized in terms of the mean stretching rate $\unicode[STIX]{x1D6EC}_{\infty }$ , as reflected by figure 7(a). Thus, while in principle the diffusion is dependent upon the entire stretching history, it is predominantly governed by the arrival time $t_{n}$ of the material strip at position $z_{n}$ and the elongation $\unicode[STIX]{x1D70C}_{n}$ . The travel time distribution over a single pore is given by the Pareto distribution (4.18). The distribution of the arrival times $t_{n}$ is given by the self-convolution of the Pareto distribution, which results in a Landau distribution. In appendix B we show explicitly that the bulk of the Landau distribution can be well approximated by the Moyal distribution.

Figure 7.

(a) Comparison of direct computation (black) and analytic approximation (6.14) (blue) for operational time $\unicode[STIX]{x1D70F}_{n}$ with pore number $n$ over four different realizations of the CTRW model. (b) Comparison of the PDF of $\ln \unicode[STIX]{x1D70F}$ based upon the analytic approximation (6.16) (solid lines), the PDF of $\ln \unicode[STIX]{x1D70F}$ from the CTRW model for (circles) $n=20$ , (pentagons) $n=40$ , (diamonds) $n=60$ , (squares) $n=80$ and (triangles) $n=160$ .

The PDF of $\ln (\unicode[STIX]{x1D70F}_{n})=\ln (t_{n})+2\ln (\unicode[STIX]{x1D70C}_{n})-\ln (n\unicode[STIX]{x1D6EC}_{\infty })$ is well approximated by replacing $\ln (t_{n})$ with its average,

where $\unicode[STIX]{x1D6FE}$ is the Euler–Mascheroni constant, $a=0.7413$ and $b=0.0064$ (appendix B). Due to the central limit theorem, $\ln (\unicode[STIX]{x1D70C}_{n})$ converges with $n$ towards the Gaussian distribution (4.16) with the mean $n\unicode[STIX]{x1D6EC}_{\infty }/2$ and variance $n\unicode[STIX]{x1D70E}^{2}/2$ , and so the PDF of $\unicode[STIX]{x1D709}_{n}:=\ln (\unicode[STIX]{x1D70F}_{n})$ also converges with $n$ towards the Gaussian distribution

with mean and variance

Figure 7(b) shows that for $\unicode[STIX]{x1D70F}>10^{3}$ the approximation (6.16) exhibits very good agreement with the PDF of $\unicode[STIX]{x1D709}_{n}$ obtained from evaluating $\unicode[STIX]{x1D70F}(t_{n})$ according to (6.8) with $\unicode[STIX]{x1D70C}(t_{n})$ given by the CTRW (5.1). The mismatch at short times is expected due to the asymptotic nature of the central limit theorem.

7 Chaotic mixing in 3D random porous media

The approximation (6.16) for $\unicode[STIX]{x1D70F}_{n}$ provides an accurate solution to the two-step CTRW (5.1), which along with the distribution (4.16) for $\unicode[STIX]{x1D70C}_{n}$ quantifies the evolution of the concentration distribution, mixing and dilution in the 3D open porous network prior to the onset of coalescence. In the following we apply this solution to determine the evolution of the mixing scale, concentration PDF, maximum concentration, scalar variance and the onset of coalescence in the 3D random porous network.