# Fluid Deformation and Mixing in Porous Media as Drivers for Chemical and Biological Processes

Tanguy Le Borgne, Joris Heyman

► To cite this version:

Tanguy Le Borgne, Joris Heyman. Fluid Deformation and Mixing in Porous Media as Drivers for Chemical and Biological Processes. Annual Review of Fluid Mechanics, 2026, 58 (1), pp.443-471. 10.1146/annurev-fluid-112723-051940 . insu-05472952

HAL Id: insu-05472952

https://insu.hal.science/insu-05472952v1

Submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

*Annual Review of Fluid Mechanics*

# Fluid Deformation and Mixing in Porous Media as Drivers for Chemical and Biological Processes

Tanguy Le Borgne<fnref n="1" />,<fnref n="2" /> and Joris Heyman<fnref n="1" />

ANNUAL
REVIEWS

<figure><img src="image_3.png" /></figure>

www.annualreviews.org

• Download figures

• Navigate cited references

• Keyword search

• Explore related articles

• Share via email or social media

Annu. Rev. Fluid Mech. 2026. 58:443–71

First published as a Review in Advance on
October 3, 2025

The *Annual Review of Fluid Mechanics* is online at
fluid.annualreviews.org

https://doi.org/10.1146/annurev-fluid-112723-051940

Copyright © 2026 by the author(s). This work is
licensed under a Creative Commons Attribution 4.0
International License, which permits unrestricted
use, distribution, and reproduction in any medium,
provided the original author and source are credited.
See credit lines of images or other third-party
material in this article for license information.

<figure><img src="image_4.png" /></figure>

Keywords

porous media, mixing, transport, dispersion, chaotic flows, chemical
reactions, microbial processes

Abstract

Porous media flows are generally viewed as inefficient mixers, where solutes may be dispersed yet poorly mixed, making mixing a critical limiting factor for a wide range of processes. The complexity and opacity of porous structures have long made these dynamics difficult to observe. With emerging experimental techniques, concepts and models of mixing in porous media are rapidly evolving. Recent advances link mixing dynamics to fluid deformation arising in flow through porous materials. Unlike diffusion and dispersion, which only dissipate chemical gradients, fluid shear and stretching amplify and sustain them. This review explores the role of fluid deformation in governing mixing, chemical reactions, and biological processes in porous media. We begin by highlighting key experimental observations that have improved our understanding of mixing in these systems. We then examine the fundamental concepts, models, and open questions surrounding fluid deformation and mixing in porous media, emphasizing their dependence on material structure, heterogeneity, dimensionality, and transient flow phenomena, as well as their interaction with chemical and biological processes.

<footnote n="1">¹Géosciences Rennes, UMR 6118 CNRS, Université de Rennes, Rennes, France; email: tanguy.le-borgne@univ-rennes.fr, joris.heyman@univ-rennes.fr</footnote>

<footnote n="2">²Nord Center, University of Oslo, Oslo, Norway</footnote>
# 1. INTRODUCTION

Fluid flow through porous materials is ubiquitous across a broad range of natural and engineered systems, including soils, rocks, biological tissues, membranes, chromatographic columns, and batteries. The large specific surface area of porous media provides a particularly favorable environment for diverse chemical reactions and biological activity. Typically occurring at low Reynolds numbers, porous media flows are shaped by fluctuations in pore aperture and connectivity, leading to a broad distribution of velocities and complex transport and mixing dynamics. One of the most striking features of transport in porous media is the emergence of non-Fickian dispersion dynamics, characterized by power-law distributions of solute travel times. The study of such anomalous transport phenomena has been the focus of extensive research over the past decades, leading to successful descriptions such as continuous time random walk models (Berkowitz et al. 2006, Neuman & Tartakovsky 2009).

However, it is now well-established that dispersion models alone fail to accurately capture mixing, as solutes can be widely dispersed yet remain poorly mixed (Dentz et al. 2011, 2023). Mixing in porous media arises from the combined effects of advection and diffusion, which drive the evolution of a scalar field toward a more homogeneous state. It plays a central role in geological, environmental, and engineering processes (Dentz et al. 2011, Rolle & Le Borgne 2019, Valocchi et al. 2019) by regulating chemical reactions, biological activity, and the dilution of transported chemical species. Velocity gradients accelerate solute mixing by stretching fluid elements, which amplifies concentration gradients and enhances the resulting diffusive fluxes. Although the kinematics of fluid deformation have been studied extensively in chaotic and turbulent flows (Batchelor et al. 1959, Ottino 1990, Villermaux 2019), the fundamental principles governing fluid deformation in porous media—and their implications for mixing, chemical reactions, and biological processes—are only beginning to be understood.

The complexity and opacity of porous structures have long hindered direct observation of fluid deformation kinematics and mixing dynamics in these systems, along with their effects on chemical and biological processes. However, emerging experimental techniques are rapidly transforming our understanding, leading to new theoretical and modeling approaches. We review key concepts and open questions related to fluid deformation kinematics and mixing in porous media, focusing on the role of pore structure, heterogeneity, dimensionality, and transient flow phenomena. We then examine how fluid deformation controls mixing in porous media and explore its broader impact on chemical and biological processes.

# 2. MIXING MECHANISMS

## 2.1. Advection, Diffusion, and Dispersion

Solute mixing in porous media emerges from the interplay between advection, mechanical dispersion, and molecular diffusion. At the pore scale, creeping laminar flows are governed by the Stokes equation, which expresses the pressure loss due to viscous dissipation as $\nabla p = \mu \nabla^2 \mathbf{v}$, where $p$ is the pressure, $\mu$ is the dynamic viscosity, and $\mathbf{v}$ is the velocity. Fluctuations in the velocity of the fluid arise from interactions with the no-slip boundary conditions on the solid pore boundaries, variations in the pore aperture, and connectivity of the pore space. These factors generally give rise to broad power-law distributions of Eulerian velocities, $p(\mathbf{v}) \sim \mathbf{v}^{\beta-2}$, with $\mathbf{v} \equiv |\mathbf{v}|$ and $\beta$ a critical exponent that controls the emergence of anomalous transport, as discussed below (Puyguiraud et al. 2021). Molecular diffusion of solutes obeys Fick's law, with diffusive flux $\mathbf{j} = -D_m \nabla C$, with $D_m$ the diffusivity and $C$ the scalar concentration. The Péclet number, $Pe = d \bar{v}/D_m$, compares the characteristic diffusion time with the mean advection time at a characteristic pore scale $d$. At scales significantly larger than the characteristic grain size, the flow is described by Darcy's law, which
expresses the pressure gradient due to viscous effects as $∇p = -μφ/k · v$, where $v$ is the average fluid velocity, $φ$ is the porosity, and $k$ is the permeability. At this scale, velocity distributions are also generally broad due to large spatial fluctuations in permeability, which are often modeled as log-normal distributions (Le Borgne et al. 2015). The fluctuations of velocity induce the dispersion of solutes. Depending on velocity statistics, the dispersion can be Fickian (typically for power-law exponents $β > 2$) or anomalous (for $β ≤ 2$) (Dentz et al. 2011). In the Fickian regime, hydrodynamic dispersion is classically modeled as an effective diffusion process, with characteristic spreading distances $σ_x = (σ_x^‖, σ_x^⊥)$ growing as $σ_x ∼ \sqrt{αvt}$, with $α = (α_‖, α_⊥)$ the dispersivities in longitudinal and transverse directions with respect to the mean flow. Experimental measurements in porous media columns suggest that $α_‖ ∼ d$, while $α_⊥ ∼ 0.1d$ (Delgado 2007).

## 2.2. Stretching-Enhanced Diffusion

As a solute plume moves through a porous medium (Figure 1), it is stretched by pore-scale velocity gradients into a collection of elongated structures that take the form of lamellae. Pore-scale velocity gradients emerge due to the no-slip and no-penetrating boundary conditions imposed by the solid matrix (Figure 2a). In laminar flows, velocity gradients have a characteristic scale given by the pore size $d$ and a typical magnitude of $\bar{v}/d$. At the pore scale, solute strips undergo smooth deformations driven by these velocity gradients. At larger scales, velocity gradients can also exist due to heterogeneities in the porosity or permeability of the solid matrix. At the Darcy scale, permeability fields are often assumed to follow log-normal distributions, characterized by a characteristic mean, variance, and correlation length (Le Borgne et al. 2015) (Figure 2c).

<figure><img src="image_5.png" /><figcaption>b Spreading scale</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>c Mixing scale</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>d Concentration variance</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>e Concentration gradient</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 1</figcaption></figure>

Spreading and mixing of a 3D solute plume in a porous bead pack, imaged with refractive index matching and laser-induced fluorescence (Heyman et al. 2020). (a) The solute is continuously injected in a steady background flow (from left to right). Schematic evolution of (b) the spreading scale, (c) the mixing scale $s$, (d) the concentration variance $σ_c^2$, and (e) the concentration gradient $∇C$, for the case of a chaotic (red) and nonchaotic (orange) flow.
a Shear driven by no-slip boundary conditions at the pore scale

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

b Shear driven by macroscopic boundary conditions at the Darcy scale

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

c Shear driven by permeability heterogeneity at the Darcy scale

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Figure 2</figcaption></figure>

Shear deformation from pore to Darcy scales. (a) Experimental imaging of a solute front in a 2D porous micromodel: (top) a concentration field, with concentrations ranging from high (yellow) to low (blue), and (bottom) concentration gradients, with gradients ranging from high (green) to low (black), highlighting the formation of thin lamellae resulting from shear deformation. Panel adapted from Borgman et al. (2023) (CC BY-NC-ND 4.0). (b) Shear induced by a hydraulic head gradient applied at the upper boundary of a homogeneous porous medium. (Top) A numerical simulation shows the transport and dilution of a solute plume (concentrations ranging from black to gray) injected at the surface. (Bottom) Mixing of a solute dye (blue) is imaged in a porous media tank. The solute is injected at the surface as three pulses, which are transported along streamlines with different velocities toward the outlet on the upper left corner. Top subpanel adapted from Bandopadhyay et al. (2018). Bottom subpanel adapted with permission from Jules et al. (2021). (c) Heterogeneous shear leading to algebraic deformation in 2D permeability fields composed of spatially variable permeabilities, with (top) a simulated flow field (low velocities in blue and high velocities in red) and (bottom) the transport of a solute plume injected at the left boundary, shown at two times (early time in red, later time in blue). Panel adapted from Le Borgne et al. (2015).

## 2.3. Dispersion Versus Mixing

As illustrated by recent 3D imaging of solute plumes transported in porous media (Figure 1), solute transport in porous media involves two coupled processes, dispersion and mixing. Dispersion describes the spatial spreading of the solute plumes and is quantified by the second spatial moment of the plume (Figure 1b) or the tracer breakthrough curves. Mixing refers to the reduction of solute concentration heterogeneity due to the combined effects of plume stretching and diffusion (Dentz et al. 2011, 2023). Dispersion and mixing are often correlated, since greater dispersion enhances plume stretching, which in turn promotes mixing. However, they can also vary independently, depending on the pore structure and the balance between advective and diffusive processes. Despite this fundamental distinction, dispersion and mixing are frequently conflated, leading to significant discrepancies between theoretical predictions and experimental observations (Dentz et al. 2011).

The concentration probability density function (PDF) provides a complete description of the mixing state of the plume (Le Borgne et al. 2015, Villermaux 2019). The second moment of the distribution, the concentration variance $\sigma_c^2 = \langle (c - \bar{c})^2 \rangle$, is a practical measure of the degree of scalar heterogeneity, with $\sigma_c^2 \to 0$ for homogeneous concentration fields (Figure 1d). Mixing can also be characterized by entropy-based metrics, such as the dilution index (Kitanidis 1994), or by the evolution of the mixing scale (Le Borgne et al. 2011), which characterizes the largest scale below which the concentration field is homogeneous (Figure 1c). Other metrics focus on concentration gradients that drive diffusive fluxes. These include the scalar dissipation rate (Le Borgne et al.
2010, Jiménez-Martínez et al. 2015), the average concentration gradient within a plume or front (Borgman et al. 2023) (Figure 1e), and the statistics of concentration increments across spatial scales, which describe intermittency in the scalar field (Le Borgne et al. 2017). Finally, mixing can also be characterized through its impact on reaction rates in mixing-limited chemical processes (Gramling et al. 2002, Werth et al. 2006, de Anna et al. 2014, Rolle & Le Borgne 2019, Valocchi et al. 2019, Sanquer et al. 2024), as discussed in more detail in Section 5.

