A Primer on the Dynamical Systems Approach to Transport in Porous Media

Since the pioneering work of Darcy (1856) regarding the quantitative relationship between fluid pressure gradient and flux through soil columns, interest in the fluid dynamics of porous media systems has grown rapidly. These dynamics govern the transport, reaction and biological activity of solutes, colloids and microorganisms in a wide range of porous materials ranging from geological media to biological tissues and engineered structures, and now underpin many vital industries, including biotechnology, energy, agriculture and water supply. All porous media are characterised by complex pore-scale architectures and many also exhibit significant heterogeneities from pore scales to much larger (e.g. regional) length scales, resulting in complex flow and transport dynamics not resolved in Darcy’s original work. Due to these inherent complexities, the conceptual progress of porous media science has not been smooth. Progress has relied on irregular advances in technology that have enabled greater precision and resolution in observation, which have then provided the data that seeded the next conceptual advance. Today, researchers have unparalleled access to pore-scale data of fluid flow, three-dimensional (3D) geometry and surface chemistry, yielding an almost overwhelmingly detailed and rich picture of the process complexity inherent in porous media. The next conceptual advance is imminent.

Figure 1 presents an example of the kinds of complex experimental data now available and awaiting interpretation. The figure shows the results of a standard push-pull experiment at low Reynolds number where the concentration distribution of a fully developed solute plume near the injection point is shown at different times $\tilde{t}$ during the pull phase. The width of the injected plume shown at $\tilde{t}=0$ has similar dimensions to those of the glass beads used to form the randomly packed porous medium. In the absence of molecular diffusion, reversibility of the push-pull flow would result in the collapse of the solute plume during the pull phase back to the original injection profile shown at $\tilde{t}=0$. Hence the solute plume shown in Figure 1 represents a residual “diffusive signature” resulting from the irreversible interactions between fluid deformation history (noting that net fluid deformation is zero) and solute diffusion. After completion of the pull phase, at $\tilde{t}=25$, this diffusive signature is observed as an unexpectedly extensive and striated distribution. This distribution indicates the presence of profound anisotropic fluid transport phenomena at the pore scale which are unseen in less-resolved investigations, and cannot be described via a conventional solute dispersion framework. Novel tools and techniques are required to understand, characterise and quantify such transport phenomena.

At the same time, it has also been discovered that quite unexpected transport phenomena can present in naturally occurring flows at the Darcy scale, i.e. well beyond the pore scale. Using a conventional poroelastic Darcian flow model, Wu et al. (2019) showed that simply by increasing the compressibility $\mathcal{C}$ of the porous medium a range of transport structures of increasing complexity can be generated. Figure 2 depicts the motion of tracer particles (black points and lines) at integer multiples of the tidal forcing period $P$ in a two-dimensional (2D) model of a compressible heterogeneous aquifer subject to both tidal boundary forcing and a background regional flow. As shown in Figure 2a, for small values of $\mathcal{C}$ tracer particles sampled at integer multiples of $P$ move along curved paths from inland (right) toward the coastal boundary (left). With increasing compressibility (Figure 2b, c), elliptic (E) and hyperbolic (H) periodic points arise in the aquifer. The local fluid motion (over one flow period) around the elliptic point (E) is rotational, while fluid stretching and compression (indicated by the stable (red) and unstable (blue) trajectories) occurs local to the hyperbolic point (H). Associated with this bifurcation in the transport structure of the flow is the appearance of a closed region (where non-diffusive tracer particles are trapped indefinitely) around the elliptic point, indicating a qualitative change in the aquifer transport dynamics. Above a critical compressibility (Figure 2d), the hyperbolic point (H) and the associated open region bifurcates into a chaotic region (indicated by the scattered points) and a regular region, while the closed region persists as a stable “island” with the elliptic point (E) at its centre. This chaotic region exhibits complex solute transport dynamics, including stochastic residence time signatures and accelerated mixing. The above transport structures are not easy to predict, detect or understand using standard techniques, yet they exert significant control on solute transport.

