## Abstract

In nature, dissipative fluxes of fluid, heat and/or reacting species couple to each other and may also couple to deformation of a surrounding porous matrix. We use the well-known analogy of Hele–Shaw flow to Darcy flow to make a model porous medium with porosity proportional to local cell height. Time- and space-varying fluid injection from multiple source/sink wells lets us create many different kinds of chaotic flows and chemical concentration patterns. Results of an initial time-dependent potential flow model illustrate that this is a partially open flow, in which parts of the material transported by the flow remain in the cell forever and parts pass through with residence time and exit time distributions that have self-similar features in the control parameter space of the stirring. We derive analytically the existence boundary in stirring control parameter space between where isolated fluid regions can and cannot remain forever in the open flow. Experiments confirm the predictions.

## 1. Introduction

The effects of subsurface chemical and biological transport are difficult to measure, model and interpret. Reasons include inherent inaccessibility, a vast span of time-scales—from microseconds of chemical reaction to millenia of matrix deformation—an irregular and unknown pore geometry spanning lengthscales of micrometres to kilometres and the fact that fluid transport itself—through reaction, erosion, deposition and biological growth—continually changes the pore structure. Indeed three fields—fluid velocity, chemical/biological species concentration and matrix deformation—are coupled and coevolve, yet this coevolution is rarely, if ever, investigated. We want to create a model that couples velocity, concentration and deformation fields, in order to study their coevolution, yet is still amenable to experiment, visualization, computation and theory. This paper describes some first steps towards that goal by showing results from a model and experiment of advection in stirred porous media.

But let us be clear. This work makes *no* claim to model in detail any geological situation. Simulating this behaviour in all its complex splendor is not possible, laboratory realizations even less so. On the other hand, it is easy to get lost in the complexities of individual cases and to forget that there might be underlying, unifying structure. We aim first to strip away all but the essentials. The experiments and models described will not, and are not meant to, reproduce or predict transport during any particular geological process. Rather they seek to build intuition of what is possible or likely, elucidate fundamental physical mechanisms, constrain eventual more elaborate computations, suggest improved field measurement strategies, suggest strategies for improving geoengineered transport and aid in interpreting field measurements.

Our jumping off point is the well-established equivalency between Darcy flow in a porous medium and Hele–Shaw (HS) flow in a thin gap between rigid plates (Homsy 1987). The mathematical connection is that both are potential flows where fluid velocity *u* is driven by a pressure gradient such that the pressure *P* acts as a velocity potential. The permeability *K* in the Darcy equation is equivalent to the square of the gap distance *b* in the HS equation. The, to our knowledge, new idea we propose and that we are working towards is to combine (i) a stirred HS cell with (ii) one metal plate and one glass plate and (iii) a corrosive fluid, say an acid. The stirred fluid velocity field distributes the acid concentration in a pattern determined by the advection–diffusion solution of the flow for a given ratio of advection to diffusion time-scales (Rothstein *et al.* 1999; Lester *et al.* 2008). The acid etches the metal surface of the HS cell changing the height locally. The cell height distribution—equivalent to a permeability distribution—feeds back to the velocity field. In this way, velocity, concentration and permeability fields are coupled and coevolve. In this paper, only advection in the stirred HS flow will be considered.

To motivate briefly our specific experimental and theoretical model consider figure 1*a*, which shows flow through rock from an injection well to an extraction well (Prommer & Stuyfzand 2005). This could just as well have been a region of natural recharge and discharge (Anderson 2005). We will make several abstractions to project natural field flows into a laboratory setting. Our first abstraction is to notice that figure 1*a* suggests a dipole, the flow between a source and a sink of fluid. A mathematical model of a dipole flow is
1.1where *F* is the complex velocity potential for flow in the complex plane between source and sink singularities on the real axis at ±1. The real part of ∇*F* is the velocity potential and the imaginary part is the streamfunction (Lamb 1932). Figure 1*b* shows contours of the dipole velocity potential in colour and of the streamfunction as dark lines. The next abstraction, leading to accessible experiments, takes the Schwarz–Christoffel transformation
1.2to map the upper half of the complex plane, the singularities and the potential of equation (1.1) onto the disc of unit radius (figure 1*c*). The lower half plane is mapped to the outside of the disc, but that need not concern us except to note that the disc boundary is a separating streamline; it is not rigid. Figure 1*d* shows the dipole vector velocity field in the disc.