Due to fluid incompressibility, stretching in one direction is accompanied by compression in the perpendicular direction. This compression reduces the width of the lamellar structures that form the support of the solute plume (Figure 1c) and increases concentration gradients and thus diffusive fluxes (Figure 1e), resulting in enhanced mixing (Figure 1d). The coupling between fluid deformation and molecular diffusion is not captured by classical dispersion theories, which primarily account for gradient dissipation. Hence, the coupling between stretching and diffusion gives rise to effective mixing laws that differ fundamentally from those described by dispersion models. Stretching-enhanced diffusion has been described and quantified within the lamella framework (Villermaux 2012; Le Borgne et al. 2013, 2015; Villermaux 2019), also called the Lagrangian stretching framework, diffusive strip, or diffuselet methods (Balkovsky & Fouxon 1999; Meunier & Villermaux 2010, 2022). The lamella framework, originally developed for mixing in turbulent and chaotic flows, is presented in more detail in the following sections.

3. FLUID DEFORMATION KINEMATICS

The kinematics of fluid deformation in porous media flows and their effect on mixing have only been investigated relatively recently (Villermaux 2012, Le Borgne et al. 2013, Lester et al. 2013). Different methods for measuring fluid deformation in porous media are highlighted in the sidebar titled Experimental Measurement of Fluid Deformation in Porous Media.

Fluid deformation processes are described through the evolution of the deformation gradient tensor F(t|x₀) (Ottino 1990, Lester et al. 2018), which characterizes how an infinitesimal fluid parcel initially located at x₀ deforms with time. The evolution of F is given by

$$
\frac{d\mathbf{F}(t|x_0)}{dt} = \epsilon(t|x_0) \cdot \mathbf{F}(t|x_0), \qquad 1.
$$

with the initial condition **F**(t = 0|x₀) = **I**, and the Lagrangian deformation rate tensor is

$$
\epsilon_{ij}(t|\mathbf{x}_0) \equiv \frac{\partial v_i}{\partial x_j} [\mathbf{X}(t|\mathbf{x}_0)], \qquad 2.
$$

where **X**(t|x₀) is the position at time t of the fluid parcel initially located at x₀, and ∂vᵢ/∂xⱼ is the velocity gradient tensor. In other words, fluid deformation results from the integration of the history of velocity gradients experienced along a Lagrangian trajectory. Elongation of an elementary material line of initial size ℓ₀ by the flow is given by

$$
\rho(t) = \frac{\|\mathbf{F}(t) \cdot \mathbf{\ell}_0\|}{\|\mathbf{\ell}_0\|}. \qquad 3.
$$

The stretching rate is defined by

$$
\gamma(t) = \frac{1}{\rho} \frac{d\rho}{dt}. \qquad 4.
$$

A key challenge lies in deterministically or statistically describing the evolution of γ or ρ in porous media flows. In the following, we review recent observations and conceptual advances in understanding the primary modes of fluid deformation in porous media.
# EXPERIMENTAL MEASUREMENT OF FLUID DEFORMATION IN POROUS MEDIA

## Particle Image Velocimetry

Stretching rates can be estimated from 3D velocity fields obtained by particle image velocimetry in porous media (Souzy et al. 2020). The velocity gradient tensor (Equation 1) and pair separation statistics (Equation 3) can be numerically computed via particle tracking based on experimentally measured velocity fields. However, elongation is highly sensitive to both the continuity of the velocity field and the accuracy of its gradients (Turuban et al. 2019). As a result, experimental noise in the velocity measurements can potentially lead to an overestimation of stretching rates.

## Mixing Interface Deformation

Continuous injection of a scalar dye produces a plume whose perimeter grows exponentially with the distance of advection, at a rate $log(ρ)/t ≈ 3λ/2$ (Equation 14), enabling the estimation of λ. However, as discussed in Section 4.2, filament aggregation limits the expansion of the mixing interface, which can result in an underestimation of stretching rates (Velásquez-Parra et al. 2024). To address this limitation, Heyman et al. (2020) have manually measured the transverse deformation of the plume, keeping track of the existence of aggregated lamellae.

## Mixing Statistics

Based on the prediction of mixing theories for the evolution of concentration statistics, it is possible to indirectly estimate stretching rates. In the lamella theory, the decay of scalar variance follows $σ_c^2 ∼ exp(-λt/2)$ (Equation 19), providing a means to quantify stretching before aggregation becomes significant (Heyman et al. 2021). Beyond this stage, aggregation theory predicts that the tail exponent of the concentration distribution increases with the number of aggregation events, which in turn depends on the stretching rate. Using this method, Kree & Villermaux (2017) measured a stretching rate $λ = 0.2 t_a^{-1}$ in bead packs, close to mixing interface deformation estimates (Heyman et al. 2020).

### 3.1. Shear Deformation

As an initially flat solute front invades a porous medium, it deforms as it flows through the grains (Borgman et al. 2023) (Figure 2a). In porous media, the simplest mode of fluid deformation is shear, which arises due to the no-slip condition at the solid boundaries (Puyguiraud et al. 2020, Yoon et al. 2021). In simple shear flows, the deformation rate tensor has a unique constant off-diagonal component. Integration of Equation 1 leads to a deformation gradient tensor with diagonal components equal to unity and an off-diagonal term proportional to t. Using Equation 3, it follows that $ρ ~ t$ and $γ ~ 1/t$ (see the Supplemental Appendix, section 1). The temporal decay of γ is a direct consequence of the progressive alignment of material lines with flow direction when the velocity gradients are uniform and steady, leading to a relatively smaller area exposed to shear deformation as time progresses. Shear-enhanced mixing has been investigated experimentally by Souzy et al. (2018) and successfully described by the lamella framework (Villermaux 2019). Experiments in 2D porous media demonstrated the impact of shear deformation on the dynamics of concentration gradients (Borgman et al. 2023) and the dilution and reaction rates (de Anna et al. 2014).

Above the pore scale, shear deformations can also arise due to nonuniform flows forming at macroscopic scales. In streambeds and hillslopes, topography-driven hydraulic head gradients generate nested streamlines (Figure 2b), which induce shear deformation and enhance the mixing of solute plumes (Bandopadhyay et al. 2018). Similarly, fluid shear also develops at freshwater-saltwater interfaces, driving mixing and reaction in coastal aquifers (De Vriendt et al. 2022). At the

**Supplemental Material** >
Darcy scale, the spatial heterogeneity of permeability produces nonuniform flows, with velocity varying over orders of magnitude in short distances (Le Borgne et al. 2015). Solute plumes transported through heterogeneous permeability fields undergo a sequence of shear deformations (Figure 2c), resulting in an algebraically growing elongation $\rho \sim t^{\alpha}$, still with a decaying stretching rate $\gamma \sim 1/t$ (Le Borgne et al. 2015). The exponent $\alpha$ increases with the heterogeneity of $k$, taking values from 0.5 (sublinear) to 2 (superlinear). Dentz et al. (2016) developed a theoretical description of elongation in 2D heterogeneous Darcy flows in the framework of coupled continuous-time random walks. They found that $\alpha$ is controlled by the statistics of low velocities, which contribute proportionally more to deformation than high velocities (see the Supplemental Appendix, section 3). For a classical power-law distribution of low velocities, $p(v) \sim v^{\beta-2}$, the elongation evolves as

$$ 
\begin{align*}
\rho &\sim t^{1/2} \quad \text{for } \beta > 4, \\
\rho &\sim t^{2/\beta} \quad \text{for } 2 < \beta < 4, \\
\rho &\sim t^{3-\beta} \quad \text{for } 1 < \beta < 2.
\end{align*} 
\quad 5.
$$

Hence, for low heterogeneity levels ($\beta > 4$), successive random shears along streamlines lead to a classical diffusive growth of the elongation $\rho \sim t^{\zeta}$ with $\zeta = 1/2$ (Equation 5). In contrast, for stronger levels of heterogeneities ($\beta < 4$), the amplification of shear deformation in low-velocity zones (see the Supplemental Appendix, section 3) leads to a superdiffusive growth, $\zeta > 1/2$ (Equation 5). Since the deformed line is confined within a characteristic longitudinal dispersion length $\sigma_x^{\parallel}$, the elongation may be related to the dispersion dynamics by $\rho \sim (\sigma_x^{\parallel})^{d_f}$, where $d_f$ is the fractal dimension of the line (Le Borgne et al. 2015). As both the elongation and the dispersion scale evolve as power laws of time, this implies that material lines in 2D heterogeneous permeability fields have stable fractal dimension, with $d_f$ increasing with the level of heterogeneity (see the Supplemental Appendix, section 4), a result qualitatively consistent with numerical simulations (Le Borgne et al. 2015). Recent advances further discuss the relationship between line stretching and non-Fickian dispersion in random flows (Lester & Dentz 2025).

## 3.2. Chaotic Stretching and Folding

While the stretching rate $\gamma$ decays in steady 2D flows, it can be sustained by successive reorientations of the shear directions (Hinch 1999). This occurs in unsteady 2D flows or in steady 3D flows, leading to chaotic mixing (Ottino 1990). Due to the inherent topological complexity of 3D porous structures, chaotic mixing has been shown to emerge as a ubiquitous feature across natural and engineered porous media (Lester et al. 2013, 2025). The resulting exponential divergence of nearby fluid trajectories is characterized by an asymptotically constant stretching rate, $\gamma(t) \to \lambda$, called the Lyapunov exponent. For a historical perspective on experimental advances that have enabled the observation and understanding of chaotic flow in steady flows through porous media, readers are referred to the sidebar titled A Brief History of the Experimental Demonstration of the Chaotic Nature of Flow at the Pore Scale. At the pore scale, chaotic advection arises from the topological complexity inherent to steady Stokes flow through 3D porous structures. It has been described for two types of architectures commonly found in natural and engineered porous materials: pore networks and granular media (Lester et al. 2025). In pore networks, the solid phase surrounds the fluid, which flows through enclosed conduits (Figure 3a). In granular media, the fluid surrounds and flows around solid inclusions (Figure 3b). As discussed in the following, both structures were found to spontaneously generate chaotic advection in a steady Stokes flow, although through different mechanisms (Lester et al. 2025).
A BRIEF HISTORY OF THE EXPERIMENTAL DEMONSTRATION OF THE CHAOTIC NATURE OF FLOW AT THE PORE SCALE

In 1989, Hulin & Plona (1989, p. 1345) proposed an insightful phenological description of the evolution of two nearby fluid particles at the pore scale, leading to irreversible mixing. They argued that “the initial flow tube gets divided into many very thin filaments. Very quickly the filaments get so thin that even a small amount of molecular diffusion is sufficient to have a particle hop to a nearby filament corresponding to a completely different flow history and coming from a different pore.” They deduced that the distance at which this process leads to irreversible mixing is a logarithmic function of Pe, which is the typical scaling arising in chaotic and turbulent flows (Villermaux 2019). However, at that time, there was no experimental means to confirm this hypothesis. Although chaotic advection in steady laminar flows was long documented in the context of static mixers (Ottino 1990, Stroock et al. 2002), its ubiquity in porous media has only been recognized more recently. In a theoretical study titled “Is Chaotic Advection Inherent to Porous Media Flows?” Lester et al. (2013) laid the foundation for recent advances in understanding chaotic mixing at the pore scale. They found that the branching and merging of flow tubes are the key to trigger chaotic mixing (Figure 3a). The simulations of Turuban et al. (2018, 2019) of Stokes flow in periodic bead packing confirmed the presence of chaotic orbits in granular media, albeit arising from a mechanism different from that observed in pore networks (Figure 3b). Experimental confirmations of the emergence of chaotic pore-scale mixing were recently obtained by advanced 3D imaging techniques (Kree & Villermaux 2017, Heyman et al. 2020, Souzy et al. 2020), more than 30 years after the first intuition of Hulin and Plona.

**3.2.1. Pore networks.** In pore networks, branching and merging of the flow at the pore junctions lead to exponential fluid deformations (Lester et al. 2013). Since the fluid is confined in closed conduits (Figure 3a, subpanel i), transverse deformation is the result of two successive steps. First, at pore branches, the flow decelerates because the total transverse flow area increases. Compression and stretching of the fluid occur in the longitudinal and transverse directions, respectively. Then, at the pore junctions, the flow accelerates again due to the constriction, which leads to stretching and compression in the longitudinal and transverse directions, respectively (Figure 3a, subpanel ii). This sequential process results in an exponential elongation of transverse interfaces when the pores recombine at an angle different from that of their branching. Reorientation occurs due to the natural tortuosity of conducts (Figure 3a, subpanel iii). The Baker map (Ottino 1990) is recovered for a perpendicular reorientation in between pores. The Lyapunov exponent is then equal to $\lambda = \log_2 t_a^{-1}$, with $t_a = x_p/v$ the advection time, where $x_p$ is the longitudinal distance between branches and junctions. In natural media, the randomness of pore reorientation angles would lead to a lower average Lyapunov, $\lambda \approx 0.12 t_a^{-1}$ (Lester et al. 2013, 2016b).

