# Temporal Evolution of Solute Dispersion in Three-Dimensional Porous Rocks

Alexandre Puyguiraud, Philippe Gouze, Marco Dentz

► To cite this version:

Alexandre Puyguiraud, Philippe Gouze, Marco Dentz. Temporal Evolution of Solute Dispersion in Three-Dimensional Porous Rocks. Water Resources Research, 2024, 60 (3), pp.e2023WR035916. 10.1029/2023wr035916 . insu-04532125

HAL Id: insu-04532125

https://insu.hal.science/insu-04532125

Submitted on 4 Apr 2024

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 Attribution 4.0 International License
Water Resources Research®

RESEARCH ARTICLE

10.1029/2023WR035916

**Key Points:**

* Pore-scale simulations of temporal evolution of solute dispersion in three-dimensional porous rocks

* The impact of medium structure and diffusion on solute transport is studied in terms of effective and ensemble dispersion coefficients

* The evolution of dispersion coefficients indicates how particles sample pore-scale velocities and how this depends on the medium structure

**Correspondence to:**

M. Dentz,
marco.dentz@csic.es

**Citation:**

Puyguiraud, A., Gouze, P., & Dentz, M. (2024). Temporal evolution of solute dispersion in three-dimensional porous rocks. *Water Resources Research*, *60*, e2023WR035916. https://doi.org/10.1029/2023WR035916

Received 24 JULY 2023

Accepted 6 MAR 2024

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

Temporal Evolution of Solute Dispersion in Three-Dimensional Porous Rocks

Alexandre Puyguiraud¹ ID, Philippe Gouze², and Marco Dentz³ ID

¹University of Rennes, CNRS, Géosciences Rennes, UMR 6118, Université de Rennes, Rennes, France, ²Géosciences Montpellier, CNRS-INSU - Montpellier University, Montpellier, France, ³Spanish National Research Council (IDAEA-CSIC), Barcelona, Spain

**Abstract** We study the temporal evolution of solute dispersion in three-dimensional porous rocks of different heterogeneity and pore structure. To this end, we perform direct numerical simulations of pore-scale flow and transport in a sand pack, which exhibits mild heterogeneity, and a Berea sandstone, which is characterized by strong heterogeneity as measured by the variance of the logarithm of the flow velocity. The impact of medium structure and pore-scale mass transfer mechanisms is probed by effective and ensemble dispersion coefficients. The former is a measure for the typical width of a plume, while the latter for the deformation, that is, the spread of a mixing front. Both dispersion coefficients evolve from the molecular diffusion coefficients toward a common finite asymptotic value. Their behavior is governed by the interplay between diffusion, pore-scale velocity fluctuations and medium structure, which determine the characteristic mass transfer time scales. We find distinctly different dispersion behaviors in the two media, which can be traced back to how particles sample pore-scale velocity variability and how this depends on the medium structure. These are key elements for the upscaling of transport, mixing and reaction from the pore to the continuum scale.

# 1. Introduction

The transport of solutes in porous media is driven by the phenomenon of dispersion, which results from the interplay between advective spreading and diffusion. The former is triggered by the spatial variability of the fluid speed which is controlled by the geometry of the connected pore network (Alim et al., 2017; Datta et al., 2013; Puyguiraud et al., 2021; Valocchi et al., 2018) while the latter is ubiquitously controlled by the concentration gradients. The heterogeneity of the porous medium, which determines the flow distribution is therefore a primary parameter that controls dispersion from the pre-asymptotic to the Fickian regimes (Dentz et al., 2004; Sherman et al., 2021). Transport in porous media is considered in many fields of academic and industrial applications from material science, engineering and medicine to groundwater hydrology, environmental technologies and petroleum engineering, and at many scales from microfluidic applications to groundwater management. Beside being necessary for understanding and predicting the spreading of chemicals such as pollutants or bionutrients, modeling dispersion is required also to understand and predict solute-solute and solute-minerals reactions that can produce new solute species and trigger mineral dissolution and precipitation features, for instance.

Dispersion in porous media has been extensively studied from the pore to the regional scale for decades (Dagan, 1990; Dentz et al., 2023; Gelhar & Axness, 1983; Saffman, 1959; Whitaker, 1967). Here we focus on hydrodynamic dispersion due to velocity fluctuations caused by the heterogeneity of the pore space. A main challenge concerns how continuum-scale solute transport can be modeled by macroscopic parameters, such as the dispersion coefficient, that can be inferred experimentally, by using direct pore scale simulations or upscaling methods such as volume averaging or stochastic modeling (Ahmadi et al., 1998; Bijeljic & Blunt, 2006; Brenner, 1980; Koch & Brady, 1985; Le Borgne et al., 2011; Puyguiraud et al., 2021; Scheven, 2013; Souzy et al., 2020). Dispersion and its pre-asymptotic behavior have been analyzed in terms of breakthrough curves, the time evolution of the spatial variance of concentration or particle distributions, or directly from particle velocities, using experiments and direct numerical pore scale simulations (Bijeljic et al., 2004; Gouze, Puyguiraud, Porcher, & Dentz, 2021; Gouze, Puyguiraud, Roubinet, & Dentz, 2021; Hulin & Plona, 1989; Khrapitchev & Callaghan, 2003; Puyguiraud et al., 2021; Sole-Mari et al., 2022). These studies show that the pore structure shapes the evolution of dispersion during the pre-asymptotic regime and then determines the asymptotic value.
Water Resources Research®

RESEARCH ARTICLE

10.1029/2023WR035916

**Key Points:**

* Pore-scale simulations of temporal evolution of solute dispersion in three-dimensional porous rocks

* The impact of medium structure and diffusion on solute transport is studied in terms of effective and ensemble dispersion coefficients

* The evolution of dispersion coefficients indicates how particles sample pore-scale velocities and how this depends on the medium structure

**Correspondence to:**

M. Dentz,
marco.dentz@csic.es

**Citation:**

Puyguiraud, A., Gouze, P., & Dentz, M. (2024). Temporal evolution of solute dispersion in three-dimensional porous rocks. *Water Resources Research*, *60*, e2023WR035916. https://doi.org/10.1029/2023WR035916

Received 24 JULY 2023

Accepted 6 MAR 2024

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

Temporal Evolution of Solute Dispersion in Three-Dimensional Porous Rocks

Alexandre Puyguiraud¹ ID, Philippe Gouze², and Marco Dentz³ ID

¹University of Rennes, CNRS, Géosciences Rennes, UMR 6118, Université de Rennes, Rennes, France, ²Géosciences Montpellier, CNRS-INSU - Montpellier University, Montpellier, France, ³Spanish National Research Council (IDAEA-CSIC), Barcelona, Spain

**Abstract** We study the temporal evolution of solute dispersion in three-dimensional porous rocks of different heterogeneity and pore structure. To this end, we perform direct numerical simulations of pore-scale flow and transport in a sand pack, which exhibits mild heterogeneity, and a Berea sandstone, which is characterized by strong heterogeneity as measured by the variance of the logarithm of the flow velocity. The impact of medium structure and pore-scale mass transfer mechanisms is probed by effective and ensemble dispersion coefficients. The former is a measure for the typical width of a plume, while the latter for the deformation, that is, the spread of a mixing front. Both dispersion coefficients evolve from the molecular diffusion coefficients toward a common finite asymptotic value. Their behavior is governed by the interplay between diffusion, pore-scale velocity fluctuations and medium structure, which determine the characteristic mass transfer time scales. We find distinctly different dispersion behaviors in the two media, which can be traced back to how particles sample pore-scale velocity variability and how this depends on the medium structure. These are key elements for the upscaling of transport, mixing and reaction from the pore to the continuum scale.

# 1. Introduction

The transport of solutes in porous media is driven by the phenomenon of dispersion, which results from the interplay between advective spreading and diffusion. The former is triggered by the spatial variability of the fluid speed which is controlled by the geometry of the connected pore network (Alim et al., 2017; Datta et al., 2013; Puyguiraud et al., 2021; Valocchi et al., 2018) while the latter is ubiquitously controlled by the concentration gradients. The heterogeneity of the porous medium, which determines the flow distribution is therefore a primary parameter that controls dispersion from the pre-asymptotic to the Fickian regimes (Dentz et al., 2004; Sherman et al., 2021). Transport in porous media is considered in many fields of academic and industrial applications from material science, engineering and medicine to groundwater hydrology, environmental technologies and petroleum engineering, and at many scales from microfluidic applications to groundwater management. Beside being necessary for understanding and predicting the spreading of chemicals such as pollutants or bionutrients, modeling dispersion is required also to understand and predict solute-solute and solute-minerals reactions that can produce new solute species and trigger mineral dissolution and precipitation features, for instance.

Dispersion in porous media has been extensively studied from the pore to the regional scale for decades (Dagan, 1990; Dentz et al., 2023; Gelhar & Axness, 1983; Saffman, 1959; Whitaker, 1967). Here we focus on hydrodynamic dispersion due to velocity fluctuations caused by the heterogeneity of the pore space. A main challenge concerns how continuum-scale solute transport can be modeled by macroscopic parameters, such as the dispersion coefficient, that can be inferred experimentally, by using direct pore scale simulations or upscaling methods such as volume averaging or stochastic modeling (Ahmadi et al., 1998; Bijeljic & Blunt, 2006; Brenner, 1980; Koch & Brady, 1985; Le Borgne et al., 2011; Puyguiraud et al., 2021; Scheven, 2013; Souzy et al., 2020). Dispersion and its pre-asymptotic behavior have been analyzed in terms of breakthrough curves, the time evolution of the spatial variance of concentration or particle distributions, or directly from particle velocities, using experiments and direct numerical pore scale simulations (Bijeljic et al., 2004; Gouze, Puyguiraud, Porcher, & Dentz, 2021; Gouze, Puyguiraud, Roubinet, & Dentz, 2021; Hulin & Plona, 1989; Khrapitchev & Callaghan, 2003; Puyguiraud et al., 2021; Sole-Mari et al., 2022). These studies show that the pore structure shapes the evolution of dispersion during the pre-asymptotic regime and then determines the asymptotic value.
Hulin and Plona (1989) and Khrapitchev and Callaghan (2003) study the reversibility of pore-scale dispersion upon flow reversal, which addresses the issue of under which conditions hydrodynamic dispersion describes solute mixing or advective solute spreading. As mentioned above, the fundamental mechanisms of hydrodynamic dispersion are pore-scale velocity fluctuations and diffusion. The former mechanism is reversible in the Stokes regime, which holds for typical applications in groundwater resources. Irreversibility, or actual solute mixing is induced by the interaction of spatial velocity fluctuations and molecular diffusion (Dentz et al., 2023). Consider for example, a solute that evolves from an extended areal source. At early times, the solute front deforms due to velocity variability within the source distribution, which leads to a complex concentration distribution, which nevertheless is partially reversible. Hydrodynamic dispersion coefficients that are defined in terms of the spatial variance of the global solute distribution, measure at pre-asymptotic this advective spreading rather than actual solute mixing. This issue was recognized by Kitanidis (1988) in the context of solute dispersion in heterogeneous porous formations, and Bouchaud and Georges (1990) in the context of random walks in quenched disordered media. These authors propose to define dispersion coefficients from the second-centered moments of the solute or particle distributions that evolve from a point-like initial condition. In the absence of local scale dispersion or molecular diffusion, these dispersion coefficients are exactly zero. In the following, we refer to this concept as effective dispersion. The dispersion concept based on the spatial variance of the solute concentration evolving from an extended areal source, is called ensemble dispersion because it describes the dispersion of the ensemble of partial plumes (Green functions) of which the solute distribution is composed. As outlined above, at pre-asymptotic times ensemble dispersion measures advective solute spreading. In fact, it measures the center of mass fluctuations of the partial plumes that evolve from the point injections that constitute the spatially extended initial distribution (Bouchaud & Georges, 1990). Thus, the two different dispersion coefficients and their temporal evolution highlight different aspects of pore-scale transport, and allow to study the impact of medium heterogeneity and small scale mass transfer mechanisms on large scale dispersion.

The evolution of effective and ensemble dispersion coefficients have been used to quantify the impact of spatial heterogeneity on solute mixing and spreading in Darcy-scale heterogeneous porous media (Attinger et al., 1999; de Barros & Dentz, 2016; De Barros et al., 2015; Dentz & de Barros, 2015; Dentz et al., 2000; Fiori, 2001; Fiori & Dagan, 2000; Vanderborght, 2001). Numerical pore-scale studies have mostly focused on the temporal behavior of longitudinal ensemble dispersion coefficients (Bijeljic & Blunt, 2006; Bijeljic et al., 2011; Gouze, Puyguiraud, Porcher, & Dentz, 2021; Puyguiraud et al., 2019, 2021; Sole-Mari et al., 2022), that is, on solute spreading, whose evolution is dominated by pore-scale velocity variability.

Here, we use the longitudinal effective and ensemble dispersion coefficients to investigate pore-scale mixing and spreading, and to shed new light on the role of the medium structure and pore-scale diffusion for continuum scale solute transport. To this end, we perform three-dimensional direct numerical simulations of pore-scale flow and solute transport in a sand-pack medium and in a Berea sandstone of distinctly different heterogeneity levels, which can be measured, for instance, by the variance of the logarithm of the flow velocity distribution. While the time behaviors of both the effective and ensemble dispersion coefficients are governed by the interplay between pore-scale diffusion and velocity fluctuations, they evolve on different characteristic time scales, which reflect the relevant pore-scale mass transfer and heterogeneity mechanisms. These findings provide new insights into the pore-scale processes that govern continuum scale solute dispersion and mixing.