From the above examples it is clear that conventional Darcy- and pore-scale tools and techniques (that are so well suited to many practical problems) are unlikely to be of much use in understanding and quantifying the transport and mixing dynamics of these more complex flows. A conceptual advance is required to push the science further, and this requires consideration of the Lagrangian kinematics of these flows to elucidate how fluid particle trajectories are organised by the underlying transport structure. This underlying transport structure governs processes ranging from solute mixing and pollutant dispersion, to chemical reactions and biological activity. The natural language of these kinematics is couched in terms of dynamical systems (chaos) theory, and in this paper we introduce the key concepts in this field required to understand solute transport in porous media using a dynamical systems approach. It is important to note that the dynamical systems approach to understanding the Lagrangian kinematics of porous media flows is not limited to understanding chaotic particle trajectories (i.e. “chaotic advection”), but, as shall be shown, this approach is also useful in understanding transport at a deeper level in regular (non-chaotic) systems, or systems such as shown in Figure 2(d) that exhibit both regular and chaotic regions. This primer on the dynamical systems approach is deliberately pitched to provide an overview of the field rather than a detailed exposition, and the interested reader is directed to the referenced works for greater detail. As we will see in the following sections, dynamical systems theory brings with it a rich set of concepts and tools relevant to modern porous media flow, and it already has a (relatively short) history in the subject.

Taking a broad (and deliberately highly selective) overview of porous media research, we may form a timeline of key concepts and breakthroughs in transport in porous media, depicted in Figure 3. The timeline shows the recent application of dynamical systems concepts to the field, in green. The dynamical systems entries are disproportionately over-represented in the timeline as each of the other research themes have accrued many thousands of articles over the decades while there are still only a handful of dynamical systems articles in the porous media literature. However we contend that the new tools and techniques provided by dynamical systems theory may facilitate rapid advances in our understanding of flow and transport in porous systems, for the following reasons.

It is now well recognized that several classes of porous media flows give rise to transport phenomena that are more complex than steady 2D flow, and so cannot be fully understood via conventional analysis (Weeks and Sposito, 1998). Select examples of complex steady Darcy flows include helical motion in steady anisotropic 3D Darcy flow (Bakker and Hemker, 2004; Cirpka et al., 2012, 2015; Chiogna et al., 2014; Ye et al., 2015), steady 3D Darcy flow in an otherwise homogeneous porous medium with spheroidal inclusions (Dato et al., 2016; Di Dato et al., 2016; Janković et al., 2003, 2009), while complex unsteady Darcy flows have been observed in both homogeneous (Bagtzoglou and Oates, 2007; Metcalfe et al., 2010b; Trefry et al., 2012b; Mays and Neupauer, 2012; Piscopo et al., 2013; Sather et al., 2022) and hetereogeneous (Trefry et al., 2019; Cho et al., 2014) media. Similarly, complex steady 3D flows have been shown to arise at the pore-scale in both open porous networks (Lester et al., 2013a, 2014b) and granular media (Kree and Villermaux, 2017; Turuban et al., 2019; Souzy et al., 2020; Heyman et al., 2020). Although the complexity of these flows arises via different mechanisms—ranging from anisotropy of the conductivity field to transient forcing to complexity of the pore-scale architecture—these flows are unified in that their fluid particle trajectories are organised in a more complicated manner than the lamellar streamlines of all steady 2D flows, and so the Lagrangian kinematics of these flows cannot be understood solely from the Eulerian velocity field. Similarly, solute mixing in all of these flows is characterised by highly striated, ramified solute distributions that cannot be fully understood from the Eulerian velocity field.

Although the Lagrangian frame has an extensive history of application in the study of transport in porous media, use of Lagrangian coordinates alone is insufficient to understand transport. Indeed, several of the studies outlined above have not explicitly resolved the Lagrangian kinematics of these flows, which is concerned with the organisation of fluid trajectories and the attendant deformation of fluid elements (see Trefry et al. (2019) for an example of organised kinematics). In this primer, we advocate for a Lagrangian coordinate description of particle trajectories, but the important next step to understand is that particle trajectories (because they form a dynamical system) organize into coherent structures that can be understood by analyzing special sets of trajectories using tools developed in dynamical systems theory. Hence the dynamical systems theory approach to transport in porous media provides a convenient framework for detecting and classifying kinematic features in the Lagrangian frame that can control fluid mixing, segregation and discharge in porous media. These kinds of transport phenomena describe complex Lagrangian structures (termed Lagrangian Coherent Structures (LCSs)) which engender potentially profound impacts on solute migration and reaction (Toroczkai et al., 1998; Tél et al., 2005; Valocchi et al., 2019). Coupled with appropriate stochastic methods (such as random walks and non-local transport theories), the Lagrangian kinematics can be used to develop quantitative predictions of solute transport, mixing and reaction that resolve the complex dynamics observed in Figures 1 and 2.