**3.2.2. Granular media.** In granular media, fluid flow is not constrained by closed conducts, allowing fluid deformation to occur simultaneously in the three directions of space (Figure 3b, subpanel i). Stagnation points naturally arise on the solid boundary of objects (Oxaal et al. 1994), due to the no-slip condition (Figure 3b, subpanel ii). In convex regions, stagnation points are generally nodes that do not impart significant transverse stretching. In concave regions, stagnation points are characterized by directions of net stretching and compression, which give rise to 2D manifolds that attract and repel fluid trajectories, respectively. These manifolds interact in the fluid domain, leading to chaotic advection (Turuban et al. 2019).

Unlike pore networks, granular media do not have well-defined pore junctions or branches. Instead, the network of contact points forms the backbone of fluid stretching and folding, leading to exponential elongation (Heyman et al. 2020). At these contact points, which are stagnation points for the flow, concavity dictates distinct directions of stretching and compression. Upstream of a
a Pore networks

i Closed pores

<figure><img src="image_16.png" /></figure>

ii Stretching, reorientation, compression

<figure><img src="image_17.png" /></figure>

iii Transverse lamellar patterns

<figure><img src="image_18.png" /></figure>

b Granular media

i Open pores

<figure><img src="image_19.png" /></figure>

ii Stretching and folding at grain contact

<figure><img src="image_20.png" /></figure>

iii Transverse lamellar patterns

<figure><img src="image_21.png" /><figcaption>Figure 3</figcaption></figure>

Mechanisms producing chaotic fluid deformation in (a) 3D porous media in pore networks and (b) granular media. (a, i) Pore networks, such as consolidated rocks, are made of closed conduits. (ii) Branches and junctions control the fluid deformation via a sequential process involving stretching, reorientation, and compression. (iii) Transverse cross section of material lines after 0 (left) and 10 (right) pore lengths. Subpanel adapted from Lester et al. (2016a). (b, i) Granular media, such as bead packs, are made of touching solids and form open cavities in between. (ii) The mechanism driving fluid deformation in bead packs involves stretching and folding at grain contact. $x_c$ is the mean distance between interacting contacts and $\delta$ the stretching factor at each contact. Stagnation points are denoted with a green asterisk. (iii) Transverse image of a tracer dye after 10 grain diameters (left) and particle tracking simulation in bead packs after the same distance (right). Subpanel adapted from Heyman et al. (2020).

contact point, the fluid decelerates, resulting in net compression along the longitudinal direction. In the transverse plane, the fluid stretches around the grains in one direction while compressing in the other by a factor $\delta$ (Lester et al. 2025) (Figure 3b, subpanel ii). Downstream of the contact point, the flow behavior is inverted: The fluid accelerates longitudinally while undergoing compression and elongation in the transverse plane by the same factor. Due to the singular nature of the contact points and the spherical geometry of the grains, stable and unstable flow manifolds always intersect orthogonally at the contact point. Flow asymmetry upstream and downstream of the contact point induces a perturbation sufficient to form a kink on material lines (Figure 3b, subpanels ii,iii). Compression and elongation in the flow-transverse plane further lead to the folding and doubling of the length of material lines. Heyman et al. (2020) proposed that the Lyapunov exponent can thus be modeled as $\lambda = (d \log 2/x_c)t_a^{-1}$, where $x_c/d$ is the average length per bead diameter $d$ of a streamline that joins two interacting bead contacts and $t_a$ is the advection time over one bead diameter ($t_a = d/\bar{v}$, with $\bar{v}$ the mean advection velocity magnitude). In random loose packings, $x_c$ was found to be $x_c \approx 3.5d$ such that $\lambda \approx 0.2 t_a^{-1}$. This corresponds to a stirring efficiency—defined as the ratio of fluid stretching rate to mean strain rate—of a few percent (Heyman et al. 2020). The value of $\lambda$ in random bead packs is still a matter of debate, with estimates ranging from $0.13 t_a^{-1}$ to $0.47 t_a^{-1}$ (see the Supplemental Appendix, section 2). The development of chaotic mixing at the pore scale in granular media was investigated numerically in periodic

Supplemental Material >
body-centered bead packing (Turuban et al. 2019). In such basic packing, a steady Stokes flow becomes fully chaotic when the mean flow orientation does not align with crystal symmetries. The Lyapunov exponent varies between 0 (for symmetric flow) and a maximum of $0.28 t_a^{-1}$ when varying the flow orientation (Lester et al. 2025).

### 3.2.3. Fluid deformation at the Darcy scale and in fractured media.
The existence of exponential fluid deformation in a steady flow at the Darcy scale is still debated. Isotropic permeability fields do not produce the helicity necessary for chaotic advection to develop in 3D steady flows, even when considering heterogeneous systems (Lester et al. 2023), but anisotropic or discontinuous fields do (Ye et al. 2015, Lester et al. 2024). Thus, chaotic advection is in principle possible at large scales in steady flows through natural discontinuous and anisotropic geological media. Fractured media—comprising quasi-2D fracture planes intersecting in three dimensions—represent highly discontinuous systems (Hyman & Jiménez-Martínez 2018, Davy et al. 2024) that may give rise to complex fluid deformation behaviors, which remain to be characterized.

### 3.2.4. Inertial, unstable, and intermittent flows.
The stretching dynamics discussed above occur in steady flows and are thus dictated essentially by the porous media structure and dimensionality. In many systems, temporal fluctuations of fluid flow can be present due to tidal or seasonal forcing, river level fluctuations, or rainfall events, for instance. At the Darcy scale, such fluctuations coupled with medium heterogeneity enhance dispersion (de Dreuzy et al. 2012), fluid stretching (Trefry & Lester 2020), and mixing (Trefry et al. 2012, Neupauer et al. 2014, Pool et al. 2016). Fluctuations can also develop as a result of nonlinear fluid dynamics, instabilities, and coupling with transport processes. Multiphase flows drive transient flow dynamics due to the interplay between fluid flow and capillary forces, leading to enhanced dispersion and reactions (Mathiesen et al. 2023, Zhang et al. 2025). Other instabilities arise due to density-driven flows (Hidalgo et al. 2015, De Wit 2020), viscous fingering (Jha et al. 2011), elastic relaxation of polymers (Browne & Datta 2024), and flow bifurcations at intermediate Reynolds numbers (Yang et al. 2024). Establishing the deformation kinematics that characterize these phenomena to predict the resulting mixing and reaction rates is an outstanding question that could be explored with the concepts and methods presented here.

### 3.2.5. Analogies with other chaotic flows.
Chaotic advection arises with a minimum of 3 degrees of freedom (Ottino 1990), that is, either in steady flows in 3D porous media or in 2D unsteady flows. In a steady laminar flow through porous media, exponential deformation and scalar gradients appear essentially in the plane normal to the mean flow (Lester et al. 2018) (Figure 4c). Thus, mixing in porous media is not far from 2D time-dependent smooth chaotic flows (Figure 4b), such as the sine flow (Ottino 1990). Scalar mixing by such an archetypal flow was thoroughly studied (Haynes & Vanneste 2005), and the resulting mixing dynamics share common properties with those observed in porous media (Heyman et al. 2021). Despite the large difference in the Reynolds number, scalar turbulence also shows striking resemblance to mixing in porous media (Figure 4a) in the so-called Batchelor regime, where the flow is smooth at the scale of scalar gradients (the Batchelor scale). This occurs at high Schmidt numbers ($Sc = v/D_m$, with $v$ the kinematic fluid viscosity), for example, in the case of water jets (Villermaux & Duplat 2003). Some aspects of the modeling framework presented here thus draw inspiration from turbulence models (Kraichnan 1974, Villermaux 2019). However, there are also key differences. Unlike turbulent and sine flows, which can occur in unbounded domains, porous media flows are inherently confined by solid boundaries, which may significantly impact mixing dynamics (Gouillart et al. 2007, Hamada et al. 2020). Furthermore, velocity fluctuations along streamlines— which drive mixing and dispersion—are primarily controlled by the static structure of the solid
a Turbulent flow

<figure><img src="image_22.png" /></figure>

b Chaotic sine flow

<figure><img src="image_23.png" /></figure>

c Porous media flow

<figure><img src="image_24.png" /><figcaption>Figure 4</figcaption></figure>

Mixing of a tracer dye in three different chaotic flows. The fluid flows from left to right and the tracer appears as black. ( a ) Instantaneous snapshot of mixing of a black solute in a turbulent duct flow ( Re ∼ 10 $^{4}$, Sc ∼ 10 $^{3}$). Panel adapted from Le Borgne et al. (2017), who used experimental data from Villermaux & Duplat (2003). ( b ) Numerical simulation of the 2D time varying sine flow for Pe ∼ 10 $^{4}$ (Heyman et al. 2024). A cross section in the space–time domain ( y, vt ) is shown. ( c ) Scalar plume forming in a steady flow through a random bead packing imaged with refractive index matching at Pe ∼ 10 $^{4}$ (Heyman et al. 2020). Although fundamentally different, all these flows produce thin solute lamellae, below the characteristic velocity length scale, with similar mixing dynamics. Stretching and compression of lamellae accelerate their diffusive mixing ( black becomes gray ) and lead to a homogenization of the mixture. Compression also favors the aggregation of nearby lamellae.

matrix, rather than by dynamic internal flow features such as turbulent eddies. This structural control may simplify the transformation from Eulerian to Lagrangian frames, offering promising opportunities to directly relate mixing behavior to the properties of the porous medium.

4. MIXING DYNAMICS

In this section, we discuss the role of fluid stretching in driving solute mixing in the context of
flows through porous media. In the sidebar titled Pore-Scale Imaging of Mixing in Porous Media,
we highlight experimental methods to quantify solute mixing at the pore scale.

**4.1. The Lamella Representation**

A solute plume transported in a heterogeneous flow field deforms into distinct thin filamentous sheets primarily aligned with the flow direction (Figure 4). As long as they do not overlap, their mixing rate is a direct function of the interplay between their elongation and diffusion (Villermaux 2012, Le Borgne et al. 2015, Lester et al. 2016b, Souzy et al. 2018). At early times and at high Péclet numbers, the solute plume can thus be viewed as a collection of independent lamellae, which deform under the stretching action of velocity gradients and diffuse due to molecular agitation (Figures 4 and 1). Along the main compression direction, each lamella follows a 1D compression–diffusion equation in a Lagrangian coordinate system that moves with the flow (Ranz 1979). This 1D approximation significantly simplifies the derivation of analytical solutions for quantifying mixing rates, as discussed below.
# PORE-SCALE IMAGING OF MIXING IN POROUS MEDIA

## 2D Micromodels

Experimental investigations with 2D micromodels (Willingham et al. 2008, de Anna et al. 2014, Jiménez-Martínez et al. 2015, Borgman et al. 2023) have made significant contributions to advancing the understanding of pore-scale mixing dynamics. However, 2D flows do not exhibit chaotic advection inherent to 3D porous media.

## 3D Tomography

Different 3D imaging techniques have been used to observe transport processes in porous media, including X-ray tomography (Hasan et al. 2020, Velásquez-Parra et al. 2024), magnetic resonance imaging (Markale et al. 2021), and refractive index matching (Kree & Villermaux 2017, Heyman et al. 2020, Souzy et al. 2020). X-ray and magnetic resonance allow the exploration of a broad variety of materials, but the contrast and resolution of scalar gradients or reaction rates can be a limitation for observing mixing processes. In contrast, optical tomography methods based on refractive index matching are limited to transparent materials (glass, plastics, or resins) but provide high-resolution images of fluorescent tracer concentrations, gradients, and reactive fluorescent probes (Guilbert et al. 2021, Sanquer et al. 2024).

## Push-Pull Laser-Sheet Imaging

Heyman et al. (2021) developed an indirect method using laser-sheet transverse imaging of a fluorescent plume advected outside an opaque porous medium. When applied in a push-pull flow protocol, which consists of reversing the flow direction after a steady push phase, this approach isolates the irreversible component of mixing induced by chaotic advection. It offers lower noise compared to tomographic techniques while enabling the investigation of opaque materials.

### 4.1.1. Compression-diffusion equation.
We consider a solute lamella of initial length $L_0$ and width $s_0$ and current elongation $\rho = L/L_0$. We define the Lagrangian frame $(\xi, \eta)$ centered on a solute lamella, where $\xi$ is the coordinate along the elongation direction and $\eta$ is the coordinate perpendicular to it. Close to a solute lamella (in two dimensions) or sheet (in three dimensions), the velocity field can be approximated as $v_\xi = \gamma\xi$ and $v_\eta = -\gamma\eta$ (Villermaux 2019). By neglecting the concentration gradients in the elongation direction $\xi$ compared to those that develop in the compression direction $\eta$, the advection diffusion equation can be approximated as a 1D compression-diffusion equation:

$$ \frac{\partial c}{\partial t} = \gamma \eta \frac{\partial c}{\partial \eta} + D_m \frac{\partial^2 c}{\partial \eta^2}. \qquad 6. $$