The last abstraction adds time dependence to the velocity field. Little work exists on mixing and transport in potential flows. Jones & Aref (1988) made an early study of a pulsed source–sink flow in the open plane. Their motivation was mainly to contrast with the (still) dominant choice of viscous flows, in order to demonstrate the purely kinematic nature of chaotic advection. Potential flow mixing was revived in one of the earliest designs of a micromixer (Evans *et al.* 1997) and later in a pulsed flow design of a DNA chip with improved hybridization time (Stremler & Cola 2006). A completely steady flow would be inadequate to capture the effects on transport of the natural time dependence of volume flow rate and direction. In order to expand the model, we add stirring controlled by two parameters for variety and flexibility. Call the dipole of equation (1.1)—or any multi-well potential flow—the base flow *ϕ*. Let *ϕ* operate for a time *τ*, after which the flow rotates through an angle *Θ* and runs again for *τ*; let this sequence repeat for a given (*τ*,*Θ*) and measure its transport properties in a similar manner to previous work on periodically reoriented flows (Metcalfe *et al.* 2006; Speetjens *et al.* 2006; Lester *et al.*2008, 2009). This stirred potential flow we will call the *rotated potential mixing*, or RPM, flow, and we will find that it leads to chaotic advection and transport in the disc.

## 2. Potential flows

Darcy’s Law for flow in porous media is (Homsy 1987)
2.1along with the incompressibility condition ∇⋅*u*=0, where *μ* is fluid viscosity and *K* is the permeability of the porous media. Velocity and pressure are averaged over volumes large enough to contain many pores, yet small compared with the entire domain. Equation (2.1) is a local relation that holds in a region where *K*, though possibly a tensor, is constant. Shortly after Darcy formulated equation (2.1), others recognized it as a potential flow and sought steady-state flow solutions from potential theory for geometries relevant to hydrology (Slichter 1899), and such solutions are still relevant and sought (Townley & Trefry 2000).

The HS flow occurs in the gap between rigid, parallel plates when the gap size *b* is much smaller than the lateral extent *R* of the plates. If we take *x* and *y* as the horizontal coordinates and *z* as the vertical coordinate across the plate gap, then the steady incompressible HS velocity field is given by the viscosity–pressure balance as (Lamb 1932)
2.2where ∇*P* is the pressure gradient in the horizontal direction. We neglect any gravity contribution. In *each* horizontal plane across the gap equation (2.2) is a potential flow with the pressure as the potential; the vertical velocity profile in the gap is parabolic, but the vertical velocity variation can be neglected for a thin enough gap (Lamb 1932). The no-slip condition at boundaries holds, but the boundary layers at the plates are thin enough to be negligible. The boundary layer thickness at horizontal walls is 𝒪(*b*), so domains with large *R* can also be neglected. It is customary to average equation (2.2) over *z* to obtain the two-dimensional HS flow equation
2.3Comparison of equations (2.3) and (2.1) shows the formal equivalence of Darcy and HS flow with the quantity *b*^{2}/12 playing the role of permeability. This equivalence has frequently been exploited in the petroleum literature for experimental visualization of transport, front propagation and fluid–fluid displacement in porous media (Homsy 1987).