One the first authors to use a dynamical systems approach to understand flow and transport in porous media was Sposito (1994), who recognised that steady Darcy flow in heterogeneous, isotropic porous media imposes strong constraints on the Lagrangian kinematics, precluding chaotic motion. Although this class of porous medium flow is not complex in the manner described above, only explicit consideration of the Lagrangian kinematics uncovers these properties. Conversely, early conceptual work by Ottino (1989); Jones and Aref (1988), Sposito (2001, 2006) and others (Lester et al., 2009, 2010b; Metcalfe et al., 2010b; Trefry et al., 2012b; Mays and Neupauer, 2012; Lester et al., 2016a) has led to the prediction, observation and engineering of complex Lagrangian kinematics in saturated porous media at the Darcy (Zhang et al., 2009; Metcalfe et al., 2010a) and field (Cho et al., 2019) scales. Complex Lagrangian kinematics have also been predicted to occur in natural environments (Trefry et al., 2019; Wu et al., 2019), e.g. compressible Darcian systems where engineered pumping activity is absent. Dynamical systems concepts have also been used to link the topology of pore-scale architectures with pore-scale solute mixing and macroscopic transport phenomena (Lester et al., 2013c, 2014b, 2016b), resulting in fundamental upscaling behaviours conditioned on the details of pore-scale Lagrangian kinematics. These mechanisms have recently been confirmed via direct experimental observations of pore-scale mixing (Souzy et al., 2020; Heyman et al., 2021, 2020), leading to the highly striated solute distributions such as those shown in Figure 1. Dynamical systems approaches also shed new light on other fluid-borne phenomena (Tél et al., 2005) and thus Lagrangian kinematics underpin the transport, spreading, clustering and reactivity of solutes, colloids and microorganisms in porous media. Hence it is natural to seek new learnings from the dynamical systems approach to transport in porous media.

But what is dynamical systems theory, what can it tell us about the fundamental characteristics of flow, and how can it be applied to porous media systems? We attempt to answer these questions in the remainder of this paper. Our approach is to provide clear and concise definitions of key dynamical terms and concepts, illustrated by examples from the recent research literature, so that the reader is armed with the basic tools and knowledge to consider future applications of this powerful theory in the context of porous media.

The dynamical systems approach to fluid transport, often termed “chaotic advection” (although the concepts may also be applied to non-chaotic systems), is the science at the intersection of fluid mechanics and nonlinear dynamical systems (Aref et al., 2017). In this case the phase space (coordinates of interest) of the dynamical system is the physical space that the fluid occupies, whether at the pore scale or the Darcy scale, and either in a two-dimensional approximation or the full three-dimensional space.

Hence chaotic advection is concerned with the trajectories of the dynamical system given by the collective motion of all fluid particles (where “particle” refers to a conceptual massless tracer particle) given by the Lagrangian kinematics of the flow. The main point in this section is that particular structures in the flow — collections of points, lines, surfaces — organize the Lagrangian transport of the entire flow, and transport can only be truly understood by uncovering the Lagrangian kinematics.

These Lagrangian coherent structures (LCS), such as those shown in Figure 2, form the “skeleton” of the flow that organises global transport (MacKay, 1994), even in the presence of random processes such as diffusion. And to understand or design transport in porous media flows (be it of heat, of chemical or biological species, of matter, in fact, of all interactions mediated by the flow) it is important and useful to find, classify, and elucidate these LCS. Broadly interpreted these Lagrangian structures consist of attracting regions, repelling regions, hold up regions, and the distinguished pathways that connect one critical region to another.

In fluids, including porous media flows from pore to regional scales, the Lagrangian trajectories are the trajectories of fluid particles. From its initial location the trajectory of every fluid particle is given by the kinematic equation

here written in 2D with the fluid particle location $(x,y)$ moved passively by the fluid velocity vector field $(V_{x},V_{y})$. Deceptively simple, the kinematic equation can be taken as the elementary definition of velocity or as defining a nonlinear dynamical system in which a given velocity field generates the Lagrangian trajectories of the fluid tracer particles. Hence these trajectories provide a coordinate transform between the Eulerian ($x,y$) and Lagrangian ($x_{0}$, $y_{0}$) frames. In keeping with the idea of building from simplest concepts first, the kinematic equation (1) has no diffusion or local dispersion.