This 1D-lamella approximation is valid for $\rho \ge 1$, when lamellae are significantly elongated and thus compressed in their transverse direction. This equation can be transformed into a diffusion equation using the dimensionless variables (Kraichnan 1974, Ranz 1979)

$$ \xi = \frac{\eta \rho}{s_0}, \quad \tau = \frac{D_m}{s_0^2} \int_0^t \rho^2 dt'. \qquad 7. $$

The variables $\xi$ and $\tau$ form a stretched spatiotemporal domain that evolves conjointly with fluid deformations. As discussed above, the elongation $\rho$ grows algebraically in shear flows and exponentially in chaotic flows. This leads to the following elementary solutions for a pulse and a front. Considering a lamella of solute introduced as a pulse with an initial profile $c = c_0 \exp(\eta^2/s_0^2)$, the
concentration profile across the lamella is

$$
c(\eta, t) = c_{\max}(t) \exp\left(-\frac{\eta^2}{s(t)^2}\right), \qquad 8.
$$

a Gaussian distribution of width and maximum respectively equal to

$$
s(t) = s_0 \frac{\sqrt{1 + 4\tau(t)}}{\rho(t)}, \quad c_{\max}(t) = \frac{c_0}{\sqrt{1 + 4\tau(t)}}. \qquad 9.
$$

The corresponding equation for a front of initial profile $c = c_0(1 + \erf(\eta/s_0))/2$ is

$$
c(\eta, t) = \frac{c_0}{2} \left( 1 + \erf \left( \frac{\eta}{s(t)} \right) \right). \qquad 10.
$$

The resulting concentration gradients across a lamella and a front scale as $\nabla c \sim c_{\max}/s$ and $\nabla c \sim c_0/s$, respectively. This framework hence predicts that concentration profiles across stretched elements follow similar shapes as for diffusion but with enhanced concentration gradients as a result of stretching.

**4.1.2. Mixing time and mixing scales.** As discussed above, the stretched time $\tau$ is the proper clock to describe the dilution of individual lamellae or a scalar front stirred in a flow. Note that in the absence of deformation, $\rho = 1$, the stretching time is proportional to the natural time ($\tau = (D_m/s_0^2)t$) and one recovers a simple diffusive process. For $\rho > 1$, fluid deformation is effective and mixing is accelerated by the stirring action of the flow as $\tau$ grows faster than linear with time (Equation 7): $\tau = (D_m/s_0^2)\nabla v^2(t^3/3)$ for a simple shear and $\tau = (D_m/(2\lambda s_0^2))(\exp(2\lambda t) - 1)$ for exponential stretching. At early times, $\tau < 1$, the evolution of the mixing scale is dominated by compression. It decays as $s = s_0/\rho$ (Equation 9), leading to a corresponding enhancement of concentration gradients. At $\tau = 1$, diffusive expansion balances compression, leading to a stabilization of the mixing scale. For $\tau > 1$, the mixing scale can either increase due to the diffusive expansion for decaying stretching rates or remain constant for constant stretching rates. This regime is characterized by an accelerated mixing as the maximum concentration of solute lamellae decays as $c_{\max} = c_0/\sqrt{1+4\tau}$ (Equation 9). The transition time from the compression-dominated regime to the mixing-dominated regime at $\tau \sim 1$ is the mixing time. In a shear-dominated flow, the mixing time and corresponding mixing scale are, respectively,

$$
t_m = \nabla v^{-1} \left( \nabla v \frac{s_0^2}{D_m} \right)^{1/3}, \quad s(t_m) = s_0 \left( \frac{s_0^2 \nabla v}{D_m} \right)^{-1/3} . \qquad 11.
$$

In chaotic flows, the mixing time and corresponding mixing scale are, respectively,

$$
t_m = \frac{1}{2\lambda} \log \left( \frac{s_0^2}{2s_B} \right) + 1, \quad s(t \ge t_m) = s_B = \sqrt{\frac{2D_m}{\lambda}}. \qquad 12.
$$

The scale $s_B$ is called the Batchelor scale, first identified in turbulent flows (Batchelor et al. 1959). In random bead packs, the Lyapunov exponent has been measured to be (Heyman et al. 2020) $\lambda \approx 0.2 t_a^{-1}$, such that, in principle, the Batchelor scale is smaller than the grain diameter $d$ for $Pe > 5$. This value can be considered as a lower bound for the incomplete mixing domain, where scalar gradients exist at the pore scale. Since diffusion in a bounded domain is faster than in an open domain (Hamada et al. 2020), the transition to incomplete mixing may occur at larger Péclet numbers. When $t \gg t_m$, the concentration of a solute lamella in a chaotic flow (Equation 9) decays as $c_{\max}(t) \rightarrow c_0 s_0 / s_B \cdot \rho^{-1}(t)$, showing the close relationship between mixing and the deformation kinematics.
### 4.1.3. Shear-driven versus chaotic mixing.
In 3D porous media, fluid deformation involves a complex coupling between longitudinal shear and transverse stretching (Lester et al. 2018). While shear deformation plays a central role in a range of systems (Figure 2), we focus here on chaotic advection, which is ubiquitous at the pore scale. The presented framework remains broadly applicable when shear dominates the deformation process, by substituting exponential with algebraic stretching (Le Borgne et al. 2015).

### 4.1.4. Mixing statistics in chaotic flows.
The variability of stretching rates plays an important role in the evolution of mixing statistics. It can be quantified by the variance of the stretching rates, $\sigma_{\lambda}^2$. In the limit of chaotic flows uncorrelated in time, Kraichnan (1974) has shown that $\lambda/\sigma_{\lambda}^2 = n_d/2$, where $n_d$ is the dimension of the space where deformations take place. Since stretching in 3D steady flows is comparable to stretching in 2D unsteady flows, we expect $n_d = 2$ and

$$ \sigma_{\lambda}^2 \approx \lambda. \qquad 13. $$

In porous media flows, this equality remains to be verified. In the case of time-correlated flows, lamellae tend to align with the persistent direction of stretching and $\sigma_{\lambda}^2 < \lambda$ (Meunier & Villermaux 2022). As fluid deformation often obeys an uncorrelated multiplicative process, the distribution of elongation of infinitesimal fluid elements $\rho$ is generally log-normal with parameters $\log \rho = \lambda t$ and $\sigma_{\log \rho}^2 = \sigma_{\lambda}^2 t$ (Lester et al. 2016b). Using Equation 13, the mean elongation follows $\langle \rho \rangle \approx \exp((3\lambda/2)t)$ and the growth of a material line of initial length $\ell_0$ in the porous domain is

$$ L(t) = \ell_0 \langle \rho \rangle \approx \ell_0 \exp((3\lambda/2)t). \qquad 14. $$

Equation 14 was used by Heyman et al. (2020) to infer the value of $\lambda$ from experimental images. The elongation sampled along the material line, a random variable called $\rho_L$, follows the statistics of the PDF of $\rho$, weighted by $\rho$ since elongated lamellae contribute proportionally more to the statistics sampled along material lines (Meunier & Villermaux 2010), which is also log-normal but with parameters $\log(\rho_L) = 2\lambda t$ and $\sigma_{\log(\rho_L)}^2 = \lambda t$. After $t_m$, the concentration along the material line is $c_{\max} = (s_0c_0/s_B)\rho_L^{-1}$, which follows a log-normal distribution with parameters $\log c_{\max} = -2\lambda t$ and $\sigma_{\log c_{\max}}^2 = \lambda t$ and a PDF

$$ p(\epsilon_{\max}) = \frac{1}{c_{\max}\sqrt{2\pi}\lambda t} \exp\left(-\frac{(\log c_{\max} + 2\lambda t - \log(s_0c_0/s_B))^2}{2\lambda t}\right). \qquad 15. $$

Note that as $t \to \infty$, this equation can be approximated as

$$ p(\epsilon_{\max}) \sim c_{\max}^{-3} \qquad 16. $$

over a large range of $c_{\max}$. This scaling was reported in 2D time-varying chaotic flows (Fereday & Haynes 2004) but not yet in porous media. The statistics of the concentration field $c$ can be obtained by convolving the PDF of $c_{\max}$ with the Gaussian lamella (Equation 8) with width $s(t) = s_B$. The PDF of $c$ has a diffusive tail $\sim c^{-1}$ for small concentrations $\log[s_B/(c_0s_0)] \ll -2\lambda t$ and follows the PDF of $c_{\max}$ at higher concentrations (Le Borgne et al. 2015). Since the PDF of $c_{\max}$ is log-normal (Equation 15), the moment of order $p$ follows

$$ \langle c_{\max}^p \rangle = (c_0 s_0 / s_B)^p \exp(-\gamma_p t), \qquad 17. $$

with $\gamma_p = p(2\lambda - p\lambda/2)$ when $p \le 2$. In particular, $\gamma_1 = 3\lambda/2$ and $\gamma_2 = 2\lambda$. As noted by Balkovsky & Fouxon (1999), higher moments ($p \ge 2$) are entirely controlled by the fraction of lamellae that have been the least stretched. Using Equation 16, the second moment is $\int c_{\max}^2 P_{\epsilon_{\max}} \sim \int c_{\max}^{-1}$, which is not integrable. In fact, the log-normal PDF should be truncated to $\rho \ge 1$ ($c_{\max} \le c_0s_0/s_B$), giving $\gamma_p = 2\lambda$ for all $p \ge 2$ (Heyman et al. 2021). In porous media, the evolution of scalar variance is
therefore highly sensitive to the existence of weakly stretched flow regions. These are also regions of high lamella curvature, where the 1D approximation (Equation 6) breaks (Thiffeault 2004). The spatial moments of c in a transverse plane of area A are given by integrating the lamella concentration across the stretched material line of length L:

$$
\langle c^p \rangle = A^{-1} L(t) \langle c_{\text{max}}^p \rangle s_B. \qquad 18.
$$

The mean and mean-squared concentrations follow

$$
\langle c \rangle = c_0 s_0 \ell_0 / A, \quad \langle c^2 \rangle = \frac{(c_0 s_0)^2 \ell_0}{A s_B} \exp(-(\lambda/2)t). \qquad 19.
$$

For conservative transport, 〈 c 〉 is constant as expected, while 〈 c 2 〉 decays in time due to the stretch-
ing and dilution of lamellae. Such estimates have been compared with dilution experiments in
bead packs (Heyman et al. 2020). Note that the lamella theory predicts that the decay exponent is
independent of Pe, as experimentally confirmed by Heyman et al. (2021) for Pe ≥ 10 3 .

**4.1.5. Limits of the independent lamella representation.** The 1D compression–diffusion equation (Equation 6) has been compared to experimental data (Souzy et al. 2018) and numerical simulations (Meunier & Villermaux 2010), showing accurate predictions for the evolution of concentration fields. Although the independent lamella representation captures essential aspects of the interplay between stretching and diffusion, it reaches fundamental limitations under certain conditions—for example, when lamellae deviate from 1D structures near cusps (Meunier & Villermaux 2010) or when interactions with solid boundaries become significant at low Péclet numbers (Hamada et al. 2020). Another limitation arises in the aggregation regime, where lamellae begin to overlap. Since the lamella elongation grows exponentially in a finite domain, solute lamellae necessarily overlap as soon as their width ceases to decay by compression. This occurs shortly after the mixing time when diffusion balances compression. In this regime, Heyman et al. (2021) observed that the PDF of solute concentrations takes an exponential shape at large concentrations (Figure 5d). The PDF is time independent in the rescaled variable c/σc, and the decay exponents of the concentration moments follow γp ∼ p, as expected for an exponential distribution. This behavior deviates from the prediction of the independent lamella model (Equation 15) and suggests the influence of lamella aggregation or the presence of weakly stretched lamellae with more complex 2D or 3D structures, as discussed below.

**4.2. Aggregation Dynamics**

After the mixing time, solute lamellae overlap as they keep elongating in a confined domain while maintaining a constant width s<sub>B</sub>. From geometrical constraints, the mean number of overlaps between lamellae is ⟨n(t)⟩ = L(t)s<sub>B</sub>/A(t), where L is the total length of material lines (Equation 14) and A is the fluid area occupied by the scalar plume. In chaotic flows, the number of overlaps is thus expected to scale as Pe<sup>1/2</sup> due to the scaling of s<sub>B</sub> (Equation 12). Due to the linearity of the advection–diffusion equation, the aggregated concentration field can be modeled as the superposition of individual lamellae obeying Equation 6. The aggregated concentration c<sub>agg</sub>, formed by the superposition of n stretched lamellae of maximum concentrations c<sub>max,i</sub>, is c<sub>agg</sub>(t) ≈ Σ<sub>i=1</sub><sup>n(t)</sup> c<sub>max,i</sub>(t). The statistics of c<sub>agg</sub> depend on the existence or absence of correlations between the components c<sub>max,i</sub> and between their mean values and the number of aggregations n.