Homsy (1987) notes the two cases where the Darcy–HS equivalence breaks down. (i) For miscible fluids with concentration gradients, Taylor dispersion can occur due to the velocity parabola across *b*; in real porous media, dispersion affects both longitudinal and transverse dispersion coefficients, but in HS it affects only longitudinal coefficients (though an evolving two-dimensional height field would recover transverse dispersion). (ii) Equivalence fails severely for immiscible (particularly gas–fluid multiphase) fluids due to interfacial forces being able to drive flow through pores: HS contains no analogous physics. There is no question that HS is a simplification of flow in rock, but to the extent that Darcy flow is a reasonable approximation, then, avoiding the known areas of breakdown, HS flow is also a reasonable starting approximation.

Note that the dipole is merely the simplest configuration of RPM flow in a disc. As many singularities as desired can be added to equation (1.1), with the only physical constraint being that to conserve mass the coefficients of the terms must sum to zero. For instance, singularities placed on the real axis map to the disc boundary, and those placed in the upper half plane map to the disc’s interior. As time dependence is simulated by switching among symmetric copies of the base flow, any base flow configuration and any of its symmetries can be used to stir and drive transport within the disc. Importantly, different symmetries will affect the transport properties, particularly the presence and bifurcation structure of islands, which are regions of the flow with restricted or non-existent fluid exchange with other regions of the flow. In fact, there are an infinite number of combinations of base flows and ways to combine them; many more than can be usefully explored without additional physical input as to what one is trying to accomplish. Also, note that the domain need not be a disc: any domain that can be mapped to the upper half plane can be used without substantial change to the analysis.