The paper is organized as follows. The methodology used to calculate flow and transport and measure dispersion are presented in Section 2. Section 3 provides a detailed study of the temporal evolution of the center of mass positions and the longitudinal effective and ensemble dispersion coefficients for the sand pack and Berea samples. It discusses the different time regimes and the implications for our understanding of the impact of pore-scale heterogeneity and mass transfer mechanisms on continuum scale transport.

## 2. Methodology

### 2.1. Pore-Scale Flow and Transport

Flow in three-dimensional porous media, described as dual solid-void structures, is described by the Stokes equation together with the continuity equation (Leal, 2007),

$$ \nabla^2 \mathbf{u}(\mathbf{x}) = -\frac{1}{\mu} \nabla p(\mathbf{x}), \qquad \nabla \cdot \mathbf{u}(\mathbf{x}) = 0, \qquad (1) $$
Hulin and Plona (1989) and Khrapitchev and Callaghan (2003) study the reversibility of pore-scale dispersion upon flow reversal, which addresses the issue of under which conditions hydrodynamic dispersion describes solute mixing or advective solute spreading. As mentioned above, the fundamental mechanisms of hydrodynamic dispersion are pore-scale velocity fluctuations and diffusion. The former mechanism is reversible in the Stokes regime, which holds for typical applications in groundwater resources. Irreversibility, or actual solute mixing is induced by the interaction of spatial velocity fluctuations and molecular diffusion (Dentz et al., 2023). Consider for example, a solute that evolves from an extended areal source. At early times, the solute front deforms due to velocity variability within the source distribution, which leads to a complex concentration distribution, which nevertheless is partially reversible. Hydrodynamic dispersion coefficients that are defined in terms of the spatial variance of the global solute distribution, measure at pre-asymptotic this advective spreading rather than actual solute mixing. This issue was recognized by Kitanidis (1988) in the context of solute dispersion in heterogeneous porous formations, and Bouchaud and Georges (1990) in the context of random walks in quenched disordered media. These authors propose to define dispersion coefficients from the second-centered moments of the solute or particle distributions that evolve from a point-like initial condition. In the absence of local scale dispersion or molecular diffusion, these dispersion coefficients are exactly zero. In the following, we refer to this concept as effective dispersion. The dispersion concept based on the spatial variance of the solute concentration evolving from an extended areal source, is called ensemble dispersion because it describes the dispersion of the ensemble of partial plumes (Green functions) of which the solute distribution is composed. As outlined above, at pre-asymptotic times ensemble dispersion measures advective solute spreading. In fact, it measures the center of mass fluctuations of the partial plumes that evolve from the point injections that constitute the spatially extended initial distribution (Bouchaud & Georges, 1990). Thus, the two different dispersion coefficients and their temporal evolution highlight different aspects of pore-scale transport, and allow to study the impact of medium heterogeneity and small scale mass transfer mechanisms on large scale dispersion.

The evolution of effective and ensemble dispersion coefficients have been used to quantify the impact of spatial heterogeneity on solute mixing and spreading in Darcy-scale heterogeneous porous media (Attinger et al., 1999; de Barros & Dentz, 2016; De Barros et al., 2015; Dentz & de Barros, 2015; Dentz et al., 2000; Fiori, 2001; Fiori & Dagan, 2000; Vanderborght, 2001). Numerical pore-scale studies have mostly focused on the temporal behavior of longitudinal ensemble dispersion coefficients (Bijeljic & Blunt, 2006; Bijeljic et al., 2011; Gouze, Puyguiraud, Porcher, & Dentz, 2021; Puyguiraud et al., 2019, 2021; Sole-Mari et al., 2022), that is, on solute spreading, whose evolution is dominated by pore-scale velocity variability.

Here, we use the longitudinal effective and ensemble dispersion coefficients to investigate pore-scale mixing and spreading, and to shed new light on the role of the medium structure and pore-scale diffusion for continuum scale solute transport. To this end, we perform three-dimensional direct numerical simulations of pore-scale flow and solute transport in a sand-pack medium and in a Berea sandstone of distinctly different heterogeneity levels, which can be measured, for instance, by the variance of the logarithm of the flow velocity distribution. While the time behaviors of both the effective and ensemble dispersion coefficients are governed by the interplay between pore-scale diffusion and velocity fluctuations, they evolve on different characteristic time scales, which reflect the relevant pore-scale mass transfer and heterogeneity mechanisms. These findings provide new insights into the pore-scale processes that govern continuum scale solute dispersion and mixing.

The paper is organized as follows. The methodology used to calculate flow and transport and measure dispersion are presented in Section 2. Section 3 provides a detailed study of the temporal evolution of the center of mass positions and the longitudinal effective and ensemble dispersion coefficients for the sand pack and Berea samples. It discusses the different time regimes and the implications for our understanding of the impact of pore-scale heterogeneity and mass transfer mechanisms on continuum scale transport.

## 2. Methodology

### 2.1. Pore-Scale Flow and Transport

Flow in three-dimensional porous media, described as dual solid-void structures, is described by the Stokes equation together with the continuity equation (Leal, 2007),

$$ \nabla^2 \mathbf{u}(\mathbf{x}) = -\frac{1}{\mu} \nabla p(\mathbf{x}), \qquad \nabla \cdot \mathbf{u}(\mathbf{x}) = 0, \qquad (1) $$
where μ is the dynamic viscosity, **u**(**x**) is the Eulerian velocity and p(**x**) is the fluid pressure at position **x** = (x₁, x₂, x₃). Here, flow is driven by the macroscopic pressure gradient, which is aligned with the x-axis of the coordinated system. Zero-flux boundary conditions are set at the solid-void interface and at the lateral domain boundaries.

Transport of solutes is described by the advection-diffusion equation (ADE) for the solute concentration c(**x**, t)

$$ \frac{\partial c(\mathbf{x}, t)}{\partial t} + \nabla \cdot [\mathbf{u}(\mathbf{x}) - D\nabla] c(\mathbf{x}, t) = 0, \qquad (2) $$

where c(**x**, t) is the solute concentration at position **x** and time t, and D is the molecular diffusion coefficient. The advection-diffusion Equation 2 is equivalent to the Langevin equation (Risken, 1996)

$$ \frac{d\mathbf{x}(t)}{dt} = \mathbf{u}[\mathbf{x}(t)] + \sqrt{2D\xi(t)}, \qquad (3) $$