**4.2.1. Random aggregation.** The simplest scenario is to consider the aggregated lamellae to be independent and to assume that their mean concentration is independent from the number of aggregations. Then, the statistics of c<sub>agg</sub> result from an *n*-convolution of the single lamella
a Longitudinal mixing and dispersion

b Transverse mixing and dispersion

c Scalar variance decay

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

d Random aggregation in concentration statistics

<figure><img src="image_28.png" /></figure>

e Self-similarity in concentration statistics

<figure><img src="image_29.png" /><figcaption>Figure 5</figcaption></figure>

Mixing by chaotic advection at the pore scale. (a) Conservative mixing experiment of a green dye in a uniform flow (bottom to top) in a 3D column made of hydrogel beads and imaged with a laser sheet in the longitudinal plane. Panel adapted with permission from Kree & Villermaux (2017). (b) Transverse cross section of a scalar plume obtained by the push-pull experimental technique (see the sidebar titled Pore-Scale Imaging of Mixing in Porous Media) in a similar bead packing (Heyman et al. 2021) showing the heterogeneity of the concentration field below the velocity scale and the transverse dispersion of the plume. (c) Evolution of the transverse scalar variance as a function of the pull time normalized by advection time in the push-pull experiment. The slope can be related to the mean stretching rate $\lambda$ through the lamella framework. Panel adapted from Heyman et al. (2021). (d) Transverse scalar concentration probability density function (PDF) as a function of distance in the experiment of Kree & Villermaux (2017), showing the expected Gamma-shape PDF (lines) that would result from an aggregation process and the exponential increase of the number of aggregations. Panel adapted with permission from Kree & Villermaux (2017). (e) Self-similar form of the push-pull scalar PDF when rescaled with the standard deviation, suggesting exponential tails in the PDF. Panel reproduced from Heyman et al. (2021).

PDF (Equation 15). Through this convolution process, the concentration PDF tends to a Gamma distribution (Villermaux & Duplat 2003), with moments

$$ \langle c_{\text{agg}} \rangle = \langle n(t) \rangle \langle c_{\text{max}}(t) \rangle = \ell_0 c_0 s_0 / A, \qquad 20. $$

$$ \sigma_{c_{\text{agg}}}^2 = \langle c_{\text{agg}} \rangle^2 / \langle n(t) \rangle \sim L^{-1}(t). \qquad 21. $$

Le Borgne et al. (2015) used the random aggregation model to predict the asymptotic dilution of a scalar plume at the Darcy scale (Figure 2c). Kree & Villermaux (2017) measured the number $n$ of convolutions needed to describe the experimental concentration PDF using a Gamma distribution
(Figure 5a,d). They found that it grows exponentially with time consistently with an exponential elongation resulting from chaotic mixing.

**4.2.2. Correlated aggregation.** Simulations of the aggregation process have been achieved via Lagrangian methods for tracking the evolution of the material line and Eulerian reconstruction of the concentration field on a fixed grid (Meunier & Villermaux 2010). Using this approach, Heyman et al. (2024) have shown that, unlike in turbulent flows, the assumption of randomness does not hold below the velocity correlation length scale in smooth chaotic flows. This is because fluid incompressibility compensates for the elongation of material lines through transverse compression, thereby accelerating the aggregation of neighboring lamellae. Thus, lamellae with high $\rho$ (i.e., low $c_{max}$) exhibit higher aggregation levels $n$ than those with low $\rho$, leading to a correlated aggregation process. Correlations lead to a variance scaling of the form $\sigma_{agg}^2 \sim L^{-\xi}(t)$, where $\xi \in [0, 1]$ is a flow-dependent exponent that varies inversely with the fractal dimension of advected material lines, $d_{\perp,1}$ (see the **Supplemental Appendix**, section 3), and thus reflects the heterogeneity of stretching. In other words, correlations in the aggregation process globally retard the homogenization of solute plumes, compared to the random case (Equation 21). In fact, the scalar decay rate is mostly insensitive to aggregation because weakly stretched filaments, which contain high concentrations, are the least aggregated. Therefore, in this regime, concentration variance can be approximated as $\sigma_{agg}^2 \sim \sigma_c^2 \sim \exp(-(\lambda/2)t)$ (Equation 19).

**4.2.3. Role of dispersion.** In porous media, the domain $A(t)$ where lamellae mix and aggregate can be limited by the presence of solid boundaries and by the slow dispersive spreading of the plume (Figure 1a). For a continuous injection into a 3D porous medium, the mixing area grows transversely as $A(t) \sim (\sigma_x^{\perp})^2 \sim \alpha_{\perp}vt$ (Figure 5a,b), while the length of the material line grows exponentially with distance. Since the filament length and the number of aggregations $n$ grow exponentially (Kree & Villermaux 2017) within an algebraically expanding dispersive area $A$, one might expect that scalar fluctuations would disappear exponentially fast within the plume. However, the continuous input of a fresh scalar field at the plume fringes can potentially sustain scalar fluctuations over time. This likely leads to a spatially variable aggregation number $n$ within the plume, with reduced aggregation near the edges compared to the center (Le Borgne et al. 2015). Although the persistence of scalar fluctuations remains an open question, the resulting concentration gradients tend to become negligible compared to the macroscopic gradient across the plume. Experimental observations (de Anna et al. 2014, Borgman et al. 2023, Sanquer et al. 2024) suggest a transition to a dispersion-controlled regime at late times. The corresponding time may be estimated by equating the macroscopic gradient $\nabla c_M \sim 1/\sqrt{\alpha_{\perp}vt}$ to the microscopic gradient $\nabla c_m \sim \exp(-\lambda t)/s_B$ (Sanquer et al. 2024). The confirmation of this transition and the characterization of the corresponding asymptotic regime require further investigation.

# 5. CONTROL OF FLUID DEFORMATION AND MIXING ON CHEMICAL AND BIOLOGICAL PROCESSES

We now discuss the consequences of the fluid deformation and mixing dynamics for chemical and biological processes in porous media (Figure 6). Fluid deformation influences reactive processes not only by enhancing solute mixing (Rolle & Le Borgne 2019, Valocchi et al. 2019) but also by shaping the spatial organization and patterns of chemical and biological species and particles within the fluid space (Tél et al. 2005, Neufeld & Hernández-Garcia 2009, Aref et al. 2017). We begin by reviewing key findings and open questions on how fluid deformation affects biomolecular reactions, which serve as prototypical models of mixing-driven processes. We then examine its role in fluid-solid interactions—such as adsorption, redox reactions, complexation,
<figure><img src="image_30.png" /><figcaption>b 2D steady reactive front</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>c 3D steady reactive front</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>d 2D transient front</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>e Fluid-solid reaction</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 6</figcaption></figure>

Impact of mixing on reactions at the pore scale. (a) Typical profile of a mixing front between two reactants *A* and *B* in the case of an instantaneous reaction (*continuous lines*), a slow reaction (*dashed lines*), and an intermediate regime (*dotted lines*). (b) Imaging of reaction rates via chemiluminescence close to a stagnation point in a micromodel in the absence and presence of 2D porous media. Panel adapted from Izumoto et al. (2023). (c) 3D refractive index imaging of a mixing-induced reaction front in a uniform flow through a porous bead pack. Panel adapted from Sanquer et al. (2024). (d) Imaging of the reaction rate produced in a sheared front through a 2D porous micromodel. Panel adapted from de Anna et al. (2014). (e) Numerical simulation of the solid–fluid adsorption reaction in a nonchaotic flow (*left*) and a chaotic flow (*right*) through a periodic bead packing. The color scale represents the concentration of the reactive solute at the same time. Panel adapted from Aquino et al. (2023).

and dissolution-precipitation—as well as in biological processes, including coexistence, preda-
tion, competition, nutrient uptake, biofilm growth, and chemotaxis. Other noteworthy processes
not addressed here—such as autocatalytic and oscillatory reactions (Izumoto et al. 2024), diffusio-
phoresis (Deseigne et al. 2014), and colloid filtration (Ouellette et al. 2008)—are also potentially
strongly influenced by the kinematics of fluid deformation in porous media. Investigating these
processes through the lens of the principles discussed below represents a promising direction for
future research.

**5.1. Mixing-Driven Reactions in the Fluid Phase**

Following experimental protocols used to study mixing-driven reactions (Figure 6a–d), we fo-
cus on bimolecular irreversible reactions A + B → C, although the approach can be generalized
to more complex reactive systems (Rolle & Le Borgne 2019, Valocchi et al. 2019). Consider a
reactive front between two species A and B, which is stretched by a heterogeneous velocity field
(Figure 6a). In the stretching regime described above, this problem can be formulated in the
independent lamella framework as a 1D compression-diffusion-reaction equation,

$$
\frac{\partial C_A}{\partial t} = \gamma(t)\eta \frac{\partial C_A}{\partial \eta} + D_m \frac{\partial^2 C_A}{\partial \eta^2} - k C_A C_B, \qquad 22.
$$

$$
\frac{\partial C_B}{\partial t} = \gamma(t)\eta \frac{\partial C_B}{\partial \eta} + D_m \frac{\partial^2 C_B}{\partial \eta^2} - k C_A C_B, \qquad 23.
$$

$$
\frac{\partial C_C}{\partial t} = \gamma(t)\eta \frac{\partial C_C}{\partial \eta} + D_m \frac{\partial^2 C_C}{\partial \eta^2} + k C_A C_B, \qquad 24.
$$

where C_A, C_B, and C_C are, respectively, the concentrations of species A, B, and C. Two dimensionless numbers can be defined for these equations: the Damköhler number Da = kC₀L₀²/D_m, with L₀ and C₀, respectively, a characteristic length scale and concentration associated with the problem, and a Péclet number associated with fluid deformation processes Pe_γ(t) = γ(t)L₀²/D_m. In the case of a shear flow, Pe_γ decreases with time as Pe_γ(t) ~ t⁻¹, while in the case of a chaotic flow, it remains constant. We denote R = kC_A C_B the local reaction rate and I = ∫Rdη the reaction intensity per unit length of the front. We also define s_r as the width of the reactive front, which can be estimated from the second spatial moment of the reaction rate across the mixing front (Bandopadhyay et al. 2017). Note that I ≈ R_s. The quantity F = C_A - C_B is conserved and thus obeys the diffusion equation in the rescaled time and space coordinate system (τ, ξ) attached to the front, with the solution given by Equation 10. Since the length of the interface increases as ρ(t), the mass of product C created per unit of time is dm_C/dt ~ ρI. The effective properties of these reactive fronts follow two distinct regimes: the kinetic-limited regime, where the characteristic reaction time is larger than the characteristic mixing time, and the mixing-limited regime in the opposite scenario.

**5.1.1. Kinetic-limited regime.** In this regime, the reaction is slow compared to mixing processes (Da ≪ Peγ), such that the concentrations of reactants, C_A and C_B, can be approximated as conservative profiles, obeying the diffusion equation in rescaled coordinates (Equation 6). Therefore, the reaction width is approximated by the conservative mixing width s and the reaction rate is derived from the product of concentrations R = kC_A C_B, where C_A and C_B are approximated by the conservative concentrations (see the **Supplemental Appendix**, section 5). This regime develops at times smaller than the reaction time t_r = 1/(kC₀). Under sustained stretching rates, where s → s_B, this regime is sustained at any time if the compression rate is greater than the reaction rate γ⁻¹ ≪ t_r (Izumoto et al. 2023).

**5.1.2. Mixing-limited regime.** The opposite limit is reached when the reaction is fast compared to mixing (Da ≫ Pe<sub>γ</sub>). In the limit of instantaneous reactions, the interpenetration of reactants does not occur: C<sub>A</sub>C<sub>B</sub> = 0 everywhere except at the front, where F = C<sub>A</sub> − C<sub>B</sub> = 0. In that case, the reaction width is infinitesimal and the reaction rate per unit length of the front is equal to the diffusive flux at the position of the front (de Anna et al. 2014), I ~ |∂F/∂η| ~ C<sub>0</sub>/s(t), where s(t) is given by Equation 9. The total mass of product created per unit time is thus equal to the length of the front, corresponding to the isocontour at F = 0, multiplied by the reaction intensity I. These results have been experimentally validated using fast chemiluminescence reactions in 2D porous micromodels (de Anna et al. 2014) (Figure 6d). The assumption of fast mixing-limited reactions allows the derivation of approximate solutions for reactive transport dynamics from conservative transport simulations (Sole-Mari et al. 2022, Shafabakhsh et al. 2024). Because reactants cannot interpenetrate, the local concentration of the product C<sub>e</sub> is estimated directly from the concentration of the conserved quantity F, C<sub>C</sub> = (1 − |F|)/2, where we have assumed
equal initial concentrations $C_A = C_B = 1$ (Shafabakhsh et al. 2024). For the more general case of mixing-limited reactions with finite Damköhler numbers, approximate solutions have been proposed, based on a linearization of the conservative profile F near the front origin and neglecting nonlinear reaction terms (Gálfi & Rácz 1988, Larralde et al. 1992). The scaling laws for $s_r$, $R$, and $I$ in this regime can be obtained from the following considerations (Gálfi & Rácz 1988): (a) The intensity of the reaction is equal to the diffusion flux and therefore independent of the reaction kinetics, $I \sim D_m \partial C_A / \partial \eta \sim D_m C_0 / s$; (b) the maximum reaction rate is $R = k C_A C_B \sim k C_A^2$ for a symmetric front; and (c) the typical value of $C_A$ within a reactive zone of length $s_r$ is $C_A \sim s_r \partial C_A / \partial \eta \sim s_r C_0 / s$. Combining (b) and (c) leads to $R \sim k C_0^2 s_r^2 / s^2$. Using $I = s_r R$ and (a), we obtain the equality $k C_0^2 s_r^3 / s^2 \sim D_m C_0 / s$, which can be used to derive $s_r$. The resulting reactive scale, maximum reaction rate, reaction intensity per unit front length, and mass of product are