Throughout this primer we focus mainly on the Lagrangian kinematics (which are purely advective), but we note that these kinematics govern fluid deformation (i.e. the stretching, shearing, and folding of fluid elements) (Lester et al., 2018), and solute transport in the presence of diffusion or local dispersion is controlled by the deformation history of fluid elements (Villermaux, 2019). Hence such additional physics can be readily incorporated, but knowledge and quantification of the Lagrangian kinematics is essential to correctly condition the random fluid stretching processes that govern the advection-diffusion process.

Equation (1) is deceptive in that it can generate complex trajectories even for simple velocity fields. An example will illustrate a simple case and expose selected hidden LCS. Referring to Figure 4, our example is a 2D, time-periodic Stokes flow in a box of fluid. This is amongst the simplest settings in which to observe LCS. One sidewall steadily slides for half the period and drives a simple vortex flow as in Figure 4(a); during the other half period, the opposing wall steadily slides in the opposite direction (Figure 4(b)). A video of a similar experiment is also available (Jana et al., 2015).

In Figure 4(c) two blobs of fluid are marked red or blue, and Figures 4(d)–4(f) show the advection-only evolution of the blobs after the indicated number of periods. The evolution of the two blobs is markedly different because their deformation is being controlled by different types of LCS. The blue blob was deliberately put into a hold-up region and undergoes a shear-like deformation that yields (at most) linear stretching, and every three periods it returns close to its initial condition. In 2D, time-periodic flows material starting in a hold-up region is truly cut off from the rest of the flow; these regions are called elliptic islands because the material inside the island rotates around an elliptic point (as exhibited in Figure 2(c,d)), and the net deformation around these periodic points is a simple rotation. In 3D or non-periodic flows hold-up regions may not be completely cut off from the rest of the flow, but are still delineated by surfaces of minimal flux across them that keeps material close for long periods of time.

The red blob, on the other hand, is rapidly deformed into a highly ramified and extended filament. The red blob was deliberately placed near a hyperbolic point (shown in Figure 4(g)), where the net fluid deformation over one flow period consists of repelling (stretching) in two directions and attracting (compressing) in the two orthogonal directions (as shown in Figure 2(c,d))). In this simple example these hyperbolic points encapsulate both the attracting and repelling LCS; moreover, in this case the distinguished pathway between attracting and repelling regions begins and ends at the same point, which is denoted a homoclinic connection. In general in 2D periodic flows heteroclinic connections, where the repelling and attracting points are different, are typical. Here the hyperbolic point and its homoclinic connection strongly deforms the red blob through repeated stretching and folding. This in turn generates exponential growth of the red blob’s interfacial area, even though the blob remains confined inside a finite volume of fluid.

That is why the shape of the red blob becomes more and more ramified with smaller and smaller intervals of unmarked fluid between the red filaments. Even for fluid particles that are initially arbitrarily close together it cannot be predicted how far apart they will end up. This destruction of long term predictability of the relative location of fluid particles is why such flows are called chaotic, even though the underlying velocity field is deterministic (Ottino, 1989).

That the chaotic part of a flow generates small-scale structure in the spatial distribution of material (the red blob) is why chaotic advection is associated with mixing. The chaotic part of the flow quickly evolves any smooth initial distribution into a complex pattern of filaments or sheets, depending on the dimensionality of the system, which converges at an exponential rate to a fractal structure. When the characteristic length scale of the blob becomes small enough, this structure is then smoothed by diffusion, manifesting as accelerated mixing. This rapid reduction of length scale and these augmented transport characteristics are purely effects of the simple kinematic equation (1).

Figure 4(g) illustrates several ways to bring out the Lagrangian structure that are commonly used in the dynamical systems literature. The black dots are a “stroboscopic map” of tracer particles similar to that used in Figure 2. Take a small set of initial conditions (a line of particles at $x=0.5$ in Figure 4(g)), and for each fluid particle evolve its trajectory and place a dot to mark the particle’s location after every period. Continue this process for many periods and the resulting collection of dots defines a stroboscopic map that visualises the long term Lagrangian structure. For instance, in Figure 4(g) the stroboscopic map clearly partitions the fluid region into three elliptical island regions and a chaotic “sea”, the dense jumble of particle locations seemingly without structure. This common technique is also often called a Poincaré section (Aref et al., 2017).