For later use we show some transport properties of the steady dipole flow, i.e. without reorientation. For the non-reoriented dipole flow, figure 2 shows the exit time, *t*_{e}, frequency distribution. This is calculated by randomly placing around 3000 points in the disc and running the map *ϕ*(*τ*=1,*Θ*=0). The average area per point is *π*/3000≈10^{−3}. The mean residence time is given by the areal flow rate *Q* divided by the disc area: *t*_{mr}=*Q*/*π*. All points exit the disc, and the maximum *t*_{e} estimates the basin exchange time as *t*_{ex}=382.86. The inset plots the initial points in the disc coloured by according to the scale in the figure. With no rotation, the spatial pattern of exiting is an orderly progression with points near the sink exiting first and those starting in the low-velocity regions near the boundary and the source exiting last.

RPM flow is meant to be a first, simple model, and it needs to have added to it scalar transport and changes to *ϕ* brought about by porosity changes. Nonetheless, just by adding reorientation time dependence to *ϕ* we can report interesting results from even this simplest case.

## 3. Advection in the partially open RPM flow

The potential flow *ϕ*(*τ*) in the disc of duration *τ* followed by a rotation *R*(*Θ*) defines a map *M*=*Rϕ*. Over the control parameter space (*τ*,*Θ*) of the RPM flow, advection in part of the domain is closed and in other parts is open. Closed advection means that fluid does not leave the domain; it remains in the disc forever. In open advection, fluid exits the domain in finite time. For example, two orbits^{1} of *M* are shown in figure 3. The large crosses are the two initial conditions. The thick (thin) lines denote the orbit for the *ϕ* (*R*) part of the map. The grey orbit leaves through the sink, and we have chosen to reinject it at the source (the only instance in this paper where the source and sink are connected), while the black orbit never exits the disc but circulates forever in an island. This is remarkable because, for *ϕ* by itself, all points in the disc exit after a finite area exchange time *t*_{ex}.

Dynamically closed flows are much more widely studied than dynamically open flows.^{2} Closed advection is governed largely by the symmetries of the base flow and boundary conditions (Ottino *et al.* 1992; Speetjens *et al.* 2006). Open advection is governed largely by filamentary unstable manifolds (Tél *et al.* 2005). As filamentary manifolds may be less familiar, we note that filamentary manifolds are fractal sets and the unstable manifolds of chaotic saddles; they form both barriers to transport and concentrating templates for chemical and biological activity. As the RPM flow exhibits transport characteristics of both closed and open flows, we may call it a *partially open flow*. We first discuss transport in the closed part of the RPM flow and then the open part.

### (a) Closed advection: islands in the stream

Transport in stirred flow is generically governed by stable and unstable manifolds of hyperbolic periodic points of the flow. Points approach stable manifolds as they move forward in time and approach unstable manifolds as they move backward in time. Periodic elliptic points are important because they have disconnected (possibly large) island regions around them that sequester or hinder (through cantori) material transport and reaction activity (Ottino *et al.* 1992; Speetjens *et al.* 2006). Islands are regions of rotation without much deformation that have Kolmogorov–Arnold–Moser (KAM) surfaces separating them from the rest of the flow. In three dimensions islands are tubes. Many books and reviews (Wiggins & Ottino 2004) give precise definitions and mathematical details, and several papers (Metcalfe *et al.* 2006) show experimental examples. At *τ*=0, the map *M*=*Rϕ* has a single elliptic point at the origin for all *Θ* (except for the trivial case of *Θ*=0) and the associated island takes up the entire disc. As *τ* rises from zero, where does the origin elliptic point go, how does it bifurcate along the way and how does the area of the associated island change? Do new elliptic points and islands arise? As islands form the regions of closed advection, the answers to these questions determine where closed advection exists in the RPM flow and if it is practically important.

The map *M* moves points in the disc forward in time. The backward-in-time map is the inverse *M*^{−1}=*ϕ*^{−1}*R*^{−1}⋅*R*^{−1}=*R*(−*Θ*), and *ϕ*^{−1}=*R*(−*π*)*ϕ*(*τ*)*R*(*π*)=*ϕ*(−*τ*). Because plane rotations commute, we can write *M*^{−1}=*R*(−*π*)[*ϕR*(−*Θ*)]*R*(*π*) and define a new map
3.1Using *J*, *M*=*R*(*Θ*)*JR*(*Θ*). The *n*-iterated forward and backward maps are
3.2
3.3
The map *J* controls the closed transport in the RPM flow. Note that *ϕ* need not be the dipole flow—it could be any flow in the disc—but that we have assumed that *Θ* and *τ* are fixed at each iteration in deriving the *n*-iterate maps. This assumption can be removed by accepting that the *n*-iterated maps no longer collapse to simple forms and allowing *Θ* and *τ* to change at every iterate.

What advection topology is induced by *J* over the parameter space (*τ*,*Θ*)? *M*^{n} has an axis of symmetry at an angle *α*=−*Θ*/2 that controls the location and evolution of periodic points (Speetjens *et al.* 2006). A symmetry axis immediately says that all period one (P1) points must lie on the symmetry axis. Consider a line of initial conditions advected by *J*⋅*R*^{−1} that rotates the line counter-clockwise without deformation, then *ϕ* advects the line towards the sink, stretching it because the velocity field is not radially uniform. The intersection of the initial line and its first iterate gives the P1 points. As *τ* is increased above zero for any given *Θ*, the origin elliptic point must move along the symmetry axis *α*. The elliptic point’s location in polar coordinates is (*r**,*α*), where *r** is to be determined. Moreover, the rotation of the matrix in brackets in equation (3.2) says that only *Θ*>0 is unique, and the double periodicity says that only *Θ*∈[−*π*,*π*] is unique. Together this means we need only to consider *Θ*∈[0,*π*] because patterns at all other values of *Θ* are reflections. This symmetry also applies to orbits of the open flow. However, if the disc dynamical system is closed by taking points that leave the sink and reinjecting them at the source, then points that never leave the disc obey the above symmetries, but reinjected points obey symmetries modified by exactly how they are reinjected. Jones & Aref (1988) explored several reinjection possibilities, but we will not pursue this point here, except to note that there are engineering problems in porous media where the inlets and outlets are connected for recirculation and where getting the reinjection model right could be crucial. Finally, not all of the reflections may hold for asymmetric base flows.