$$s_r \sim L_0 Da^{-1/3} (s/L_0)^{1/3}, \quad R \sim \frac{C_0 D_m}{L_0^2} Da^{1/3} (s/L_0)^{-4/3}, \quad I \sim \frac{C_0 D_m}{s}, \quad m_c = \ell_0 \int_0^t dt \rho I.$$

Thus, they are directly related to the conservative, time-dependent mixing scale $s(t)$ that results from the balance between compression and diffusion (Equation 9). The resulting reactive properties for uniform flows, as well as shear and exponential deformations, are summarized in the **Supplemental Appendix**, section 5. These scalings have been verified from numerical simulations of reactive fronts in shear flows (Bandopadhyay et al. 2017), from imaging of a compressed reactive front at a stagnation point in a 2D microfluidic model (Izumoto et al. 2023) (Figure 6b), and from 3D imaging of reactive fronts in chaotic flows at the pore scale (Sanquer et al. 2024) in bead packs (Figure 6c).

### 5.1.3. Limits and open questions.
The reactive front properties discussed above provide a first-order framework to understand the interplay between fluid stretching, mixing, and reaction. However, these predictions neglect some important mechanisms that need to be investigated further. One such mechanism is the interaction between different front elements through aggregation (Sanquer et al. 2024), which occurs after the mixing time, as discussed above for conservative mixing. In this regime, the length of the mixing interface stabilizes due to the competing effects of aggregation, which shortens it, and stretching, which continuously elongates it (Villermaux 2018, Farhat et al. 2025). Different hypotheses on reaction statistics in this regime are discussed by Guilbert & Villermaux (2021). Beyond the simple bimolecular reactions discussed above, other key processes requiring further investigation include the interplay between mixing and nonlinear reactions (Hubert et al. 2020, Le Traon et al. 2021) and the coupling of mixing with precipitation/dissolution dynamics and patterns (Schuszter & Brau 2016, Sharma et al. 2023), both of which are expected to yield new effective reaction kinetics.

## 5.2. Solid-Fluid Reactions

Because of their large specific surface, porous media host a large range of chemical and biological reactions at the interface between fluid and solid phases, including adsorption, complexation, redox reactions, and precipitation or dissolution. Reaction kinetics are often modeled assuming well-mixed (uniform) pore-scale reactant concentrations. However, transport limitations can lead to large deviations from well-mixed estimates by reducing the access of solutes to reactive surfaces (Battiato & Tartakovskiy 2011). The comparison of chaotic and nonchaotic flows in crystalline porous media (Aquino et al. 2023) shows that pore-scale chaotic mixing significantly increases reaction efficiency compared to nonchaotic flows (Figure 6e). Stagnation points and hyperbolic manifolds, which generate chaotic mixing, accelerate transport to/from the surface, disrupting the chemical boundary layer (Turuban et al. 2019).
The link between mixing dynamics on the pore scale and fluid-solid reactions was described by Aquino et al. (2023) by combining the continuous-time random walk framework with the restart theory. Consider a reactive transport problem where a transported solute species is consumed by a reaction that occurs on solid surfaces in a porous medium of grain size *d*. The ratio of characteristic diffusion to reaction time at the solid-fluid interface is characterized by the Damköhler number $Da = t_D/t_r$, with the diffusion time $t_D = d^2/D_m$ and the reaction time $t_r = (kc_sL)^{-1}$, where $D_m$ is the diffusion coefficient, *k* is the reaction rate, and $c_s$ is the surface concentration of the reactant at the solid surface. The effective asymptotic reaction rate is defined as $k_e = |dM/dt|/M$ when $t \to \infty$, with *M* the mass of the solute in the pore space, and the reaction efficiency is $\epsilon = k_e/k_e^{\text{wm}}$ with $k_e^{\text{wm}} = dAV^{-1}t_r^{-1}$ the effective reaction rate under well-mixed conditions, where *A* is the reactive surface area per unit volution *V*. Based on the restart theory, the reaction on the surface is assumed to be limited by a fluid-mixing process in the pore space characterized by a characteristic renewal time $t_{re}$. For cases such as Poiseuille flows or nonchaotic flows, where mass transfer to the reactive surface is driven by diffusion, the reaction efficiency is a function of $Da$ but independent of $Pe$. In chaotic flows, mixing becomes faster than diffusion, with $t_{re} < t_D$ (Figure 6e), and the chaotic restart predicts that the reaction efficiency is (Aquino et al. 2023)

$$ \epsilon = (1 + Da\sqrt{2t_{re}/t_D})^{-1}. \qquad 25. $$

When the bulk of the flow is fully chaotic, the stirring rate is proportional to the mean velocity. Hence, the characteristic timescale for mixing in the bulk scales as $t_{re} \sim \lambda^{-1}Pe^{-1}$, with $\lambda$ the Lyapunov exponent, and the Péclet number is $Pe = dv/D_m$. Inserting this renewal time in Equation 25 leads to a scaling $\epsilon \sim Da^{-1}(\lambda t_a)^{1/2}Pe^{1/2}$ for $Pe > 1$. However, in porous media, the no-slip condition interferes with chaotic mixing close to walls (Gouillart et al. 2007) because low velocities near the interface become the limiting factor for advective exploration. The corresponding timescale can be determined by estimating the characteristic length scale $l_r$ of the depletion layer in the vicinity of the reactive surface (Aquino et al. 2023). The characteristic time to exit from a layer of size $l_r$ from the surface by advection is $t_a \sim d(\nabla v l_r)^{-1}$, with $\nabla v$ the average velocity gradient at the pore scale. The characteristic diffusion time to transport reactant to the surface from a distance $l_r$ is $t_D \sim l_r^2/D$. Equating these two characteristic times leads to the characteristic size of the reactive layer $l_r \sim (dD/\nabla v)^{-1/3} \sim Pe^{-1/3}$, which is consistent with Lévêque’s (1928) scaling for mass or thermal transfer in shear flows. The characteristic timescale for mixing across this layer is thus $t_{re} \sim l_r^2/D \sim Pe^{-2/3}$. Inserting this restart time in Equation 25 leads to $\epsilon \sim Da^{-1}Pe^{1/3}$ for $Pe > 1$. This scaling was verified from numerical simulations of reactive transport with fluid-solid reaction in porous media (Aquino et al. 2023), confirming that the effective solid-fluid reaction rate is governed by the interplay between chaotic mixing in the pore space and shear close to the walls. In summary, the fluid-solid reaction efficiencies expected for different porous media flows are $\epsilon \sim 1$ (kinetics limited), $\epsilon \sim Da^{-1}$ [nonchaotic (e.g., Poiseuille) flow], $\epsilon \sim Da^{-1}(\lambda t_a)^{1/2}Pe^{1/2}$ (chaotic flow with slip), and $\epsilon \sim Da^{-1}Pe^{1/3}$ (chaotic flow with no slip) (Aquino et al. 2023).

## 5.3. Microbial Processes

As for the reactive processes described above, we argue that fluid deformation and mixing should play a major role in triggering or sustaining microbial dynamics in porous media. Microbial processes have traditionally been studied under well-mixed or diffusion-dominated conditions, assumed to prevail at the microscale. However, the discovery that chaotic advection in porous media can sustain the existence of heterogeneous chemical landscapes on the pore scale challenges this paradigm. Theoretically, chaotic flows are known to strongly affect biological dynamics (Tél et al. 2005). However, their consequences for microbial processes in porous media are so far largely unknown.
**5.3.1. Mixing fronts as microbiological hot spots.** Microorganisms get energy for their metabolism from redox reactions, which occur through a transfer of electrons from an electron donor (e.g., dissolved minerals) to an electron acceptor (e.g., oxygen). Electron donors and acceptors are often spatially segregated, as they are produced from different sources (e.g., dissolved oxygen is abundant in shallow water while dissolved minerals are concentrated in deep anoxic water in soils and rocks). Hence, a large part of the redox reactions that drive microbial activity occur in the mixing fronts between fluids of different chemical compositions (McClain et al. 2003). Field observations have shown that fluid mixing-driven microbial hot spots occur at groundwater-surface water interfaces (Stegen et al. 2016) and in deep fractured rock environments (Bochet et al. 2020). The dynamics of reactive fronts discussed above are therefore directly relevant to microbial processes, with the added complexity that reactive transport equations must be coupled with models describing bacterial transport and growth. Addressing this challenge would offer the potential to establish a relationship between mixing dynamics and biomass production in porous media.

**5.3.2. Microbial processes on solid surfaces sustained by fluid mixing.** Bacteria often establish surface-attached communities, forming biofilms that provide numerous advantages (Hall-Stoodley et al. 2004). In this sessile state, transport processes within the pore space play a key role in their environment, facilitating the renewal of nutrients, as well as the transport of metabolic products and signaling molecules (Nguyen et al. 2021). As a result, mixing dynamics and their effects on fluxes at fluid-solid interfaces should shape both microbial habitats and resource availability. The scaling laws described above apply to first-order surface reaction kinetics, which may serve as a crude approximation for microbial uptake. However, bacterial metabolism is generally nonlinear and is often modeled using Michaelis–Menten dynamics. The interaction between mixing processes and such nonlinear reactions should give rise to new effective kinetics (Hubert et al. 2020), which could be established by extending the framework presented here.

**5.3.3. Coexistence, competition, and predation.** Beyond mixing, fluid deformation also governs the spatial distribution of microorganisms by interacting with traits such as motility, size, chirality, and shape anisotropy (Lauga 2016). Theoretical and computational studies suggest that chaotic flows can concentrate self-propelled particles along fractal manifolds (Neufeld & Hernández-Garcia 2009), though the implications of this phenomenon for bacterial dispersion, mixing, and attachment are not known. In contrast to well-mixed systems—where competitive exclusion typically favors the most optimally adapted species—the spatial heterogeneity generated by chaotic advection may promote species coexistence by enabling distinct populations to occupy and exploit separate ecological niches (Károlyi et al. 2000, Tél et al. 2005). Understanding how fluid deformation kinematics in porous media influence coexistence, competition, and predation hence remains a key open question for future research.

**5.3.4. Chemotaxis.** A key interaction between mixing dynamics and microbial processes arises from the ability of bacteria to sense and respond to chemical gradients—a process known as chemotaxis (Keller & Segel 1971, Taylor & Stocker 2012, Bhattacharjee et al. 2021, Gao et al. 2023). While diffusion and dispersion processes can only dissipate chemical gradients, fluid stretching provides a mechanism to amplify and sustain them, as discussed in Section 4.1.1. Shear deformation has been shown to enhance chemotaxis in the presence of shear deformation within 2D porous media (de Anna et al. 2021). Extending this understanding to chaotic mixing in 3D porous media—and its impact on chemotactic behavior—remains a key open question that could be effectively addressed using the framework presented here. Furthermore, under fluid flow in porous media, rheotaxis effects induced by the interaction of fluid shear near solid surfaces and bacterial physical characteristics, including finite size, chirality, and anisotropic shape (Marcos
et al. 2012, Lauga 2016), may interfere with chemotaxis. Determining the relative importance of chemotaxis and rheotaxis in porous media remains an unresolved and significant research question.

# 6. CONCLUSIONS AND PERSPECTIVES

This review synthesizes recent advances in experimental observations and theoretical descriptions of mixing dynamics and mixing-driven processes in porous media from the perspective of fluid deformation. We review experimental challenges, assess current understanding, identify key open questions, and highlight the strengths and limitations of this framework compared to traditional dispersion-based models.