The blue and red lines in Figure 4(g) are structures termed respectively stable manifolds and unstable manifolds, which are associated with the central hyperbolic periodic point (Ottino, 1989). The stable manifold is the most repelling material curve, while the unstable manifold is the most attracting material curve. While the stable and unstable manifolds are time inverses of each other, we focus on the red unstable manifold and note how much it resembles the deformed red blob, particularly the deformed red blob after 6 periods (Figure 4(f)). This is no coincidence. The unstable manifold is a fractal structure that fills the chaotic sea, and material starting near the central hyperbolic point is stretched along the unstable manifold and eventually converges to the unstable manifold. The unstable manifolds of other hyperbolic points, if they exist, interweave together to form the chaotic seas. The unstable manifold is the fundamental mixing pattern (or “mixing template”, see Aref et al. (2017)). Technically to generate chaotic advection requires transverse intersection of the unstable and stable manifolds, and manifold interaction is in fact the generic mechanism behind—and topological condition for—chaos in both steady and time-periodic flows.

So far, we have given a very brief overview through one particular example of how Lagrangian critical regions, consisting of attracting, repelling, and hold-up regions along with the distinguished connections between critical regions, organize the Lagrangian transport of passive fluid particles. An immediate question is then where does this kind of transport occur in porous media flows? The answer is wherever there are sufficient degrees of freedom (DOFs) in the fluid velocity field, be that as an upscaled continuum flow or in pore spaces, chaotic advection is the norm rather than the exception. The sufficient DOFs for chaos in continuous systems is three, meaning that 3D steady or 2D unsteady flows can possibly exhibit chaotic advection (Speetjens et al., 2021), whereas steady 2D flows are intrinsically non-chaotic. Thus we must reconsider if steady 2D flow approximations can adequately represent many porous media flows at both the pore and Darcy scales.

Consider first the Stokes flow in pore spaces discussed in section 3. Here the 3D pore spaces intrinsically have many hyperbolic points in the skin friction field of the solid matrix surface (Lester et al., 2013c) that generically generate chaotic fluid trajectories. Suggestive of this generic chaos is the similarity of the stretched and folded solute plume in the experiment of figure 7(c) with the stretched and folded filaments of Figure 4. However, in 3D periodic or non-periodic flows not all of the organizing LCSs are as well characterised mathematically or experimentally as they are in 2D; there is a richer set of of Lagrangian structures and connections between them that one can investigate and use (Haller, 2015; Speetjens et al., 2021).

These mechanisms also apply to evolution of the diffusive dye blob shown in Figure 1, which is most clearly elucidated by first considering the push phase through the bead pack as a complex steady 3D flow. Unlike the time-periodic 2D flow considered in Figure 4, there do not exist periodic points and lines in this flow, but there still exist stable and unstable manifolds (attracting and repelling LCS) that organise transport and impart exponential stretching of fluid elements. The diffusive dye shadows the unstable manifold of the push flow, and diffuses transverse to this structure in a filamentous manner whilst being stretched exponentially along the backbone of the unstable manifold. Eventually these filaments merge together via diffusion, leading to a well-mixed state. However, for the push-pull experiment, the flow is reversed before this state is reached. Due to reversal the stable and unstable manifolds are exchanged, and the dye blob now evolves along the stable manifolds of the push flow, leading to the highly striated structures shown in Figure 1 for different reversal times $\bar{t}$. While the original dye blob would be recovered upon flow reversal for non-diffusive dyes, for diffusive dyes, the irreversible nature of diffusion means that the resultant dye structure contains artifacts from the stretching history of the blob, and so provides a diffusive snapshot of the stable manifolds of the push flow.

Similarly, the Poincare sections shown in Figure 2 clearly illustrate the LCS of the time-periodic poroelastic Darcy flow, where elliptic and hyperbolic points, elliptic islands and stable and unstable manifolds are shown. In this case, Lagrangian transport is trivial for small values of the medium compressibility $\mathcal{C}$ (Figure 2(a)), but with increasing compressibility a topological bifurcation occurs, leading to formation of elliptic and hyperbolic points (Figure 2(b,c)). Under a further increase of the medium compressibility (Figure 2(d)), the stable and unstable manifolds form a heteroclinic connection (not shown), and a chaotic sea forms around the elliptic island. These LCS govern the transport of both non-diffusive and diffusive species in this time-periodic poroelastic Darcy flow, and can lead to trapping of solutes, regions of accelerated mixing and transport that cannot be observed without elucidation of the LCS. In general, the tools and techniques associated with chaotic advection are necessary to understand solute transport in a wide range of porous media applications.