where ξ(t) is a Gaussian white noise with mean ⟨ξᵢ⟩ = 0 and covariance ⟨ξⱼ(t)ξₖ(t)⟩ = δⱼₖδ(t - t'); δⱼₖ is the Kronecker delta.

The characteristic pore length ℓ₀, the mean velocity magnitude ⟨v⟩ = ⟨|v(**x|)⟩ and the diffusion coefficient D set the advection time τᵥ = ℓ₀/⟨v⟩ and the characteristic diffusion time τₔ = ℓ₀²/D. The two time scales are compared by the Péclet number Pe = τₔ/τᵥ = ⟨v⟩ℓ₀/D.

## 2.2. Dispersion Coefficients

We analyze pore-scale solute dispersion and the underlying heterogeneity and transport mechanisms in terms of the temporal evolution of suitably defined longitudinal dispersion coefficients. To this end, we consider the concentration distribution c(**x**, t) that evolves from the normalized plane source

$$ c(\mathbf{x}, t = 0) = \rho(\mathbf{x}) = \phi^{-1} \delta(x_1) \frac{\mathbb{I}(\mathbf{x} \in \Omega_f)}{wh}, \qquad (4) $$

where Ω_f denotes the fluid domain. We choose a uniform injection across a plane, which mimics the mixing interface between two initially segregated constituents. The indicator function I(·) is one if its argument is true and zero else, w and h denote the width and height of the medium and φ is porosity. The injection plane is large enough such that

$$ \int_{\Omega} d\mathbf{x}\rho(\mathbf{x}) = \phi, \qquad (5) $$

where Ω denotes the bulk domain, that is, the union of fluid domain and solid domain. The solute distribution can be decomposed into partial plumes (Green functions) g(**x**, t|**x'**) that satisfy Equation 2 for the initial conditions

$$ g(\mathbf{x}, t = 0|\mathbf{x}') = \delta(\mathbf{x} - \mathbf{x}')\mathbb{I}(\mathbf{x}' \in \Omega_f). \qquad (6) $$

The concentration distribution c(**x**, t) can be written in terms of g(**x**, t|**x'**) as

$$ c(\mathbf{x}, t) = \int_{\Omega} d\mathbf{x}' \rho(\mathbf{x}') g(\mathbf{x}, t|\mathbf{x}'). \qquad (7) $$

### 2.2.1. Effective Dispersion

We define the longitudinal effective dispersion coefficient in terms of the second centered moments of the Green function g(**x**, t|**x'**)
where $\mu$ is the dynamic viscosity, $\mathbf{u}(\mathbf{x})$ is the Eulerian velocity and $p(\mathbf{x})$ is the fluid pressure at position $\mathbf{x} = (x_1, x_2, x_3)$. Here, flow is driven by the macroscopic pressure gradient, which is aligned with the x-axis of the coordinated system. Zero-flux boundary conditions are set at the solid-void interface and at the lateral domain boundaries.

Transport of solutes is described by the advection-diffusion equation (ADE) for the solute concentration $c(\mathbf{x}, t)$

$$ \frac{\partial c(\mathbf{x}, t)}{\partial t} + \nabla \cdot [\mathbf{u}(\mathbf{x}) - D\nabla] c(\mathbf{x}, t) = 0, \quad (2) $$

where $c(\mathbf{x}, t)$ is the solute concentration at position $\mathbf{x}$ and time $t$, and $D$ is the molecular diffusion coefficient. The advection-diffusion Equation 2 is equivalent to the Langevin equation (Risken, 1996)

$$ \frac{d\mathbf{x}(t)}{dt} = \mathbf{u}[\mathbf{x}(t)] + \sqrt{2D\xi(t)}, \quad (3) $$

where $\xi(t)$ is a Gaussian white noise with mean $\langle \xi_i \rangle = 0$ and covariance $\langle \xi_j(t)\xi_k(t') \rangle = \delta_{jk}\delta(t-t')$; $\delta_{jk}$ is the Kronecker delta.

The characteristic pore length $\ell_0$, the mean velocity magnitude $\langle v \rangle = \langle |v(\mathbf{x})| \rangle$ and the diffusion coefficient $D$ set the advection time $\tau_v = \ell_0/\langle v \rangle$ and the characteristic diffusion time $\tau_D = \ell_0^2/D$. The two time scales are compared by the Péclet number $Pe = \tau_D/\tau_v = \langle v \rangle \ell_0/D$.

## 2.2. Dispersion Coefficients

We analyze pore-scale solute dispersion and the underlying heterogeneity and transport mechanisms in terms of the temporal evolution of suitably defined longitudinal dispersion coefficients. To this end, we consider the concentration distribution $c(\mathbf{x}, t)$ that evolves from the normalized plane source

$$ c(\mathbf{x}, t=0) = \rho(\mathbf{x}) = \phi^{-1} \delta(x_1) \frac{\mathbb{I}(\mathbf{x} \in \Omega_f)}{wh}, \quad (4) $$

where $\Omega_f$ denotes the fluid domain. We choose a uniform injection across a plane, which mimics the mixing interface between two initially segregated constituents. The indicator function $\mathbb{I}(\cdot)$ is one if its argument is true and zero else, $w$ and $h$ denote the width and height of the medium and $\phi$ is porosity. The injection plane is large enough such that

$$ \int_{\Omega} d\mathbf{x} \rho(\mathbf{x}) = \phi, \quad (5) $$

where $\Omega$ denotes the bulk domain, that is, the union of fluid domain and solid domain. The solute distribution can be decomposed into partial plumes (Green functions) $g(\mathbf{x}, t|\mathbf{x}')$ that satisfy Equation 2 for the initial conditions

$$ g(\mathbf{x}, t=0|\mathbf{x}') = \delta(\mathbf{x} - \mathbf{x}') \mathbb{I}(\mathbf{x}' \in \Omega_f). \quad (6) $$

The concentration distribution $c(\mathbf{x}, t)$ can be written in terms of $g(\mathbf{x}, t|\mathbf{x}')$ as

$$ c(\mathbf{x}, t) = \int_{\Omega} d\mathbf{x}' \rho(\mathbf{x}') g(\mathbf{x}, t|\mathbf{x}'). \quad (7) $$

### 2.2.1. Effective Dispersion

We define the longitudinal effective dispersion coefficient in terms of the second centered moments of the Green function $g(\mathbf{x}, t|\mathbf{x}')$
$$ \kappa_L(t; \mathbf{x}') = m_2(t; \mathbf{x}') - m_1(t; \mathbf{x}')^2, \qquad (8) $$

which are measures for the spatial extension of the Green function. The ith longitudinal moment of g(x, t|x') is defined by

$$ m_i(t; \mathbf{x}') = \int d\mathbf{x} x_i^i g(\mathbf{x}, t|\mathbf{x}'). \qquad (9) $$

The first moments $m_1(t; \mathbf{x}')$ determine the center of mass position of $g(\mathbf{x}, t|\mathbf{x}')$. The average of $\kappa_L(t; \mathbf{x}')$ over all Green functions defines the effective second centered moment

$$ \kappa_L^{\text{eff}}(t) = \int d\mathbf{x}' \rho(\mathbf{x}') \kappa_L(t; \mathbf{x}'). \qquad (10) $$

It is a measure for the average streamwise width of the Green function. Half the temporal rate of growth of $\kappa_L^{\text{eff}}(t)$ defines the effective dispersion coefficients

$$ D_L^{\text{eff}}(t) = \frac{1}{2} \frac{d}{dt} \kappa_L^{\text{eff}}(t). \qquad (11) $$

Thus, effective dispersion characterizes the typical spreading behavior for a solute evolving from a point-like injection.

### 2.2.2. Ensemble Dispersion

The longitudinal ensemble dispersion coefficients are defined in terms of the moments of the global concentration $c(\mathbf{x}, t)$,

$$ m_i(t) = \int d\mathbf{x}^i c(\mathbf{x}, t) = \int d\mathbf{x}' \rho(\mathbf{x}') m_i(t; \mathbf{x}'). \qquad (12) $$

As per the second equality sign, the moments can be seen as averages over the ensemble of Green functions and as such are named the ensemble moments in the following. The second centered ensemble moment is defined by

$$ \kappa_L^{\text{ens}}(t) = m_2(t) - m_1(t)^2. \qquad (13) $$

It is a measure for the streamwise spatial extension of the concentration distribution, or equivalently for the ensemble of Green functions. The temporal rate of growth of the second centered ensemble moments is measured by the ensemble dispersion coefficients

$$ D_{ij}^{\text{ens}}(t) = \frac{1}{2} \frac{d}{dt} \kappa_{ij}^{\text{ens}}(t). \qquad (14) $$

### 2.2.3. Center of Mass Fluctuations

The difference between the ensemble and effective variances,

$$ \delta \kappa_L^m(t) = \int d\mathbf{x}' \rho(\mathbf{x}') [m_1(t; \mathbf{x}') - m_1(t)]^2, \qquad (15) $$

quantifies the variance of the center of mass fluctuations of the Green functions that constitute the solute plume. Along the same lines, the difference between the ensemble and effective dispersion coefficients measures the dispersion of the center of mass positions of the Green functions
$$ \kappa_L(t; \mathbf{x}') = m_2(t; \mathbf{x}') - m_1(t; \mathbf{x}')^2, \qquad (8) $$

which are measures for the spatial extension of the Green function. The ith longitudinal moment of g(x, t|x') is defined by

$$ m_i(t; \mathbf{x}') = \int d\mathbf{x} x_i^i g(\mathbf{x}, t|\mathbf{x}'). \qquad (9) $$

The first moments $m_1(t; \mathbf{x}')$ determine the center of mass position of g(x, t|x'). The average of $\kappa_L(t; \mathbf{x}')$ over all Green functions defines the effective second centered moment

$$ \kappa_L^{\text{eff}}(t) = \int d\mathbf{x}' \rho(\mathbf{x}') \kappa_L(t; \mathbf{x}'). \qquad (10) $$

It is a measure for the average streamwise width of the Green function. Half the temporal rate of growth of $\kappa_L^{\text{eff}}(t)$ defines the effective dispersion coefficients

$$ D_L^{\text{eff}}(t) = \frac{1}{2} \frac{d}{dt} \kappa_L^{\text{eff}}(t). \qquad (11) $$

Thus, effective dispersion characterizes the typical spreading behavior for a solute evolving from a point-like injection.

### 2.2.2. Ensemble Dispersion

The longitudinal ensemble dispersion coefficients are defined in terms of the moments of the global concentration $c(\mathbf{x}, t)$,

$$ m_i(t) = \int d\mathbf{x}^i c(\mathbf{x}, t) = \int d\mathbf{x}' \rho(\mathbf{x}') m_i(t; \mathbf{x}'). \qquad (12) $$

As per the second equality sign, the moments can be seen as averages over the ensemble of Green functions and as such are named the ensemble moments in the following. The second centered ensemble moment is defined by

$$ \kappa_L^{\text{ens}}(t) = m_2(t) - m_1(t)^2. \qquad (13) $$

It is a measure for the streamwise spatial extension of the concentration distribution, or equivalently for the ensemble of Green functions. The temporal rate of growth of the second centered ensemble moments is measured by the ensemble dispersion coefficients

$$ D_{ij}^{\text{ens}}(t) = \frac{1}{2} \frac{d}{dt} \kappa_{ij}^{\text{ens}}(t). \qquad (14) $$

### 2.2.3. Center of Mass Fluctuations

The difference between the ensemble and effective variances,

$$ \delta\kappa_L^m(t) = \int d\mathbf{x}' \rho(\mathbf{x}') [m_1(t; \mathbf{x}') - m_1(t)]^2, \qquad (15) $$

quantifies the variance of the center of mass fluctuations of the Green functions that constitute the solute plume. Along the same lines, the difference between the ensemble and effective dispersion coefficients measures the dispersion of the center of mass positions of the Green functions
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Eulerian velocity pdfs for the sand pack (blue circles) and the Berea sandstone (red squares). Inlay: The three-dimensional pore geometry of (left) the sand pack sample (5 mm³) and of (right) the Berea sandstone (1 mm³). The gray and blue colors represent the pore space and the solid phase, respectively.</figcaption></figure>

$$ \delta D_{ij}^{m}(t) = \frac{1}{2} \frac{d}{dt} \delta \kappa_{ij}^{m}(t). \qquad (16) $$

## 2.3. Numerical Simulations

In the following, we describe the studied porous media, the numerical solution of the pore-scale flow problem and of the transport problem using random walk particle tracking.

### 2.3.1. Porous Media

We study two three-dimensional porous media of different complexity, (a) a Berea sandstone sample and (b) a sand pack sample, which are illustrated in Figure 1. The Berea sample displays a complex pore structure with a porosity of $\phi = 0.18$ and a permeability of $6.7 \times 10^{-12} \text{ m}^2$ (Gouze, Puyguiraud, Roubinet, & Dentz, 2021). Sandstones are more or less diagenized (altered and cemented) buried sand structures and the Berea sandstone is often considered as a typical example of high permeability sedimentary reservoir (Churcher et al., 1991). The sand pack sample has a high porosity of $\phi = 0.37$ with a more regular structure of the pore space and a higher permeability value of $2.9 \times 10^{-11} \text{ m}^2$. For both samples the characteristic pore diameter

was evaluated from image analysis to be about $4.5 \times 10^{-5} \text{ m}$ and the characteristic length is $l_0 = 1.5 \times 10^{-4} \text{ m}$ (Gouze, Puyguiraud, Roubinet, & Dentz, 2021). The sand-pack image (Sand Pack LV60C) is obtained from the Imperial College image repository (Imperial College Consortium on Pore-scale Imaging and Modelling, 2014). It is a compact packing of irregular quartz grains of variable size that is a proxy of sub-surface aquifers (Di Palma et al., 2019). The different degrees of heterogeneity of the two samples can be quantified in terms of the distribution of flow speeds as discussed in the next section.

### 2.3.2. Fluid Flow

As mentioned above, both pore geometries are based on X-Ray microtomography images. The geometries are meshed using regular hexahedron cells (voxels). This type of mesh has two major advantages. First, it perfectly fits the voxels of the X-Ray tomography images, and second, it allows for a simple and computationally efficient velocity interpolation scheme, which is required for the transport simulation based on random walk particle tracking (Mostaghimi et al., 2012). Each of the images is decomposed in $900^3$ voxels of length $l_m = 1.060 \times 10^{-6} \text{ m}$ for the Berea sandstone and $l_m = 5.001 \times 10^{-6} \text{ m}$ for the sand pack. This implies that the throat (or pore constriction) diameters are resolved in average with 9 voxels for the sand pack and 45 voxels for the Berea sample, which ensures an accurate computation of the flow field (Gouze, Puyguiraud, Porcher, & Dentz, 2021). Note that the voxel size for the sand pack is comparable to the one chosen in Sole-Mari et al. (2022). Fluid flow in the pore space is solved numerically using the SIMPLE algorithm implemented in OpenFOAM (Weller et al., 1998). Pressure boundary conditions are set at the inlet ($x = 0$) and outlet ($x = 900l_m$) of the domains. No-slip boundary conditions are prescribed at the void-solid interface and at the lateral boundaries of the domain. Note that imposing impermeable lateral boundaries impacts transverse dispersion but hardly affects longitudinal dispersion and mixing which are controlled by the velocity variation along the trajectory and transverse diffusion within a pore (or throat). Once the solver has converged, the flow velocities are extracted at the centers of the interfaces of the mesh (i.e., at the six faces of each of the regular hexahedra that form the mesh) in the normal direction to the face.

The ratio between the mean flow speed $\langle v \rangle$ and the mean flow velocity $\langle u \rangle$ in streamwise direction defines the advective tortuosity $\chi = \langle v \rangle/\langle u \rangle$. For the Berea sample, we find $\chi = 1.64$, and for the sand pack $\chi = 1.32$. Since for Stokes flow, the flow velocities scale with the pressure gradient, the flow field is determined for a unit gradient and then scaled for the Péclet scenario under consideration. For example, for $Pe = 200$, the mean flow speeds are $\langle v \rangle = 2.67 \times 10^{-3} \text{ m/s}$. The mean streamwise velocities can be obtained from the respective tortuosity values. Figure 1 shows the distribution of flow speeds for the two media, which can be used to quantify the flow heterogeneity of the two porous medium samples and highlight the differences between them (Alhashmi et al., 2016).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Eulerian velocity pdfs for the sand pack (blue circles) and the Berea sandstone (red squares). Inlay: The three-dimensional pore geometry of (left) the sand pack sample (5 mm³) and of (right) the Berea sandstone (1 mm³). The gray and blue colors represent the pore space and the solid phase, respectively.</figcaption></figure>

$$ \delta D_{ij}^{m}(t) = \frac{1}{2} \frac{d}{dt} \delta \kappa_{ij}^{m}(t). \qquad (16) $$

## 2.3. Numerical Simulations

In the following, we describe the studied porous media, the numerical solution of the pore-scale flow problem and of the transport problem using random walk particle tracking.

### 2.3.1. Porous Media

We study two three-dimensional porous media of different complexity, (a) a Berea sandstone sample and (b) a sand pack sample, which are illustrated in Figure 1. The Berea sample displays a complex pore structure with a porosity of $\phi = 0.18$ and a permeability of $6.7 \times 10^{-12} \text{ m}^2$ (Gouze, Puyguiraud, Roubinet, & Dentz, 2021). Sandstones are more or less diagenized (altered and cemented) buried sand structures and the Berea sandstone is often considered as a typical example of high permeability sedimentary reservoir (Churcher et al., 1991). The sand pack sample has a high porosity of $\phi = 0.37$ with a more regular structure of the pore space and a higher permeability value of $2.9 \times 10^{-11} \text{ m}^2$. For both samples the characteristic pore diameter

was evaluated from image analysis to be about $4.5 \times 10^{-5} \text{ m}$ and the characteristic length is $l_0 = 1.5 \times 10^{-4} \text{ m}$ (Gouze, Puyguiraud, Roubinet, & Dentz, 2021). The sand-pack image (Sand Pack LV60C) is obtained from the Imperial College image repository (Imperial College Consortium on Pore-scale Imaging and Modelling, 2014). It is a compact packing of irregular quartz grains of variable size that is a proxy of sub-surface aquifers (Di Palma et al., 2019). The different degrees of heterogeneity of the two samples can be quantified in terms of the distribution of flow speeds as discussed in the next section.

### 2.3.2. Fluid Flow

As mentioned above, both pore geometries are based on X-Ray microtomography images. The geometries are meshed using regular hexahedron cells (voxels). This type of mesh has two major advantages. First, it perfectly fits the voxels of the X-Ray tomography images, and second, it allows for a simple and computationally efficient velocity interpolation scheme, which is required for the transport simulation based on random walk particle tracking (Mostaghimi et al., 2012). Each of the images is decomposed in $900^3$ voxels of length $l_m = 1.060 \times 10^{-6} \text{ m}$ for the Berea sandstone and $l_m = 5.001 \times 10^{-6} \text{ m}$ for the sand pack. This implies that the throat (or pore constriction) diameters are resolved in average with 9 voxels for the sand pack and 45 voxels for the Berea sample, which ensures an accurate computation of the flow field (Gouze, Puyguiraud, Porcher, & Dentz, 2021). Note that the voxel size for the sand pack is comparable to the one chosen in Sole-Mari et al. (2022). Fluid flow in the pore space is solved numerically using the SIMPLE algorithm implemented in OpenFOAM (Weller et al., 1998). Pressure boundary conditions are set at the inlet ($x = 0$) and outlet ($x = 900l_m$) of the domains. No-slip boundary conditions are prescribed at the void-solid interface and at the lateral boundaries of the domain. Note that imposing impermeable lateral boundaries impacts transverse dispersion but hardly affects longitudinal dispersion and mixing which are controlled by the velocity variation along the trajectory and transverse diffusion within a pore (or throat). Once the solver has converged, the flow velocities are extracted at the centers of the interfaces of the mesh (i.e., at the six faces of each of the regular hexahedra that form the mesh) in the normal direction to the face.

The ratio between the mean flow speed $\langle v \rangle$ and the mean flow velocity $\langle u \rangle$ in streamwise direction defines the advective tortuosity $\chi = \langle v \rangle / \langle u \rangle$. For the Berea sample, we find $\chi = 1.64$, and for the sand pack $\chi = 1.32$. Since for Stokes flow, the flow velocities scale with the pressure gradient, the flow field is determined for a unit gradient and then scaled for the Péclet scenario under consideration. For example, for $Pe = 200$, the mean flow speeds are $\langle v \rangle = 2.67 \times 10^{-3} \text{ m/s}$. The mean streamwise velocities can be obtained from the respective tortuosity values. Figure 1 shows the distribution of flow speeds for the two media, which can be used to quantify the flow heterogeneity of the two porous medium samples and highlight the differences between them (Alhashmi et al., 2016).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Snapshots of the conservative simulation for the Berea sandstone for Pe = 2,000 at three different times t = 0.15τ_v, t = 0.8τ_v, and t = 5τ_v. The density of particles represents the concentration.</figcaption></figure>

The flow heterogeneity is measured by the variance $\sigma_f^2$ of the natural logarithm $f = \ln v$ of the flow speed $v$. For the Berea sandstone sample, we obtain $\sigma_f^2 = 10$, for the sand pack sample $\sigma_f^2 = 2$, that is, the Berea sample is significantly more heterogeneous.

### 2.3.3. Random Walk Particle Tracking

Solute transport is modeled using random walk particle tracking (Noetinger et al., 2016). The numerical simulation is based on the discretized version of the Langevin Equation 3,

$$x(t + \Delta t) = x(t) + u[x(t)]\Delta t + \sqrt{2D\Delta t}\xi(t), \quad (17)$$

where $\xi = (\xi_1, \xi_2, \xi_3)$. The $\xi_i$ are independent random variables that are uniformly distributed in $[-\sqrt{3}, \sqrt{3}]$. The central limit theorem ensures that the sum of these uniform random variables is Gaussian distributed with zero mean and unit variance. The particle velocities $u[x(t)]$ are interpolated from the velocities at the voxel faces using the algorithm of Mostaghimi et al. (2012), which implements a quadratic interpolation in the void voxels that are in contact with the solid and thus guarantees an accurate representation of the flow field in the vicinity of the solid-void interface. The time step is variable and chosen such that the particle displacement at a given step is shorter than or equal to the side length of a voxel. The time step varies from $\Delta t = 10^{-8}$ s at early times to get an accurate resolution of the moments to $\Delta t = 10^{-3}$ s at late times to ensure faster simulations. The diffusion coefficient is set to $D = 10^{-9}$ m²/s.

To investigate the effective and ensemble dispersion coefficients, particles are seeded at all voxels in the void space within a plane perpendicular to the mean flow direction, see Figure 2. There are $N_0 = 1.5 \times 10^5$ voxels in the injection plane for the Berea sample, and $N_0 = 3 \times 10^5$ for the sand pack. At each position, $N_p = 10^3$ particles are injected, which makes a total of $N = 1.5 \times 10^8$ particles for the Berea sandstone and $N = 3 \times 10^8$ for the sand pack. We consider this scenario for Pe = 200 and Pe = 2,000. Each partial injection of particles defines a partial plume, or Green function. The first and second moments of a partial plume originating from $x'_j$ are then obtained from the particle positions $x_k(t|x'_j)$

$$m_i(t; \mathbf{x}'_j) = \frac{1}{N_p} \sum_{k=1}^{N_p} x_{i,k}(t; \mathbf{x}'_j)^i, \quad (18)$$

where $i = 1, 2$ and the subscript $k$ denotes the $k$th particle. The effective and ensemble variances are given by

$$\kappa_L^{\text{eff}}(t) = \frac{1}{N_0} \sum_{j=1}^{N_0} \kappa_L(t; \mathbf{x}'_j), \quad (19)$$

$$\kappa_L^{\text{ens}}(t) = \frac{1}{N_0} \sum_{j=1}^{N_0} m_2(t; \mathbf{x}'_j) - \left[ \frac{1}{N_0} \sum_{j=1}^{N_0} m_1(t; \mathbf{x}'_j) \right]^2. \quad (20)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Snapshots of the conservative simulation for the Berea sandstone for Pe = 2,000 at three different times t = 0.15τ<sub>v</sub>, t = 0.8τ<sub>v</sub>, and t = 5τ<sub>v</sub>. The density of particles represents the concentration.</figcaption></figure>

The flow heterogeneity is measured by the variance σ<sub>f</sub><sup>2</sup> of the natural logarithm f = ln v of the flow speed v. For the Berea sandstone sample, we obtain σ<sub>f</sub><sup>2</sup> = 10, for the sand pack sample σ<sub>f</sub><sup>2</sup> = 2, that is, the Berea sample is significantly more heterogeneous.

2.3.3. Random Walk Particle Tracking

Solute transport is modeled using random walk particle tracking (Noetinger et al., 2016). The numerical simulation is based on the discretized version of the Langevin Equation 3,

$$
\mathbf{x}(t + \Delta t) = \mathbf{x}(t) + \mathbf{u}[\mathbf{x}(t)]\Delta t + \sqrt{2D\Delta t}\xi(t), \quad (17)
$$

where ξ = (ξ₁, ξ₂, ξ₃). The ξᵢ are independent random variables that are uniformly distributed in [-√3, √3]. The central limit theorem ensures that the sum of these uniform random variables is Gaussian distributed with zero mean and unit variance. The particle velocities **u**[**x**(*t*)] are interpolated from the velocities at the voxel faces using the algorithm of Mostaghimi et al. (2012), which implements a quadratic interpolation in the void voxels that are in contact with the solid and thus guarantees an accurate representation of the flow field in the vicinity of the solid-void interface. The time step is variable and chosen such that the particle displacement at a given step is shorter than or equal to the side length of a voxel. The time step varies from Δ*t* = 10⁻⁸ s at early times to get an accurate resolution of the moments to Δ*t* = 10⁻³ s at late times to ensure faster simulations. The diffusion coefficient is set to *D* = 10⁻⁹ m²/s.

To investigate the effective and ensemble dispersion coefficients, particles are seeded at all voxels in the void space within a plane perpendicular to the mean flow direction, see Figure 2. There are N₀ = 1.5 × 10⁵ voxels in the injection plane for the Berea sample, and N₀ = 3 × 10⁵ for the sand pack. At each position, Nₚ = 10³ particles are injected, which makes a total of N = 1.5 × 10⁸ particles for the Berea sandstone and N = 3 × 10⁸ for the sand pack. We consider this scenario for Pe = 200 and Pe = 2,000. Each partial injection of particles defines a partial plume, or Green function. The first and second moments of a partial plume originating from **x**'j are then obtained from the particle positions **x**_k (t|**x**'j)

$$
m_i(t; \mathbf{x}'_j) = \frac{1}{N_p} \sum_{k=1}^{N_p} x_{i,k}(t; \mathbf{x}'_j)^i, \quad (18)
$$

where *i* = 1, 2 and the subscript *k* denotes the *k*th particle. The effective and ensemble variances are given by

$$
\kappa_L^{\text{eff}}(t) = \frac{1}{N_0} \sum_{j=1}^{N_0} \kappa_L(t; \mathbf{x}'_j), \quad (19)
$$

$$
\kappa_L^{\text{ens}}(t) = \frac{1}{N_0} \sum_{j=1}^{N_0} m_2(t; \mathbf{x}'_j) - \left[ \frac{1}{N_0} \sum_{j=1}^{N_0} m_1(t; \mathbf{x}'_j) \right]^2 . \quad (20)
$$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Temporal evolution of the center of mass position for (a) sand pack and (b) Berea samples for Pe = 200 and Pe = 2,000. Temporal evolution of the variance of the center of mass position for (c) the sand pack and (d) Berea samples for Pe = 200 and Pe = 2,000. The dashed vertical lines denote (black) the advection time scale $\tau_v$, (yellow and orange) the respective diffusion time scales $\tau_D$.</figcaption></figure>

### 3. Dispersion Behavior

In this section, we analyze the dispersion behavior in the sand pack and Berea samples. Figure 2 displays three snapshots of the concentration distribution for the Berea sandstone at Pe = 2,000. The concentration distribution is heterogeneous and characterized by fast solute transport along preferential flow paths and retention in slow flowing regions. In the following, we discuss the evolution of the mean displacement, and the longitudinal effective and ensemble dispersion coefficients defined in Section 2.2 for the sand pack and the Berea sandstone samples. In the following figures, time is non-dimensionalized by the advection time $\tau_v$.

#### 3.1. Center of Mass

Figure 3 shows the evolution of the streamwise center of mass position $m_1(t)$ of the global solute distribution $c(x, t)$ in the top panels. The bottom panels show the rate of change $\delta D_L^m(t)$ of the variance of the center of mass positions $m_1(t|x')$ of partial plumes $g(x, t|x')$ defined by Equation 16. The center of mass of the global plume moves with the mean flow velocity $\langle u \rangle$, while the center of mass velocities of the partial plumes evolve from the velocities at the respective injection points toward the mean flow velocity. At short times $t \ll \tau_v$, that is, travel distances shorter than the average pore size, the center of mass velocities are approximately constant, which implies $m_1(t; x') = u_1(x')t$ and therefore

$$ \delta D_L^m(t) = \sigma_0^2 t, \qquad (21) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. Temporal evolution of the center of mass position for (a) sand pack and (b) Berea samples for Pe = 200 and Pe = 2,000. Temporal evolution of the variance of the center of mass position for (c) the sand pack and (d) Berea samples for Pe = 200 and Pe = 2,000. The dashed vertical lines denote (black) the advection time scale $\tau_v$, (yellow and orange) the respective diffusion time scales $\tau_D$.</figcaption></figure>

### 3. Dispersion Behavior

In this section, we analyze the dispersion behavior in the sand pack and Berea samples. Figure 2 displays three snapshots of the concentration distribution for the Berea sandstone at Pe = 2,000. The concentration distribution is heterogeneous and characterized by fast solute transport along preferential flow paths and retention in slow flowing regions. In the following, we discuss the evolution of the mean displacement, and the longitudinal effective and ensemble dispersion coefficients defined in Section 2.2 for the sand pack and the Berea sandstone samples. In the following figures, time is non-dimensionalized by the advection time $\tau_v$.

#### 3.1. Center of Mass

Figure 3 shows the evolution of the streamwise center of mass position $m_1(t)$ of the global solute distribution $c(x, t)$ in the top panels. The bottom panels show the rate of change $\delta D_L^m(t)$ of the variance of the center of mass positions $m_1(t|x')$ of partial plumes $g(x, t|x')$ defined by Equation 16. The center of mass of the global plume moves with the mean flow velocity $\langle u \rangle$, while the center of mass velocities of the partial plumes evolve from the velocities at the respective injection points toward the mean flow velocity. At short times $t \ll \tau_v$, that is, travel distances shorter than the average pore size, the center of mass velocities are approximately constant, which implies $m_1(t; x') = u_1(x')t$ and therefore

$$ \delta D_L^m(t) = \sigma_0^2 t, \qquad (21) $$
where $\sigma_0^2$ denotes the initial velocity variability. The initial particle velocities persist until the plume starts sampling the flow field by transverse diffusion across streamlines, and by advection along the streamlines. This ballistic early time regime is observed for both the sand pack and Berea samples.

### 3.1.1. Sand Pack Sample

The evolution of $\delta D_L^m(t)$ for the sand pack sample is characterized by two regimes. The early time ballistic regime, and a sharp decay after a maximum that is assumed on the advective time scale $\tau_v$. This is at first counter-intuitive because transverse diffusion is the only mechanism that makes the partial plume sample the flow heterogeneity such that the differences between the center of mass positions of different partial plumes decrease. Thus, one would expect that the relevant time scale is set by the characteristic pore length and diffusion, that is, by the diffusion time $\tau_D$. Sampling occurs indeed by diffusion in transverse direction. However, the distance $\ell_c$ to sample a new velocity depends on the flow rate because streamtubes in low velocity regions are wider than in high velocity regions. Since the flow rate is constant in a streamtube, $Q_c = \ell_c^2 \langle v \rangle$, with $Q_c$ a characteristic flow rate, the decorrelation length becomes $\ell_c = \sqrt{Q_c/\langle v \rangle}$. Thus, the time scale at which particles decorrelate is

$$ \tau_c = \frac{\ell_c^2}{D} = \frac{Q_c}{D\ell_0} \tau_v. \qquad (22) $$

From Figure 3, we observe that $\tau_c \approx \tau_v$, which means that the characteristic flow rate is $Q_c \approx D\ell_0$.

### 3.1.2. Berea Sandstone Sample

For the Berea sample, we observe three different regimes for $\delta D_L^m(t)$. The early time regime is ballistic as discussed above. The start of the second regime is marked by the advective time scale $\tau_v$ as observed for the sand pack. Here, however, $\delta D_L^m(t)$ does not assume a maximum on the advective time scale and then decays, but keeps increasing until the diffusion time $\tau_D$, where it reaches maximum and then shows a rapid decay. A possible explanation of this behavior in the second time regime is the contrast between particles that are initialized at moderate to high flow velocities, which decorrelate due to transverse velocity sampling on the one hand, and the persistence of particles in low velocity conducts on the other hand, which can give rise to the observed sub-linear increase of $\delta D_L^m(t)$. These low velocities are eliminated on the time scale $\tau_D$, which sets the maximum transition time along a conduct. In other words, transition times of particles that move at low velocities along a conduct are cut-off at the diffusion time scale (Puyguiraud et al., 2021).

In summary, the evolution of the center of mass fluctuations is marked by the advection time scale for the sand pack sample, and by the advection and diffusion time scales for the Berea sample. The fact that the intermediate regime is not present for the sand pack sample can be explained by the spatial medium structures of the two samples shown in Figure 1. The structure of the Berea sample can be seen as a connected network of conducts, while the sand pack is more a connected network of pore bodies. These differences are also reflected in the evolutions of the effective and ensemble dispersion coefficients discussed in the next section.

## 3.2. Ensemble and Effective Dispersion

Figures 4 and 5 show the evolution of the effective and ensemble dispersion coefficients for the sand pack and Berea samples. One observes a marked difference between the ensemble and effective dispersion coefficients at short and intermediate times. At early times $t < \tau_0 = D/\langle v \rangle^2 = Pe^{-1}\tau_v$, diffusion dominates over advection, and both the ensemble and effective dispersion coefficients are equal to the molecular diffusion coefficient $D$. Note that the time scale $\tau_0$ here is much smaller than the diffusion time across a pore throat. Thus, at $t < \tau_0$ the diffusion process is not hindered by the confining porous medium and therefore both dispersion coefficients are equal to the molecular diffusion coefficient.

For $\tau_0 < t < \tau_v$, advection starts dominating over diffusion. As outlined in the previous section, particles are transported at their initial velocities that persist over the characteristic length scale $\ell_0$. Thus, the ensemble dispersion coefficients evolve ballistically in this regime

$$ D_L^{\text{ens}}(t) = \sigma_0^2 t, \qquad (23) $$
where $\sigma_0^2$ denotes the initial velocity variability. The initial particle velocities persist until the plume starts sampling the flow field by transverse diffusion across streamlines, and by advection along the streamlines. This ballistic early time regime is observed for both the sand pack and Berea samples.

### 3.1.1. Sand Pack Sample

The evolution of $\delta D_L^m(t)$ for the sand pack sample is characterized by two regimes. The early time ballistic regime, and a sharp decay after a maximum that is assumed on the advective time scale $\tau_v$. This is at first counter-intuitive because transverse diffusion is the only mechanism that makes the partial plume sample the flow heterogeneity such that the differences between the center of mass positions of different partial plumes decrease. Thus, one would expect that the relevant time scale is set by the characteristic pore length and diffusion, that is, by the diffusion time $\tau_D$. Sampling occurs indeed by diffusion in transverse direction. However, the distance $\ell_c$ to sample a new velocity depends on the flow rate because streamtubes in low velocity regions are wider than in high velocity regions. Since the flow rate is constant in a streamtube, $Q_c = \ell_c^2 \langle v \rangle$, with $Q_c$ a characteristic flow rate, the decorrelation length becomes $\ell_c = \sqrt{Q_c/\langle v \rangle}$. Thus, the time scale at which particles decorrelate is

$$ \tau_c = \frac{\ell_c^2}{D} = \frac{Q_c}{D\ell_0} \tau_v. \qquad (22) $$

From Figure 3, we observe that $\tau_c \approx \tau_v$, which means that the characteristic flow rate is $Q_c \approx D\ell_0$.

### 3.1.2. Berea Sandstone Sample

For the Berea sample, we observe three different regimes for $\delta D_L^m(t)$. The early time regime is ballistic as discussed above. The start of the second regime is marked by the advective time scale $\tau_v$ as observed for the sand pack. Here, however, $\delta D_L^m(t)$ does not assume a maximum on the advective time scale and then decays, but keeps increasing until the diffusion time $\tau_D$, where it reaches maximum and then shows a rapid decay. A possible explanation of this behavior in the second time regime is the contrast between particles that are initialized at moderate to high flow velocities, which decorrelate due to transverse velocity sampling on the one hand, and the persistence of particles in low velocity conducts on the other hand, which can give rise to the observed sub-linear increase of $\delta D_L^m(t)$. These low velocities are eliminated on the time scale $\tau_D$, which sets the maximum transition time along a conduct. In other words, transition times of particles that move at low velocities along a conduct are cut-off at the diffusion time scale (Puyguiraud et al., 2021).

In summary, the evolution of the center of mass fluctuations is marked by the advection time scale for the sand pack sample, and by the advection and diffusion time scales for the Berea sample. The fact that the intermediate regime is not present for the sand pack sample can be explained by the spatial medium structures of the two samples shown in Figure 1. The structure of the Berea sample can be seen as a connected network of conducts, while the sand pack is more a connected network of pore bodies. These differences are also reflected in the evolutions of the effective and ensemble dispersion coefficients discussed in the next section.

## 3.2. Ensemble and Effective Dispersion

Figures 4 and 5 show the evolution of the effective and ensemble dispersion coefficients for the sand pack and Berea samples. One observes a marked difference between the ensemble and effective dispersion coefficients at short and intermediate times. At early times $t < \tau_0 = D/\langle v \rangle^2 = Pe^{-1}\tau_v$, diffusion dominates over advection, and both the ensemble and effective dispersion coefficients are equal to the molecular diffusion coefficient $D$. Note that the time scale $\tau_0$ here is much smaller than the diffusion time across a pore throat. Thus, at $t < \tau_0$ the diffusion process is not hindered by the confining porous medium and therefore both dispersion coefficients are equal to the molecular diffusion coefficient.

For $\tau_0 < t < \tau_v$, advection starts dominating over diffusion. As outlined in the previous section, particles are transported at their initial velocities that persist over the characteristic length scale $\ell_0$. Thus, the ensemble dispersion coefficients evolve ballistically in this regime

$$ D_L^{\text{ens}}(t) = \sigma_0^2 t, \qquad (23) $$
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Dispersion coefficients for the sand pack sample. Top panels: (a) Ensemble dispersion coefficients for (red solid line) Pe = 2,000 and (orange solid line) Pe = 200, and (b) corresponding effective dispersion coefficients. The vertical dashed lines denote the corresponding diffusion time scale $\tau_D = \tau_v Pe$. The dashed lines in panel (a) show expression (26) for (purple) $\alpha = 3.1$ and (blue) $\alpha = 3.9$. The dotted lines show expression (25) for (purple) $\beta = 13.1$ and (blue) $\beta = 16.3$. Bottom panels: (Black solid lines) Ensemble and (blue solid lines) effective dispersion coefficients for (c) Pe = 200 and (d) Pe = 2,000. The vertical black dashed lines denote the decorrelation time scale $\tau_c = \tau_v$, the blue dashed lines denotes the respective diffusion time scales. The horizontal dash-dotted lines denote the asymptotic short time and long time values.</figcaption></figure>

where $\sigma_0^2$ is the initial velocity variance. It behaves in the same way as $\Delta D_L^m(t)$, see Equation 21.

This effect of the center of mass fluctuations between partial plumes is removed by the definition of the effective dispersion coefficients as the average dispersion coefficient of the partial plumes. For $\tau_0 < t < \tau_v$, a partial plume is translated by its initial velocity. As its size increases by diffusion, the plume gets sheared by the transverse velocity contrast. Therefore, the effective dispersion coefficients $D_L^{eff}(t)$ first remain at the value of the molecular diffusion coefficient and then increase steeply due to shear dispersion. Figures 4b and 5b show that the increase of the effective dispersion coefficients occurs for high Pe at earlier non-dimensional times than for low Pe. This indicates that the shear rate does not scale linearly with $\langle u \rangle$. In fact, a typical shear rate can be written as

$$ \gamma = \frac{\langle v \rangle}{\ell_{\gamma}}, \qquad (24) $$

where $\ell_{\gamma}$ is the scale of transverse velocity contrast. The latter is proportional to the typical streamtube size. That is, as $\ell_{\gamma}^2 \langle v \rangle = \text{constant}$, we have $\ell_{\gamma} \sim \langle v \rangle^{-1/2}$. The characteristic shear length scale decreases with increasing flow rate, and thus the shear rate scales as $\gamma \sim \langle u \rangle^{3/2}$. The characteristic shear time scale is then $\tau_{\gamma} = \gamma^{-1} \propto \tau_v / \langle v \rangle^{1/2}$. This dependence explains the differences in the time behaviors of the effective dispersion coefficients for different Pe.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dispersion coefficients for the sand pack sample. Top panels: (a) Ensemble dispersion coefficients for (red solid line) Pe = 2,000 and (orange solid line) Pe = 200, and (b) corresponding effective dispersion coefficients. The vertical dashed lines denote the corresponding diffusion time scale $\tau_D = \tau_v Pe$. The dashed lines in panel (a) show expression (26) for (purple) $\alpha = 3.1$ and (blue) $\alpha = 3.9$. The dotted lines show expression (25) for (purple) $\beta = 13.1$ and (blue) $\beta = 16.3$. Bottom panels: (Black solid lines) Ensemble and (blue solid lines) effective dispersion coefficients for (c) Pe = 200 and (d) Pe = 2,000. The vertical black dashed lines denote the decorrelation time scale $\tau_c = \tau_v$, the blue dashed lines denotes the respective diffusion time scales. The horizontal dash-dotted lines denote the asymptotic short time and long time values.</figcaption></figure>

where $\sigma_0^2$ is the initial velocity variance. It behaves in the same way as $\Delta D_L^m(t)$, see Equation 21.

This effect of the center of mass fluctuations between partial plumes is removed by the definition of the effective dispersion coefficients as the average dispersion coefficient of the partial plumes. For $\tau_0 < t < \tau_v$, a partial plume is translated by its initial velocity. As its size increases by diffusion, the plume gets sheared by the transverse velocity contrast. Therefore, the effective dispersion coefficients $D_L^{eff}(t)$ first remain at the value of the molecular diffusion coefficient and then increase steeply due to shear dispersion. Figures 4b and 5b show that the increase of the effective dispersion coefficients occurs for high Pe at earlier non-dimensional times than for low Pe. This indicates that the shear rate does not scale linearly with $\langle u \rangle$. In fact, a typical shear rate can be written as

$$ \gamma = \frac{\langle v \rangle}{\ell_{\gamma}}, \qquad (24) $$

where $\ell_{\gamma}$ is the scale of transverse velocity contrast. The latter is proportional to the typical streamtube size. That is, as $\ell_{\gamma}^2 \langle v \rangle = \text{constant}$, we have $\ell_{\gamma} \sim \langle v \rangle^{-1/2}$. The characteristic shear length scale decreases with increasing flow rate, and thus the shear rate scales as $\gamma \sim \langle u \rangle^{3/2}$. The characteristic shear time scale is then $\tau_{\gamma} = \gamma^{-1} \propto \tau_v / \langle v \rangle^{1/2}$. This dependence explains the differences in the time behaviors of the effective dispersion coefficients for different Pe.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Dispersion coefficients for the Berea sandstone sample. Top panels: (a) Ensemble dispersion coefficients for (red solid line) Pe = 2,000 and (orange solid line) Pe = 200, and (b) corresponding effective dispersion coefficients. The vertical dashed lines denote the corresponding diffusion time scale $\tau_D = \tau_v Pe$. Bottom panels: (Black solid lines) Ensemble and (blue solid lines) effective dispersion coefficients for (c) Pe = 200 and (d) Pe = 2,000. The vertical black dashed lines denote the decorrelation time scale $\tau_c = \tau_v$, the blue dashed lines denote the respective diffusion time scales. The horizontal dash-dotted lines denote the asymptotic short time and long time values.</figcaption></figure>

The early time ballistic and shear dispersion behaviors for $t < \tau_v$ are observed for both the sand pack and Berea samples. For $t > \tau_v$ the dispersion behaviors are different as discussed in the following.

### 3.2.1. Sand Pack Sample

Figures 4a–4d show the evolution of the ensemble and effective dispersion coefficients for the sand pack sample. For times $t > \tau_v$, that is for mean travel distances larger than the average pore size, particles start sampling different flow velocities along their trajectories, and the ballistic behavior for the ensemble dispersion coefficients breaks down, see Figure 4a.

For purely advective transport, the ensemble dispersion coefficients continue growing non-linearly with time, which can be traced back to the broad distribution of transition time across pores (Puyguiraud et al., 2019). At finite Pe, the ensemble dispersion coefficients first follow the purely advective behavior and eventually cross over toward their asymptotic value on the time scale. The effective dispersion coefficients shown in Figure 4 cross over toward their asymptotic values, also on the time scale $\tau_v$. As shown in Figures 4c and 4d, they converge with $D_L^{\text{ens}}(t)$.

As mentioned in Section 3.1, these behaviors are at first sight counter-intuitive because we expect the deviation from the purely advective behavior observed for $D_L^{\text{ens}}(t)$ and the convergence of $D_L^{\text{eff}}(t)$ toward $D_L^{\text{ens}}(t)$ to be governed by diffusion. For ensemble dispersion, diffusion is the mechanism that decorrelates subsequent (low
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Dispersion coefficients for the Berea sandstone sample. Top panels: (a) Ensemble dispersion coefficients for (red solid line) Pe = 2,000 and (orange solid line) Pe = 200, and (b) corresponding effective dispersion coefficients. The vertical dashed lines denote the corresponding diffusion time scale $\tau_D = \tau_v Pe$. Bottom panels: (Black solid lines) Ensemble and (blue solid lines) effective dispersion coefficients for (c) Pe = 200 and (d) Pe = 2,000. The vertical black dashed lines denote the decorrelation time scale $\tau_c = \tau_v$, the blue dashed lines denote the respective diffusion time scales. The horizontal dash-dotted lines denote the asymptotic short time and long time values.</figcaption></figure>

The early time ballistic and shear dispersion behaviors for $t < \tau_v$ are observed for both the sand pack and Berea samples. For $t > \tau_v$ the dispersion behaviors are different as discussed in the following.

### 3.2.1. Sand Pack Sample

Figures 4a–4d show the evolution of the ensemble and effective dispersion coefficients for the sand pack sample. For times $t > \tau_v$, that is for mean travel distances larger than the average pore size, particles start sampling different flow velocities along their trajectories, and the ballistic behavior for the ensemble dispersion coefficients breaks down, see Figure 4a.

For purely advective transport, the ensemble dispersion coefficients continue growing non-linearly with time, which can be traced back to the broad distribution of transition time across pores (Puyguiraud et al., 2019). At finite Pe, the ensemble dispersion coefficients first follow the purely advective behavior and eventually cross over toward their asymptotic value on the time scale. The effective dispersion coefficients shown in Figure 4 cross over toward their asymptotic values, also on the time scale $\tau_v$. As shown in Figures 4c and 4d, they converge with $D_L^{\text{ens}}(t)$.

As mentioned in Section 3.1, these behaviors are at first sight counter-intuitive because we expect the deviation from the purely advective behavior observed for $D_L^{\text{ens}}(t)$ and the convergence of $D_L^{\text{eff}}(t)$ toward $D_L^{\text{ens}}(t)$ to be governed by diffusion. For ensemble dispersion, diffusion is the mechanism that decorrelates subsequent (low
velocities in time and thus leads to the separation of $D_L^{\text{ens}}(t)$ from the (anomalous) purely advective behavior. Similarly, the mechanism by which the effective dispersion coefficients converge toward the ensemble dispersion coefficients is due to decorrelation of the velocities of the particles that start from the same point, which is due to diffusion in transverse direction. Thus one would expect that the dispersion coefficients evolve on the diffusion time scale $\tau_D$.

As discussed in Section 3.1.1, the decorrelation mechanism is indeed transverse diffusion across a length scale that is related to a typical streamtube width. Thus, the decorrelation time $\tau_c$ is given by Equation 22, which is proportional to $\tau_v$. This observation explains the temporal evolution of the ensemble and effective dispersion coefficients on the time scale $\tau_v$.

Sole-Mari et al. (2022) study the longitudinal ensemble dispersion coefficients for transport in granular porous media. These authors find that the evolution time-scale depends on $Pe$ and varies between the diffusion time scale for $Pe$ smaller than a critical $Pe_c$ and the advection time scale for $Pe \gg Pe_c$. The estimate of these authors is based on the heuristic expression

$$D_L^{\text{ens}}(t) = D + \delta D_L^{\infty} \left[ 1 - \exp\left(-\sqrt{t/\beta\tau_v}\right) \right], \quad (25)$$

which we slightly modified such that $D_L^{\text{ens}}(t=0) = D$ is equal to the molecular diffusion coefficient. The asymptotic mechanical dispersion coefficient is denoted by $\delta D_L^{\infty}$ and $\beta$ is a fitting parameter, which Sole-Mari et al. (2022) relate to $Pe$. Fits are performed using non-linear least squares. We find $\beta \approx 16$ for $Pe = 2,000$ and $\beta \approx 13$ for $Pe = 200$. As shown in Figure 4a, this expression can describe the crossover behavior from ballistic to asymptotic, but is not able to reproduce the ballistic behavior observed at times $t \ll \tau_v$ because it behaves as $D_L^{\text{ens}}(t) \propto \sqrt{t}$ in this time regime. Here, we use the following heuristic expression

$$D_L^{\text{ens}}(t) = D + \frac{\sigma_0^2}{\langle v \rangle} \ell_0(v) \alpha [1 - \exp(-t/\alpha\tau_v)], \quad (26)$$

to fit the data for the longitudinal ensemble dispersion coefficients with $\alpha$ a fitting parameter. At short times $t \ll \tau_v$, Equation 26 shows the correct ballistic behavior as $D_L^{\text{ens}}(t) = D + \sigma_0^2 t$. Figure 4a shows the best fits with $\alpha = 3.1$ for $Pe = 200$ and $\alpha = 3.9$ for $Pe = 2,000$. Note however, that unlike Equation 25, Equation 26 describes too fast a crossover from the ballistic to the asymptotic behavior. The effective dispersion coefficients show a more complex behavior and cannot be described by these simple regressions.

### 3.2.2. Berea Sandstone Sample

Figures 5a–5d show the evolution of the ensemble and effective dispersion coefficients for the Berea sandstone sample. As seen in Figure 5a, the initial ballistic behavior for the ensemble dispersion coefficients breaks down on the time scale $\tau_v$ when particles start sampling different flow velocities along their trajectories. For purely advective transport, we observe anomalous dispersion characterized by a super-linear growth of the ensemble dispersion coefficients, which can be traced back to broad distributions of advective particle transition times (Puyguiraud et al., 2019). Unlike for the sand pack, here the crossover toward the constant asymptotic long time values occurs on the diffusion time scale $\tau_D$. As discussed in Section 3.1.2, the temporal decorrelation of low velocities is due to diffusion along pore channels with the characteristic time scale $\tau_D$ (Puyguiraud et al., 2021). Also the crossover of the effective dispersion coefficient toward the asymptotic value shown in Figure 5b occurs on the time scale $\tau_D$.

The convergence of the effective and ensemble dispersion coefficients shown in Figures 5c and 5d on the other hand, occurs on the decorrelation time scale $\tau_c$, see Equation 22. This time scale is set by transverse diffusion across streamtubes, which is the mechanism by which particles that originate from the same initial position start decorrelating and sampling different flow velocities. Thus, for $t > \tau_c$ particle trajectories become independent, particles “forget” that they originate from the same initial condition, and dispersion is enhanced due to the increasing difference in travel distance between different particles. This is the same mechanism that causes longitudinal ensemble dispersion, and therefore the two dispersion coefficients converge. These complex
velocities in time and thus leads to the separation of $D_L^{\text{ens}}(t)$ from the (anomalous) purely advective behavior. Similarly, the mechanism by which the effective dispersion coefficients converge toward the ensemble dispersion coefficients is due to decorrelation of the velocities of the particles that start from the same point, which is due to diffusion in transverse direction. Thus one would expect that the dispersion coefficients evolve on the diffusion time scale $\tau_D$.

As discussed in Section 3.1.1, the decorrelation mechanism is indeed transverse diffusion across a length scale that is related to a typical streamtube width. Thus, the decorrelation time $\tau_c$ is given by Equation 22, which is proportional to $\tau_v$. This observation explains the temporal evolution of the ensemble and effective dispersion coefficients on the time scale $\tau_v$.

Sole-Mari et al. (2022) study the longitudinal ensemble dispersion coefficients for transport in granular porous media. These authors find that the evolution time-scale depends on $Pe$ and varies between the diffusion time scale for $Pe$ smaller than a critical $Pe_c$ and the advection time scale for $Pe \gg Pe_c$. The estimate of these authors is based on the heuristic expression

$$D_L^{\text{ens}}(t) = D + \delta D_L^{\infty} \left[ 1 - \exp\left(-\sqrt{t/\beta\tau_v}\right) \right], \quad (25)$$

which we slightly modified such that $D_L^{\text{ens}}(t=0) = D$ is equal to the molecular diffusion coefficient. The asymptotic mechanical dispersion coefficient is denoted by $\delta D_L^{\infty}$ and $\beta$ is a fitting parameter, which Sole-Mari et al. (2022) relate to $Pe$. Fits are performed using non-linear least squares. We find $\beta \approx 16$ for $Pe = 2,000$ and $\beta \approx 13$ for $Pe = 200$. As shown in Figure 4a, this expression can describe the crossover behavior from ballistic to asymptotic, but is not able to reproduce the ballistic behavior observed at times $t \ll \tau_v$ because it behaves as $D_L^{\text{ens}}(t) \propto \sqrt{t}$ in this time regime. Here, we use the following heuristic expression

$$D_L^{\text{ens}}(t) = D + \frac{\sigma_0^2}{\langle v \rangle} \ell_0(v) \alpha [1 - \exp(-t/\alpha\tau_v)], \quad (26)$$

to fit the data for the longitudinal ensemble dispersion coefficients with $\alpha$ a fitting parameter. At short times $t \ll \tau_v$, Equation 26 shows the correct ballistic behavior as $D_L^{\text{ens}}(t) = D + \sigma_0^2 t$. Figure 4a shows the best fits with $\alpha = 3.1$ for $Pe = 200$ and $\alpha = 3.9$ for $Pe = 2,000$. Note however, that unlike Equation 25, Equation 26 describes too fast a crossover from the ballistic to the asymptotic behavior. The effective dispersion coefficients show a more complex behavior and cannot be described by these simple regressions.

### 3.2.2. Berea Sandstone Sample

Figures 5a–5d show the evolution of the ensemble and effective dispersion coefficients for the Berea sandstone sample. As seen in Figure 5a, the initial ballistic behavior for the ensemble dispersion coefficients breaks down on the time scale $\tau_v$ when particles start sampling different flow velocities along their trajectories. For purely advective transport, we observe anomalous dispersion characterized by a super-linear growth of the ensemble dispersion coefficients, which can be traced back to broad distributions of advective particle transition times (Puyguiraud et al., 2019). Unlike for the sand pack, here the crossover toward the constant asymptotic long time values occurs on the diffusion time scale $\tau_D$. As discussed in Section 3.1.2, the temporal decorrelation of low velocities is due to diffusion along pore channels with the characteristic time scale $\tau_D$ (Puyguiraud et al., 2021). Also the crossover of the effective dispersion coefficient toward the asymptotic value shown in Figure 5b occurs on the time scale $\tau_D$.

The convergence of the effective and ensemble dispersion coefficients shown in Figures 5c and 5d on the other hand, occurs on the decorrelation time scale $\tau_c$, see Equation 22. This time scale is set by transverse diffusion across streamtubes, which is the mechanism by which particles that originate from the same initial position start decorrelating and sampling different flow velocities. Thus, for $t > \tau_c$ particle trajectories become independent, particles “forget” that they originate from the same initial condition, and dispersion is enhanced due to the increasing difference in travel distance between different particles. This is the same mechanism that causes longitudinal ensemble dispersion, and therefore the two dispersion coefficients converge. These complex
behaviors cannot be captured by simple heuristic expressions as Equations 25 and 26, but require a more in-depth
theoretical analysis, which is beyond the scope of this paper.

4. Conclusions

We investigate solute dispersion in three-dimensional porous rocks using detailed numerical simulations of pore-scale flow and transport. We consider a sand-like medium, and a Berea sandstone sample. The two media have quite distinct pore structures, which manifests in distinct pore-scale flow variability. The latter is quantified by the distribution of Eulerian flow speeds. The degree of flow heterogeneity is measured by the variance of the logarithm of the flow speed, which is significantly higher for the Berea sample than for the sand pack sample.

We use longitudinal effective and ensemble dispersion coefficients to probe pore-scale mixing and spreading behaviors and shed light on the impact of medium structure and diffusive mass transfer on solute and particle transport. Effective dispersion coefficients are defined in terms of the spatial average of the second-centered moments of the partial plumes (Green functions) that constitute the global solute distribution. Thus, they can be seen as a measure for the typical width of a mixing front. Ensemble dispersion coefficients on the other hand are defined in terms of the second centered moments of the global solute plume and as such are a measure for mixing front deformation due to the flow variability within the initial plume. The fundamental mechanisms that cause hydrodynamic dispersion are pore-scale flow variability and molecular diffusion, which govern the evolution of both the effective and ensemble dispersion coefficients. Their detailed time behaviors and the time scale on which they evolve provide information on the salient small scale mass transfer and heterogeneity mechanisms.

The early time behavior of the ensemble coefficient is ballistic as a result of the spatial persistence of flow velocities in the initial plume. The effective coefficients on the other hand are significantly smaller than their ensemble counterparts. Their early time evolution is dominated by shear dispersion, which results from the velocity gradients within the partial plumes, whose lateral extent initially increases by diffusion. The two dispersion coefficients start converging when the lateral extent of the partial plumes is large enough for the efficient sampling of the flow heterogeneity, and it is here, where dispersion in the sand pack and Berea sandstone behaves differently. For the sand pack, the evolution of effective dispersion is marked by the characteristic diffusion time across a streamtube, which sets the time for both convergence to ensemble dispersion and its asymptotic behavior. For the Berea sandstone, this time scale marks the time for convergence of effective and ensemble dispersion, which, however, still evolve non-linearly with time until they retrieve their asymptotic long time value on the time scale for diffusion over a typical pore length. The evolution of solute dispersion reflects the medium structure, which determines the microscopic mass transfer mechanisms. In fact, we hypothesize that the observed behaviors originate in the network-like medium structure in case of the Berea sample, and the strong connectivity of pores in the sand pack. While the behavior of ensemble dispersion can be captured by travel-time based approaches like the continuous time random walk in terms of flow variability and medium structure, it is still elusive how to quantify effective dispersion in these terms. Nevertheless, these findings give insight into how particles sample pore-scale velocity variability and how it depends on the medium structure. These are key elements for the upscaling of transport, mixing and reaction from the pore to the continuum scale.

Data Availability Statement

The simulation data displayed in the figures are deposited at the DIGITAL.CSIC open repository and are available
in Puyguiraud et al. (2023).

Acknowledgments

This project has received funding from the
European Union's Horizon 2020 research
and innovation programme under the
Marie Skłodowska-Curie grant agreement
No 899546. MD gratefully acknowledges
the support of the Spanish Ministry of
Science and Innovation through the project
HydroPore (PID2019-106887GB-C31).

References

Ahmadi, A., Quintard, M., & Whitaker, S. (1998). Transport in chemically and mechanically heterogeneous porous media. *Advances in Water Resources*, 22(1), 59–86. https://doi.org/10.1016/s0309-1708(97)00032-8

Alhashmi, Z., Blunt, M., & Bijeljic, B. (2016). The impact of pore structure heterogeneity, transport, and reaction conditions on fluid–fluid reaction rate studied on images of pore space. *Transport in Porous Media*, 115(2), 215–237. https://doi.org/10.1007/s11242-016-0758-z

Alim, K., Parsa, S., Weitz, D. A., & Brenner, M. P. (2017). Local pore size correlations determine flow distributions in porous media. *Physical Review Letters*, 119(14), 144501. https://doi.org/10.1103/PhysRevLett.119.144501

Attinger, S., Dentz, M., Kinzelbach, H., & Kinzelbach, W. (1999). Temporal behavior of a solute cloud in a chemically heterogeneous porous medium. *Journal of Fluid Mechanics*, 386, 77–104. https://doi.org/10.1017/s0022112099004334

Bijeljic, B., & Blunt, M. J. (2006). Pore-scale modeling and continuous time random walk analysis of dispersion in porous media. *Water Resources Research*, 42(1), W01202. https://doi.org/10.1029/2005WR004578
behaviors cannot be captured by simple heuristic expressions as Equations 25 and 26, but require a more in-depth
theoretical analysis, which is beyond the scope of this paper.

4. Conclusions

We investigate solute dispersion in three-dimensional porous rocks using detailed numerical simulations of pore-scale flow and transport. We consider a sand-like medium, and a Berea sandstone sample. The two media have quite distinct pore structures, which manifests in distinct pore-scale flow variability. The latter is quantified by the distribution of Eulerian flow speeds. The degree of flow heterogeneity is measured by the variance of the logarithm of the flow speed, which is significantly higher for the Berea sample than for the sand pack sample.

We use longitudinal effective and ensemble dispersion coefficients to probe pore-scale mixing and spreading behaviors and shed light on the impact of medium structure and diffusive mass transfer on solute and particle transport. Effective dispersion coefficients are defined in terms of the spatial average of the second-centered moments of the partial plumes (Green functions) that constitute the global solute distribution. Thus, they can be seen as a measure for the typical width of a mixing front. Ensemble dispersion coefficients on the other hand are defined in terms of the second centered moments of the global solute plume and as such are a measure for mixing front deformation due to the flow variability within the initial plume. The fundamental mechanisms that cause hydrodynamic dispersion are pore-scale flow variability and molecular diffusion, which govern the evolution of both the effective and ensemble dispersion coefficients. Their detailed time behaviors and the time scale on which they evolve provide information on the salient small scale mass transfer and heterogeneity mechanisms.

The early time behavior of the ensemble coefficient is ballistic as a result of the spatial persistence of flow velocities in the initial plume. The effective coefficients on the other hand are significantly smaller than their ensemble counterparts. Their early time evolution is dominated by shear dispersion, which results from the velocity gradients within the partial plumes, whose lateral extent initially increases by diffusion. The two dispersion coefficients start converging when the lateral extent of the partial plumes is large enough for the efficient sampling of the flow heterogeneity, and it is here, where dispersion in the sand pack and Berea sandstone behaves differently. For the sand pack, the evolution of effective dispersion is marked by the characteristic diffusion time across a streamtube, which sets the time for both convergence to ensemble dispersion and its asymptotic behavior. For the Berea sandstone, this time scale marks the time for convergence of effective and ensemble dispersion, which, however, still evolve non-linearly with time until they retrieve their asymptotic long time value on the time scale for diffusion over a typical pore length. The evolution of solute dispersion reflects the medium structure, which determines the microscopic mass transfer mechanisms. In fact, we hypothesize that the observed behaviors originate in the network-like medium structure in case of the Berea sample, and the strong connectivity of pores in the sand pack. While the behavior of ensemble dispersion can be captured by travel-time based approaches like the continuous time random walk in terms of flow variability and medium structure, it is still elusive how to quantify effective dispersion in these terms. Nevertheless, these findings give insight into how particles sample pore-scale velocity variability and how it depends on the medium structure. These are key elements for the upscaling of transport, mixing and reaction from the pore to the continuum scale.

Data Availability Statement

The simulation data displayed in the figures are deposited at the DIGITAL.CSIC open repository and are available
in Puyguiraud et al. (2023).

Acknowledgments

This project has received funding from the
European Union's Horizon 2020 research
and innovation programme under the
Marie Skłodowska-Curie grant agreement
No 899546. MD gratefully acknowledges
the support of the Spanish Ministry of
Science and Innovation through the project
HydroPore (PID2019-106887GB-C31).

References

Ahmadi, A., Quintard, M., & Whitaker, S. (1998). Transport in chemically and mechanically heterogeneous porous media. *Advances in Water Resources*, 22(1), 59–86. https://doi.org/10.1016/s0309-1708(97)00032-8

Alhashmi, Z., Blunt, M., & Bijeljic, B. (2016). The impact of pore structure heterogeneity, transport, and reaction conditions on fluid–fluid reaction rate studied on images of pore space. *Transport in Porous Media*, 115(2), 215–237. https://doi.org/10.1007/s11242-016-0758-z

Alim, K., Parsa, S., Weitz, D. A., & Brenner, M. P. (2017). Local pore size correlations determine flow distributions in porous media. *Physical Review Letters*, 119(14), 144501. https://doi.org/10.1103/PhysRevLett.119.144501

Attinger, S., Dentz, M., Kinzelbach, H., & Kinzelbach, W. (1999). Temporal behavior of a solute cloud in a chemically heterogeneous porous medium. *Journal of Fluid Mechanics*, 386, 77–104. https://doi.org/10.1017/s0022112099004334

Bijeljic, B., & Blunt, M. J. (2006). Pore-scale modeling and continuous time random walk analysis of dispersion in porous media. *Water Resources Research*, 42(1), W01202. https://doi.org/10.1029/2005WR004578
Bijeljic, B., Mostaghimi, P., & Blunt, M. J. (2011). Signature of non-Fickian solute transport in complex heterogeneous porous media. *Physical Review Letters, 107*(20), 204502. https://doi.org/10.1103/physrevlett.107.204502

Bijeljic, B., Muggeridge, A. H., & Blunt, M. J. (2004). Pore-scale modeling of longitudinal dispersion. *Water Resources Research, 40*(11), W11501. https://doi.org/10.1029/2004wr003567

Bouchaud, J.-P., & Georges, A. (1990). Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications. *Physics Reports, 195*(4-5), 127-293. https://doi.org/10.1016/0370-1573(90)90099-n

Brenner, H. (1980). Dispersion resulting from flow through spatially periodic porous media. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 297*, 81-133.

Churcher, P., French, P., Shaw, J., & Schramm, L. (1991). Rock properties of Berea sandstone, Baker dolomite, and Indiana limestone. In *SPE International Symposium on Oilfield Chemistry*.

Dagan, G. (1990). Transport in heterogeneous porous formations: Spatial moments, ergodicity, and effective dispersion. *Water Resources Research, 26*(6), 1281-1290. https://doi.org/10.1029/wr026i006p01281

Datta, S. S., Chiang, H., Ramakrishnan, T. S., & Weitz, D. A. (2013). Spatial fluctuations of fluid velocities in flow through a three-dimensional porous medium. *Physical Review Letters, 111*(6), 064501. https://doi.org/10.1103/PhysRevLett.111.064501

de Barros, F. P. J., & Dentz, M. (2016). Pictures of blockscale transport: Effective versus ensemble dispersion and its uncertainty. *Advances in Water Resources, 91*, 11-22. https://doi.org/10.1016/j.advwatres.2016.03.004

De Barros, F. P. J., Fiori, A., Boso, F., & Bellin, A. (2015). A theoretical framework for modeling dilution enhancement of non-reactive solutes in heterogeneous porous media. *Journal of Contaminant Hydrology, 175*, 72-83. https://doi.org/10.1016/j.jconhyd.2015.01.004

Dentz, M., Cortis, A., Scher, H., & Berkowitz, B. (2004). Time behavior of solute transport in heterogeneous media: Transition from anomalous to normal transport. *Advances in Water Resources, 27*(2), 155-173. https://doi.org/10.1016/j.advwatres.2003.11.002

Dentz, M., & de Barros, F. (2015). Mixing-scale dependent dispersion for transport in heterogeneous flows. *Journal of Fluid Mechanics, 777*, 178-195. https://doi.org/10.1017/jfm.2015.351

Dentz, M., Hidalgo, J. J., & Lester, D. (2023). Mixing in porous media: Concepts and approaches across scales. *Transport in Porous Media, 146*(1-2), 5-53. https://doi.org/10.1007/s11242-022-01852-x

Dentz, M., Kinzelbach, H., Attinger, S., & Kinzelbach, W. (2000). Temporal behavior of a solute cloud in a heterogeneous porous medium: 1. Point-like injection. *Water Resources Research, 36*(12), 3591-3604. https://doi.org/10.1029/2000WR900162

Di Palma, P. R., Guyennon, N., Parmigiani, A., Huber, C., Heße, F., & Romano, E. (2019). Impact of synthetic porous medium geometric properties on solute transport using direct 3D pore-scale simulations. *Geofluids, 2019*, 1-13. https://doi.org/10.1155/2019/6810467

Fiori, A. (2001). On the influence of local dispersion in solute transport through formations with evolving scales of heterogeneity. *Water Resources Research, 37*(2), 235-242. https://doi.org/10.1029/2000wr900245

Fiori, A., & Dagan, G. (2000). Concentration fluctuations in aquifer transport: A rigorous first-order solution and applications. *Journal of Contaminant Hydrology, 45*(1-2), 139-163. https://doi.org/10.1016/s0169-7722(00)00123-6

Gelhar, L. W., & Axness, C. L. (1983). Three-dimensional stochastic analysis of macrodispersion in aquifers. *Water Resources Research, 19*(1), 161-180. https://doi.org/10.1029/wr019i001p00161

Gouze, P., Puyguiraud, A., Porcher, T., & Dentz, M. (2021). Modeling longitudinal dispersion in variable porosity porous media: Control of velocity distribution and microstructures. *Frontiers in Water, 3*, 766338. https://doi.org/10.3389/frwa.2021.766338

Gouze, P., Puyguiraud, A., Roubinet, D., & Dentz, M. (2021). Pore-scale transport in rocks of different complexity modeled by random walk methods. *Transport in Porous Media, 146*(1-2), 139-158. https://doi.org/10.1007/s11242-021-01675-2

Hulin, J. P., & Plona, T. J. (1989). "Echo" tracer dispersion in porous media. *Physics of Fluids A: Fluid Dynamics, 1*(8), 1341-1347. https://doi.org/10.1063/1.857309

Imperial College Consortium on Pore-scale Imaging and Modelling. (2014). LV60C sandpack. (Technical Report). https://doi.org/10.6084/m9.fishshare.1189272.v1

Khrapitchev, A. A., & Callaghan, P. T. (2003). Reversible and irreversible dispersion in a porous medium. *Physics of Fluids, 15*(9), 2649-2660. https://doi.org/10.1063/1.1596914

Kitanidis, P. K. (1988). Prediction by the method of moments of transport in a heterogeneous formation. *Journal of Hydrology, 102*(1-4), 453-473. https://doi.org/10.1016/0022-1694(88)90111-4

Koch, D. L., & Brady, J. F. (1985). Dispersion in fixed beds. *Journal of Fluid Mechanics, 154*, 399-427. https://doi.org/10.1017/s0022112085001598

Leal, L. G. (2007). *Advanced transport phenomena: Fluid mechanics and convective transport processes*. Cambridge University Press. https://doi.org/10.1017/CBO9780511800245

Le Borgne, T., Bolster, D., Dentz, M., Anna, P., & Tartakovsky, A. (2011). Effective pore-scale dispersion upscaling with a correlated continuous time random walk approach. *Water Resources Research, 47*(12), W12538. https://doi.org/10.1029/2011WR010457

Mostaghimi, P., Bijeljic, B., & Blunt, M. J. (2012). Simulation of flow and dispersion on pore-space images. *SPE Journal, 17*(04), 1131-1141. https://doi.org/10.2118/135261-PA

Noetinger, B., Roubinet, D., Russian, A., Le Borgne, T., Delay, F., Dentz, M., et al. (2016). Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale. *Transport in Porous Media, 115*(2), 1-41. https://doi.org/10.1007/s11242-016-0693-z

Puyguiraud, A., Gouze, P., & Dentz, M. (2019). Upscaling of anomalous pore-scale dispersion. *Transport in Porous Media, 128*(2), 837-855. https://doi.org/10.1007/s11242-019-01273-3

Puyguiraud, A., Gouze, P., & Dentz, M. (2021). Pore-scale mixing and the evolution of hydrodynamic dispersion in porous media. *Physical Review Letters, 126*(16), 164501. https://doi.org/10.1103/PhysRevLett.126.164501

Puyguiraud, A., Gouze, P., & Dentz, M. (2023). Temporal evolution of solute dispersion in three-dimensional porous rocks [Dataset]. *DIGITAL CSIC*. https://doi.org/10.20350/digitalCSIC/15442

Risken, H. (1996). *The Fokker-Planck equation*. Springer Heidelberg.

Saffman, P. (1959). A theory of dispersion in a porous medium. *Journal of Fluid Mechanics, 6*(03), 321-349. https://doi.org/10.1017/s0022112059000672

Scheven, U. (2013). Pore-scale mixing and transverse dispersivity of randomly packed monodisperse spheres. *Physical Review Letters, 110*(21), 214504. https://doi.org/10.1103/physrevlett.110.214504

Sherman, T., Engdahl, N. B., Porta, G., & Bolster, D. (2021). A review of spatial Markov models for predicting pre-asymptotic and anomalous transport in porous and fractured media. *Journal of Contaminant Hydrology, 236*, 103734. https://doi.org/10.1016/j.jconhyd.2020.103734

Sole-Mari, G., Bolster, D., & Fernandez-Garcia, D. (2022). A closer look: High-resolution pore-scale simulations of solute transport and mixing through porous media columns. *Transport in Porous Media, 146*(1-2), 1-27. https://doi.org/10.1007/s11242-021-01721-z
Bijeljic, B., Mostaghimi, P., & Blunt, M. J. (2011). Signature of non-Fickian solute transport in complex heterogeneous porous media. *Physical Review Letters*, **107**(20), 204502. https://doi.org/10.1103/physrevlett.107.204502

Bijeljic, B., Muggeridge, A. H., & Blunt, M. J. (2004). Pore-scale modeling of longitudinal dispersion. *Water Resources Research*, **40**(11), W11501. https://doi.org/10.1029/2004wr003567

Bouchaud, J.-P., & Georges, A. (1990). Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications. *Physics Reports*, **195**(4-5), 127-293. https://doi.org/10.1016/0370-1573(90)90099-n

Brenner, H. (1980). Dispersion resulting from flow through spatially periodic porous media. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences*, **297**, 81-133.

Churcher, P., French, P., Shaw, J., & Schramm, L. (1991). Rock properties of Berea sandstone, Baker dolomite, and Indiana limestone. In *SPE International Symposium on Oilfield Chemistry*.

Dagan, G. (1990). Transport in heterogeneous porous formations: Spatial moments, ergodicity, and effective dispersion. *Water Resources Research*, **26**(6), 1281-1290. https://doi.org/10.1029/wr026i006p01281

Datta, S. S., Chiang, H., Ramakrishnan, T. S., & Weitz, D. A. (2013). Spatial fluctuations of fluid velocities in flow through a three-dimensional porous medium. *Physical Review Letters*, **111**(6), 064501. https://doi.org/10.1103/PhysRevLett.111.064501

de Barros, F. P. J., & Dentz, M. (2016). Pictures of blockscale transport: Effective versus ensemble dispersion and its uncertainty. *Advances in Water Resources*, **91**, 11-22. https://doi.org/10.1016/j.advwatres.2016.03.004

De Barros, F. P. J., Fiori, A., Boso, F., & Bellin, A. (2015). A theoretical framework for modeling dilution enhancement of non-reactive solutes in heterogeneous porous media. *Journal of Contaminant Hydrology*, **175**, 72-83. https://doi.org/10.1016/j.jconhyd.2015.01.004

Dentz, M., Cortis, A., Scher, H., & Berkowitz, B. (2004). Time behavior of solute transport in heterogeneous media: Transition from anomalous to normal transport. *Advances in Water Resources*, **27**(2), 155-173. https://doi.org/10.1016/j.advwatres.2003.11.002

Dentz, M., & de Barros, F. (2015). Mixing-scale dependent dispersion for transport in heterogeneous flows. *Journal of Fluid Mechanics*, **777**, 178-195. https://doi.org/10.1017/jfm.2015.351

Dentz, M., Hidalgo, J. J., & Lester, D. (2023). Mixing in porous media: Concepts and approaches across scales. *Transport in Porous Media*, **146**(1-2), 5-53. https://doi.org/10.1007/s11242-022-01852-x

Dentz, M., Kinzelbach, H., Attinger, S., & Kinzelbach, W. (2000). Temporal behavior of a solute cloud in a heterogeneous porous medium: 1. Point-like injection. *Water Resources Research*, **36**(12), 3591-3604. https://doi.org/10.1029/2000WR900162

Di Palma, P. R., Guyennon, N., Parmigiani, A., Huber, C., Heße, F., & Romano, E. (2019). Impact of synthetic porous medium geometric properties on solute transport using direct 3D pore-scale simulations. *Geofluids*, **2019**, 1-13. https://doi.org/10.1155/2019/6810467

Fiori, A. (2001). On the influence of local dispersion in solute transport through formations with evolving scales of heterogeneity. *Water Resources Research*, **37**(2), 235-242. https://doi.org/10.1029/2000wr900245

Fiori, A., & Dagan, G. (2000). Concentration fluctuations in aquifer transport: A rigorous first-order solution and applications. *Journal of Contaminant Hydrology*, **45**(1-2), 139-163. https://doi.org/10.1016/s0169-7722(00)00123-6

Gelhar, L. W., & Axness, C. L. (1983). Three-dimensional stochastic analysis of macrodispersion in aquifers. *Water Resources Research*, **19**(1), 161-180. https://doi.org/10.1029/wr019i001p00161

Gouze, P., Puyguiraud, A., Porcher, T., & Dentz, M. (2021). Modeling longitudinal dispersion in variable porosity porous media: Control of velocity distribution and microstructures. *Frontiers in Water*, **3**, 766338. https://doi.org/10.3389/frwa.2021.766338

Gouze, P., Puyguiraud, A., Roubinet, D., & Dentz, M. (2021). Pore-scale transport in rocks of different complexity modeled by random walk methods. *Transport in Porous Media*, **146**(1-2), 139-158. https://doi.org/10.1007/s11242-021-01675-2

Hulin, J. P., & Plona, T. J. (1989). "Echo" tracer dispersion in porous media. *Physics of Fluids A: Fluid Dynamics*, **1**(8), 1341-1347. https://doi.org/10.1063/1.857309

Imperial College Consortium on Pore-scale Imaging and Modelling. (2014). LV60C sandpack. (Technical Report). https://doi.org/10.6084/m9.fishshare.1189272.v1

Khrapitchev, A. A., & Callaghan, P. T. (2003). Reversible and irreversible dispersion in a porous medium. *Physics of Fluids*, **15**(9), 2649-2660. https://doi.org/10.1063/1.1596914

Kitanidis, P. K. (1988). Prediction by the method of moments of transport in a heterogeneous formation. *Journal of Hydrology*, **102**(1-4), 453-473. https://doi.org/10.1016/0022-1694(88)90111-4

Koch, D. L., & Brady, J. F. (1985). Dispersion in fixed beds. *Journal of Fluid Mechanics*, **154**, 399-427. https://doi.org/10.1017/s0022112085001598

Leal, L. G. (2007). *Advanced transport phenomena: Fluid mechanics and convective transport processes*. Cambridge University Press. https://doi.org/10.1017/CBO9780511800245

Le Borgne, T., Bolster, D., Dentz, M., Anna, P., & Tartakovsky, A. (2011). Effective pore-scale dispersion upscaling with a correlated continuous time random walk approach. *Water Resources Research*, **47**(12), W12538. https://doi.org/10.1029/2011WR010457

Mostaghimi, P., Bijeljic, B., & Blunt, M. J. (2012). Simulation of flow and dispersion on pore-space images. *SPE Journal*, **17**(04), 1131-1141. https://doi.org/10.2118/135261-PA

Noetinger, B., Roubinet, D., Russian, A., Le Borgne, T., Delay, F., Dentz, M., et al. (2016). Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale. *Transport in Porous Media*, **115**(2), 1-41. https://doi.org/10.1007/s11242-016-0693-z

Puyguiraud, A., Gouze, P., & Dentz, M. (2019). Upscaling of anomalous pore-scale dispersion. *Transport in Porous Media*, **128**(2), 837-855. https://doi.org/10.1007/s11242-019-01273-3

Puyguiraud, A., Gouze, P., & Dentz, M. (2021). Pore-scale mixing and the evolution of hydrodynamic dispersion in porous media. *Physical Review Letters*, **126**(16), 164501. https://doi.org/10.1103/PhysRevLett.126.164501

Puyguiraud, A., Gouze, P., & Dentz, M. (2023). Temporal evolution of solute dispersion in three-dimensional porous rocks [Dataset]. *DIGITAL CSIC*. https://doi.org/10.20350/digitalCSIC/15442

Risken, H. (1996). *The Fokker-Planck equation*. Springer Heidelberg.

Saffman, P. (1959). A theory of dispersion in a porous medium. *Journal of Fluid Mechanics*, **6**(03), 321-349. https://doi.org/10.1017/s0022112059000672

Scheven, U. (2013). Pore-scale mixing and transverse dispersivity of randomly packed monodisperse spheres. *Physical Review Letters*, **110**(21), 214504. https://doi.org/10.1103/physrevlett.110.214504

Sherman, T., Engdahl, N. B., Porta, G., & Bolster, D. (2021). A review of spatial Markov models for predicting pre-asymptotic and anomalous transport in porous and fractured media. *Journal of Contaminant Hydrology*, **236**, 103734. https://doi.org/10.1016/j.jconhyd.2020.103734

Sole-Mari, G., Bolster, D., & Fernandez-Garcia, D. (2022). A closer look: High-resolution pore-scale simulations of solute transport and mixing through porous media columns. *Transport in Porous Media*, **146**(1-2), 1-27. https://doi.org/10.1007/s11242-021-01721-z
Souzy, M., Lhuissier, H., Méheust, Y., Borgne, T. L., & Metzger, B. (2020). Velocity distributions, dispersion and stretching in three-dimensional porous media. *Journal of Fluid Mechanics*, **891**, A16. https://doi.org/10.1017/jfm.2020.113

Valocchi, A. J., Bolster, D., & Werth, C. J. (2018). Mixing-limited reactions in porous media. *Transport in Porous Media*, **130**(1), 157–182.
https://doi.org/10.1007/s11242-018-1204-1

Vanderborght, J. (2001). Concentration variance and spatial covariance in second-order stationary heterogeneous conductivity fields. *Water Resources Research*, **37**(7), 1893–1912. https://doi.org/10.1029/2001wr900009

Weller, H. G., Tabor, G., Jasak, H., & Fureby, C. (1998). A tensorial approach to computational continuum mechanics using object-oriented techniques. *Computers in Physics*, **12**(6), 620–631. https://doi.org/10.1063/1.168744

Whitaker, S. (1967). Diffusion and dispersion in porous media. *AIChE Journal*, **13**(3), 420–427. https://doi.org/10.1002/aic.690130308
Souzy, M., Lhuissier, H., Méheust, Y., Borgne, T. L., & Metzger, B. (2020). Velocity distributions, dispersion and stretching in three-dimensional porous media. *Journal of Fluid Mechanics*, **891**, A16. https://doi.org/10.1017/jfm.2020.113

Valocchi, A. J., Bolster, D., & Werth, C. J. (2018). Mixing-limited reactions in porous media. *Transport in Porous Media*, **130**(1), 157–182.
https://doi.org/10.1007/s11242-018-1204-1

Vanderborght, J. (2001). Concentration variance and spatial covariance in second-order stationary heterogeneous conductivity fields. *Water Resources Research*, **37**(7), 1893–1912. https://doi.org/10.1029/2001wr900009

Weller, H. G., Tabor, G., Jasak, H., & Fureby, C. (1998). A tensorial approach to computational continuum mechanics using object-oriented techniques. *Computers in Physics*, **12**(6), 620–631. https://doi.org/10.1063/1.168744

Whitaker, S. (1967). Diffusion and dispersion in porous media. *AIChE Journal*, **13**(3), 420–427. https://doi.org/10.1002/aic.690130308