To find *r** take the symmetry axis (*r*,*α*) and apply . Where the symmetry axis and its first iterate intersect is defined as the location of a P1 point. As the dipole flow has an exact expression for the time an orbit spends on any given streamline, we can derive an expression for *r** as an implicit function of *τ* and *Θ*. Where there is no solution, there is no intersection, and the original P1 elliptic point no longer exists. In an earlier paper (Metcalfe *et al.* 2008), we gave an approximation to this *island existence boundary*; now, we have derived the exact expression for the island existence boundary in (*τ*,*Θ*) space as (Lester & Metcalfe 2008)
3.4which is shown in figure 4 as the curve rising from *τ*=0. A selection of Poincaré maps illustrates island existence and the symmetry effects. As the RPM flow is open (points can and do leave through the sink), these are non-standard Poincaré maps that show only the closed parts of the flow, and the Poincaré maps of the disc are either centred on the parameter point to which they belong or are displaced along a line whose end is on the parameter point to which they belong. Maps show the last 60 of 100 iterates of *J*. The thin lines mark the symmetry axis and the initial conditions of the map. The rotation of the symmetry axis with 2*Θ* and the movement of the P1 elliptic point towards the disc boundary with increasing *τ* is clear. For values above the boundary curve all points exit after three to five iterates of *J*. To the right of the point labelled *Θ*_{crit} is an example of something the symmetry analysis does not reveal. Within the curves labelled *τ*_{u} and *τ*_{l}, which are computed numerically, the elliptic point undergoes a homoclinic bifurcation to a hyperbolic point, but retains a KAM boundary to form a non-elliptic isolated mixing region (Bresler *et al.* 1997).

### (b) Open advection

Tél *et al.* (2005) review reaction and transport in open chaotic flows. By open Tél *et al.* mean a flow that is relatively uniform far upstream and far downstream from a region of chaotic interaction or stirring. This allows treatment of transport in the flow as a chaotic scattering problem. A key theoretical concept is exit time distributions from which one can derive exit rates. The orbits of points with the longest exit times can be followed to approximate the filamentary manifolds controlling transport in the chaotic scattering region. On the other hand, the chemical engineering world would treat reaction and transport in the RPM flow in the same manner as a throughflow reactor and would generate residence time distributions (Nauman & Buffham 1983). We show an example of each approach for advection in the RPM flow, emphasizing somewhat the parametric variation, i.e. how residence or exit times change over the entirety of the control parameter space (*τ*,*Θ*) of the flow.

Figure 5 shows an example of an exit time distribution. We initialize by randomly placing about 3000 points in the disc, then iterate these points with *M* and record the time it takes for each to exit through the sink. There is no reinjection, and after 10*t*_{ex} we declare the remaining points to be in an island. (Island size variation from *t*_{ex} to 16*t*_{ex} is ±2 points.) The frequency distribution shows only points that did exit. Several things are noteworthy. The exit order has some simplicity. Points in a scallop shape near the sink leave, similarly to the unreoriented flow. Then, new points are rotated near the sink and they leave in their turn. However, points in the boundary between scallops have markedly longer exit times. Instead of exiting with their neighbours, these veins flow around the island first before exiting. They are 2–4% of the total and account for the toe in the frequency distribution above 0.

For the unreoriented flow *ϕ*, a point released at the source moves down the central streamline and exits through the sink after a time *t*_{r0}=133.342. Figure 6 shows the frequency distribution of residence times *t*_{r} of a point released at the source over a grid of about 1 million points in (*τ*,*Θ*) for the RPM flow. The inset plots the frequency on a log scale showing the distribution’s power-law tails. A fit of the tails—grey lines in the figure—gives exponents of 0.222±0.006 for the short-time tail and −0.275±0.002 for the long-time tail. Both exponents are close to, but statistically different from, 1/4. Stirring expands the residence time distribution by three orders of magnitude both above and below the unstirred residence time. Figure 7 shows the (*τ*,*Θ*) parameter space colour coded by . There is clear structure, including what seems to be self-similar resonance ‘tongues’ in the lower left corner of the space. Neither the power laws in the frequency distribution nor the self-similar parameter space structure currently have a solid theoretical understanding.

