1. Introduction
Flow and solute transport in heterogeneous porous and fractured media is determined by the spatial variability of hydraulic conductivity, which may vary over orders of magnitude (Bear Reference Bear1972). Steady flows through such media are characterised by organisation of the flow along the characteristic length scales of the medium, and broad distributions of flow speeds and thus residence times of transported scalars. These features affect both the transport and deformation of fluid elements, and the mixing and dispersion of passive and reactive scalars. They lead to non-Fickian large-scale transport properties, which manifest in broad distributions of first-passage times and nonlinear evolutions of dispersion.
Although closely related, dispersion, stirring and mixing are distinct physical processes that control solute transport and dilution (Villermaux Reference Villermaux2019). In porous or fractured media, the medium itself stirs the fluid as it permeates through the medium driven by a pressure gradient. At the pore scale, the medium stirs the fluid by the flow variability in the complex pore space (Lester et al. Reference Lester, Metcalfe and Trefry2013), which has been investigated experimentally (Souzy et al. Reference Souzy, Lhuissier, Méheust, Le Borgne and Metzger2020; Heyman et al. Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020) and numerically (Lester et al. Reference Lester, Dentz and Le Borgne2016). In the three-dimensional pore space, advection is typically chaotic (Lester et al. Reference Lester, Heyman, Méheust and Le Borgne2025) and stretching exponential. This is different on the Darcy scale. First, for steady flow in isotropic porous media, advection cannot be chaotic due to topological constraints (Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2022), that is, stretching is necessarily sub-exponential. Second, the medium stirs the bulk fluid through spatial variations in permeability, which causes flow variability and generates local velocity gradients. This leads to the deformation of fluid or solute fronts, which increases concentration gradients, promotes diffusive mass transfer and thus solute mixing, chemical reactions and biological processes.
Here, we address the question of how variations in the Darcy-scale medium heterogeneity act on the dispersion and stretching of fluid elements. Although dispersion and stretching are closely related in steady porous medium flows, there also are important differences. Longitudinal dispersion is solely governed by velocity fluctuations along streamlines (Gelhar Reference Gelhar1993; Berkowitz et al. Reference Berkowitz, Cortis, Dentz and Scher2006; Dagan Reference Dagan2012), whereas transverse dispersion is controlled by the geometry of the streamlines (Lester et al. Reference Lester, Dentz, Singh and Bandopadhyay2023). Conversely, fluid stirring is controlled by velocity gradients (shear and vorticity) normal to streamlines (Dentz et al. Reference Dentz, Lester, Le Borgne and de Barros2016b ). Although they are distinct, these processes are governed by the underlying hydraulic conductivity field, hence strong links may be expected between them. It is important to point out that the dispersion of fluid elements and solute particles are equal in the limit of purely advective transport. Solute dispersion is in general dependent on the Péclet number, which measures the ratio between the characteristic diffusion and advection time scales. Here, we focus on the dispersion and stretching of fluid elements. The two processes are illustrated in figure 1 that shows both the dispersion and deformation of a material line.
Evolution of particles initially positioned perpendicularly to the mean flow. The particle positions are shown for
$t=0$
,
$t=50\tau _v$
and
$t= 100 \tau _v$
, where
$\tau_v = \ell_Y/\langle u \rangle$
with
$\ell_Y$
the characteristic heterogeneity length scale. The colouring of the particles is proportional to the accumulated stretching experienced by each particle, defined as
$\rho$
. The background shows the longitudinal velocity as a colormap from blue to red (low to high velocity zones). The variance of the log of permeability used for this figure is
$\sigma ^2_{\kern-1pt f}=1$
and the number of particles
$N_p = 10^4$
. The figure illustrates a solute occupying space from left to right separated by the propagation front. The high deformed zones of the front will be where the solute will be more prompt to mixing and thus to reactions between the invading and invaded fluids.
Figure 1. Long description
A heat map displays particle positions and longitudinal velocity in a fluid flow. The x-axis is labeled ‘x/l_Y’ and the y-axis is labeled ‘y/l_Y’. The color scale on the right ranges from blue to red, indicating low to high velocity zones. Particles are colored proportionally to their accumulated stretching, with higher stretching shown in yellow and lower stretching in purple. The particles are initially positioned perpendicular to the mean flow and their positions are shown at different times. The background color map illustrates the longitudinal velocity, with blue representing low velocity and red representing high velocity. The variance of the log of permeability is 3.0, and the number of particles is 10000.
Stochastic modelling (Rubin Reference Rubin2003) provides a framework to systematically quantify the spatial variability of hydraulic conductivity and its impact on flow and transport. In this framework, hydraulic conductivity is typically modelled as a multi-Gaussian spatial random field with log-normal distribution of conductivity point values (Gelhar Reference Gelhar1993; Dagan Reference Dagan2012). Although criticised for the lack of connectivity of high permeability zones (Gómez-Hernández & Wen Reference Gómez-Hernández and Wen1998) and the assumption of Gaussianity of the log increments (Painter Reference Painter1996), the multi-Gaussian representation is as an idealised model both simple enough to lead to general conclusions and complex enough to disclose some sufficient conditions between medium heterogeneity and dispersion. Perturbation theory (Gelhar & Axness Reference Gelhar and Axness1983) predicts the asymptotic value of longitudinal dispersion as a function of heterogeneity through the variance of the log conductivity, the characteristic correlation scale and the mean velocity. For increasing variance, however, the linear relation between the dispersion coefficient and log-conductivity variance breaks down (Fiori et al. Reference Fiori, Janković and Dagan2003; de Dreuzy et al. Reference de Dreuzy, Beaudoin and Erhel2007; Gotovac et al. Reference Gotovac, Cvetkovic and Andricevic2009). Furthermore, the asymptotic regime may be reached only at extremely large times, and transport can in general not be described in terms of the advection–dispersion equation characterised by the average flow velocity and constant macrodispersion coefficients (Adams & Gelhar Reference Adams and Gelhar1992; Levy & Berkowitz Reference Levy and Berkowitz2003; Kang et al. Reference Kang, Le Borgne, Dentz, Bour and Juanes2015). Approaches to account for non-Fickian or anomalous large-scale dispersion include spatio-temporally non-local advection–dispersion equations (Neuman & Tartakovsky Reference Neuman and Tartakovsky2009), fractional advection–dispersion equations and multi-rate mass transfer approaches (Frippiat & Holeyman Reference Frippiat and Holeyman2008), continuous time random walk (CTRW) (Berkowitz et al. Reference Berkowitz, Cortis, Dentz and Scher2006) and time domain random walk models (Noetinger et al. Reference Noetinger, Roubinet, Russian, Le Borgne, Delay, Dentz, de Dreuzy and Gouze2016). In this work we employ the CTRW as an upscaling approach that acknowledges the spatial correlation of particle speeds and the impact of broad residence time distributions on the large-scale particle dynamics (Berkowitz et al. Reference Berkowitz, Cortis, Dentz and Scher2006; Le Borgne et al. Reference Le Borgne, Dentz and Carrera2008; Dentz et al. Reference Dentz, Kang, Comolli, Le Borgne and Lester2016a ; Comolli et al. Reference Comolli, Hakoun and Dentz2019).
The deformation of elementary fluid volumes by velocity gradients is a key process for scalar mixing, chemical reactions and biological processes in fluid flows. Stretching rates are a key input into mixing theories, along with the Péclet number, that compares diffusive and advective mass transfer over a characteristic length scale (Villermaux Reference Villermaux2019). Fluid deformation has been the subject of research over decades in unsteady turbulent flows (Cocke Reference Cocke1969; Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987; Villermaux Reference Villermaux2019). Flows through porous and fractured media, albeit often steady, can be analysed in a similar way due to the complexity induced by the medium structure. One important difference is that the steady nature of the flow imposes topological constraints in that pathlines of the flow are invariant, forming fixed streamlines that govern stirring and dispersion. One consequence of this constraint is that one of the principal directions of the deformation gradient tensor must coincide with the local velocity direction (Tabor Reference Tabor1992). Additionally, the magnitude of this principal deformation is directly proportional to the local velocity magnitude, prohibiting net stretching along the flow direction and only allowing persistent stretching and contraction in the directions normal to the flow (Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018). Furthermore, for divergence-free flow, the sum of the principal stretching rates (often characterised as Lyapunov exponents) must be zero. Hence, steady three-dimensional flow can only exhibit one unique Lyapunov exponent.
In this study we focus attention on stretching and dispersion in steady two-dimensional (2-D) flows in heterogeneous porous media. For steady 2-D flows, this stretching dynamics is even further constrained by the reduction in the geometric degrees of freedom. From the arguments above, the Lyapunov exponent in these flows must be zero (Metcalfe et al. Reference Metcalfe, Lester and Trefry2023), corresponding to sub-exponential (non-chaotic) fluid stretching. These topological constraints also limit transverse dispersion to be zero in the purely advective limit (Attinger et al. Reference Attinger, Dentz and Kinzelbach2004), as invariant streamlines are confined to 2-D surfaces, hence neighbouring streamlines cannot continually separate as they are advected. Although these constraints appear to strongly limit the applicability of studies of flow and transport in 2-D porous media, they actually are directly relevant to steady 3-D flow in an important class of heterogeneous porous materials with isotropic hydraulic conductivity fields. For steady 3-D flow in this class of porous materials, the helicity density
$\mathcal{H}(\boldsymbol{x})\equiv \boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot }(\boldsymbol{\nabla } \times \boldsymbol{u}(\boldsymbol{x}))$
, where
$\boldsymbol{u}(\boldsymbol{x})$
is the Eulerian flow velocity (Moffatt Reference Moffatt1969) is identically zero (Sposito Reference Sposito1994), meaning that they are non-chaotic and their streamlines are confined to families of 2-D streamsurfaces (Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2021) within the flow. As such, the dynamics of fluid dispersion (Lester et al. Reference Lester, Dentz, Singh and Bandopadhyay2023) and stretching Lester et al. (Reference Lester, Dentz, Bandopadhyay and Le Borgne2022) in steady 3-D flow in isotropic heterogeneous porous media is very similar to that of steady 2-D flow in heterogeneous porous media.
Coupled CTRW models have been proposed (Dentz et al. Reference Dentz, Lester, Le Borgne and de Barros2016b ; Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018) to analyse deformation in 2-D and 3-D steady random flows. While these studies show that algebraic (sub-exponential) stretching in 2-D steady random flows is driven by correlations between velocity and shear, the link between these stretching rates and medium properties for heterogeneous Darcy flow is not established. Le Borgne et al. (Reference Le Borgne, Dentz and Villermaux2015) study stretching in 2-D porous media, and propose a multiplicative stochastic process to describe stretching evolution, which, however, needs to be calibrated by observations of the stretching process itself and does not reveal its relation to the flow field of medium structure, nor its relation to the dispersion dynamics.
In this paper, we aim to close these gaps by simultaneously analysing the dispersion and deformation of fluid elements in 2-D heterogeneous porous media. The analysis of particle velocities and deformation along streamlines will provide insights into the heterogeneity-induced stochastic dynamics of these processes, and identify the relevant heterogeneity mechanisms. This allows us to link the large-scale stretching and dispersion behaviours to the underlying Eulerian flow organisation and medium structure. Beyond the fundamental interest of flow and deformation in 2-D porous media, they are relevant for example, in the context of aquifers whose extension is much larger than their depth, for flow and advection-dominated transport in rough fractures and in 3-D isotropic porous media.
The paper is organised as follows. Section 2 presents the governing equations for flow, particle transport and fluid deformation and the numerical methodology. Section 3 develops the stochastic model for fluid deformation and dispersion, and § 4 discusses how the medium heterogeneity controls the dispersion and deformation behaviours observed in heterogeneous porous media in the light of this model. Section 5 provides a summary and conclusions.
2. Methodology
We introduce here the continuum-scale flow model for porous media, the kinematic equations for particle motion and fluid deformation and the observables of interest.
2.1. Flow
We consider flow through a 2-D continuum-scale heterogeneous porous medium, which is described by the Darcy equation (Bear Reference Bear1972)
\begin{align} \boldsymbol{q}(\boldsymbol{x}) = -K(\boldsymbol{x}) \boldsymbol{\nabla } h(\boldsymbol{x}), \end{align}
where
$\boldsymbol{x} = (x_1,x_2)^\top$
is the position vector,
$\boldsymbol{q}$
denotes the Darcy velocity,
$K(\boldsymbol{x})$
hydraulic conductivity and
$h(\boldsymbol{x})$
hydraulic head. Both the fluid and the solid matrix are considered incompressible such that
$\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{q} = 0$
, which yields the groundwater flow equation
\begin{align} K(\boldsymbol{x}) \boldsymbol{\nabla }^2 h(\boldsymbol{x}) + \boldsymbol{\nabla } K(\boldsymbol{x}) \boldsymbol{\cdot }\boldsymbol{\nabla } h(\boldsymbol{x}) = 0. \end{align}
The spatially variable hydraulic conductivity field
$K(\boldsymbol{x})$
is modelled as a multi-Gaussian random field (Rubin Reference Rubin2003). Its marginal probability density function (PDF), that is the distribution of the point values of
$K(\boldsymbol{x})$
, is denoted by
$p_K(k)$
. In this work, we consider the log-normal and gamma distributions
\begin{align} p_K(k) = \frac {\exp \left [- \frac {(\ln k - \mu _{\kern-1pt f})^2}{2 \sigma _{\kern-1pt f}^2} \right ]}{k \sqrt {2 \pi \sigma _{\kern-1pt f}^2}}, && p_K(k) = \left (\frac {k}{k_0}\right )^{\alpha - 1} \frac {\exp (-k/k_0)}{\varGamma (\alpha ) k_0}, \end{align}
respectively, where
$\varGamma (\alpha )$
denotes the gamma function. The log-normal PDF is characterised by the mean
$\mu _{\kern-1pt f}$
and variance
$\sigma _{\kern-1pt f}^2$
of
$f(\boldsymbol{x}) = \ln K(\boldsymbol{x})$
and the gamma distribution by the exponent
$\alpha$
and the scale
$k_0$
.
The conductivity fields
$K(\boldsymbol{x})$
are derived from the stationary multivariate Gaussian random fields
$Y(\boldsymbol{x})$
, which are characterised by a unit Gaussian marginal density, that is, they have zero mean and unit variance
$\sigma _Y^2 = \langle Y(\boldsymbol{x})^2 \rangle = 1$
. The angular brackets denote the ensemble average over all realisations of
$Y(\boldsymbol{x})$
. The covariance function of
$Y(\boldsymbol{x})$
is denoted by
\begin{align} \langle Y(\boldsymbol{x})Y(\boldsymbol{x}') \rangle = C_Y(\boldsymbol{x} - \boldsymbol{x}'). \end{align}
In this work, we employ the isotropic Gaussian covariance function
\begin{align} C_Y(\boldsymbol{x}) = \exp \left (- \frac {\boldsymbol{x}^2}{2 \ell _Y^2} \right )\! , \end{align}
where
$\ell _Y$
denotes the correlation length. The conductivity field is obtained according to the point-wise transform
\begin{align} K(\boldsymbol{x}) = P_K^{-1}(\varPhi [Y(\boldsymbol{x})]), \end{align}
where
$P_K(k)$
is the cumulative distribution function (CDF) of
$p_K(k)$
,
$P_K^{-1}(u)$
its inverse and
$\varPhi (y)$
is the CDF of the unit Gaussian PDF. For the log-normal distribution, the map is given by
\begin{align} K(\boldsymbol{x}) = \exp [f(\boldsymbol{x})], \qquad f(\boldsymbol{x}) = \mu _{\kern-1pt f} + \sigma _{\kern-1pt f} Y(\boldsymbol{x}). \end{align}
For the gamma distribution, there is no analytical expression for the map between
$Y(\boldsymbol{x})$
and
$K(\boldsymbol{x})$
and we use the inverse or Smirnov transform method (Devroye Reference Devroye1986). Note that the effective hydraulic conductivity for log-normally distributed
$K(\boldsymbol{x})$
is equal to the geometric mean conductivity
$\exp (\mu _{\kern-1pt f})$
(Keller Reference Keller1964; Rubin Reference Rubin2003). In the following we set
$\mu _{\kern-1pt f} = 0$
, which is equivalent to normalising all the log-normal
$K(\boldsymbol{x})$
fields by their effective conductivity.
We consider constant head boundary conditions such that the mean Darcy velocity
$\langle \boldsymbol{q}(\boldsymbol{x}) \rangle$
is constant and aligned with the
$x_1$
-direction of the coordinate system. The flow velocity is denoted by
$\boldsymbol u(\boldsymbol{x}) = \boldsymbol{q}(\boldsymbol{x})/\phi$
, where
$\phi$
is porosity. Porosity is here assumed to be constant and set to one, which is equivalent to rescaling time. We refer in the following to
$\boldsymbol{u}(\boldsymbol{x})$
as the Eulerian velocity and to its magnitude
$u(\boldsymbol{x}) \equiv |\boldsymbol{u}(\boldsymbol{x})|$
as the Eulerian speed. The PDF of
$u(\boldsymbol{x})$
is denoted by
$p_e(u)$
. It is obtained by spatial sampling according to
\begin{align} p_e(u) = \lim _{V \to \infty } \frac {1}{V} \int \limits _{\varOmega } \text{d} \boldsymbol x \, \delta [u - u(\boldsymbol x)|], \end{align}
with
$\varOmega$
is the sampling domain,
$V$
its volume. Due to ergodicity, spatial sampling of
$u(x)$
is equivalent to ensemble sampling, and thus,
$p_e(u) = \langle \delta [u - u(\boldsymbol x)] \rangle$
. The ratio between the mean flow speed and the mean streamwise flow velocity defines the advective tortuosity
$\chi$
(Comolli et al. Reference Comolli, Hakoun and Dentz2019)
\begin{align} \chi = \frac {\langle u(\boldsymbol{x}) \rangle }{\langle u_1(\boldsymbol{x}) \rangle }. \end{align}
The numerical methodology for the generation of the hydraulic conductivity fields and the numerical solution for the Darcy flow field
$\boldsymbol q(\boldsymbol{x})$
are detailed in Appendix B. Figure 2 illustrates a multi-Gaussian spatial random field and the distribution of streamwise velocities for a log-normal hydraulic conductivity distribution. For the numerical simulations used in §§ 3 and 4 we consider a log-normally distributed
$K(\boldsymbol{x})$
characterised by
$\sigma _{\kern-1pt f}^2 = 0.25, 1, 3, 6$
, and a gamma-distributed
$K(\boldsymbol{x})$
with
$\alpha = 0.25, 0.5, 0.75$
.
(a) Log-normal hydraulic conductivity field with
$\sigma ^2_{\kern-1pt f} = 1$
. (b) Distribution of streamwise velocities for the corresponding log-normally distributed hydraulic conductivity field.
2.2. Particle motion
The position
$\boldsymbol{x}(t,\boldsymbol{a})$
of a particle originating at
$\boldsymbol{x}(t = 0,\boldsymbol{a}) = \boldsymbol{a}$
is described by the kinematic equation
\begin{align} \frac {\text{d} \boldsymbol{x}(t, \boldsymbol {a})}{\text{d}t} = {\boldsymbol{v}}(t,\boldsymbol{a}), \end{align}
where
${\boldsymbol{v}}(t,\boldsymbol{a}) \equiv \boldsymbol {u} [\boldsymbol{x}(t,\boldsymbol {a}) ]$
is the isochronically sampled Lagrangian velocity. The distance
$s(t,\boldsymbol{a})$
a particle travels along a streamline satisfies
\begin{align} \frac {\text{d} s(t, \boldsymbol{a})}{\text{d}t} = v(t,\boldsymbol{a}), \end{align}
where
$v(t, \boldsymbol{a}) \equiv |{\boldsymbol{v}}(t,\boldsymbol{a})|$
is the isochronically sampled Lagrangian speed. Its PDF is denoted by
$p_t(v)$
. It is equal to the PDF of Eulerian flow speeds,
$p_t(v) = p_e(v)$
(Dentz et al. Reference Dentz, Kang, Comolli, Le Borgne and Lester2016a
; Comolli et al. Reference Comolli, Hakoun and Dentz2019).
We can express the advection equation (2.10) as a function of distance
$s$
along the streamline by using the map (2.11) from
$t \to s$
, which gives
\begin{align} \frac {\text{d} \boldsymbol{x}(s,\boldsymbol{a})}{\text{d}s} = \frac { {\boldsymbol{v}}(s,\boldsymbol{a})}{ v(s,\boldsymbol{a})}, \qquad \frac {\text{d} t}{\text{d}s} = \frac {1}{ v(s,\boldsymbol{a})} , \end{align}
where
$ {\boldsymbol{v}}(s)\equiv \boldsymbol{u}[ \boldsymbol{x}(s,\boldsymbol{a})]$
is the equidistant Lagrangian velocity. Its magnitude defines the equidistance Lagrangian speed
$ v(s,\boldsymbol{a}) \equiv | {\boldsymbol{v}}(s,\boldsymbol{a})|$
. Its PDF is denoted by
$p_s(v)$
. It is related to the Eulerian speed PDF
$p_e(v)$
by flux weighting (Dentz et al. Reference Dentz, Kang, Comolli, Le Borgne and Lester2016a
; Comolli et al. Reference Comolli, Hakoun and Dentz2019)
\begin{align} p_s(v) = \frac {v p_e(v)}{\langle u \rangle }. \end{align}
The initial particle positions are distributed along a line of length
$L \gg \ell _Y$
perpendicular to the mean flow, either uniformly or flux weighted. The initial particle distribution is denoted by
$\rho (\boldsymbol{a})$
. For a uniform initial distribution,
$\rho (\boldsymbol{a})=L^{-1} \delta (a_1)$
, where
$a_1$
is the first coordinate of the vector
$\boldsymbol a$
. This initial distribution corresponds to the initial speed distribution
$p_0(v) = p_e(v)$
. The numerical particle tracking methodology is detailed in Appendix B.
2.3. Fluid deformation
We focus on the evolution of the elongation of an infinitesimal material fluid element, whose length and orientation are described by the vector
$ \boldsymbol{z}(t, \boldsymbol{a}) = \boldsymbol{x}(t,\boldsymbol{a} + \delta \boldsymbol{a}) - \boldsymbol{x}(t,\boldsymbol{a})$
. For simplicity of notation, in the following we set
$\boldsymbol{z}(t) = \boldsymbol{z}(t, \boldsymbol{a})$
and drop the argument
$\boldsymbol{a}$
indicating the initial position throughout. The evolution of
$\boldsymbol z(t)$
is governed by
\begin{align} \frac {\text{d} \boldsymbol{z}(t) }{\text{d}t} = \boldsymbol{\epsilon }(t) \boldsymbol{z}(t), \end{align}
where the velocity gradient is defined by
$\boldsymbol{\epsilon }(t) = \boldsymbol{\nabla } \otimes \boldsymbol{u}[\boldsymbol{x}(t,\boldsymbol{a})]^\top \equiv \boldsymbol{\nabla } {\boldsymbol{v}}(t)^\top$
. The transpose is denoted by
$\top$
. The initial condition is given by
$\boldsymbol z(t = 0) = \delta \boldsymbol{a}$
.
Note that the evolution of the infinitesimal line element
$\boldsymbol{z}(t)$
is also related to the deformation gradient tensor
$\boldsymbol{F}(t)\equiv \text{d}\boldsymbol{x}/\text{d}\boldsymbol{a}$
that characterises the transform between the Eulerian and and Lagrangian frames as
$\boldsymbol{z}(t)=\boldsymbol{F}(t)\boldsymbol{z}(0)$
, and analogous to (2.14),
$\boldsymbol{F}(t)$
evolves as
$\dot {\boldsymbol{F}}(t)=\boldsymbol{\epsilon }(t)\boldsymbol{F}(t)$
. The Cauchy–Green tensor
$\boldsymbol{C}(t)=\boldsymbol{F}(t)^\top \boldsymbol{F}(t)$
which characterises fluid deformation plays a central role in controlling e.g. solute mixing and chemical reactions, hence stretching of the line element
$\boldsymbol{z}(t)$
is intimately related to these fluid-borne processes.
Equation (2.14) is transformed into the coordinate system aligned with the streamline along which the fluid particle is moving (Dentz et al. Reference Dentz, Lester, Le Borgne and de Barros2016b ). The transformation is performed through the rotation matrix defined by
\begin{align} \boldsymbol{A} (t) = \frac {1}{v(t)}\begin{bmatrix} v_1(t) & -v_2(t) \\[3pt] v_2(t) & v_1(t) \end{bmatrix}\!. \end{align}
We refer to these coordinates as Protean coordinates in the following (Adachi Reference Adachi1983). The strip
$\boldsymbol z(t)$
is transformed according to
\begin{align} \boldsymbol z'(t) = \boldsymbol A^\top (t) \boldsymbol z(t). \end{align}
Thus, we obtain for the evolution of
$\boldsymbol z'(t)$
\begin{align} \frac {\text{d} \boldsymbol{z}'(t)}{\text{d}t} = \boldsymbol \epsilon '(t) \boldsymbol z'(t), \end{align}
where the velocity gradient transforms as a pseudotensor
\begin{align} \boldsymbol \epsilon '(t) = \boldsymbol{A}^\top (t) \boldsymbol{\epsilon }(t) \boldsymbol{A}(t) + \frac {\text{d} \boldsymbol{A}^\top (t)}{\text{d}t} \boldsymbol{A} (t); \end{align}
see also Appendix A. In the following, we drop the primes for simplicity of notation. The transformed velocity gradient tensor has the upper triangular form
\begin{align} \boldsymbol{\epsilon } (t) = \begin{bmatrix} \epsilon _{11} (t) & \sigma (t) \\[3pt] 0 & -\epsilon _{11} (t) \end{bmatrix}\! . \end{align}
Its components are given by
\begin{align} \epsilon _{11}(t) &= \frac {1}{v(t)}\frac {\text{d} v(t)}{\text{d}t} , \\[-12pt] \nonumber \end{align}
\begin{align} \sigma (t) &= \frac {1}{v(t)^2} \left [v_1(t)^2 \sigma (t)- v_2(t)^2 \sigma (t) -4 v_1(t) v_2(t) \epsilon _{11}(t) \right ]. \\[10pt] \nonumber \end{align}
We defined the shear rate as
$\sigma (t) = \epsilon _{12}(t) + \epsilon _{21}(t)$
. Due to incompressibility,
$\epsilon _{22}(t) = -\epsilon _{11}(t)$
. Using the upper triangular form of
$\boldsymbol \epsilon (t)$
, the kinematic equation (2.17) can be integrated explicitly, which gives
\begin{align} z_1(t) &= \frac {v(t)}{ v(0)} \left [ z_1(0) + z_2(0) \int _0^t {\text{d}}t \sigma (t) \frac { v(0)^2}{ v(t)^2 } \right ] , \\[-12pt] \nonumber \end{align}
\begin{align} z_2(t) &= \frac {v(0)}{ v(t)} z_2(0) . \\[10pt] \nonumber \end{align}
These equations show that the principal stretches given by the diagonal components
$\epsilon _{11}^\prime (t)$
of
$\boldsymbol{\epsilon } (t)$
generate zero net stretching as they directly generate the fluctuating terms
$v(t)/v(0)$
and
$v(0)/v(t)$
in (2.22a
), (2.22b
), respectively. Hence, not only are the Lyapunov exponents associated with these principal stretching directions zero, but the time-averaged principal stretches are likewise. Conversely, the off-diagonal component
$\sigma ^\prime (t)$
of
$\boldsymbol{\epsilon } (t)$
(which corresponds to a fluctuating shear) leads to persistent algebraic stretching of fluid elements, as reflected by the integral in (2.22a
). Note that rotation into the streamline coordinate system (2.15) corresponds to enforcement of the topological constraints (discussed in § 1) associated with the existence of streamlines in steady flow. Here, the alignment of the coordinate frame with the streamlines recovers the zero Lyapunov exponent inherent to steady 2-D flow, whereas naive stochastic models of fluid stretching based on (2.14) require ad hoc approaches to enforce sub-exponential line growth (Le Borgne et al. Reference Le Borgne, Dentz and Villermaux2013, Reference Lester, Dentz, Bandopadhyay and Le Borgne2015). This example reinforces how flow topology controls mixing and transport in porous media.
Equations (2.22a
), (2.22b
) can be formulated in terms of equidistant sampling through the variable change
$t \to s$
. That is, we set
$\text{d}t = ds / v(s)$
in (2.22) to obtain
\begin{align} z_1(s) &= \frac { v(s)}{ v(0)} z_1(0) + v(s) v(0) z_2(0) \zeta (s), && \zeta (s) = \int _0^s {\text{d}}s \frac { \sigma (s)}{ v(s)^3 } , \\[-12pt] \nonumber \end{align}
\begin{align} z_2(s) &= \frac {v(0)}{ v(s)} z_2(0). \\[10pt] \nonumber \end{align}
The strip vector at time
$t$
is then given by
\begin{align} \boldsymbol z(t) = {\boldsymbol z}[s(t)]. \end{align}
The strip length is given by
$\ell (t) = |\boldsymbol z(t)| = \sqrt {z_1(t^2) + z_2(t)^2}$
, its initial length by
$\ell (t = 0) = \ell _0 = |\delta \boldsymbol{a}|$
. The numerical methodology for the solution of the deformation problem is detailed in Appendix B.
2.4. Observables
Figure 3 illustrates the evolution of the dispersion of a line source of fluid particles and the length and orientation of fluid elements in a heterogeneous Darcy velocity field. In the following, we define the observables to probe the dispersion and deformation behaviours. The evaluation of the observables from the detailed numerical transport and deformation simulations is detailed in Appendix B.
(a) Evolution of a set of
$10^5$
particles that are initially distributed uniformly along a line perpendicular to the mean flow direction in a log-normal conductivity field with
$\sigma ^2_{\kern-1pt f} = 1$
and
$\ell _Y = 1$
. Particle positions are shown at
$t = 20 \tau _v$
and
$t= 90 \tau _v$
illustrating longitudinal dispersion. (b) Illustration of line stretching and rotation in the same flow field as in the (a). The black solid lines show the evolution of the vector
$\boldsymbol z(t,\boldsymbol{a})$
for four different initial positions
$\boldsymbol{a}$
, obtained by numerical integration of (2.14) along the trajectories
$\boldsymbol{x}(t,\boldsymbol{a})$
shown by the light grey solid lines, which are obtained by numerical integration of (2.10) as described in Appendix B. The background is a map of the flow magnitude.
2.4.1. Dispersion
We consider the displacement mean and variance of longitudinal particle position average across different realisations of the hydraulic conductivity field
\begin{align} \mu (t) = \left \langle \int \mathrm{d}\boldsymbol{a} \, \rho (\boldsymbol{a}) x_1(t, \boldsymbol{a}) \right \rangle , && \kappa (t) = \left \langle \int \mathrm{d}\boldsymbol{a} \, \rho (\boldsymbol{a}) \left [x_1(t, \boldsymbol{a}) - \mu (t)\right ]^2 \right \rangle \!. \end{align}
The displacement mean indicates the average longitudinal position of the tortuous line shown in figure 3, and the displacement variance its dispersion.
The residence or transit times of fluid particles between the inlet and a control position are measured by the first-passage time defined by
\begin{align} \tau (x_1, \boldsymbol{a}) = \min \left [ t \mid x_1(t, \boldsymbol{a}) \geqslant x_1 \right ]. \end{align}
The distribution of first-passage times
$f(t,x_1)$
or breakthrough curves is obtained by sampling over all particles for a given realisation and averaging over all disorder realisations
\begin{align} f(t,x_1) = \left \langle \int \mathrm{d}\boldsymbol{a}\,\rho (\boldsymbol{a})\,\delta [t-t(x_1\chi )] \right \rangle \!. \end{align}
2.4.2. Lagrangian deformation
We consider the mean elongation of a material strip as a measure for the average fluid deformation
\begin{align} \lambda (t) = \langle \ell (t) \rangle . \end{align}
Furthermore, we consider the PDF of
$\ell (t)$
, which is obtained through averaging over all particles in single realisations and over the ensemble of disorder realisations, as
\begin{align} p_\ell (l) = \left \langle \int \mathrm{d}\boldsymbol{a}\,\rho (\boldsymbol{a})\,\delta [l-\ell (t,\boldsymbol{a})] \right \rangle \!. \end{align}
3. Stochastic dynamics of stretching
Based on the formulation of stretching in the Protean frame outlined in § 2.3, here, we derive a CTRW approach for fluid deformation in heterogeneous porous media that provides both a link to the dispersion process and to the flow and medium structure. We first recall the CTRW formulation for particle dispersion in § 3.1. Then we develop the model for fluid deformation based on the analysis of the statistical properties of shear rate series in § 3.2. The two stochastic models are synthesised in a single CTRW approach for stretching and dispersion in § 3.3.
3.1. Continuous time random walk for particle dispersion
The kinematic equations (2.12) for streamwise particle motion are written as (Comolli et al. Reference Comolli, Hakoun and Dentz2019; Dentz & Hyman Reference Dentz and Hyman2023; Dentz & Massoudieh Reference Dentz and Massoudieh2025)
\begin{align} \frac {\text{d} x_1(s)}{d s} = \chi ^{-1}, \qquad \frac {\text{d} t(s)}{d s} = \frac {1}{ v(s)}, \end{align}
where the streamline distance is projected onto the streamwise distance by tortuosity
$\chi$
. Following Le Borgne et al. (Reference Le Borgne, Dentz and Carrera2008) and Dentz et al. (Reference Dentz, Kang, Comolli, Le Borgne and Lester2016a
), the speed series
$\{ v(s)\}$
is modelled as a stationary Markov process. Note that, here, the particle trajectories depend on the initial condition only through the initial velocity distribution. The explicit dependence on
$\boldsymbol a$
is dropped due to the ergodicity of the stochastic process that models particle motion.
The stationary distribution of
$\{ v(s)\}$
is given by
$p_s(v)$
, which is defined in terms of the Eulerian speed distribution by (2.13). The initial distribution of particle positions
$\rho (\boldsymbol{a})$
translates in this framework into a distribution of initial velocities
$p_0(v)$
, which correspond to the velocities that the particle sampled in the injection zone. The Markov process for particle speeds is defined in terms of the normal score transform of
$ v(s)$
\begin{align} w(s) = \varPhi ^{-1}\{P_s[ v(s)]\}, \qquad v(s) = P_s^{-1}(\phi [ w(s)]), \end{align}
where
$\varPhi (w)$
is the cumulative density function (CDF) of the unit Gaussian PDF
$\phi (w)$
, and
$P_s(v)$
is the CDF of the stationary speed PDF
$p_s(v)$
. The series
$\{ w(s)\}$
of normal scores is by definition Gaussian distributed. The evolution of
$ w(s)$
is modelled by the Ornstein–Uhlenbeck process (Morales et al. Reference Morales, Dentz, Willmann and Holzner2017; Hakoun et al. Reference Hakoun, Comolli and Dentz2019)
\begin{align} \frac {\text{d} w(s)}{\text{d}s} = - w(s) \ell _c^{-1} + \sqrt {2 \ell _c^{-1}} \eta _1(s), \end{align}
which implies that it is also stationary (Doob Reference Doob1942). The random function
$ \eta _1(s)$
denotes a unit Gaussian white noise. In the following, we use the notation
$\eta _i(s)$
(
$i = 1, 2, \ldots$
) to refer to independent unit Gaussian random noises. The two
$\eta _i(s)$
and
$\eta _k(s)$
are independent for
$i \neq k$
.
This formulation is equivalent to generating the speed series by a Gaussian copula (Dentz & Hyman Reference Dentz and Hyman2023; Dentz & Massoudieh Reference Dentz and Massoudieh2025). The CTRW (3.1) based on the speed process described above represents a stochastic model for particle motion in steady random flows and has been used for the modelling of dispersion from the pore to the reservoir scales in porous and fractured media (Puyguiraud et al. Reference Puyguiraud, Gouze and Dentz2019; Comolli et al. Reference Comolli, Hakoun and Dentz2019; Hyman & Dentz Reference Hyman and Dentz2021).
In this framework, the displacement mean and variance in streamwise direction are
\begin{align} \mu (t) = \mathcal{\chi }^{-1} \langle s(t) \rangle , \qquad \kappa (t) = \mathcal{\chi }^{-2} \langle [s(t) - \langle s(t) \rangle ]^2 \rangle , \end{align}
where
$s(t) = \max [s|t(s) \leqslant t]$
. The angular brackets here denote the ensemble average over all particle realisations. The first-passage time distribution is
\begin{align} f(t, x_1) = \langle \delta [t - t(x_1 \mathcal{\chi })] \rangle , \end{align}
where
$t(s)$
is given by (3.1).
3.2. Stochastic dynamics of deformation
The deformation process, given by (2.23) depends on the speed series
$\{ v(s)\}$
and the series of shear rates
$\{ \sigma (s) \}$
. In fact, the auxiliary process
$ \zeta (s)$
defined in (2.23a
) evolves as
\begin{align} \frac {\text{d} \zeta (s)}{\text{d}s} = \frac { \sigma (s)}{v(s)^3}. \end{align}
The upscaled model for Lagrangian fluid deformation requires definition of the stochastic process for the series of shear rates
$\{ \sigma (s) \}$
, which is correlated with the particle speeds
$ v(s)$
. The average of
$\sigma (s)$
is equal to zero due to the stationarity of the velocity field
$\boldsymbol{u}(\boldsymbol{x})$
. In order to identify the relation to the speed process, we consider the mean absolute shear rate conditional to the particle speed,
$\langle |\sigma || v\rangle$
. Figure 4 shows the conditional mean sampled along particle trajectories in log-normal and gamma-distributed conductivity fields. We consistently observe a linear relation between the conditional mean absolute shear rate and velocity,
$\langle |\sigma ||v\rangle = k_c v$
, with
$k_c$
a proportionality constant. The inverse of
$k_c$
denotes a characteristic shear length scale. Thus, we write
$ \sigma (s)$
as the following multiplicative process:
\begin{align} \sigma (s) = k_c v(s) \xi (s), \end{align}
where
$ \xi (s)$
is a stochastic process with zero mean and unit average
$\langle | \xi (s)| \rangle = 1$
. Thus, the shear rate is defined entirely by the speed process
$v(s)$
and the process
$\xi (s)$
defined below. The numerical values for
$k_c$
for the log-normal and gamma-distributed conductivity fields under consideration in this work are given in table 1. Note that for the gamma distribution,
$k_c$
appears to be constant.
Values of
$k_c$
, see (3.7), for the (bottom) gamma-distributed and (top) log-normal conductivity fields.
Table 1. Long description
A table comparing values of k c for gamma-distributed and log-normal conductivity fields. The table has three columns and four rows, including headers. Column headers are sigma squared f values (0.25, 1, 3), k c values, and alpha values. Row labels are k c and alpha. Row 1: sigma squared f = 0.25, k c = 0.029, alpha = 0.75. Row 2: sigma squared f = 1, k c = 0.065, alpha = 0.5. Row 3: sigma squared f = 3, k c = 0.097, alpha = 0.25. Row 4: k c = 0.074 for all sigma squared f values.
Conditional mean absolute shear rate for (a) log-normal distribution of
$K(\boldsymbol{x})$
and (b) gamma distribution of
$K(\boldsymbol{x})$
. The dashed lines denote the linear scaling with
$v$
.
Figure 4. Long description
Two scatter plots are presented side by side. Panel A on the left shows the conditional mean absolute shear rate for a log-normal distribution of permeability. The x-axis is labeled with the variable v and spans from 10^-4 to 10^2. The y-axis is labeled with the variable <|σ’||u’> and spans from 10^-5 to 10^2. The plot includes data points in different colors representing various values of σ²f: 0.25 (red circles), 1 (blue circles), 3 (purple circles), and 6 (green circles). A dashed line denotes the linear scaling with v. Panel B on the right shows the conditional mean absolute shear rate for a gamma distribution of permeability. The x-axis is labeled with the variable v and spans from 10^-6 to 10^1. The y-axis is labeled with the variable <|σ’||u’> and spans from 10^-6 to 10^1. The plot includes data points in different colors representing various values of α: 0.75 (red squares), 0.5 (blue squares), and 0.25 (green squares). A dashed line denotes the linear scaling with v.
(a) (symbols) The PDF of
$ \xi (s)$
sampled from the numerical simulations and (solid line) the corresponding Gaussian PDF. (b) Evolution of
$ \xi (s)$
in terms of streamline distance and covariance function of
$ \xi (s)$
, which decays on the correlation scale
$\ell _c$
as shown by the red dotted line representing a function proportional to
$\exp (-s/\ell _c)$
.
The auxiliary process (3.6) for
$ \zeta (s)$
thus can be written as
\begin{align} \frac {\text{d} \zeta (s)}{\text{d}s} = \frac {k_c \xi (s)}{ v(s)^2}. \end{align}
Figure 5 shows the PDF of
$ \xi (s) = \sigma (s)/ v(s) k_c$
and illustrates a single
$\xi (s)$
-trajectory. Its PDF can be well approximated by a Gaussian distribution. Furthermore, its correlation function is approximately exponential and decorrelates on the same scale
$\ell _c$
as the speed process
$ v(s)$
. Thus, we model
$ \xi (s)$
as the Ornstein–Uhlenbeck process
\begin{align} \frac {\text{d} \xi (s)}{\text{d}s} = - \xi (s) \ell _c^{-1} + \sqrt {\pi \ell _c^{-1}} \eta _2(s). \end{align}
The initial values
$ \xi (s = 0) = \xi _0$
are distributed according to
\begin{align} p_\xi (y) = \exp \big (- y^2/\pi \big )/\pi , \end{align}
which is equal to the stationary distribution of
$ \xi (s)$
. It is Gaussian with the mean absolute value equal to one. Fluid deformation in this framework is obtained from (2.23) by using the stochastic processes for
$ v(s)$
and
$ \sigma (s)$
outlined above and projection from
$s \to t$
according to (3.1)
\begin{align} \ell (t) = |{\boldsymbol{z}}[s(t)]|. \end{align}
The PDF of strip elongations
$\ell (t)$
is obtained as
\begin{align} p_\ell (l,t) = \left \langle \delta [l - \ell (t)] \right \rangle . \end{align}
3.3. Summary
For the convenience of the reader, we summarise the CTRW model for dispersion and stretching of fluid particles. It can be written in differential form as
\begin{align} x(s+\text{d}s) &= x(s) + \frac {\text{d}s}{\chi }, \quad t(s + ds) = t(s) + \frac {\text{d}s}{ v(s)}, \quad v(s) = P_s^{-1}(\phi [ w(s)]), \\[-28pt] \nonumber \end{align}
\begin{align} w(s + \text{d}s) = w(s) - \ell _c^{-1} w(s) + \sqrt {2 \ell _c^{-1} \text{d}s} Z_1(s), \\[0pt] \nonumber \end{align}
respectively, for particle position, speed and time and
\begin{align} z_1(s) = \frac { v(s)}{ v(0)} z_1(0) + v(s) v(0) z_2(0) \zeta (s), \quad z_2(s) = \frac {v(0)}{ v(s)} z_2(0), \\[-28pt] \nonumber \end{align}
\begin{align} \zeta (s+ ds) = \zeta (s) + \frac {k_c \xi (s)}{ v(s)^2}, \quad \xi (s+ds) = \xi (s) + \ell _c^{-1} + \sqrt {\pi \ell _c^{-1} ds} Z_2(s), \\[0pt] \nonumber \end{align}
for the deformations. The
$Z_1(s)$
and
$Z_2(s)$
are independent unit Gaussian random variables. The distribution of initial values
$ \xi _0$
is given by (3.10). The initial velocity distribution
$p_0(v)$
and the distribution of the initial strip orientations are determined by the system under consideration.
Velocity distributions for (a) log-normal fields (b) gamma-distributed conductivity fields. The dashed lines in (b) denotes the scalings
$v^{\beta - 1}$
with (top to bottom)
$\beta = 0.3, 0.55, 0.85$
. Note that the power-law behaviours for the gamma fields are cutoff at a minimum velocity
$v_m$
.
Figure 6. Long description
Panel A: A scatter plot shows the probability density function (PDF) of velocity (v) for log-normal fields. The x-axis represents velocity (v) on a logarithmic scale ranging from 10^-6 to 10^2, and the y-axis represents the PDF on a logarithmic scale ranging from 10^-8 to 10^0. Different colored markers represent different variances (σ^2_f) of 0.25, 1, 3, and 6. Panel B: A scatter plot shows the PDF of velocity (v) for gamma-distributed conductivity fields. The x-axis represents velocity (v) on a logarithmic scale ranging from 10^-10 to 10^0, and the y-axis represents the PDF on a logarithmic scale ranging from 10^-5 to 10^5. Different colored markers represent different alpha (α) values of 0.25, 0.5, and 0.75. Dashed lines indicate scalings with power-law behaviors cutoff at a minimum velocity.
3.4. Discussion of model parameters
The CTRW framework for dispersion and deformation outlined above uses the Eulerian speed distribution
$p_e(v)$
, advective tortuosity
$\chi$
, correlation length
$\ell _c$
and proportionality constant
$k_c$
as inputs. The Eulerian speed distribution and advective tortuosity are flow attributes and do not depend on the Lagrangian particle dynamics. For moderately heterogeneous log-normal conductivity fields, it has been found that the speed distribution is also log-normal (Cvetkovic et al. Reference Cvetkovic, Cheng and Wen1996). For strongly heterogeneous media, the distribution of low flow speeds can be related to the statistics of low hydraulic conductivities (Fiori et al. Reference Fiori, Janković, Dagan and Cvetković2007; Hakoun et al. Reference Hakoun, Comolli and Dentz2019; Hidalgo et al. Reference Hidalgo, Neuweiler and Dentz2022). Figure 6 shows the speed PDFs resulting from the log-normal and gamma-distributed conductivity fields confirming these observations. For the gamma fields, the speed PDFs behave as
$p_e(v) \propto v^{\beta -1}$
at small speeds down to a cutoff at around
$10^{-7} \langle u \rangle$
. The values of the exponent are
$\beta = 0.3, 0.55, 0.85$
for
$\alpha = 0.25, 0.5, 0.75$
. They do not coincide with
$\alpha$
, but are consistent with the values observed in Hakoun et al. (Reference Hakoun, Comolli and Dentz2019).
The length
$\ell _c$
is the correlation distance of the streamwise Lagrangian speed
$ v(s)$
. It is related to the correlation length
$\ell _Y$
of the multi-Gaussian field
$Y(\boldsymbol{x})$
underlying the random conductivity fields
$K(\boldsymbol{x})$
. Analytical expressions for
$\ell _c$
are only known from perturbation theory for weakly heterogeneous fields (Cvetkovic et al. Reference Cvetkovic, Cheng and Wen1996), which gives
$\ell _c=8/3 \ell _Y$
. For strong heterogeneity, we currently rely on empirical expressions (Cvetkovic et al. Reference Cvetkovic, Cheng and Wen1996; Hakoun et al. Reference Hakoun, Comolli and Dentz2019).
As outlined above, the distribution of initial particle speeds and strip orientations depend on the system under consideration. In order to compare the CTRW model with the detailed numerical simulations discussed in Appendix B, the initial speed distribution is set equal to the Eulerian speed PDF,
$p_0(v) = p_e(v)$
, which corresponds to the uniform initial particle distribution. The distribution for the initial orientation is set equal to the one used for the detailed simulations in the Protean frame given in (B11).
4. How porous medium heterogeneity controls stirring
We study here how the porous medium heterogeneity controls stirring. To this end, we consider detailed numerical simulations of the stretching and dispersion behaviours in the heterogeneous hydraulic conductivity fields given in § 2.1. The medium heterogeneity is characterised by the
$K(\boldsymbol{x})$
statistic, which is mapped through the flow equation (2.2) onto the flow velocity
$\boldsymbol q(\boldsymbol{x})$
. We analyse here how the observed dispersion and stretching behaviours can be related to the medium heterogeneity in terms of the distribution of Eulerian flow velocities using the CTRW model presented in § 3. Note that the dispersion behaviours for this type of conductivity field was discussed in depth in Comolli et al. (Reference Comolli, Hakoun and Dentz2019). Thus, here we focus on stretching. We present the dispersion behaviour for completeness and to demonstrate the capability of the CTRW model of § 3 to quantify and predict both dispersion and stretching.
Top panel: (left) First-passage time distribution and (right) displacement variance for the gamma-distributed conductivity fields. The straight dashed lines in the left panel indicate the scalings
$t^{-1-\beta }$
, and in the right the scalings
$t^{2-\beta }$
. The power-law scalings are visible in the time regime between the peak time and a cutoff time scale that is given by the lower cutoff velocity in the velocity distribution shown in figure 6, which translates into a cutoff of the particle transition times at around
$10^6 \tau _v$
. Bottom panel: The same for the log-normal conductivity fields. The symbols denote the data from the detailed numerical simulations, the straight lines the CTRW model. The black dashed line in the left panel denotes the fit of the data for
$\sigma f^2 = 6$
with an inverse Gaussian distribution, the analytical solution of the advection–dispersion equation with constant velocity and longitudinal dispersion coefficient. The black dashed line in the right panel denotes the asymptotic linear scaling of the displacement variances at large times.
4.1. Dispersion
Figure 7 shows the first-passage time distributions and displacement variances for particle motion in gamma-distributed and log-normal conductivity fields and their estimates based on the CTRW model. For both scenarios, the CTRW model provides an excellent description of the direct numerical data. It reproduces all non-Fickian transport features and shows the link between the particle dispersion and the medium heterogeneity in terms of the Eulerian velocity field.
In order to further analyse the observed scalings and non-Fickian dispersion features, we first note that particle velocities are independent for separation distances larger than
$\ell _c$
and vary only weakly for distances larger than
$\ell _c$
. Thus, the equations of motion (3.13a
) of a particle can be coarse grained as
\begin{align} x_{n+1} = x_n + \frac {\ell _c}{\chi }, \qquad t_{n+1} = t_n + \tau _n, \qquad \tau _n = \frac {\ell _c}{v_n}, \end{align}
where the
$v_n$
are independent random velocities. For the uniform initial particle distribution under consideration here, the initial particle speeds
$v_0$
are distributed according to
$p_e(v)$
. For all subsequent steps, they are distributed according to
$p_s(v)$
, which is defined in terms of
$p_e(v$
by (2.13). The PDF of transition times
$\tau _n$
for
$n \gt 0$
is denoted by
$\psi (t)$
, which can be expressed in terms of the Eulerian speed distribution by
\begin{align} \psi (t) = \frac {\tau _c \ell _c}{t^3} p_e(\ell _c/t), \end{align}
where we defined
$\tau _c = \ell _c/\langle u \rangle$
. The initial transition times
$\tau _0 = \ell _c/v_0$
are given accordingly by
\begin{align} \psi _0(t) = \frac {\ell _c}{t^2} p_e(\ell _c/t). \end{align}
In this coarse-grained framework, the first-passage time of a particle at a distance
$x_n = n \ell _c/\chi$
is given by
\begin{align} t_n = \tau _0 + \sum \limits _{i = 1}^{n-1} \tau _i. \end{align}
That is, the first-passage times are the sums of independent random transition times. Their distribution is given by the convolution of
$\psi _0(t)$
and the
$n-1$
-fold convolution of
$\psi (t)$
. Thus, the scaling of the first-passage time distribution can be understood in terms of the scaling of the transition time distribution, and thus, in terms of the velocity distributions.
For the gamma fields, the Eulerian velocity distribution scales as
$p_e(v) \propto v^{\beta -1}$
with
$0 \lt \beta \lt 1$
for velocities in the interval
$v_m \ll v \ll \langle u \rangle$
, that is, between a lower cutoff velocity
$v_m$
and the mean velocity
$\langle u \rangle$
; see figure 6. Thus, the initial transition time distribution scales as
$\psi _0(t) \sim t^{-1-\beta }$
for
$\tau _c \ll t \ll \tau _m$
, while the transition time distributions for all subsequent steps scale as
$\psi (t) \sim t^{-2-\beta }$
. That is,
$\psi _0(t)$
has a stronger tail behaviour than
$\psi (t)$
, and thus it dominates the tail behaviour of the first-passage time distribution, which scales as
$f(t,x_1) \propto t^{-1-\beta }$
for
$t \ll \tau _m$
. The power-law tails are cutoff at times of around
$\tau _m = \ell _c/v_m \approx 10^6 \tau _v$
. The displacement means for the gamma fields (not shown) are evolving linearly with time and are equal to
$\mu (t) = \langle u \rangle t$
. The displacement variance evolves ballistically at short times
$t \lt \tau _v$
, and according to a power law at increasing time. The CTRW theory predicts the behaviour
$\kappa (t) \propto t^{2-\beta }$
for
$\tau _c \ll t \ll \tau _m$
(Comolli et al. Reference Comolli, Hakoun and Dentz2019); see figure 7(a).
For the log-normal conductivity fields, the first-passage time distributions do not display heavy tails for small
$\sigma _{\kern-1pt f}^2 \leqslant 1$
; see figure 7(b). For larger
$\sigma _{\kern-1pt f}^2$
, we observe strong tailing, which again can be traced back to the broad distribution of Eulerian flow speeds, which displays power-law tangents over more than one order of magnitude for intermediate values of
$v$
. The displacement mean evolves as
$\mu (t) = \langle u \rangle t$
. The displacement variance shows ballistic behaviour at early times and then crosses over toward a linear long-time scaling. For
$\sigma _{\kern-1pt f}^2 = 3$
and
$6$
, there is a broad power-law cross-over, which originates in the broad distribution of flow speeds shown in figure 6. While the displacement variance asymptotically tends toward Fickian, that is, linear scaling for all
$\sigma _{\kern-1pt f}^2$
, the first-passage time distribution may still display heavy tailed behaviour. This is illustrated in figure 7(b) for the case
$\sigma _{\kern-1pt f}^2 = 6$
. The dashed black line denotes the fit of an inverse Gaussian distribution ((11) in Kreft & Zuber (Reference Kreft and Zuber1978)) with constant velocity and longitudinal dispersion coefficient. While the bulk behaviour is well represented by the analytical solution, the tail behaviour cannot be captured.
4.2. Stretching
In the previous section, we have seen, how the medium heterogeneity impacts on particle dispersion and how this can be quantified in terms of the CTRW model. In this section, we analyse the stretching behaviours in the gamma and log-normal conductivity fields in terms of the mean elongation of a material strip and the PDF of strip elongations. In order to link the stretching behaviour to the medium heterogeneity, we use the CTRW model presented in § 3.
Mean strip elongation for (left) gamma and (right) log-normal conductivity fields. Full symbols denote the detailed numerical simulations, open symbols the CTRW model. The solid lines in the left panel denote the power-law scaling obtained from the Lévy walk model (4.8) for the stretching evolution. The asymptotic scaling is reached at around
$10^4\tau _v$
. The elongation data are multiplied by multiples of 5 in the left panel and by multiples of 5/2 in the right panel for clearer illustration. The colour code is the same as in figure 7.
Figure 8. Long description
Panel A: A line graph depicts mean strip elongation for gamma conductivity fields. The x-axis represents the normalized time (t/tau_v) on a logarithmic scale ranging from 10^-1 to 10^5. The y-axis represents the mean elongation on a logarithmic scale ranging from 10^1 to 10^7. The graph includes full symbols for detailed numerical simulations and open symbols for the CTRW model. Solid lines denote power-law scaling from the Lévy walk model. The elongation data are multiplied by multiples of 5 for clearer illustration. Panel B: A line graph depicts mean strip elongation for log-normal conductivity fields. The x-axis represents the normalized time (t/tau_v) on a logarithmic scale ranging from 10^-1 to 10^5. The y-axis represents the mean elongation on a logarithmic scale ranging from 10^1 to 10^6. The graph includes full symbols for detailed numerical simulations and open symbols for the CTRW model. The elongation data are multiplied by multiples of 5/2 for clearer illustration. The color code is consistent with figure 7.
Figure 8 shows the evolution of the mean elongation
$\lambda (t)$
for the gamma and log-normal conductivity fields obtained using Protean coordinates, as outlined in Appendix B.3. The data from the detailed numerical model are compared with the estimates from the CTRW model presented in § 3. The CTRW model is able to describe all the salient features of the evolution of the mean elongation in terms of only the Eulerian speed distribution. It provides a direct link between medium and flow heterogeneity and the stretching behaviour. In the following, the observed stretching behaviours are discussed in more detail in the light of this relation.
At early times, the evolution of
$\ell (t)$
is linear, which can be understood by considering equations (2.22), and noting that, at times
$t \ll \tau _c$
, particle speeds and shear rates can be considered constant and equal to their initial values
$v_0$
and
$\sigma _0$
. Inserting these short-time approximations into (2.22) gives
\begin{align} z_1(t) = \delta a_1 + \delta a_2 \sigma _0 t, \qquad z_2(t) = \delta a_2. \end{align}
Thus, the short time evolution of the strip elongation
$\ell (t) = \sqrt {z_1(t)^2 + z_2(t)^2}$
is given by
\begin{align} \ell (t) = \sqrt {(\delta a_1 + \delta a_2 \sigma _0 t)^2 + \delta a_2^2} = \ell _0 + \frac {\delta a_1 \delta a_2}{\ell _0} \sigma _0 t + \cdots , \end{align}
where the dots denote higher-order contributions. Hence, the early time behaviour is linear.
For increasing times, the mean elongation crosses over to a power-law behaviour as
\begin{align} \lambda (t) \propto t^{\gamma }, \end{align}
where
$\gamma$
is a scaling exponent that is specified in the following. In order to assess the long-time evolution of fluid deformation, Dentz et al. (Reference Dentz, Lester, Le Borgne and de Barros2016b
) formulated the stretching process as a coupled CTRW or Lévy walk. In fact, and as outlined in § 4.1, for separation distances larger than the correlation length
$\ell _c$
, subsequent particle speeds can be considered independent, while the particle velocities vary only weakly for distances smaller than
$\ell _c$
. In analogy to the coarse graining of the particle motion, we coarse grain the deformation process on the scale
$\ell _c$
and write (3.8) as
\begin{align} \zeta _{n+1} = \zeta _n + k_c \ell _c^{-1} \sqrt {\pi } \xi _n \tau _n^2, \qquad t_{n+1} = t_n + \tau _n. \end{align}
The
$\xi _n$
values are independent unit Gaussian random variables. The process (4.8) denotes a coupled CTRW, and for a power-law distribution of transition times, a Lévy walk (Zaburdaev et al. Reference Zaburdaev, Denisov and Klafter2015; Dentz et al. Reference Dentz, Le Borgne, Lester and de Barros2015). The coupling is expressed by the fact that the increment of the
$\zeta _n$
process is proportional to the square of the time increment
$\tau _n$
. The strip elongation after
$n$
random walk steps is then given by
\begin{align} \ell _n = \sqrt {\left [\frac {v_n}{v_0} z_1(0) + v_n z_2(0) \zeta _n\right ]^2 + \left [\frac {v_0}{v_n} z_2(0) \right ]^2}\!. \end{align}
In the limit of large
$\zeta _n$
one can expand this expression, and one obtains for the strip elongation at leading order
\begin{align} \ell _n = v_n v_0 \delta a_2 | \zeta _n|, \end{align}
where we used the fact that
$z_2(0) = \delta a_2$
. Thus, the evolution of the strip elongation can be rationalised in terms of the Lévy walk (4.8). This coarse-grained picture is valid for travel distances much larger than
$\ell _c$
, or equivalently for times much larger than
$\tau _c$
.
The power-law behaviour given by (4.7) can now be understood on the basis of the Lévy walk (4.8) and the scalings derived in Rebenshtok et al. (Reference Rebenshtok, Denisov, Hänggi and Barkai2014) and Dentz et al. (Reference Dentz, Le Borgne, Lester and de Barros2015). Recall that the distribution of transition times
$\psi (t)$
exhibits the power-law behaviour
$\psi (t) \propto t^{-2-\beta }$
for
$\tau _c \ll t \ll \tau _m$
, as outlined in § 4.1. Dentz et al. (Reference Dentz, Le Borgne, Lester and de Barros2015) show that, for
$0 \lt \beta \lt 1$
, which is the case for the gamma fields, Lévy walk theory predicts that the scaling exponent in (4.7) is equal to
$\gamma = 2 - \beta$
. The data for
$\lambda (t)$
from the detailed numerical simulations and the CTRW model agree with these theoretical scalings. For the log-normal conductivity fields, the Eulerian speed PDFs do not exhibit a clear scaling at small speeds. The scaling behaviours for the elongation in this case fall into the class defined by
$\beta \gt 3$
. In this case, Lévy walk theory predicts the scaling exponent
$\gamma = 1/2$
(Dentz et al. Reference Dentz, Le Borgne, Lester and de Barros2015). At intermediate times, power-law cross-overs can be observed that are related to power-law tangents in the speed PDF.
The PDFs of strip elongations for (a) gamma fields for
$\alpha = 0.25$
and (b) log-normal conductivity fields for
$\sigma _{\kern-1pt f}^2 = 6$
at times
$t = 1, 5, 20\tau _c$
. The symbols denote the detailed numerical simulations, the solid lines the CTRW model and the dash-dotted lines fits of the log-normal distribution for
$t = 20 \tau _c$
. The PDFs are shifted along the vertical axis by multiplication with a constant factor to avoid overlap.
The distribution of strip elongations
$p_\ell (l,t)$
plays a central role for the quantification of the concentration content of a disperse mixture (Le Borgne et al. Reference Le Borgne, Dentz and Villermaux2013, Reference Lester, Dentz, Bandopadhyay and Le Borgne2015). Regions of strong stretching lead to a fast decay of solute concentrations by diffusion, because stretching leads to a steepening of concentration gradients. This dynamics is quantified by the diffusive strip method (Villermaux Reference Villermaux2019), which allows us to map the elongation PDF directly on the PDF of concentration point values.
Figure 9 shows the PDF
$p_\ell (l,t)$
for the gamma and log-normal fields for the strongest heterogeneities with
$\alpha = 0.25$
and
$\sigma _{\kern-1pt f}^2 = 6$
at different times. The Supplementary Information provides the results for the elongation PDFs for all the analysed hydraulic conductivity fields. The CTRW model is in excellent agreement with the data from the detailed numerical simulations at all observation times. At short times,
$p_\ell (l,t)$
the PDF is sharply peaked at the initial elongation. At increasing times, it becomes skewed toward larger values of elongation following the coupled CTRW process, and the peak is smoothed
These results show how the medium heterogeneity impacts on and controls the elongation and thus the stretching behaviour. Both the mean elongation and its distribution can be linked to the medium heterogeneity via the CTRW model.
4.3. Discussion
The way porous media stir a fluid is very different from the way chaotic or turbulent flows stir a free fluid. We are considering Darcy-scale flow, that is, the bulk fluid is stirred by the spatial variability as it is permeating through the medium driven by a pressure gradient. The ensuing spatial variability of the steady flow field leads to the deformation of a material line as shown in figure 3, and the dispersion of solute fronts. Furthermore, material strips experience a series of deformation events, stretching and compression, as they are transported along the streamlines of the bulk flow, also shown in figure 3. Steady isotropic Darcy is special in the sense that chaotic flow is prohibited (Lester et al. Reference Lester, Dentz, Bandopadhyay and Le Borgne2022), and thus stretching is sub-exponential, and in general characterised by power-law stretching behaviours, as shown in the previous section. The works by Le Borgne et al. (Reference Le Borgne, Dentz and Villermaux2013) and Le Borgne et al. (Reference Le Borgne, Dentz and Villermaux2015) proposed to model the evolution of strip elongation as a multiplicative stochastic process with coefficients whose time dependence is adjusted from the elongation data so as to match the observed mean elongation behaviour. This framework predicts the PDF of elongation to be log-normal. As shown in figure 9, the PDF is not log-normal at times smaller than
$20 \tau _c$
. For
$t = 20 \tau _c$
a log-normal distribution provides a reasonable approximation of the data. However, this modelling approach is based on a stochastic dynamics that is obtained by fitting the stretching time series as a geometric Brownian motion, and therefore it does not provide a link between the stretching behaviour and medium heterogeneity
In this paper, we have taken a different route with the aim of understanding and quantifying how the medium heterogeneity stirs the bulk fluid and how it controls the stretching behaviour. This implies the formulation of the stretching problem in Protean coordinates aligned with the streamlines of the flow in § 2.3. This framework reveals the dynamics of stretching that results from the topological constraints of steady Darcy flow. It shows that stretching is dominated by the series of shear deformations of variable durations that a strip experiences as it is transported through the heterogeneous medium. The duration of shear episodes varies because the residence or transition times of a material strip vary according to the local flow velocity. This formulation allows for a first-principles derivation of the stochastic stretching model in § 3. The model belongs to the class of CTRWs, which naturally reflects the fact that flow and, in this case, stretching, are organised along the spatial structure imprinted in the hydraulic conductivity distribution. In § 4, this model is used to quantitatively link the observed dispersion and stretching behaviours to the Eulerian velocity distribution.
5. Summary and conclusions
We study the mechanisms and stochastic dynamics of stretching and dispersion of fluid elements in heterogeneous two-dimensional porous media. We consider the quantitative link between stretching and dispersion on the one hand, and between the stretching behaviour and the medium and flow properties on the other. This leads to an upscaled modelling approach and identification of the key parameters to predict both large-scale fluid deformation and dispersion, which are key inputs for mixing models in porous media (Villermaux Reference Villermaux2019).
In order to identify the stochastic dynamics of fluid deformation, we first formulate the deformation problem in Protean coordinates that are aligned along streamlines. This formulation reveals the key role played by the Lagrangian series of particle speeds and shear rates for the evolution of the elongation of a fluid material strip. We focus on particle speeds and shear rates sampled equidistantly along streamlines, which reflect the spatial organisation of the steady flow fields. Their stochastic evolutions are quantified through stationary Markov processes. Both are characterised by a single correlation length scale. Following the CTRW approach for particle dispersion, the speed process is represented by an Ornstein–Uhlenbeck process for the normal score transform of particle speeds. For the shear rate process, we find that (i) it is linearly correlated to the particle speed, (ii) that its conditional distribution is Gaussian and (iii) that it evolves on the same correlation scale as the speed process. Based on these findings, we model the stochastic evolution of shear rates also as an Ornstein–Uhlenbeck process. This CTRW formulation renders a single stochastic framework for the quantification and analysis of dispersion and fluid deformation.
The framework is employed to analyse deformation of fluid particles obtained from direct numerical simulations in gamma-distributed and log-normal hydraulic conductivity fields. Our analysis confirms existing results for dispersion in terms of first-passage time distributions and displacement moments. We observe common signatures of non-Fickian dispersion, such as a broad distribution of first-passage times and nonlinear evolutions of the displacement variance. These non-Fickian features are also reflected by the mean elongation of material strips, which is a measure of stretching. The initially linear behaviour evolves toward asymptotic power-law scalings that are determined by the behaviour of the Eulerian speed PDF. These behaviours are fully predicted by the derived CTRW stretching model in terms of the Eulerian speed PDF. The observed scalings are related analytically to the scalings of the speed PDF by noting that the coarse-grained deformation process follows a coupled CTRW, or, in the case of power-law distributions of particle transition times, a Lévy walk. In summary, the derived CTRW model predicts the stretching and dispersion behaviour in terms of the Eulerian speed PDF and thus elucidates how the medium heterogeneity controls stirring.
The understanding of the heterogeneity controls on stretching is central for the quantification and prediction of mixing and chemical reactions in porous media (Villermaux Reference Villermaux2019; Rolle & Le Borgne Reference Rolle and Le Borgne2019). Stretching of a material strip leads to strong concentration gradients in the direction perpendicular to it, which accelerate mixing and chemical reactions. In a heterogeneous porous medium, regions of strong stretching are related to reaction hotspots (Rolle & Le Borgne Reference Rolle and Le Borgne2019), which can explain dissolution and precipitation patterns in porous media (Hidalgo et al. Reference Hidalgo, Dentz, Cabeza and Carrera2015; Cil et al. Reference Cil, Xie, Packman and Buscarnera2017; Kottsova et al. Reference Kottsova, Kong, Zitha, Saar, Bruhn, Knornschild, Allaz, Archer and Brehme2025) and give insight into karst formation (De Vriendt et al. Reference De Vriendt, Pool and Dentz2020) and the spatial distribution of ore deposits, which are typically located along structures that promote mixing.
Supplementary material
Supplementary material are available at https://doi.org/10.1017/jfm.2026.12192.
Acknowledgements
K.F., J.J.H. and M.D. acknowledge funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement Number 956457 through the project CoPerMix.
Declaration of Interests
The authors report no conflict of interest.
Appendix A. Transformation of the velocity gradient tensor into Protean coordinates
In order to derive (2.18), we take the time derivative of (2.16), which is given by
\begin{align} \frac {\text{d} \boldsymbol z'}{\text{d}t} = \frac {\text{d}}{\text{d}t}\big [\boldsymbol A^{T}(t) \boldsymbol z(t)\big ] = \frac {\text{d} \boldsymbol A^{T}}{\text{d}t} \boldsymbol z(t) + \boldsymbol A^{T}(t)\frac {\text{d}\boldsymbol z(t)}{\text{d}t}. \end{align}
Using (2.14) gives
\begin{align} \frac {\text{d} \boldsymbol z'}{\text{d}t} = \frac {\text{d} \boldsymbol A^{T}}{\text{d}t} \boldsymbol z(t) + \boldsymbol A^{T}(t) \boldsymbol \epsilon (t) \boldsymbol z(t) = \left [ \frac {\text{d} \boldsymbol A^{T}}{\text{d}t} + \boldsymbol A^{T}(t) \boldsymbol \epsilon (t) \right ] \boldsymbol z(t). \end{align}
Using that
$\boldsymbol z(t) = \boldsymbol A(t) \boldsymbol z'(t)$
gives
\begin{align} \frac {\text{d} \boldsymbol z'}{\text{d}t} = \left [ \frac {\text{d} \boldsymbol A^{T}}{\text{d}t} + \boldsymbol A^{T}(t) \boldsymbol \epsilon (t) \right ] \boldsymbol A(t) \boldsymbol z'(t) = \left [ \frac {\text{d} \boldsymbol A^{T}}{\text{d}t} \boldsymbol A(t) + \boldsymbol A^{T}(t) \boldsymbol \epsilon (t) \boldsymbol A(t) \right ] \boldsymbol z'(t). \end{align}
The expression in the square brackets on the right side defines the transformed velocity gradient in (2.18).
Appendix B. Numerical simulations
B.1 Field generation and flow simulations
As outlined in § 2.1, the hydraulic conductivity fields are generated by transformation from a Gaussian-distributed multi-Gaussian random field
$Y(\boldsymbol{x})$
. The latter is generated using the fast Fourier transform following Ruan & McLaughlin (Reference Ruan and McLaughlin1998) on a two-dimensional computational domain that is discretised into
$N = 1000 \times 2000$
nodes in directions transverse and longitudinal to the mean flow. The method assures periodicity in both directions. The domain dimensions are
$100 \ell _Y \times 200 \ell _Y$
. We generate
$N_R = 200$
realisations of each random field.
The conductivity fields with a log-normal marginal distribution are generated from
$Y(\boldsymbol{x})$
according to (2.7). In the numerical simulations, we set
$\mu _{\kern-1pt f} = 0$
and probe different values for
$\sigma _{\kern-1pt f}^2$
. For the gamma fields, there is no analytical map, because there is no closed form expression for the inverse of the cumulative gamma distribution. Here, we use a gamma distribution that is truncated from below at
$k_{\textit{min}}$
. Thus, its cumulative distribution is given by
\begin{align} P_K(k) = \frac {P_{\varGamma }(k) - P_{\varGamma }(k_{\textit{min}})} {1 - P_{\varGamma }(k_{\textit{min}})}, \end{align}
where
$P_\varGamma (k)$
is the CDF of the gamma distribution
\begin{align} p_\varGamma (k) = \left (\frac {k}{k_0}\right )^{\alpha - 1} \frac {\exp (-k/k_0)}{\varGamma (\alpha ) k_0}. \end{align}
For the numerical mapping using the Smirnov transform, conductivity values are organised in logarithmically spaced bins. The values of the exponent
$\alpha$
are chosen to be between zero and one. Furthermore, we set
$k_0 = 1/\alpha$
and
$k_{\textit{min}} = 10^{-5} k_0$
so that the arithmetic mean conductivity is approximately one.
The flow equation (2.2) is solved using a finite volume method based on a two-point flux approximation, as described in Aarnes et al. (Reference Aarnes, Gimse and Lie2007). The mean flow velocity
$\langle u \rangle$
is imposed along the
$1$
-direction. The advective time is thus defined as
$\tau _v = \ell _Y/\langle u \rangle$
. In the perpendicular directions, periodic boundary conditions are imposed for hydraulic head. This scheme on a regular grid is equivalent to a finite difference scheme. It provides the hydraulic head at the cell centres and the flow velocities perpendicular to the cell faces. In order to determine the local velocity gradient the flow field is interpolated by a third-order spline for the velocity components
$u_1(\boldsymbol{x})$
and
$u_2(\boldsymbol{x})$
. Figure 2 illustrates the hydraulic conductivity distribution and the resulting Darcy flow field.
B.2 Displacement moments and first-passage times
For the numerical calculation of displacement moments and first-passage time distributions, we use the Pollock method (Pollock Reference Pollock1988) to numerically integrate the kinematic equation (2.10). Particles are injected uniformly along a straight line perpendicular to the mean flow direction. The line is located at a distance of
$5 \ell _Y$
from the left domain boundary to avoid boundary effects and represents
$90\%$
of the perpendicular length of the domain.
The longitudinal displacement mean and variance are determined from the particle positions by averaging of
$N_p = 10^4$
particles in
$N_R = 200$
medium realisations
\begin{align} \mu (t) = \frac {1}{N_R} \sum \limits _{i = 1}^{N_R} \frac {1}{N_p}\sum \limits _{j = 1}^{N_p} x_1^{(i,j)}(t), && \kappa (t) = \frac {1}{N_R} \sum \limits _{i = 1}^{N_R} \frac {1}{N_p}\sum \limits _{j = 1}^{N_p} \left [x_1^{(i,j)}(t) - \mu (t)\right ]^2, \end{align}
where
$x_1^{(i,j)}(t)$
denotes the longitudinal position of the
$j$
th particle in the
$i$
th disorder realisation. The first-passage times are determined according to (2.26) by the minimum time a particle needs to arrive at the control plane perpendicular to the mean flow direction. The PDF of first-passage times is determined by
\begin{align} f(t,x_1) = \frac {1}{N_R} \sum \limits _{i = 1}^{N_R} \frac {1}{N_p}\sum \limits _{j = 1}^{N_p} \sum _{k = 1}^{N_b} \frac {\mathbb I(t_i \leqslant \tau ^{(i,j)}(x_1) \leqslant t_{k+1})}{t_{k+1} - t_{k}}, \end{align}
where
$\tau ^{(i,j)}$
denotes the first-passage time of particle
$j$
in realisation
$i$
,
$\mathbb I(\boldsymbol{\cdot })$
is the indicator function, which is one if its argument is true and zero else, and the
$t_k$
values are the edges of the time bins to determine the numerical histogram. Due to the broad range of first-passage times, we use a logarithmic binning.
B.3 Lagrangian fluid deformation
The numerical calculation of fluid deformation is based on formulation (2.23) in the Protean coordinate frame, in combination with particle tracking in streamline distance
$s$
as expressed by (2.12). We use an Euler forward scheme to integrate the latter
\begin{align} x_1(s_k + \Delta s) &= x_1(s_k) + \frac { v_1(s_k) \Delta s}{ v(s_k)} , \\[-12pt] \nonumber \end{align}
\begin{align} x_2(s_k + \Delta s) &= x_2(s_k) + \frac { v_2(s_k) \Delta s}{ v(s_k)} , \\[-12pt] \nonumber \end{align}
\begin{align} t(s_k + \Delta s) & = t(s_k) + \frac {\Delta s}{ v(s_k)}, \\[10pt] \nonumber \end{align}
with
$s_k = k \Delta s$
and
$\Delta s = 10^{-2} \ell _Y$
. The system (2.23) is discretised accordingly. It can be written as
\begin{align} z_1(s_{k}) &= \frac { v(s_k)}{v(0)} z_1(0) + v(s_k) v(0) \zeta (s_k) , \\[-12pt] \nonumber \end{align}
\begin{align} z_2(s_{k}) &= \frac {v(0)}{ v(s_k)} z_2(0), \\[10pt] \nonumber \end{align}
where
$ \zeta (s_k)$
represents the integral term on the right side of (2.23a
). It is evaluated by the trapezoidal rule as
\begin{align} \zeta (s_k) = \frac {\Delta s}{2} \left [\frac {\sigma (0)}{v(0)^3} + \frac {\sigma (s_k)}{ v(s_k)^3} \right ] + \sum _{n = 1}^{s_{k-1}} \Delta s \frac { \sigma (s_n)}{ v(s_n)^3}. \end{align}
The initial values
$z_1(0)$
and
$z_2(0)$
are
\begin{align} z_1(0) = \ell _0 \cos (\theta ), \qquad z_2(0) = \ell _0 \sin (\theta ), \end{align}
where the angle
$\theta$
is drawn from a Gaussian distribution with mean
$\pi /2$
and variance one. The observables at a given time
$t$
are obtained by projecting step number
$k$
onto time through
$k_t = \max [k|t(s_k) \leqslant t]$
\begin{align} \boldsymbol z(t) = {\boldsymbol z}(s_{k_t}). \end{align}
Furthermore, we use a two-particle method, which numerically calculates the distance vector
$\boldsymbol z(t) = \boldsymbol{x}(t,\boldsymbol a) - \boldsymbol{x}(t,\boldsymbol a + \delta \boldsymbol a)$
between two particles that are initially spaced by
$\ell _0 = |\delta \boldsymbol a| = 10^{-5} \ell _Y$
. The particle positions are determined by the Pollock method (Pollock Reference Pollock1988). The particle pairs are placed along a line perpendicular to the mean flow direction at a distance of
$5 \ell _Y$
from the left domain boundary. The initial orientation is chosen consistent with (B11) as
\begin{align} \boldsymbol z(t = 0) = \boldsymbol A(0) {\boldsymbol z}(s = 0), \end{align}
where
$\boldsymbol A(0)$
is defined by (2.15).
t=0
t=50τv
t=100τv
τv=ℓY/⟨u⟩
ℓY
ρ
σf2=1
Np=104
σf2=1
105
σf2=1
ℓY=1
t=20τv
t=90τv
z(t,a)
a
x(t,a)
kc
K(x)
K(x)
v
ξ(s)
ξ(s)
ξ(s)
ℓc
exp(−s/ℓc)
vβ−1
β=0.3,0.55,0.85
vm
t−1−β
t2−β
106τv
σf2=6
104τv
α=0.25
σf2=6
t=1,5,20τc
t=20τc
Open access