For Darcy type flows Figure 5 is a guide. Steady Darcy flows, even in 3D, are non-chaotic unless the conductivity field is sufficiently complex. The constraint that causes this behaviour is discussed in section 4. Unsteady Darcy flows in either 2D or 3D can always be made chaotic if one has control over the time-dependent flow forcing. In natural Darcy flows, time-dependent forcing can unleash chaotic fluid trajectories, but as illustrated in Figure 2, the transition to chaotic behaviour is dependent upon the details of the system at hand.

To end this section we point out one more feature of Figure 4, our simple prototype chaotic flow: the typical Lagrangian picture consists of both hold-up regions (fluid somewhat or completely isolated from the rest of the domain) and well-mixed regions (generated by the manifold interactions of the attracting and repelling regions). Applications of porous media flows often require maximizing one or the other, i.e. isolating fluid or mixing it. It is not always easy to find the particular flow forcing sequence that will maximize fluid isolation or mixing. However, it can be done, and Figure 12 shows an experiment designed to confine fluid in a Darcy flow, while Figure 13 shows a field experiment designed to maximize mixing of emplaced fluid. Developing general tools to take a porous media flow and quickly tailor it for a particular purpose is an area of ongoing research.

The classical picture of pore-scale fluid flow and solute dispersion follows a steady 2D paradigm as shown in the top part of Figure 6, which is reproduced from Bear and Verruijt (1987). This 2D conceptual image illustrates the main features that govern solute transport, including mechanical dispersion due to (a) strong flow velocity gradients between grains, (b) flow separation around grains and (c) molecular diffusion of the solute. Although these solute transport mechanisms have been elaborated significantly from many contributions over the past few decades, it is sometimes overlooked that the 2D nature of this conceptual image itself imposes severe constraints on the admissible fluid dynamics, which in turn may limit our mental picture of solute transport. Although conceptual, these 2D topological constraints implicitly skew our understanding of the true 3D mechanisms which drive solute transport and sometimes lead us to neglect other important phenomena which may occur in 3D domains.

An important topological constraint that applies to all steady 2D flows is that streamlines cannot cross over each other, and so separating streamlines that pass either side of a grain must eventually recombine. Hence, the separation distance between neighbouring streamlines may only fluctuate as they move downstream; they cannot grow or shrink without bound. In mathematical terms, this topological constraint means that fluid elements in steady 2D flow may only grow at most algebraically with time (Dentz et al., 2016) (due to e.g. shear between grains), which in turn places limitations on the rates of solute dispersion and mixing in such flows (Le Borgne et al., 2015). This constrained streamline behaviour in steady 2D flow also limits transverse dispersion as molecular diffusion is the only mechanism by which solutes can spread transversely to the mean flow direction. Conceptually, and in the absence of molecular diffusion, solute molecules simply follow initial streamlines and so cannot spread transversely, hence mechanical transverse dispersion is zero in 2D. Experimental observation of this 2D splitting-recombination phenomenon is shown in the bottom part of Figure 6, where the solute plume in the 2D porous medium splits and recombines several times, yet the net transverse dispersion is limited and is governed by the solute molecular diffusivity.

Although it is tempting to conceptualise steady flow through 3D porous architectures as an extrusion of the dynamics in Figure 6, these topological constraints do not apply to 3D steady flows as the additional degree of freedom—the third spatial dimension—removes the restriction that a fluid streamline is bounded by its neighbours. The result is that the streamline may wander freely throughout the fluid domain, as shown in Figure 7(a): in steady 3D flows neighbouring streamlines may diverge without bound as they are advected downstream, leading to wholly new transport mechanisms (such as chaotic advection) which are not readily apparent in conventional dispersion models. This freedom also admits braiding of streamlines, depicted in Figure 7(a), that is associated with exponential growth of transverse material lines.