### (c) Experiment

To create an experimental RPM flow, we made a circular HS cell of 120 mm diameter with a *b*=120 μm gap between the bottom flat metal plate and the top flat glass plate. At the disc edge of the metal plate, we drilled 360 evenly spaced 0.6 mm diameter holes. Below the holes is a rotatable manifold that closes all the holes except for two holes 180° apart. These two openings are the inlet and outlet wells, and the openings in the manifold lead to fixed reservoirs. By stopping the flow and rotating the manifold, we can generate a dipole flow in any orientation at 1° increments. The working fluid was a 50–70% glycerin–water mixture. Further experimental details will be given elsewhere.

Figure 8 shows an example of iterating the map *M* (a repeated sequence of flow and reorientation) for *τ*=0.1 and *Θ*=*π*/4, i.e. below the island existence boundary. The top line shows the experiment, where the cell is initially filled with red fluorescent dye and clear fluid is pumped through the inlet. The next line shows simulations at the same conditions, and the bottom line gives the number of iterations of *M*. For iterations above 13, the red island persists essentially unchanged. Note what appears to be filamentary manifolds developing and thinning. Figure 9 shows the island existence boundary decorated with several experiments at the points indicated by the grey dots. Below the line islands of some size always exist, though closer to the boundary island size becomes smaller. Above the line all the dyed fluid initially in the cell eventually leaves.

## 4. Conclusions and discussion

We have described an experiment and theoretical model for advection in porous media. It is a first step in realizing a fully coupled model of advection, reaction and matrix deformation. However, with stirred advection alone there are new results. The model is a rotated potential flow in a disc with any number of source–sink wells. The model is a *partially open flow* that has characteristics of both dynamically closed and dynamically open flows. Using symmetry, we have derived the character and extent of the closed (island) parts of the flow and where in the stirring control parameter space isolated regions can exist. Experiments and computations verified the theory. In the open part of the flow exit time distributions and experiments show the filamentary manifolds predicted to control transport in the open part of the flow. Residence time distributions for fluid entering the disc at the source show power-law tails and a self-similar structure in parameter space. The next step is to add reactive transport to achieve our goal of a three-field, fully coupled model of flow in porous media.

## Acknowledgements

This work has been supported by CSIRO’s Minerals Down Under flagship and CSIRO’s Complex Systems Science initiative. Laboratory work by P. Kulkarni at CSIRO Highett was funded by a Minerals Down Under/iVEC internship and by the Research Foundation of the City University of New York. Tony Kilpatrick, Tony Swallow, Robert Stewart and Dean Harris provided technical support at various stages of building the experiment. Some of the computations used CSIRO’s Advanced Scientific Computing parallel cluster, and we are grateful to Gareth Williams for helping adapt code.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

↵1 The discontinuities in the orbits are caused by rotating the points instead of the flow. In an experiment, the orbits are continuous and the flow is rotated. In this sense figure 3—and subsequent figures—can be thought of as being viewed from a reference frame rotating with the flow.

↵2 Duct flows are closed in the plane but extruded in the

*z*direction, e.g. a pipe, and are often called open flows (Wiggins & Ottino 2004; Metcalfe*et al.*2006) because fluid exits each plane. Are duct flows open or closed in the sense used here? Insofar as being a dynamical system, duct flows are closed because in ducts time is just a parameter that is interchangeable with distance along the duct, and transverse to the extruded direction the duct flow is closed. In contrast, in a source–sink flow material really leaves the dynamical system in finite time unless it is explicitly reinjected into the dynamical system.

- © 2010 The Royal Society