Experimental, numerical, and theoretical advances have made it clear that fluid deformation plays a central role in driving mixing dynamics in porous media. Due to the opacity and complexity of these confined environments, only recently have new experimental methods enabled the analysis of fluid deformation kinematics. As a result, the governing laws and regimes of fluid deformation and resulting mixing dynamics are still being explored. As formalized in the lamella framework, mixing properties are controlled by the elongation and aggregation statistics of stretched solute lamellae, highlighting the need to quantify not only average deformation properties but also their fluctuations and spatial correlations. While fluid deformation processes are now relatively well-understood in idealized porous media, characterizing them across the diverse range of natural and engineered porous materials and transient flow regimes remains an open challenge.

The interplay between fluid deformation, mixing, and fluid-fluid or fluid-solid reactions gives rise to emergent reactive patterns and reaction kinetics. Future research should expand on this approach to explore how fluid deformation and mixing interact with nonlinear reactions—common in many biochemical systems—including adsorption, redox processes, complexation, dissolution-precipitation, and autocatalytic or oscillatory reactions.

The widespread occurrence of chaotic advection in porous media challenges prevailing assumptions about microbial processes in porous media, revealing a mechanism that sustains heterogeneous chemical landscapes at the microscale. Beyond mixing, fluid deformation kinematics also shape the spatial distributions of microorganisms by interacting with their motility, finite size, chirality, and anisotropic morphology. The interplay between fluid deformation, mixing, and microbial activity—including chemotaxis, colonization, biofilm formation, coexistence, and competition—thus opens new avenues for exploring microbial dynamics in porous environments.

## SUMMARY POINTS

1. Recent experimental breakthroughs have opened a new window on fluid deformation and mixing in porous media by allowing a direct visualization of these dynamics in 3D porous media.

2. Porous media flow properties including spatial structure, heterogeneity, dimensionality, and temporal dynamics control stretching kinematics from the pore to the Darcy scales. Deformation modes range from algebraic elongation in shear-dominated flows to exponential elongation under chaotic advection.

3. Fluid deformation is a key driver for mixing processes in porous media, as it enhances concentration gradients and maintain fluctuations at small scales, while diffusion and dispersion can only dissipate them.
4. Deformation-mediated mixing and reactions alter reaction patterns and effective kinetics compared to dispersion-based models.

5. Fluid deformation and mixing control microbial processes in porous media by shaping the chemical landscapes experienced by microorganisms and by influencing their spatial distribution through interactions with their motility, finite size, chirality, and anisotropic shape.

## FUTURE ISSUES

1. General fluid deformation laws in porous and fractured material should be established. Fluid deformation kinematics have been established only for simple porous architectures. The concepts presented here provide a framework to investigate them for different porous media structures (e.g., rocks, fractured media, fibers, membranes, batteries, biological tissues) and unsteady flows relevant to important environmental, engineering, and health applications.

2. Mixing dynamics should be upscaled from the pore to the Darcy scales. The description of mixing through the lens of fluid deformation processes has significantly advanced the characterization, understanding, and prediction of mixing in porous media, extending beyond the scope of traditional dispersion models. Encapsulating this understanding into effective transport models at the scale of practical applications is a key outstanding challenge for future research and knowledge transfer.

3. Reaction patterns and kinetics arising from fluid deformation and mixing in porous media should be determined. Linking fluid deformation and mixing in porous media with nonlinear reactions—such as adsorption, redox processes, complexation, dissolution-precipitation, and autocatalytic or oscillatory reactions—opens new avenues for exploring a range of fluid-borne chemical processes.

4. Interactions between deformation kinematics, mixing, and microbial activity in porous media should be unraveled. The interaction between fluid deformation, mixing and biological activity—spanning chemotaxis, microbial colonization, biofilm formation, coexistence, and competition—offers new research directions and may advance our understanding of microbial dynamics in porous environments.

## DISCLOSURE STATEMENT

The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.

## ACKNOWLEDGMENTS

The authors have received funding from the European Research Council [grant 101042466 (CHORUS)] and from the European program H2020-MSCA-ITN [grant 956457 (COPERMIX)]. We are also grateful to Daniel Lester and Marco Dentz for their advice on the manuscript and to Howard Stone for reviewing the manuscript. We thank the students, postdocs, and collaborators who have made significant contributions to the development of the
framework presented here—through experiments, simulations, and theoretical work—including T. Aquino, A. Bandopadhyay, O. Borgman, P. de Anna, K. de Vriendt, M. Dentz, S. Izumoto, J. Jimenez-Martinez, D. Lester, G. Linga, Y. Méheust, B. Metzger, H. Sanquer, P. Shafabakhsh, M. Souzy, R. Turuban, and E. Villermaux.

## LITERATURE CITED

Aquino T, Le Borgne T, Heyman J. 2023. Fluid-solid reaction in porous media as a chaotic restart process. *Phys. Rev. Lett.* 130(26):264001

Aref H, Blake JR, Budišić M, Cardoso SS, Cartwright JH, et al. 2017. Frontiers of chaotic advection. *Rev. Mod. Phys.* 89(2):025007

Balkovsky E, Fouxon A. 1999. Universal long-time properties of Lagrangian statistics in the Batchelor regime and their application to the passive scalar problem. *Phys. Rev. E* 60(4):4164-74

Bandopadhyay A, Davy P, Le Borgne T. 2018. Shear flows accelerate mixing dynamics in hyporheic zones and hillslopes. *Geophys. Res. Lett.* 45(21):11659-68

Bandopadhyay A, Le Borgne T, Méheust Y, Dentz M. 2017. Enhanced reaction kinetics and reactive mixing scale dynamics in mixing fronts under shear flow for arbitrary Damköhler numbers. *Adv. Water Resour.* 100:78-95

Batchelor GK, Howells ID, Townsend AA. 1959. Small-scale variation of convected quantities like temperature in turbulent fluid part 2. The case of large conductivity. *J. Fluid Mech.* 5(1):134-39

Battiato I, Tartakovskiy D. 2011. Applicability regimes for macroscopic models of reactive transport in porous media. *J. Contam. Hydrol.* 120:18-26

Berkowitz B, Cortis A, Dentz M, Scher H. 2006. Modeling non-Fickian transport in geological formations as a continuous time random walk. *Rev. Geophys.* 44(2):RG2003

Bhattacharjee T, Amchin DB, Ott JA, Kratz F, Datta SS. 2021. Chemotactic migration of bacteria in porous media. *Biophys. J.* 120(16):3483-97

Bochet O, Bethencourt L, Dufresne A, Farasin J, Pédrot M, et al. 2020. Iron-oxidizer hotspots formed by intermittent oxic-anoxic fluid mixing in fractured rocks. *Nat. Geosci.* 13(2):149-55

Borgman O, Turuban R, Géraud B, Le Borgne T, Méheust Y. 2023. Solute front shear and coalescence control concentration gradient dynamics in porous micromodel. *Geophys. Res. Lett.* 50(5):e2022GL101407

Browne CA, Datta SS. 2024. Harnessing elastic instabilities for enhanced mixing and reaction kinetics in porous media. *PNAS* 121(29):e2320962121

Davy P, Le Goc R, Darcel C, Pinier B, Selroos JO, Le Borgne T. 2024. Structural and hydrodynamic controls on fluid travel time distributions across fracture networks. *PNAS* 121(47):e2414901121

de Anna P, Jimenez-Martinez J, Tabuteau H, Turuban R, Le Borgne T, et al. 2014. Mixing and reaction kinetics in porous media: an experimental pore scale quantification. *Environ. Sci. Technol.* 48(1):508-16

de Anna P, Pahlavan AA, Yawata Y, Stocker R, Juanes R. 2021. Chemotaxis under flow disorder shapes microbial dispersion in porous media. *Nat. Phys.* 17(1):68-73

de Dreuzy J-R, Carrera J, Dentz M, Le Borgne T. 2012. Asymptotic dispersion for two-dimensional highly heterogeneous permeability fields under temporally fluctuating flow. *Water Resour. Res.* 48(1):W01532

De Vriendt K, Le Borgne T, Pool M, Dentz M. 2022. Subsurface mixing dynamics across the salt-freshwater interface. *Geophys. Res. Lett.* 49(7):e2021GL097094

De Wit A. 2020. Chemo-hydrodynamic patterns and instabilities. *Annu. Rev. Fluid Mech.* 52:531-55

Delgado JMPQ. 2007. Longitudinal and transverse dispersion in porous media. *Chem. Eng. Res. Des.* 85(9):1245-52

Dentz M, Hidalgo JJ, Lester D. 2023. Mixing in porous media: concepts and approaches across scales. *Transp. Porous Media* 146(1):5-53

Dentz M, Le Borgne T, Englert A, Bijeljic B. 2011. Mixing, spreading and reaction in heterogeneous media: a brief review. *J. Contam. Hydrol.* 120-121(C):1-17

Dentz M, Lester DR, Le Borgne T, de Barros FPJ. 2016. Coupled continuous-time random walks for fluid stretching in two-dimensional heterogeneous media. *Phys. Rev. E* 94(6):061102
Deseigne J, Cottin-Bizonne C, Stroock AD, Bocquet L, Ybert C. 2014. How a "pinch of salt" can tune chaotic mixing of colloidal suspensions. *Soft Matter* 10(27):4795–99

Farhat S, Bolster D, Sole-Mari G. 2025. Linking mixing interface deformation to concentration gradients in porous media. *Phys. Rev. Fluids* 10(2):024501

Fereday DR, Haynes PH. 2004. Scalar decay in two-dimensional chaotic advection and Batchelor-regime turbulence. *Phys. Fluids* 16(12):4359–70

Gálfi L, Rácz Z. 1988. Properties of the reaction front in an $A + B \rightarrow C$ type reaction-diffusion process. *Phys. Rev. A* 38(6):3151–54

Gao B, Wang X, Ford RM. 2023. Chemotaxis along local chemical gradients enhanced bacteria dispersion and PAH bioavailability in a heterogenous porous medium. *Sci. Total Environ.* 859:160004

Gouillart E, Kuncio N, Dauchot O, Dubrulle B, Roux S, Thiffeault JL. 2007. Walls inhibit chaotic mixing. *Phys. Rev. Lett.* 99(11):114501

Gramling CM, Harvey CF, Meigs LC. 2002. Reactive transport in porous media: a comparison of model prediction with laboratory visualization. *Environ. Sci. Technol.* 36(11):2508–14

Guilbert E, Almarcha C, Villermaux E. 2021. Chemical reaction for mixing studies. *Phys. Rev. Fluids* 6(11):114501

Guilbert E, Villermaux E. 2021. Chemical reactions rectify mixtures composition. *Phys. Rev. Fluids* 6(11):L112501

Hall-Stoodley L, Costerton JW, Stoodley P. 2004. Bacterial biofilms: from the natural environment to infectious diseases. *Nat. Rev. Microbiol.* 2(2):95–108

Hamada M, Cueto-Felgueroso L, de Anna P. 2020. Diffusion limited mixing in confined media. *Phys. Rev. Fluids* 5(12):124502

Hasan S, Niasar V, Karadimitriou NK, Godinho JR, Vo NT, et al. 2020. Direct characterization of solute transport in unsaturated porous media using fast X-ray synchrotron microtomography. *PNAS* 117(38):23443–49

Haynes PH, Vanneste J. 2005. What controls the decay of passive scalars in smooth flows? *Phys. Fluids* 17(9):097103

Heyman J, Lester D, Le Borgne T. 2021. Scalar signatures of chaotic mixing in porous media. *Phys. Rev. Lett.* 126(3):034505

Heyman J, Lester DR, Turuban R, Méheyst Y, Le Borgne T. 2020. Stretching and folding sustain microscale chemical gradients in porous media. *PNAS* 117(24):13359–65

Heyman J, Villermaux E, Davy P, Le Borgne T. 2024. Mixing as a correlated aggregation process. *J. Fluid Mech.* 992:A6

Hidalgo JJ, Dentz M, Cabeza Y, Carrera J. 2015. Dissolution patterns and mixing dynamics in unstable reactive flow. *Geophys. Res. Lett.* 42(15):6357–64

Hinch EJ. 1999. Mixing: turbulence and chaos—an introduction. In *Mixing: Chaos and Turbulence*, ed. H Chaté, E Villermaux, J-M Chomaz. Springer

Hubert A, Aquino T, Tabuteau H, Méheyst Y, Le Borgne T. 2020. Enhanced and non-monotonic effective kinetics of solute pulses under Michaelis–Menten reactions. *Adv. Water Resour.* 146:103739

Hulin JP, Plona TJ. 1989. "Echo" tracer dispersion in porous media. *Phys. Fluids A* 1(8):1341–47

Hyman JD, Jiménez-Martnez J. 2018. Dispersion and mixing in three-dimensional discrete fracture networks: nonlinear interplay between structural and hydraulic heterogeneity. *Water Resour. Res.* 54(5):3243–58

Izumoto S, Escala DM, Heyman J, Le Borgne T, De Wit A. 2024. Control of chemical waves by fluid stretching and compression. *Phys. Rev. Lett.* 133(21):218001