The origin of chaotic advection in 3D pore-scale flows can be traced back to stagnation points that arise on grains or at pore boundaries in all porous media (Lester et al., 2013b), and fluid deformation local to these points generates persistent exponential stretching of fluid elements in the bulk. Conversely, in steady 2D flows, such local fluid stretching at stagnation points is cancelled out at downstream reattachment points. The persistent exponential stretching in 3D leads to the highly striated and filamentous solute distributions transverse to the mean flow direction shown in Figure 7(c). Flow separation and persistent exponential stretching of material elements is wholly responsible for the complex steady solute plume structures observed in 3D pore-scale flows. As these dynamics are inadmissible in steady 2D flow, they can be overlooked in conventional modelling approaches.

It is important to note that exponential fluid stretching qualitatively augments solute mixing, as these mixing dynamics are singular in the presence of exponential fluid stretching (Cerbelli et al., 2017) in that the rate of mixing in the limit of vanishing diffusivity is not the same as that for purely advective flow. Conversely, sub-exponential stretching (which is characteristic of regular fluid flow) exhibits non-singular mixing rates which limit to zero with vanishing molecular diffusivity. Hence, the impact of chaotic advection on solute mixing is profound.

The impact of pore-scale chaotic advection extends beyond solute mixing, and also acts to augment longitudinal dispersion (Lester et al., 2014a, c), chemical and biological activity (Tél et al., 2005) and transport of colloids and micro-organisms (Haller, 2011). As the Péclet number $Pe$ (which characterises the relative timescales of advection to diffusion) is typically small in pore-scale flows, it is commonly assumed that diffusion is rapid compared to advective transport and so chaotic advection has minimal impact. However, Heyman et al. (2020) showed that exponential fluid stretching results in incomplete solute mixing at the pore scale even for small Péclet numbers ($Pe>5$): chaotic advection acts to augment diffusive mixing and transport across a wide range of pore-scale systems.

Heyman et al. (2021) has recently developed a novel method to characterise pore-scale chaotic mixing in opaque porous media based on the push-pull flow shown in Figure 1. This method has been used to recover the extent of chaotic mixing in model systems such as glass beads, as well as to estimate the extent of chaotic mixing in natural systems such as randomly packed gravel. This technique has the potential to characterise chaotic mixing across a broad range of natural porous media, ranging from granular matter to biological pore networks. Biological tissues may be regarded as porous networks that can host chaotic mixing, with applications in a range of contexts including living brain tissue (Reichold et al., 2009) (Figure 8a), vascular networks (Ghaffarizadeh et al., 2015) (Figure 8b) and lung alveoli (Tsuda et al., 2002, 2011). Conversely, Turuban et al. (2018, 2019) show that highly ordered and engineered porous materials such as crystalline lattices of uniform spheres can exhibit a broad range of pore-scale mixing dynamics ranging from non-chaotic to chaotic mixing that is stronger than that exhibited by any known random system. These mixing dynamics strongly depend upon the orientation of the crystalline axes with respect to the mean flow direction. Zhao (2012) have used disordered engineered porous metal foams (Figure 9) to develop compact heat exchangers that accelerate heat transfer by generating chaotic mixing in the pore space. These advances illustrate the potential for chaotic mixing principles to assist the development of novel porous materials exhibiting tuneable mixing and transport properties, with myriad industrial and medical applications.

In previous sections we have seen how dynamical systems approaches allow us to detect, classify and quantify the various kinematic structures inherent in a wide array of physical/biological systems and across dimensionalities. These outcomes are not only conceptual — they also provide us with information that can be applied to old and new problems for benefit. Here we discuss in greater depth two examples where the study of Lagrangian kinematics can led to advances in understanding and quantification of fundamental transport and mixing processes, and optimisation of engineering design and performance for porous systems.

Complex transport dynamics (such as chaotic advection) abounds in porous media flow. It is inherent to almost all pore-scale flows, and arises in many scenarios at the Darcy and regional scales. Complex Lagrangian kinematics differ from conventional mixing dynamics in that fluid deformation may be exponential in time, and the underlying dynamics often cannot be fully understood in the Eulerian frame. More than just a quantitative change, exponential fluid stretching profoundly augments solute mixing and transport, chemical reactions and biological activity in porous media. Specifically, unlike sub-exponential fluid stretching, exponential fluid stretching leads to singular solute mixing (Cerbelli et al., 2017) in that the limit of vanishing diffusivity does not correspond (qualitatively or quantitatively) to the zero diffusivity case. Hence the impacts of chaotic advection are typically not captured by conventional methods of modelling transport in porous media, and so unfamiliar tools and techniques are required to visualise, understand, quantify and, in engineering applications, ultimately control and exploit this phenomenon.

Hitherto the dynamical systems approach has not been widely appreciated within the community of practitioners of flow and transport in porous media, and this article unabashedly hopes to address that. This is not to say that chaotic advection, or more precisely transport mediated by Lagrangian coherent structures, is a magic bullet that solves all outstanding problems the field. On the contrary, we merely state that students and practitioners should be aware of this transport modality when circumstances warrant, and include appropriate techniques amongst the arsenal of conventional methods (stochastic methods, upscaling techniques, uncertainty analysis) to understand and analyse transport in porous media. However, these circumstances are not rare, rather they are common in porous media flows. As chaotic dynamics are the norm rather than exception in classical systems, chaos arises in most situations except for those where kinematic or topological constraints explicitly forbid it. For historical reasons the mental picture cultivated for porous media transport excluded complex Lagrangian kinematics, but now the dynamical systems approach to transport should be in your toolkit.

One important conceptual tool is the advective template of the flow which is given by the evolution of non-diffusive fluid tracer particles, and so explicitly resolves the system Lagrangian kinematics, associated coherent structures and topology. For any particular application of porous media flow, this advective template necessarily provides an incomplete description of the system as almost always there will be additional physical phenomena such as dispersion, diffusion, chemical reactions and biological activity, all of which interact nonlinearly with the advective template. We argue, however, that for every application knowledge of the underlying advective template is in and of itself useful and informative as this forms the organising structure upon which all additional fluid-borne phenomena play out. It is for this reason that we have focused this primer mainly on the advective template, yet nearly everywhere else in the porous media flow literature it is neglected. We note that many additional tools and techniques applicable to porous media can be found in the dynamical systems (chaos) literature in general and chaotic advection in particular.

Although deliberately limited in scope, several key messages are conveyed in this primer. First, chaotic advection is common to many porous media applications across all length scales. Second, these complex Lagrangian kinematics qualitatively augment associated fluid-borne phenomena, neither of which can be fully understood without utilisation of appropriate tools and techniques. Third, there exist easily identifiable kinematic constraints in particular applications that explicitly preclude chaotic advection, and chaos is expected when these constraints are absent. Fourth, fundamental studies of chaotic advection point toward development of novel classes of porous hierarchical materials (from pore to Darcy scales) with tuneable transport properties. Fifth, saturated porous media may now be regarded as a reaction vessel whose mixing state can be manipulated to promote reactivity or segregation according to purpose. Finally, many natural systems exhibit chaotic advection and complex Lagrangian transport due to the unsteady forcings inherent to these systems.

A new way of looking at a problem opens up myriad new and interesting research questions. Key knowledge gaps in the field include: (i) development and application of methods for characterisation of chaotic advection in all real (i.e. not model) natural and engineered porous materials, (ii) how the kinematics of steady 3D flow at the pore-scale influence macroscopic fluid mixing and chemical transport in saturated porous media, (iii) how to develop upscaled models of dispersion and dilution which capture these processes, (iv) design rules for engineered porous materials with tuneable transport properties, (v) development of topologically consistent upscaling methods for Darcy flows that do not introduce spurious kinematics, (vi) merging of chaotic advection methods with stochastic methods and conventional approaches to transport in porous media. This topological approach to fluid physics at the pore scale is a relatively new area and opportunities abound for further development. Another knowledge gap is how these concepts may be extended to real porous architectures via the development of stochastic models for the evolution of the fluid deformation tensor conditioned on tomographic imaging of the pore space and computational fluid dynamics simulations. Finally, more experiments are needed looking for the effects of chaotic advection across porous media applications.

Our goal with this paper is to implant in the reader the chaotic way of thinking and to generate some familiarity with the requisite mathematical apparatus. By this route we hope to open up new ways of designing porous media reaction technology and new understandings of ecological and biological porous transport. Let us know if we succeed.

The authors declare that no funds, grants, or other support were received during the preparation of this manuscript. The authors have no relevant financial or non-financial interests to disclose. The manuscript was written and drafted by all authors with equivalent contributions. All authors read and approved the final manuscript. No data sets were generated or analysed during the current study.