Izumoto S, Heyman J, Huisman JA, De Vriendt K, Soulaine C, et al. 2023. Enhanced mixing and reaction in converging flows: theory and pore-scale imaging. *Water Resour. Res.* 59(8):e2023WR034749

Jha B, Cueto-Felgueroso L, Juanes R. 2011. Fluid mixing from viscous fingering. *Phys. Rev. Lett.* 106(19):194502

Jiménez-Martínez J, Anna PD, Tabuteau H, Turuban R, Borgne TL, Méheyst Y. 2015. Pore-scale mechanisms for the enhancement of mixing in unsaturated porous media and implications for chemical reactions. *Geophys. Res. Lett.* 42(13):5316–24

Jules V, Lajeunesse E, Devauchelle O, Guérin A, Jaupart C, Lagrée PY. 2021. Flow and residence time in a two-dimensional aquifer recharged by rainfall. *J. Fluid Mech.* 917:A13
Kalda J. 2000. Simple model of intermittent passive scalar turbulence. *Phys. Rev. Lett.* 84(3):471-74

Károlyi G, Péntek Á, Scheuring I, Tél T, Toroczkai Z. 2000. Chaotic flow: the physics of species coexistence. *PNAS* 97(25):13661-65

Keller EF, Segel LA. 1971. Model for chemotaxis. *J. Theor. Biol.* 30(2):225-34

Kitanidis PK. 1994. The concept of the dilution index. *Water Resour. Res.* 30(7):2011-26

Kraichnan RH. 1974. Convection of a passive scalar by a quasi-uniform random straining field. *J. Fluid Mech.* 64(4):737-62

Kree M, Villermaux E. 2017. Scalar mixtures in porous media. *Phys. Rev. Fluids* 2(10):104502

Larralde H, Araujo M, Havlin S, Stanley HE. 1992. Reaction front for A + B → C diffusion-reaction systems with initially separated reactants. *Phys. Rev. A* 46(2):855-59

Lauga E. 2016. Bacterial hydrodynamics. *Annu. Rev. Fluid Mech.* 48:105-30

Le Borgne T, Dentz M, Bolster D, Carrera J, de Dreuzy J-R, Davy P. 2010. Non-Fickian mixing: temporal evolution of the scalar dissipation rate in heterogeneous porous media. *Adv. Water Resour.* 33(12):1468-75

Le Borgne T, Dentz M, Davy P, Bolster D, Carrera J, et al. 2011. Persistence of incomplete mixing: a key to anomalous transport. *Phys. Rev. E* 84(1):015301

Le Borgne T, Dentz M, Villermaux E. 2013. Stretching, coalescence, and mixing in porous media. *Phys. Rev. Lett.* 110(20):204501

Le Borgne T, Dentz M, Villermaux E. 2015. The lamellar description of mixing in porous media. *J. Fluid Mech.* 770:458-98

Le Borgne T, Huck PD, Dentz M, Villermaux E. 2017. Scalar gradients in stirred mixtures and the deconstruction of random fields. *J. Fluid Mech.* 812:578-610

Le Traon C, Aquino T, Bouchez C, Maher K, Le Borgne T. 2021. Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction. *Geochim. Cosmochim. Acta* 306:189-209

Lester D, Dentz M. 2025. Line stretching in random flows. Preprint, arXiv:2504.17982v1 [physics.flu-dyn]

Lester DR, Dentz M, Le Borgne T. 2016a. Chaotic mixing in three-dimensional porous media. *J. Fluid Mech.* 803:144-74

Lester DR, Dentz M, Le Borgne T, de Barros FP. 2018. Fluid deformation in random steady three-dimensional flow. *J. Fluid Mech.* 855:770-803

Lester DR, Dentz M, Singh P, Bandopadhyay A. 2023. Under what conditions does transverse macrodispersion exist in groundwater flow? *Water Resour. Res.* 59(3):e2022WR033059

Lester DR, Heyman J, Méheust Y, Le Borgne T. 2025. A unified theory of pore-scale chaotic advection. *J. Fluid Mech.* 1017:A13

Lester DR, Metcalfe G, Trefry MG. 2013. Is chaotic advection inherent to porous media flow? *Phys. Rev. Lett.* 111(17):174101

Lester DR, Trefry MG, Metcalfe G. 2016b. Chaotic advection at the pore scale: mechanisms, upscaling and implications for macroscopic transport. *Adv. Water Resour.* 97:175-92

Lester DR, Trefry MG, Metcalfe G, Dentz M. 2024. Is chaotic advection inherent to heterogeneous Darcy flow? Preprint, arXiv:2412.05419v1 [physics.flu-dyn]

Lévêque A. 1928. Les lois de la transmission de chaleur par convection. *Ann. Mines* 13:201-99, 305-62, 381-415

Marcos, Fu HC, Powers TR, Stocker R. 2012. Bacterial rheotaxis. *PNAS* 109(13):4780-85

Markale I, Cimmarusti GM, Britton MM, Jiménez-Martnez J. 2021. Phase saturation control on mixing-driven reactions in 3D porous media. *Environ. Sci. Technol.* 55(13):8742-52

Mathiesen J, Linga G, Misztal M, Renard F, Le Borgne T. 2023. Dynamic fluid connectivity controls solute dispersion in multiphase porous media flow. *Geophys. Res. Lett.* 50(16):e2023GL105233

McClain ME, Boyer EW, Dent CL, Gergel SE, Grimm NB, et al. 2003. Biogeochemical hot spots and hot moments at the interface of terrestrial and aquatic ecosystems. *Ecosystems* 6:301-12

Meunier P, Villermaux E. 2010. The diffusive strip method for scalar mixing in two dimensions. *J. Fluid Mech.* 662:134-72

Meunier P, Villermaux E. 2022. The diffuselet concept for scalar mixing. *J. Fluid Mech.* 951:A33

Neufeld Z, Hernández-Garcia E. 2009. *Chemical and Biological Processes in Fluid Flows: A Dynamical Systems Approach*. Imperial College Press
Neuman SP, Tartakovsky DM. 2009. Perspective on theories of non-Fickian transport in heterogeneous media. *Adv. Water Resour.* 32(5):670–80

Neupauer RM, Meiss JD, Mays DC. 2014. Chaotic advection and reaction during engineered injection and extraction in heterogeneous porous media. *Water Resour. Res.* 50(2):1433–47

Nguyen J, Lara-Gutiérrez J, Stocker R. 2021. Environmental fluctuations and their effects on microbial communities, populations and individuals. *FEMS Microbiol. Rev.* 45(4):fuaa068

Ottino JM. 1990. Mixing, chaotic advection, and turbulence. *Annu. Rev. Fluid Mech.* 22:207–54

Ouellette NT, O'Malley PJJ, Gollub JP. 2008. Transport of finite-sized particles in chaotic flow. *Phys. Rev. Lett.* 101(17):174504

Oxaal U, Flekkøy EG, Feder J. 1994. Irreversible dispersion at a stagnation point: experiments and lattice Boltzmann simulations. *Phys. Rev. Lett.* 72(22):3514–17

Pool M, Dentz M, Post VE. 2016. Transient forcing effects on mixing of two fluids for a stable stratification. *Water Resour. Res.* 52(9):7178–97

Puyguiraud A, Gouze P, Dentz M. 2021. Pore-scale mixing and the evolution of hydrodynamic dispersion in porous media. *Phys. Rev. Lett.* 126(16):164501

Puyguiraud A, Perez LJ, Hidalgo JJ, Dentz M. 2020. Effective dispersion coefficients for the upscaling of pore-scale mixing and reaction. *Adv. Water Resour.* 146:103782

Ranz WE. 1979. Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows. *AlChE J.* 25(1):41–47

Rolle M, Le Borgne T. 2019. Mixing and reactive fronts in the subsurface. *Rev. Mineral. Geochem.* 85(1):111–42

Sanquer H, Heyman J, Hanna K, Le Borgne T. 2024. Microscale chaotic mixing as a driver for chemical reactions in porous media. *Environ. Sci. Technol.* 58(20):8899–908

Schuszter G, Brau F. 2016. Calcium carbonate mineralization in a confined geometry. *Environ. Sci. Technol.* 3(4):156–59

Shafabakhsh P, Le Borgne T, Renard F, Linga G. 2024. Resolving pore-scale concentration gradients for transverse mixing and reaction in porous media. *Adv. Water Resour.* 192:104791

Sharma RP, Deng J, Kang PK, Szymczak P. 2023. Effects of mixing at pore intersections on large-scale dissolution patterns and solute transport. *Geophys. Res. Lett.* 50(21):e2023GL105183

Sole-Mari G, Bolster D, Fernández-Garcia D. 2022. A closer look: high-resolution pore-scale simulations of solute transport and mixing through porous media columns. *Transp. Porous Media* 146:85–11

Souzy M, Lhuissier H, Méheust Y, Le Borgne T, Metzger B. 2020. Velocity distributions, dispersion and stretching in three-dimensional porous media. *J. Fluid Mech.* 891:A16

Souzy M, Zaier I, Lhuissier H, Le Borgne T, Metzger B. 2018. Mixing lamellae in a shear flow. *J. Fluid Mech.* 838:R3

Stegen JC, Fredrickson JK, Wilkins MJ, Konopka AE, Nelson WC, et al. 2016. Groundwater-surface water mixing shifts ecological assembly processes and stimulates organic carbon turnover. *Nat. Commun.* 7(1):11237

Stroock AD, Dertinger SK, Ajdari A, Mezic I, Stone HA, Whitesides GM. 2002. Chaotic mixer for microchannels. *Science* 295(5555):647–51

Taylor JR, Stocker R. 2012. Trade-offs of chemotactic foraging in turbulent water. *Science* 338(6107):675–79

Tél T, de Moura A, Grebogi C, Károlyi G. 2005. Chemical and biological activity in open flows: a dynamical system approach. *Phys. Rep.* 413(2–3):91–196

Thiffeault JL. 2004. Stretching and curvature of material lines in chaotic flows. *Physica D* 198(3):169–81

Trefry M, Lester D. 2020. Lagrangian complexity persists with multimodal flow forcing in compressible porous systems. *Transp. Porous Media* 135(3):555–86

Trefry MG, Lester DR, Metcalfe G, Ord A, Regenauer-Lieb K. 2012. Toward enhanced subsurface intervention methods using chaotic advection. *J. Contam. Hydrol.* 127(1–4):15–29

Turuban R, Lester DR, Heyman J, Borgne TL, Méheust Y. 2019. Chaotic mixing in crystalline granular media. *J. Fluid Mech.* 871:562–94

Turuban R, Lester DR, Le Borgne T, Méheust Y. 2018. Space-group symmetries generate chaotic fluid advection in crystalline granular media. *Phys. Rev. Lett.* 120(2):024501

Valocchi AJ, Bolster D, Werth CJ. 2019. Mixing-limited reactions in porous media. *Transp. Porous Media* 130(1):157–82
Velásquez-Parra A, Marone F, Griffa M, Jiménez-Martnez J. 2024. Chaotic transport of solutes in unsaturated porous media. *Environ. Sci. Technol.* 58(28):12643-52

Villermaux E. 2012. Mixing by porous media. *C. R. Mech.* 340(11-12):933-43

Villermaux E. 2018. On shapes and forms: population balance dynamics of corrugated stirred fronts. *C. R. Phys.* 19(5):306-15

Villermaux E. 2019. Mixing versus stirring. *Annu. Rev. Fluid Mech.* 51:245-73

Villermaux E, Duplat J. 2003. Mixing as an aggregation process. *Phys. Rev. Lett.* 91(18):184501

Werth CJ, Cirpka OA, Grathwohl P. 2006. Enhanced mixing and reaction through flow focusing in heterogeneous porous media. *Water Resour. Res.* 42(12):W12414

Willingham TW, Werth CJ, Valocchi AJ. 2008. Evaluation of the effects of porous media structure on mixing-controlled reactions using pore-scale modeling and micromodel experiments. *Environ. Sci. Technol.* 42(9):3185-93

Yang W, Chen MA, Lee SH, Kang PK. 2024. Fluid inertia controls mineral precipitation and clogging in pore to network-scale flows. *PNAS* 121(28):e2401318121

Ye Y, Chiogna G, Cirpka OA, Grathwohl P, Rolle M. 2015. Experimental evidence of helical flow in porous media. *Phys. Rev. Lett.* 115(19):194502

Yoon S, Dentz M, Kang PK. 2021. Optimal fluid stretching for mixing-limited reactions in rough channel flows. *J. Fluid Mech.* 916:A45

Zhang X, Dou Z, Hamada M, de Anna P, Jimenez-Martinez J. 2025. Enhanced reaction kinetics in stationary two-phase flow through porous media. *Environ. Sci. Technol.* 59(2):1334-43