# Flow Kinematics in Three-Dimensional Porous Media of Varying Pore Size Distribution Using Smoothed Particle Hydrodynamics

Fan Chen, Junfeng Sun, Antoine Wautier, Mathieu Souzy

► To cite this version:

Fan Chen, Junfeng Sun, Antoine Wautier, Mathieu Souzy. Flow Kinematics in Three-Dimensional Porous Media of Varying Pore Size Distribution Using Smoothed Particle Hydrodynamics. Water Resources Research, 2025, 61 (10), 10.1029/2025WR040413. http://dx.doi.org/10.1029/2025WR040413

HAL Id: hal-05301729

https://hal.inrae.fr/hal-05301729v1

Submitted on 7 Oct 2025

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/2025WR040413

**Key Points:**

* The flow kinematics and advection-induced dispersion process in 3D porous media of spherical packing are studied via SPH

* Velocity distribution and the dispersion process are found quantitatively align with experimental observation in quasi-monodisperse case

* Increasing polydispersity results in narrower velocity distribution, restricted dispersion and decrease in streamline tortuosity

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

F. Chen,
fan.chen@inrae.fr

**Citation:**

Chen, F., Sun, J., Wautier, A., & Souzy, M. (2025). Flow kinematics in three-dimensional porous media of varying pore size distribution using smoothed particle hydrodynamics. *Water Resources Research, 61*, e2025WR040413. https://doi.org/10.1029/2025WR040413

Received 5 MAR 2025
Accepted 18 SEP 2025

**Author Contributions:**

**Conceptualization:** Fan Chen,

Antoine Wautier, Mathieu Souzy

**Formal analysis:** Fan Chen, Junfeng Sun,

Antoine Wautier, Mathieu Souzy

**Funding acquisition:** Mathieu Souzy

**Investigation:** Fan Chen, Junfeng Sun,

Antoine Wautier, Mathieu Souzy

**Methodology:** Fan Chen,

Antoine Wautier, Mathieu Souzy

**Resources:** Antoine Wautier

**Software:** Fan Chen, Junfeng Sun

**Supervision:** Antoine Wautier,

Mathieu Souzy

**Validation:** Junfeng Sun

**Visualization:** Fan Chen, Junfeng Sun

© 2025. The Author(s).
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>

Flow Kinematics in Three-Dimensional Porous Media of
Varying Pore Size Distribution Using Smoothed Particle
Hydrodynamics

Fan Chen¹ ID, Junfeng Sun² ID, Antoine Wautier¹ ID, and Mathieu Souzy¹ ID

¹Aix Marseille University, INRAE, RECOVER, Aix-en-Provence, France, ²Key Laboratory of Transportation Tunnel Engineering of Ministry of Education, Southwest Jiaotong University, Chengdu, China

**Abstract** The effect of pore size distribution on the flow kinematics and transport properties within a three-dimensional porous medium is investigated through numerical simulations using the Smoothed Particle Hydrodynamics (SPH) method. The method is first validated for a model porous medium within a monodisperse random spherical packing, for which the velocity distribution of the fluid flowing through the pores (i.e., the interstitial fluid velocity) and the dispersion process are found to be in both qualitative and quantitative agreement with previous experimental results. When varying the pore size distribution of the porous medium by using polydisperse beads (of different diameters), the interstitial fluid velocity distributions get narrower, and the streamlines' tortuosity decreases. This is interpreted as a result of the narrower pore size distribution reported for polydisperse microstructures. Although the dispersion process remains qualitatively the same among the investigated microstructures, with an initial ballistic trend followed by a transient seemingly anomalous regime and eventually a Fickian regime, the transverse dispersion process is found to be quantitatively reduced for polydisperse microstructure (i.e., with a narrower pore size distribution), consistently with the reported decrease in streamlines' tortuosity.

**Plain Language Summary** Water often flows through soils, rocks, and filters in ways that are much more complex than simple straight channels. Inside these porous materials, the pores have different shapes and sizes, which can strongly influence how water and dissolved substances (such as pollutants, nutrients, or medicines) move and mix. In this study, we use a computer simulation method named Smoothed Particle Hydrodynamics (SPH) to look closely at how fluid flows inside three-dimensional porous media made of spherical beads. First, we check that the method gives reliable results of transport properties by comparing it to laboratory experiments. Then, we study how changing the bead size distribution—going from nearly uniform beads to mixtures of different bead sizes—changes the flow behavior. We find that when the pore size distribution becomes narrower (which happens when the beads are more varied in size), the fluid paths become less tortuous and the range of velocities becomes smaller. As a consequence, the spreading of a solute cloud (dispersion) is reduced. These results help us better understand how microscopic structures control transport in porous media, which is important in many natural and engineering scenarios, including groundwater contamination, soil remediation, carbon storage, and drug delivery in biological tissues.

1. Introduction

Solute transport in porous media is of fundamental importance in a broad range of scientific fields and applications, including contaminant migration, improved oil recovery (Wang et al., 2023), the long-term fate of nuclear waste, geological storage of carbon dioxide, soil remediation (Dentz et al., 2016), drug delivery and nutrient transport in brain or plant tissues (Peyrounette et al., 2018), or heat exchangers, filters and catalytic processes (Anitha et al., 2022) which also rely on transport through porous media (Dentz et al., 2011; Le Borgne et al., 2015; Puyguiraud et al., 2021). Due to the inherent danger of contamination or the essence of soluble chemical purification in water supply, characterizing and modelling solute transport, that is, where does it spread and how long does it take to mix, is of crucial importance within three dimensional porous media. Modelling these systems requires an understanding of how the solid microstructure dictates how the fluid within the pores (from here on named the interstitial fluid) flows.

Since porous media usually exhibit a tortuous microstructure (Adeaga et al., 2023) with a wide size distribution of intricate pores, the characterization of the flow becomes complex as the transport mechanism result from the
**Writing – original draft:** Fan Chen,
Junfeng Sun, Antoine Wautier,
Mathieu Souzy

**Writing – review & editing:** Fan Chen,
Antoine Wautier, Mathieu Souzy

coupling between the flow and the solid interfaces (Skibitzke et al., 1963). As a result, the transport processes are fundamentally controlled by the distribution of pore-scale interstitial fluid velocities, which is itself driven by the variations in local flow rates owing to heterogeneous distribution of porosity (Dentz et al., 2018, 2022; Morales et al., 2017; Villermaux, 2019). For geomaterials at large scale, the presence of preferential paths often channelizes the flow along high-permeability zones (Koltermann & Gorelick, 1996). On the one hand, these preferential paths promote efficient dispersion by allowing the solute to efficiently spread while spatially exploring the high permeability regions. On the other hand, whether the (dispersed) solute is efficiently mixed depends on how the fluid material lines are elongated by the flow, which relies on the stretching history set by both the high and low permeability regions of the medium (Le Borgne et al., 2015). Considering that the range of permeability of various natural geologic porous media can span more than 13 orders of magnitude (Freeze & Cherry, 1979), in a highly heterogeneous medium, both preferential flow paths and low-permeability zones exist (Feehley et al., 2000; Rosenbom et al., 2014), resulting in a wide variance of hydraulic conductivity.

The presence of these low permeability regions also results in low flow velocity regions which are essential to characterize the dispersion process. Note that there are crucial distinctions between dispersion and mixing, and that these two processes are fundamentally entangled and complementary to characterize transport properties. Let's consider a blob of dye which is injected through a given medium. Dispersion addresses the question of how the blob explores its surrounding environment: what is the typical size of the “envelope” (area in 2D, volume in 3D) within which the blob of dye is contained? On the other hand, one may wonder how long does it take to mix the injected blob of dye, that is, how long does it take for the dye concentration probability distribution to get homogeneous? To address such question one has to characterize the mixing process, resulting from the coupling between the history of stretching of fluid material lines (which is also a purely advective process sets by the kinematic properties of the flow field) and the effect of molecular diffusion: mixing is stretching-enhanced molecular diffusion (Villermaux, 2019). Typically, in a porous medium, the blob of dye will rapidly form a lamellae-like topology (Le Borgne et al., 2015; Ottino et al., 1979), with elongated and folded lamellae (sheets in 3D) contained within the aforementioned “dispersion envelope.” The concentration field can be seen as a set of lamellae, whose concentration depends on the interaction between local stretching and diffusion. Compression perpendicular to the stretching direction enhances molecular diffusive smoothing by steepening the scalar gradient. Since the effect of molecular diffusion can be a posteriori accounted for using a dedicated framework for modelling mixing, such as the lamellar description which has already been extensively developed and used in porous media (Le Borgne et al., 2015; Meunier & Villermaux, 2010; Ranz, 1979; Villermaux, 2012, 2019), characterizing the kinematics of the flow is crucial to assess the transport properties, and understanding how the flow kinematic is controlled by the microstructure remains an open question. While stretching coupled with molecular diffusion describes when mixing of the lamellae will occur, the dispersion process describes where mixing will occur by prescribing the envelope in which the elongated lamellae are contained. Thus, understanding the effect of varying the pore size distribution over the flow kinematic in porous media is of paramount importance to characterize mixing in geomaterials.

A strong dispersion process generally leads to efficient mixing by promoting strong elongation of the transported solute. Consequently, characterizing dispersion has been a crucial question over the last decades, in particular in geomaterials. The low flow velocity regions in the low permeability zones induce what has been commonly referred to “non-Fickian” or “anomalous” advective dispersion of pollutants or solute (Levy & Berkowitz, 2003; Souzy et al., 2020; Stephen & Simpson, 1987; Tanguy et al., 2010). The dispersion process is usually characterized through the temporal evolution of the mean squared displacement of fluid solute $\sigma^2 \propto t^a$. While Fickian dispersion relates to the case where $\alpha = 1$, non-Fickian dispersion can be either sub-diffusive ($\alpha < 1$), or super-diffusive ($\alpha > 1$). Multiple physical theories and numerical models of non-Fickian transport have been developed in the past 30 years (Guo et al., 2021), among which stochastic models such as the Continuous Time Random Walk Method (CTRW) (Anna et al., 2013; Souzy et al., 2020). The dispersion process using CTRW fundamentally depends on the probability distribution of interstitial fluid velocities (Lista, 2023; Souzy et al., 2020). In particular, the shape of the distribution over the low-velocity range has been shown to set the anomalous behavior of the dispersion process (Berkowitz et al., 2006; Souzy et al., 2020). While the link between the dispersion process and the velocity distribution received peculiar attention and remains an active topic nowadays, the shape of the velocity distributions remains an open question in particular within the low-velocity range. Using a dynamic-range velocity measurement technique combined with index matching, (Souzy et al., 2020) reported a flat distribution at low velocities, in agreement with numerical simulations (Dentz et al., 2018). On the other hand,
using similar optical index matching techniques, (Carrel et al., 2018; Morales et al., 2017) reported in their experiments velocities exhibiting power-law distributions within the low-velocity range. Deciphering the effect of the microstructure and the pore size distribution on these results is therefore crucial to bring light to this open debate and could therefore benefit from reliable numerical simulations.

Based on this observation, the current paper aims at providing a better understanding of how the microstructure of porous medium at the pore scale affects transport properties. Specifically, several key research questions are addressed: how do pore-scale velocity distributions emerge from the geometry of sphere packings, and to what extent can numerical methods quantitatively reproduce experimentally reported flow statistics? What is the influence of poly-dispersity of the granular packing—that is, increasing variance in grain size within a packing—on the velocity distributions and streamline tortuosity? Does it modify the dispersion process, by promoting or impeding the transient anomalous dispersion regime?

To address such questions, we aim at using numerical simulations benchmarked on previous experimental results in a specific configuration, before varying parameters such as the pore size distribution. Geomaterials commonly exhibit small pore size which constrains the flow. Consequently, the typical Reynolds number in the porous medium is much smaller than 1, a signature of negligible inertial effects. Within such framework, various numerical methods adopting the Stokes equation have been developed and employed to simulate this solid-fluid coupled phenomenon. Various numerical approaches consist of representing a porous medium by the Discrete Element Method (DEM) and coupling with a spatially resolved fluid scheme adopting Finite Volume Method (FVM) (Moukalled et al., 2015), Finite Element Method (FEM) (Knight et al., 2020), Pore Finite Volume network (PFV) (Catalano et al., 2014), or Smoothed Particle Hydrodynamics (SPH). Unlike FEM or FVM using Eulerian framework, the grain-scale SPH has a natural mesh-free Lagrangian perspective, making the tracing of the fluid dispersion dynamics becomes straightforward. This is a valuable asset which make SPH well suited for investigating dispersion throughout numerical simulations. In contrast, tracing the fluid phase dispersion using FVM and FEM may require specific post-processing algorithms based on their Eulerian velocity fields outputs. Furthermore, although such supplementary asset is not used in this study, conventional Computational Fluid Dynamics (CFD) methods rely on grid resolution and may face challenges of frequent re-meshing in case of irregular boundary interfaces or changing geometries due to fluid-driven granular material movement (Yueting et al., 2023). Recently, SPH has been implemented in varied large-scale engineering application scenarios (Canelas et al., 2016; Manenti et al., 2019) and validated for modeling incompressible flows at low Reynolds numbers at grain-scale, in spite of the potential occurrence of small pressure fluctuations near curved boundaries (Morris et al., 1997).

In the current work, we aim to investigate the effect of varying the pore size distribution of a porous medium (by varying the bead size within a granular packing) over the interstitial fluid velocity distribution and the dispersion process, through numerical simulations performed using SPH. The numerical SPH method is first validated by comparing the distributions of the velocity magnitude and individual components, and the dispersion process, with previous experiments from Souzy et al. (2020) realized using optical refractive index matching and particle tracking techniques to measure the three-dimensional velocity field in an isotropic porous medium composed of randomly packed solid spheres of the same diameter, at the pore scale (typically few hundreds of beads composing the medium). These experiments were used to show that the dispersion of non-diffusive tracers nonetheless becomes Fickian beyond a time set by the smallest effective velocity of the tracers, and also to reconstruct the stretching histories of fluid material lines, experimentally evidencing the chaotic nature of advection within three-dimensional porous media. After validating the SPH method over this high resolution experimental data set, this numerical tool is used in the current study to explore the influence of microstructure configurations induced by varying the particle size distribution (PSD), and thus the pore size distribution.

The paper is organized as follows: first the brief framework of SPH is introduced in Section 2. In Section 3, the numerical setups for simulating the flow in a porous medium composed of a random monodisperse spherical packing are introduced, mimicking a configuration similar to the experiments by Souzy et al. (2020). Velocity distributions, dispersion process, and streamlines' tortuosity are analyzed and compared in detail throughout Section 4 for this model microstructure configuration. In Section 5, simulations are then performed to investigate and discuss the effect of polydispersity by conducting a series of simulations in porous media with varied pore size distribution.
## 2. SPH: Numerical Methodology

As stated above, SPH presents various assets and features which makes it a particularly interesting numerical approach for investigating flow in complex porous microstructures where interactions between fluid and solid phases are crucial when the pore-scale dynamics influence the macroscopic flow properties. Using SPH approach, adaptive resolution of the flow is performed, which facilitates handling large deformations and complex boundary interactions without complex grid remeshing protocols. A brief introduction of the applied SPH algorithm in open-source software DualSPHysics and its boundary conditions is given below (Domínguez et al., 2022).

### 2.1. Governing Equations of SPH

SPH is a Lagrangian mesh-less method that discretizes continuum fluid phase using smooth points or particles. Such an approach is commonly applied in a broad spectrum of engineering domains, spanning from free-surface flows (Lind et al., 2015; Lyu et al., 2021; Mokos et al., 2016; Monaghan, 1994; Shao et al., 2012) to human cardiac function simulations (Zhang et al., 2021). The discretized Navier-Stokes equations are locally integrated according to the physical properties of neighbouring particles (Crespo et al., 2014). Density, pressure, and velocity fields of a fluid particle at position $\mathbf{r}$ are estimated by an integral representation of the smooth kernel function $W(\mathbf{r} - \mathbf{r}', h)$:

$$f(\mathbf{r}) = \int_{\Omega} f(\mathbf{r}') W(\mathbf{r} - \mathbf{r}', h) d\mathbf{r}' \quad (1)$$

$$W(r, h) = \frac{1}{\pi h^3} \begin{cases} \frac{3}{4}q^3 - \frac{3}{2}q^2 + 1 & 0 \le q < 1 \\ \frac{1}{4}(2-q)^3 & 1 \le q < 2 \\ 0 & q \ge 2 \end{cases} \quad (2)$$

$$q = |\mathbf{r} - \mathbf{r}'|/h \quad (3)$$

where $\Omega$ represents the support domain of a particle with position vector $\mathbf{r}$; $\mathbf{r}'$ denotes a neighbouring particle position vector in $\Omega$; $q$ is the distance between these two neighbouring particles, normalized by the smoothing length $h$ which defines the spatial radius over which the properties of a particle are smoothed or averaged with the neighbouring particles, that is, the characteristic length scale of the influencing area.

The mass and momentum conservation in the SPH framework are described by Domínguez et al. (2022):

$$\frac{d\rho_a}{dt} = \sum_b m_b (\mathbf{v}_a - \mathbf{v}_b) \cdot \nabla_a W(\mathbf{r}_{ab}, h) + \Phi_a \quad (4)$$

where a density diffusion term $\Phi_a$ is introduced in the continuity equation from Diego and Colagrossi (2009) to eliminate spurious numerical noise in the pressure field fluctuations:

$$\Phi_a = 2\delta hc_0 \sum_b \frac{m_b (\rho_b - \rho_a)}{\rho_b} \frac{(\mathbf{r}_a - \mathbf{r}_b) \cdot \nabla_a W(\mathbf{r}_a - \mathbf{r}_b, h)}{|\mathbf{r}_a - \mathbf{r}_b|^2} \quad (5)$$

where $a$ and $b$ in the subscripts denote fluid particles; $m$, $\rho$, and $\mathbf{v}$ refer respectively to the mass, density, and velocity of fluid particle; $\nabla_a$ is the kernel gradient, computed based on the coordinates of particle $a$. Note that $\delta$ is a constant model parameter which is commonly set to $\delta = 0.1$ (Li et al., 2020; Zhan et al., 2019). The final form of the momentum equation incorporates the pressure term, laminar viscosity, and external forces given as:

$$\frac{d\mathbf{v}_a}{dt} = -\sum_b m_b \left( \frac{P_b + P_a}{\rho_b \cdot \rho_a} \right) \nabla_a W_{ab} + \sum_b m_b \left[ \frac{4\nu_0 \mathbf{r}_{ab} \cdot \nabla_a W_{ab}}{(\rho_a + \rho_b)(r_{ab}^2 + \eta^2)} \right] \mathbf{v}_{ab} + \mathbf{g} \quad (6)$$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Explanatory diagram of the ghost particle method in the fluid domain at the SPH-DEM interface, next to a wall (on the left) and a corner (on the right).</figcaption></figure>

where $v_0$ is the kinematic viscosity (typically $10^{-6}$ for water); $\mathbf{r}_{ab} = \mathbf{r}_a - \mathbf{r}_b$ is the relative distance between particles a and b; $\mathbf{v}_{ab} = \mathbf{v}_a - \mathbf{v}_b$ is the relative velocity vector. $\nabla_a W_{ab}$ is the gradient of the smoothing kernel evaluated at particle a due to particle b. $\eta^2$ is a small regularization parameter to avoid singularities in the denominator. g represents the gravitational acceleration.

## 2.2. SPH Boundary Conditions

As simulations aim to investigate the interstitial flow within the porous medium, the static assemblies of solid grains are generated using Discrete Element Method (DEM), while the DEM inter-grain contact algorithm are turned-off during the simulation process. The DEM particles' surfaces are considered as rigid boundaries, using integrated stresses on the surfaces to compute interaction forces exerted from the DEM solid particles to the SPH fluid phase. DEM solids are discretized as SPH boundary particles, maintaining the same numerical resolution as the fluid phase.

To simulate the fluid-particle boundary interactions, the ghost particle method (Bouscasse et al., 2013), extrapolated pressure method (Canelas et al., 2016), and dynamic boundary condition (Xiong et al., 2023) are employed, aligning with the fluid-boundary treatment. Figure 1 shows a ghost particle mirrored into the fluid domain by the DEM boundary particle interface. Subsequently, first-order consistent SPH spatial interpolation is conducted on the surrounding fluid particles of the ghost node to determine its density and gradient. Then, a linear extrapolation is applied to derive the density of the boundary particle, which yields:

$$ \rho_b = \rho_g + (\mathbf{r}_b - \mathbf{r}_g) \cdot (\nabla \rho_g) \quad (7) $$

where $\mathbf{r}_b$ and $\mathbf{r}_g$ represent the position of the boundary particle and its corresponding ghost particle, respectively, and $\langle \nabla \rho_g \rangle$ is the corrected SPH gradient at the ghost node. The ghost particle velocity is thus determined through a Shepard filter sum:

$$ \mathbf{u}_g = \frac{\sum_{g'} \mathbf{u}_{g'}' W V_{g'}}{\sum_{g'} W V_{g'}} \quad (8) $$

where $g'$ accounts for the neighbouring particles of g. The no-slip boundary condition at the interface is then set by tuning the velocity of the boundary particle to $\mathbf{u}_g$ with its direction reversed. A more detailed description on the outline framework of the SPH algorithm can be found in Domínguez et al. (2022).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Numerically simulated SPH domain and 3D porous medium microstructure generated using DEM. (b) Fluid boundary conditions and (c) SPH particle spacing $d_p$.</figcaption></figure>

## 3. Simulating 3D Porous Media

Before using the model to explore the variation of pore size distribution in porous material of increasing microstructure complexity, benchmark configurations are used to validate the model. In addition to a first validation using the reference benchmark configuration of a single sphere water-entry experiment (Aristoff et al., 2010; Xu et al., 2020) (see details in Supporting Information S1), results of the SPH simulations are compared to experimental results within a three-dimensional porous medium composed of a random pack of quasi-monodisperse spheres, imposing a constant Stokes flow.

### 3.1. Parameters and Boundary Conditions

This first simulation in a 3D porous medium is performed using a microstructure consisting of a random pack of around 470 spheres, as shown in Figure 2. The solid volume fraction is 0.58, that is, a porosity $\epsilon = 0.42$. Such a configuration intends to mimic experiments by Souzy et al. (2020) performed over a volume $5d \times 5d \times 9d$ contained in the bulk of a large ($12d \times 12d \times 100d$) column of random pack of $d = 2$ mm diameter beads. Using optical refractive index matching and particle tracking techniques, the authors measured the velocity field to investigate the dispersion and stretching processes, experimentally evidencing the chaotic nature of advection within the three-dimensional porous medium.

To prevent crystallization within the microstructure, simulated spheres are taken from a narrow bi-disperse size distribution (i.e., a bimodal distribution). This is performed by imposing that the sphere mass distribution is bifurcated between two distinct sizes: half of the total mass is allocated to spheres with a diameter of 70 mm and the other half to particles with a diameter of 80 mm. In the following, the mass median diameter $D_{50\%}^w = 75$ mm (also often referred to as MMD) will be used as the characteristic length scale of the system. The introduction of such a degree of dispersity with such a small size difference between spheres size results in a quasi-monodispersed packing with its corresponding microstructure yet prohibiting crystallization effect, as shown in Figure 2. Note that from herein, such configuration is labeled as PSD1 (Particle Size Distribution 1), to
<table><caption>**Table 1**

Numerical Parameters Used in the SPH Model</caption><thead><tr><td>Parameter</td><td>Unit</td><td>Value</td></tr></thead><tbody><tr><td>Mass density</td><td>kg/m³</td><td>1,000</td></tr><tr><td>Kinematic viscosity</td><td>m²/s</td><td>1</td></tr><tr><td>Density diffusion value (Φ<sub>a</sub>)</td><td>-</td><td>0.1</td></tr><tr><td>Initial particle spacing (d<sub>p</sub>)</td><td>mm</td><td>2</td></tr></tbody></table>

distinguish it from the other microstructure in poly-dispersed granular packings that will be further introduced and discussed in Section 5.

As depicted in Figure 2, the SPH fluid region dimensions are $4D_{50\%}^w \times 4D_{50\%}^w \times 13.3D_{50\%}^w$. To prevent rigid boundary effects where porosity would increase at the particle-boundary interface, the generated DEM domain is larger than the SPH fluid region, thus allowing for a homogeneous porosity within the SPH fluid region of interest. The lateral boundaries of the SPH region are treated as non-slip boundary. Additionally, to prohibit fluid inflow fluctuations at the entrance of the DEM region, the

SPH domain is elongated in the longitudinal flow direction at both ends. Such a method is implemented to get stabilized and homogeneous in-flow and out-flow within the porous medium.

The interstitial fluid, characterized by a density $\rho = 1000 \text{ kg/m}^3$, flows in the longitudinal direction (x-axis) by imposing a stable flow rate $Q = 3.6 \times 10^{-5} \text{ m}^3/\text{s}$ using fixed-velocity Dynamic Boundary Condition (DBC) as documented in DualSPHysics (Crespo et al., 2014), resulting in an average velocity magnitude within the porous medium $\langle U \rangle = 3.24 \times 10^{-3} \text{ m/s}$. To ensure a laminar flow with a Reynolds number $Re = \frac{\langle U \rangle}{\nu} D_{50\%}^w = 2.43 \times 10^{-4} \ll 1$, the kinematic viscosity of the fluid in the simulation is tuned to $\nu = 1 \text{ m}^2/\text{s}$. The simulations aiming at solely investigating the effect of the advection field (i.e., the flow kinematics) and the resulting dispersion process, molecular diffusion effect (i.e thermal/Brownian agitation) are not implemented. Note that this is consistent with (Souzy et al., 2017), as the Péclet number in the experiments was large enough ($>10^6$) such that molecular diffusion was negligible compared to the effect of the advection. To ensure accurate simulation precision while limiting the computational cost at the same time, the initial particle spacing is set such that $d_p = 2 \text{ mm}$ (thus $D_{50\%}/d_p = 37.5$) to provide a high resolution of fluid flow at the irregular particle surface. Overall, a total number of 5,204,463 SPH fluid particles are simulated, and the detailed parameters are summarized within Table 1.

## 3.2. Initial Transient Flow

When modeling low Reynolds numbers incompressible flows at grain scale, SPH induces small-scale transient pressure fluctuations nearby curved boundaries (Morris et al., 1997). This results in an overall transient flow at the beginning of a simulation. Such fluctuations however rapidly vanish and the flow eventually reach an almost steady state. As experiments were performed applying a steady flow by imposing a constant flow rate, simulations are performed until the flow temporal fluctuations get negligible. To satisfy such criterion quantitatively, the probability distribution function of the local velocity difference $\Delta U$ over a temporal lag $t/\tau$ is computed and plotted in Figure 3a, where $\langle U \rangle$ denotes the average velocity and $\tau = D_{50\%}/\langle U \rangle$ is the Stokes time, i.e. The typical time it takes for a fluid element to travel the characteristic length scale. As expected, although asymmetrical initially due to transient flow, the PDF of velocity differences quickly converge to a symmetrical narrow peaked distribution centered around 0, a signature that the flow gets steady over the simulated domain. This steady state regime can also be contemplated by computing the temporal evolution of the standard deviation of $\Delta U$, see

<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) The probability distribution function of the temporal velocity variations $P(\Delta U/\langle U \rangle)$ converges to a peaked distribution after a short initial transient; (b) Temporal evolution of the corresponding standard deviation of $P(\Delta U/\langle U \rangle)$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Internal cross-section of the model's 3D velocity field and probability distribution function of velocity field components along (b) the flow direction, (c) the transverse directions, and (d) of the velocity magnitude.</figcaption></figure>

Figure 3b. For $t/\tau > 1.2$, the flow field becomes steady and only negligible remaining temporal velocity fluctuations occur. In the following, every result is taken in the steady state regime for $t/\tau > 1.2$.

# 4. Results and Discussion for a Monodisperse Configuration

In this section, numerical simulations of the flow kinematics for the PSD1 microstructure mimicking experiments performed by Souzy et al. (2020) in a monodisperse configuration are compared with experimental measurements.

## 4.1. Velocity Field Distribution

The 3D velocity field is shown together with the porous medium microstructure for the configuration PSD1 in Figure 4a, highlighting heterogeneous velocity distributions with low velocity and high velocity regions. The distribution of the velocity components ($U_x, U_y, U_z$) and magnitude $U = \sqrt{U_x^2 + U_y^2 + U_z^2}$ are presented in Figures 4b–4d. Note that the velocity distributions obtained through the numerical simulation are similar to both previous experiments (Carrel et al., 2018; Datta et al., 2013; Holzner et al., 2015; Souzy et al., 2020) and previous numerical simulations (Dentz et al., 2018): the distribution of the streamwise component ($U_x$) is skewed towards positive values and has an exponential tail, revealing prevailing fluid flow alignment with preset flow orientation.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Distribution of streamlines orientation $P(\cos \theta)$ relative to the longitudinal direction, and average velocity magnitude for a specified orientation $\langle U \rangle_{\cos \theta}$. (b) Trajectories of non-diffusive SPH particles during the process of being advected numerically in the 3D porous medium (PSD1).</figcaption></figure>

The tail of the distribution for large $U_x$ is well fitted by an exponential $P(U_x/\langle U \rangle) = e^{-K*(U_x/\langle U \rangle)}$ with $K = 0.41$. Interestingly, the distribution also exhibits values $U_x < 0$, indicating the presence of back-flow regions as already reported (Beguin et al., 2013; Carrel et al., 2018; Dentz et al., 2018; Holzner et al., 2015; Souzy et al., 2020), although with a noticeably weak probability. The distributions of the transverse velocities ($U_y$, $U_z$) are symmetric as expected, and have exponential shapes with an estimated declining slope of transverse velocity $K \approx 0.55-0.6$. Last, the distribution of the velocity magnitude also exhibits an exponential tail for high velocity magnitudes. Interestingly, the distribution seems to reach a plateau at low velocity, in agreement with dynamic range velocity measurements reported by Souzy et al. (2020). This is particularly important, because as will be further discussed in the next section, it is the characteristics of the velocity distribution at small velocities which control the long-term dispersion behaviour (Bijeljic et al., 2011). Interestingly, the distribution of velocity magnitude obtained from the numerical simulations exhibits a higher probability of large velocity magnitude events when compared to experimental measurements. Note that contrary to the experiments where the fluid freely exits and enters from the lateral sides of the experimentally investigated domain, there is a no-slip boundary condition on the lateral boundaries of the SPH simulations. We expect this difference in the boundary condition to result in the promotion of slightly higher-velocity regions compared to experiments. It could also result from technical limitations in Souzy et al. (2020) which makes it difficult to measure experimentally velocities larger than $5\langle U \rangle$ within their setup.

In addition, Figure 5a presents the statistics of the streamline local orientations $P(\cos \theta)$ and the mean velocity magnitude for a given orientation $\langle U \rangle_{\cos \theta}/\langle U \rangle$, where $\theta$ defines the streamline orientation with regard to the longitudinal flow direction (x-axis) (see inset of Figure 5a). Interestingly, the probability distribution function of velocity orientation $P(\cos \theta)$ follows a similar exponential fitting trend as reported experimentally in Souzy et al. (2020) as $P(\cos \theta) \propto e^{a \cos \theta}$, with a similar reported value $a \approx 4$. This trend simply reflects that streamlines are strongly aligned with the longitudinal direction, as exemplified over the streamlines shown in Figure 5b. Besides, results show that the magnitude of the velocity is also strongly correlated with its orientation. The average projected velocity magnitude for a given orientation $\theta$ increases when streamlines are more aligned with the flow, following the exponential trend $\langle U \rangle_{\cos \theta}/\langle U \rangle \propto e^{b \cos \theta}$, where $b \approx 1.5$, also in quantitative agreement with experiments.

## 4.2. Fluid Dispersion Law

Dispersion is a purely advective process which characterizes how the flow explores its surrounding. Typically, when a cloud of non-diffusive passive tracers is injected within a flow, the dispersion process will characterize the evolution of the “envelope” containing all the tracers. Note that non-diffusive passive tracers always flow following the streamline on which they are initially placed, as in the absence of molecular diffusion (no random Brownian motion due to thermal agitation) for point tracers there is no mechanism allowing them to “jump” from their streamline onto a neighbouring streamline. Although experimentally such quantification relies on injecting a cloud of non-Brownian passive particles (i.e., few microns size particles which are too large to experience thermal
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Successive images of a cloud of SPH particles highlighted in red, which progressively disperse while flowing through the microstructure PSD1 at different time of advection from (a) to (d). Flow direction is from left to right.</figcaption></figure>

agitation), and monitoring the spatial variance of the passive tracers while flowing through the medium, accessing such information is trivialised thanks to the use of SPH simulations. By recording the spatial coordinates of each individual SPH particle, the dispersion can thus be quantitatively characterized as each SPH particle is equivalent to a passive tracer which would flow on its local streamline throughout the medium. Figure 6 shows an example of a cloud of SPH particles (tagged in red) initially contained within a spherical space of diameter $40d_p$ and progressively dispersing in the surrounding medium, as each of these particles is passively advected. As can be noted in Figure 6, after $t/\tau = 4.8$, some SPH particles have already reached the boundary of the simulated domain. To allow a long-time advection ($t/\tau \approx 100$), each SPH particle that reaches the rigid boundary of the fluid domain is re-injected into the upward stream at a location where the fluid velocity is the closest (both in magnitude and orientation) to its current velocity. Note that although such re-injection procedure may promote velocity decorrelation, the no-slip condition imposed on the lateral boundaries of the fluid domain and the important length of the fluid domain in the longitudinal direction ($13.3D_{50\%}^w$) result in very low frequency of such re-injection events, such that it was found to have no quantitative effect on the dispersion statistics when compared to analysis performed (over shorter extent of time) without implementing re-injection protocol. Hence, the re-injection procedure allows to explore dispersion over much larger time-scale. Results are gathered from 5 million passive SPH particles for $t/\tau \le 96$, and plotted on Figure 7.

The dispersion process is quantified by measuring the variance of all the particles' positions, both in the longitudinal and transverse directions, that is, $\sigma_i^2 = \langle l_i^2 \rangle - \langle l_i \rangle^2$, where $l_i$ represents the displacement along the x, y, or
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Successive positions of a cloud of SPH particles flowing through the microstructure PSD1 with their respective spatial distributions along the flow direction in short (a) and long (b) time-scale. The cloud progressively spreads and the spatial distribution gets wider as the tracers are advected. (c) Corresponding variances of the position of the tracers $\sigma_i^2 = \langle t_i^2 \rangle - \langle t \rangle^2$ with $i = x, y$ or $z$ from simulation (solid line) and experiments (dashed line) by Souzy et al. (2020).</figcaption></figure>

z axis. As depicted in Figure 7, after a short ballistic regime, where $\sigma_x^2 \propto t^2$ for $t/\tau \lesssim 1$, the variance becomes seemingly anomalous for $1 \lesssim t/\tau \lesssim 10$. It can be interpreted as super-diffusive in the stream-wise direction ($\sigma_x^2 \propto t^{1.4}$) and sub-diffusive $\sigma_y^2, \sigma_z^2 \propto t^{0.8}$ along the transverse orientation (y, z). More importantly, for $t/\tau > 10$, both longitudinal and transverse dispersion eventually relax to a Fickian regime ($\propto t^{1.0}$) with variances following
$$ \sigma_i^2 = 2Dt, \qquad (9) $$

where D is the effective dispersion coefficient (in m²·s⁻¹).

Note that similar results were reported in the transient regime by Kang et al. (2014) over a similar overall range of t/τ, and also by Souzy et al. (2020) whose experimental results are reported on Figure 7c for comparison. These were obtained by monitoring the evolution of a blob of small non-diffusive particles, after injecting it in the bulk of a porous medium made of a random monodisperse bead packing. The particles were non-diffusive fluorescent tracers (PMMA B-particles from MF-Rhodamine) with a diameter of 3.23 µm. Once the injection is completed, a constant flow rate is imposed, and the evolution of the blob is captured within one plane illuminated by a laser sheet. Similarly to what is done here using numerical advection, the dispersion process was quantified (through image processing) by measuring the spatial variance of the tracer locations both in the longitudinal and the transverse directions.

The ballistic dispersion regime lasts typically until t ~ τ, which reveals that the velocity of the tracers remains correlated over a typical length D<sup>w</sup><sub>50%</sub>. Beyond, velocities in the longitudinal direction are essentially positive and the (super-diffusive) transient behaviour is controlled by the occurrence of the smallest velocities, which significantly delay the longitudinal transport of some of the tracers as observable on Figure 7b. However, this transient apparently anomalous regime eventually relaxes to a Fickian regime. This could appear as a strike contradiction with the theoretical work by Berkowitz et al. (2006), as they have shown that a velocity distribution exhibiting a plateau at low velocity should necessarily lead to long term anomalous dispersion behaviour (with σ<sub>x</sub><sup>2</sup> ∝ t log(t)): the presence of a plateau at low velocity in the distribution means that there are low velocity regions in the flow such that any passively advected tracer which would initially explore such region will remain with its velocity correlated for a duration which is inversely proportional to its initial velocity. The dispersion process becomes Fickian once every advected tracer is decorrelated from its initial value, any flow exhibiting an asymptotic plateau at low velocity will result in non-Fickian anomalous transport.

However here, as illustrated on Figure 4d, although reaching an apparent plateau at low velocities, it is not an asymptotic plateau as the lowest values experienced by SPH particles is U<sub>min</sub> ~ ⟨U⟩/10. Such limitation arises from the fact that SPH particle velocities are computed using kernel functions that average over a finite size, which prevents exploring lower velocities in the very-close neighbouring region of solid boundaries. Consequently, and following previous works (Bijeljic et al., 2013; Dentz et al., 2018; Souzy et al., 2020), this sets a finite minimal velocity which has a crucial impact on the longitudinal dispersion, precisely because the dispersion anomaly is controlled by the low velocities. The Fickian regime is reached once the velocities of all the tracers become decorrelated, thus after a typical time t<sub>F</sub>/τ ≈ 20. This is not only consistent but also in excellent quantitative agreement with Continuous Time Random Walk performed by Souzy et al. (2020) based on a similar exponential distribution of velocities but discarding all velocities below U<sub>min</sub> ~ ⟨U⟩/10, highlighting how such approach is robust to describe how the velocity distribution prescribes the dispersion process. Note that although this could be considered as a limitation of the SPH numerical method as it prohibits the SPH particles to experience velocities U < ⟨U⟩/10, flows in real porous medium also exhibit similar finite minimal velocity: if one considers a passive tracer injected within a porous medium, (a) either the tracer is large enough (≥ few microns) such that it is non-Brownian, in which case due to steric finite size effect it cannot explore the slowest streamlines in the close vicinity of the solid interfaces, (b) either it is smaller but then thermal agitation will “kick” it out of the low velocity streamlines such that it cannot remain there for extended duration, which in turn leads to the same result as if it would have prohibited tracers from exploring the low velocity region of the distribution, and which promote lagrangian velocity decorrelation, thus resulting here also in a Fickian long term dispersion behaviour (Souzy et al., 2020). Consequently, although with different fundamental origins, both SPH simulations and real-life situations end up with a range of low velocities which cannot be explored by advected tracers. The quantitative correspondence of such finite minimal velocity between experiments and the SPH simulations makes the use of the SPH method particularly well suited for investigating and characterizing the dispersion process in such systems.

Contrastingly, the reported sub-diffusive transient dispersion process in the transverse direction is a consequence of the tortuosity of the streamlines inside the three dimensional porous medium, which induces an alternance
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Average total streamline elongation *L* and average projected length of the streamlines on both the longitudinal flow direction *L*<sub>x</sub> and the transverse direction *L*<sub>y</sub>. (b) The distribution of tortuosity *P*(*T* = *L*/*L*<sub>x</sub>) follows a Gaussian distribution which gets narrower over time.</figcaption></figure>

between positive and negative transverse velocities with similar magnitudes, which leads to a transient sub-diffusive regime before eventually reaching the Fickian regime after the same typical time $t_F$.

4.3. Flow Streamline Tortuosity

This tortuosity of the streamlines reflects the flow complexity within a porous medium. As such, it can be used as a proxy to quantify the meandering of the streamlines, as exemplified on Figure 5 or in the inset of Figure 8a showing a typical portion of a streamline. While the trajectory is mainly directed towards the flow longitudinal direction, it alterns in the transverse directions between positive and negative displacements. Thus, the total elongation of its paths *L* is the cumulative sum of not only its longitudinal displacements along *x*, but also the transverse displacements along *y* and *z*. Figure 8a shows the evolution of the average streamline total elongation ⟨*L*⟩ normalized by *D*<sub>50%</sub><sup>w</sup>, and the respective projection on both the longitudinal flow direction ⟨*L*<sub>x</sub>⟩ and the transverse direction ⟨*L*<sub>y</sub>⟩. For sake of simplicity, and due to the transverse isotropy, only ⟨*L*<sub>z</sub>⟩ is plotted. As expected, ⟨*L*<sub>x</sub>/D<sub>50%</sub><sup>w</sup>⟩ grows linearly with a slope of 1. Noticeably, ⟨*L*⟩ also exhibits a linear trend, although with an increased slope, reflecting the streamlines entanglement. One can thus characterize the tortuosity of the flow by computing the tortuosity index $T = L/L_x$, the ratio between the total streamline elongation and the displacement in the longitudinal direction. For PSD1 the tortuosity index ⟨*T*⟩ = 1.33. In addition, $L_z$ fluctuates around 0 for 0 ≤ *t*/*T* ≤ 100, a signature of the transverse isotropy of the system.

Figure 8b shows distributions of tortuosity computed by analyzing the temporal evolution of the trajectories over the 5,204,463 advected SPH fluid particles. Two remarks should be made here: (a) the apparent unsteadiness of the tortuosity distributions is directly attributed to the temporal advection approach used here; as discussed before, the flow exhibits very low velocities regions, such that SPH advected particles within such regions remains there for long durations. Consequently for short times $t/\tau$ the velocity of each fluid particle is not yet decorrelated from its initial value, resulting in partial streamline exploration, thus temporally affecting our measurements of $P(T)$ until each advected SPH particle explores a large enough portion of the local streamline on which it is transported. Streamline tortuosity being a static property dictated by the local microstructure and imposed flow conditions, in the following it will be characterized using the average value ⟨*T*⟩ for large enough $t/\tau$. (b) Over the explored range of time $t/\tau \le 96$, $P(T)$ is well-fitted by a normal distribution, which gets narrower over time. For PSD1 and $t/\tau = 96$, the standard deviation of the best Gaussian curve fit plotted in Figure 8b is $\sigma_T = 0.026$. This reveals that some streamlines are more tortuous than others. Whether this is an artifact feature resulting from the way tortuosity is computed by performing Lagrangian advection over a set of random streamlines over a *finite* time duration remains an open question, as one can argue that for infinitely long streamlines, the random microstructure would asymptotically lead to a Dirac distribution of tortuosity. Addressing such an intriguing question would, however, require measuring streamline tortuosity using a non-Lagrangian approach and over flow exploring much further in the longitudinal flow direction, which is out of the scope of the current work.
<table><caption>**Table 2**
Simulated Microstructures Ranging From Quasi-Monodispersed (PSD1) to Widely Bidispersed (PSD4) Granular Packings</caption><thead><tr><th>Parameter</th><th>Unit</th><th>PSD1</th><th>PSD2</th><th>PSD3</th><th>PSD4</th></tr></thead><tbody><tr><td>Small spheres size (50% in mass)</td><td>mm</td><td>80</td><td>50</td><td>40</td><td>30</td></tr><tr><td>Large spheres size (50% in mass)</td><td>mm</td><td>70</td><td>100</td><td>110</td><td>120</td></tr><tr><td>Size ratio</td><td></td><td>1.14</td><td>2.00</td><td>2.75</td><td>4.00</td></tr></tbody></table>

# 5. Effect of Varying the Polydispersity on the Flow Kinematics

## 5.1. Simulated Microstructures

After investigating the case of a quasi-monodisperse microstructure and comparing it to experimental measurements, we now turn to the effect of microstructure characterization made of granular packings of progressive increase poly-dispersity over the flow kinematics. This is achieved by investigating four bi-modal size distributions, where the mass median

diameter is maintained constant, $D_{50\%}^w = 75$ mm. Similarly to what was aforementioned for PSD1, the mass distribution is bifurcated between two distinct sizes: half of the total mass is allocated to spheres with a larger diameter and the other half to particles with a smaller size. Including PSD1, the ratios of the two different diameters in all the PSDs are 1.14, 2.00, 2.75, and 4.00 from quasi-monodisperse (PSD1, low polydispersity) to widely bidisperse (PSD4, high polydispersity) as given in Table 2. The choice of such microstructures allows configurations displaying a constant porosity $\epsilon = 0.42$ (thus a solid fraction $\phi = 0.58$). The generated numerical packings are displayed in Figure 9a.

The local shortest distance from any point in the void space to the nearest solid boundary ($L_s$) is computed for each point in the void space using a maximal inscribed sphere method (Chen et al., 2023), and can be used as a proxy to characterize the local pore size. Figure 9b shows that the probability distributions $P(L_s)$ for the various microstructures are found to be always peaked around $0.05D_{50\%}^w$, and that the average distance to a solid interface $\langle L_s \rangle \approx D_{50\%}^w/10$.

It is worth noting that the definition of the size of a pore varies within the literature; among the various approaches, one consists in defining the local pore size from the radius of the maximal inscribed sphere one could place on the skeleton of the void space. The skeleton could be obtained by extracting the centerline of the flow in the pore space while preserving the flow topology (Alim et al., 2017). Following such an approach, the corresponding pore size distributions $P(pore\ size/D_{50\%}^w)$ for the various investigated PSDs are plotted (in green) on Figure 9b. This reveals that although exhibiting a slight quantitative shift with the distribution of the normalized local distance to a solid interface aforementioned, both approaches exhibit similar trend with increasing the polydispersity of the medium, supporting the fact that the distribution of $P(L_s/D_{50\%}^w)$ is a good proxy of the pore size distribution.

Interestingly, the distribution is wider for the quasi-monodisperse granular packing (PSD1), suggesting that increasing polydispersity of PSD results in a narrower distribution of the size of the pores composing the porous medium network. This result can be interpreted as a consequence of the fact that for widely bi-dispersed configuration such as PSD4, most of the pores are composed of void spaces between mostly small beads whose diameter is $<D_{50\%}^w$, while for the quasi-monodisperse configuration PSD1, pores are composed of larger void spaces between beads whose diameter is $\approx D_{50\%}^w$. This can also be quantitatively contemplated by estimating for PSD1 (respectively PSD2, PSD3, PSD4) the ratio between the standard deviation of the pore size distribution

<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Various bi-dispersed microstructures with a size ratio ranging from 1.14 (PSD1) to 4 (PSD4). (b) Corresponding distributions of normalised local distance to a solid interface $P(L_s)$ and skeleton pore size.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Distributions of the velocity components for various simulated microstructure: (a) $U_x$ along the mean flow direction and (b) $U_z$ along the transverse direction. (c) Distribution of the velocity magnitude.</figcaption></figure>

and its mean value, $\frac{\sigma_{Poresize}}{\langle Poresize \rangle} = 0.562$ (respectively 0.522, 0.518, 0.464), revealing the narrowing of the relative width of the pore size distribution when increasing the polydispersity.

## 5.2. Narrowing the Velocity Distributions

Figure 10 shows the distributions of velocity components and velocity magnitudes for the various investigated microstructures. Several comments can be made:

1. The distributions of the velocity component $U_x$ along the flow direction is skewed toward positive values with an exponential tail for large values, and reveals the persistence of a small portion of regions of backflow for all the investigated microstructures,

2. In the directions transverse to the flow, the velocity distributions are symmetrical and consistently exhibit exponential tails,

3. The distribution of velocity magnitude are found to be flat at low velocity $P(U) \propto U^0$ for all the investigated microstructure,

4. When increasing polydispersity, the distributions of both the individual components $U_x$, $U_y$, $U_z$ and the velocity magnitude $U = \sqrt{U_x^2 + U_y^2 + U_z^2}$ gets narrower. Specifically, the wider the bi-dispersion size ratio is, the faster the tail of the longitudinal velocity distribution decreases following $P(U_x) \propto e^{-K*(U_x/\langle U \rangle)}$ with values of $K$ ranging from 0.41 for PSD1 up to 0.77 for PSD4. Similarly, the slope of the exponential tails of the transverse velocities distributions increases from 0.60 for PSD1 to 1.25 for PSD4.

This is interpreted as a direct consequence of the pore size distribution which promotes larger velocities through larger pores for the quasi-monodisperse microstructure PSD1, while in contrast, the widely bi-dispersed microstructure PSD4 exhibits a narrower pore size distribution, resulting in a narrower distribution of velocity. Hence, this suggests a correlation between the velocity magnitude and the local distance to a solid interface. Such correlation can also be contemplated through the velocity magnitude density maps of Figure 11: although results exhibit a strong scattering as the distribution of velocity magnitudes for a given local distance to a solid interface is quite large, an overall trend emerges, where larger pore size results in larger velocity magnitudes. Most interestingly, these density maps clearly show that: (a) the correlation between local distance to a solid interface and velocity magnitude is robust among the investigated PSDs, (b) an increase in microstructure polydispersity results both in a decrease in the pore sizes and of the velocity magnitudes compared to a quasi-monodisperse configuration, as moving from PSD1 to PSD4 results in a progressive shrinkage of the density map on both axis.

Note that although counter-intuitive, such result is consistent considering that when increasing polydispersity, the presence of smaller solid beads filling the gaps between larger ones creates more continuous and uniform flow channels, as it involves a larger amount of solid beads, which in turn induce that the streamlines develop through a larger amount of smaller pores. Such structures prohibit strong local velocity variations, leading to a narrower velocity distribution. The aforementioned results highlight the role of the solid microstructure over the flow kinematic and suggest that polydispersity promotes flow velocity homogeneity by restricting large velocity fluctuations.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Density maps of the velocity magnitude as a function of the local distance to a solid interface $L_s$ for the configurations (a) PSD1, (b) PSD2, (c) PSD3, and (d) PSD4, reveal similarities in the density when increasing polydispersity, but also shrinkage on both axis.</figcaption></figure>

5.3. Decreasing the Dispersion Process

Following the methodology introduced in Section 4.2, the dispersion process is investigated among the various microstructures by measuring the spatial variance of a cloud of SPH particles advected within each

<figure><img src="image_14.png" /><figcaption>**Figure 12.** (a) Spatial variances of a cloud of SPH particles advected in the flow for various microstructures. (b) PDF of streamline tortuosity at $t/\tau = 96$ for various microstructures.</figcaption></figure>
microstructure, both in the longitudinal and transverse directions, as shown on Figure 12a. Independently of the considered microstructure, after an initial ballistic regime, the variance always exhibits a transient trend before eventually relaxing to a Fickian regime for $t/\tau \gtrsim 20$. Consistently with the narrower distributions of velocity reported above, the variance for increased polydispersity configurations exhibits similar slopes but with lower quantitative values reflecting a decrease in the effective dispersion coefficient (from Equation 9), in both the longitudinal and transverse directions. Note that this is expected as the size of the dispersion envelope is set by the differences in the trajectories experienced by the fastest fluid elements and the slowest ones within the flow. Consequently, narrower velocity distributions result in a decrease of the dispersion process, both in the longitudinal and the transverse directions as reported on Figure 12a.

## 5.4. Decrease in Streamline Tortuosity

Finally, the distribution of streamlines tortuosity $P(T)$ (obtained following the procedure explained in Section 4.3) is displayed for the various investigated microstructures on Figure 12b. Interestingly, results reveal that tortuosity remains normally distributed among the investigated microstructures and that tortuosity distributions shift to lower values when increasing polydispersity. Although the absolute value of the standard deviation remains almost constant $\sigma_T \approx 0.0271$, the relative standard deviation $\sigma_T/T$ increases as the polydispersity increases. This suggests that flow streamlines are exhibiting stronger meandering discrepancies for polydisperse microstructures but are also simultaneously, on average, more aligned with the longitudinal flow direction. This is consistent with the decrease of the spatial variance magnitude in the transverse direction reported in Section 5.3: in the extreme case of a flow exhibiting an absence of transverse dispersion process, the streamlines tortuosity $T$ would be unity as all streamlines would be aligned with the flow direction.

## 5.5. Porous Medium Permeability

Although simulations are performed in porous medium maintaining constant both the porosity $\epsilon = 0.42$ and the mass median diameter $D_{50\%}^w = 75$ mm, the resulting permeability $K$ of the investigated porous medium displays strong variation, as one can expect from the modification of both the pore size distributions and streamline tortuosity $T$. Following (Schulz et al., 2019) and incorporating the tortuosity in the Kozeny-Carman equation yields:

$$K = \frac{\epsilon^3}{c \langle T \rangle^2 A^2 (1 - \epsilon)^2} \quad (10)$$

where $c = 4.8$ is the Kozeny constant for porous medium made of an assembly of perfect spheres (Ozgumus et al., 2014) and $A$ is the specific surface area for a distribution of spheres. $A$ is defined as the ratio of the fluid-solid interfacial area to the solid volume (Nguyen & Indraratna, 2020), such that for a sphere packing $A = 6 \sum \frac{p_i^w}{d_i}$ where $p_i^w$ is the mass probability of a particle of diameter $d_i$.

Using Equation 10 for the various PSDs with the average tortuosity $\langle T \rangle$ reported above, this yields respectively to permeability $K = 3.9, 2.2, 1.3$ and $0.7$ cm$^2$. Interestingly, although a decrease in tortuosity from PSD1 to PSD4 is expected to result in an increase in permeability, a widely bidispersed microstructure exhibits a larger ratio between small and large spheres, resulting not only in a narrower pore size distribution but also in an increase of the specific surface area $A$, which in turn decreases the permeability $K$. Consequently, although flow simulations where performed imposing a constant flow rate, it's interesting to note that the hydraulic conductivity of the widely bidisperse configuration PSD4 is $\approx 5$ times smaller than that of the quasi monodisperse PSD1 configuration.

# 6. Conclusion

Using the existing SPH approach to simulate fluid flow in various microstructures, this work investigates statistically the flow kinematics and the advection-induced dispersion process in three-dimensional porous medium made of a random pack of spheres. By imposing a constant flow rate resulting in a steady laminar flow, results have been first analyzed in comparison with previous experimental observations regarding velocity distributions, dispersion process and streamline tortuosity for a quasi-monodisperse case. Results are found to be both in
qualitative and quantitative agreement with experiments, highlighting how SPH approach are efficient for simulating such complex configurations and capture critical flow features inherent to three dimensional porous media. The method is then extended to investigate the effect of increased polydispersity over the flow kinematic and dispersion process by simulating various microstructures. Results could be summarized as follows:

1. Distributions of velocity magnitude are found to be always flat at low velocities and exhibit an exponential tail for large velocity magnitudes, in agreement with numerical results in Dentz et al. (2018) and experiments from Souzy et al. (2020). Velocity components in the transverse direction are symmetrically distributed with exponential tails, and the velocity distribution in the longitudinal direction is skewed towards positive values with an exponential tail.

2. Increasing polydispersity for three-dimensional porous medium made of random packs of spheres results in narrower pore size distributions and, thus, narrower velocity distributions.

3. The dispersion process remains qualitatively the same independently of the degree of polydispersity. The variance of the spatial position of an advected cloud of passive tracers follows an initial ballistic regime before eventually relaxing to a Fickian regime after a transient, which is set by the minimal velocity experienced by the considered particles advected within the flow. By narrowing velocity distributions, increasing polydispersity also restricts quantitatively the dispersion process in both longitudinal and transverse directions.

4. When increasing polydispersity, tortuosity remains normally distributed, suggesting that although relative discrepancies in tortuosity increase, the average streamlines are less tortuous in the overfilled polydisperse configurations.

These results highlight how the microstructure complexity in three-dimensional porous medium sets the flow kinematics and transport processes characteristics within such media. The existing SPH approach developped by Domínguez et al. (2022) appears to be a reliable and promising numerical tool for future investigations aiming at exploring porous medium of increased microstructure complexity, such as porous medium made of a random pack of non-spherical beads, the presence of local cracks, or the dynamic modification of the preferential flow streamlines by flowing a particulate suspension within the medium which would result in local clogging of some pores.

Another important perspective of the future advances in solute transport studies will benefit from integrating process-based models with high-resolution geo-hydro-physical tomographic surveys (Busato et al., 2019; Yeh & Liu, 2000). Such approaches have the potential to move beyond the use of effective parameters that obscure the underlying complexity, and instead provide spatially explicit insights into heterogeneity and mixing dynamics in porous media (Su et al., 2020). By linking imaging techniques such as hydraulic tomography and time-lapse cross-hole electrical resistivity tomography with numerical modeling, one may reconstruct spatially variable property fields and their dynamics (e.g., exchange-flow pathways and mixing zones), thereby grounding solute transport interpretations for better resolved structures and processes. These aforementioned investigations would be particularly relevant for the diversity of applications ranging from contaminant migration, filters and catalytic processes, geological storage of carbon dioxide, soil remediation, or drug delivery and nutrient transport in brain or plant tissues.

**Global Research Statement**

This work is in cooperation with Assistant Professor Hao Xiong from Shenzhen University, China, and the authors appreciate his support on the SPH numerical modeling tool used in this research.

**Conflict of Interest**

The authors declare no conflicts of interest relevant to this study.

**Data Availability Statement**

The open-source code DualSPHysics [Software] was used in this work (Crespo et al., 2014), available at https://dual.sphysics.org/downloads/. The experimental data for validation and comparison in this work are available from Souzy et al. (2020) [Dataset].
Acknowledgments

The authors would like to sincerely thank
the three anonymous reviewers for their
constructive comments and insightful
suggestions, which have greatly improved
the quality of this paper. Fan Chen
acknowledges the China Scholarship
Council (CSC) for the support (No.
20220700024). This research was funded,
in whole or in part, by l'Agence Nationale
de la Recherche (ANR), project ANR-23-
CE51-0024-01. For the purpose of open
access, the author has applied a CC-BY
public copyright licence to any Author
Accepted Manuscript (AAM) version
arising from this submission.

References

Adeaga, O., Adegbola, A., & Mudashiru, L. (2023). Simulation of fluid flow in subsurface porous media: An introductory overview (Vol. 7, p. 7). https://doi.org/10.9790/1813-0704026672

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

Anitha, S., Pichumani, M., & Thomas, T. (2022). Physical and mathematical modelling of fluid and heat transport phenomena in porous media. *Engineering Materials, 66*–688. https://doi.org/10.1007/978-3-030-85397-6_21

Anna, P., Borgne, T., Dentz, M., Tartakovsky, A., Bolster, D., & Davy, P. (2013). Flow intermittency, dispersion, and correlated continuous time random walks in porous media. *Physical Review Letters, 110*(18), 184502. https://doi.org/10.1103/PhysRevLett.110.184502

Aristoff, J., Truscott, T., Techet, A., & Bush, J. (2010). The water entry of decelerating spheres. *Physics of Fluids, 22*(3), 032102. https://doi.org/10.1063/1.3309454

Beguin, R., Philippe, P., & Faure, Y.-H. (2013). Pore-scale flow measurements at the interface between a sandy layer and a model porous medium: Application to statistical modeling of contact erosion. *Journal of Hydraulic Engineering, 139*(1), 1–11. https://doi.org/10.1061/(asce)hy.1943-7900.0000641

Berkowitz, B., Cortis, A., Dentz, M., & Scher, H. (2006). Modeling non-Fickian transport in geological formation as a continuous time random walk. *Reviews of Geophysics, 44*(2). https://doi.org/10.1029/2005RG000178

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., Raeini, A., Mostaghimi, P., & Blunt, M. J. (2013). Predictions of non-Fickian solute transport in different classes of porous media using direct simulation on pore-scale images. *Physical Review, 87*(1), 013011. https://doi.org/10.1103/physreve.87.013011

Bouscasse, B., Colagrossi, A., Marrone, S., & Antuono, M. (2013). Nonlinear water wave interaction with floating bodies in SPH. *Journal of Fluids and Structures, 42*, 112–129. https://doi.org/10.1016/j.jfluidstructs.2013.05.010

Busato, L., Boaga, J., Perri, M. T., Majone, B., Bellin, A., & Cassiani, G. (2019). Hydrogeophysical characterization and monitoring of the hyporheic and riparian zones: The vermigliana creek case study. *Science of the Total Environment, 648*, 1105–1120. https://doi.org/10.1016/j.scitotenv.2018.08.179

Canelas, R., Crespo, A., Domínguez, J., Ferreira, R., & Gesteira, M. (2016). SPH-DCDEM model for arbitrary geometries in free surface solid-fluid flows. *Computer Physics Communications, 202*, 131–140. https://doi.org/10.1016/j.cpc.2016.01.006

Carrel, M., Morales, V. L., Dentz, M., Derlon, N., Morgenroth, E., & Holzner, M. (2018). Pore-scale hydrodynamics in a progressively bioclogged three-dimensional porous medium: 3-d particle tracking experiments and stochastic transport modeling. *Water Resources Research, 54*(3), 2183–2198. https://doi.org/10.1002/2017WR021726

Catalano, E., Chareyre, B., & Barthélémy, E. (2014). Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects. *International Journal for Numerical and Analytical Methods in Geomechanics, 38*(1), 51–71. https://doi.org/10.1002/nag.2198

Chen, F., Jiang, S., Xiong, H., Yin, Z.-Y., & Chen, X. (2023). Micro pore analysis of suffusion in filter layer using tri-layer CFD-DEM model. *Computers and Geotechnics, 156*, 105303. https://doi.org/10.1016/j.compgeo.2023.105303

Crespo, A., Domínguez, J., Rogers, B., Gesteira, M., Longshaw, S., Canelas, R., et al. (2014). DualSPHysics: Open-source parallel CFD solver based on Smoothed Particle Hydrodynamics (SPH). *Computer Physics Communications. https://doi.org/10.1016/j.cpc.2014.10.004*

Datta, S. S., Chiang, H., Ramakrishnan, T., & 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

Dentz, M., Borgne, T., Englert, A., & Bijeljic, B. (2011). Mixing, spreading and reaction in heterogeneous media: A brief review. *Journal of Contaminant Hydrology, 120–121*, 1–17. https://doi.org/10.1016/j.jconhyd.2010.05.002

Dentz, M., Borgne, T. L., Lester, D. R., & Barros, a. F. P. J. d. (2016). Mixing in groundwater. In *The Handbook of Groundwater Engineering* (3rd ed.). CRC Press. (Num Pages: 30).

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

Dentz, M., Icardi, M., & Hidalgo, J. J. (2018). Mechanisms of dispersion in a porous medium. *Journal of Fluid Mechanics, 841*, 851–882. (Publisher: Cambridge University Press). https://doi.org/10.1017/jfm.2018.120

Diego, M., & Colagrossi, A. (2009). A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH. *Computer Physics Communications, 180*(6), 861–872. https://doi.org/10.1016/j.cpc.2008.12.004

Domínguez, J. M., Fourtakas, G., Altomare, C., Canelas, R. B., Tafuni, A., García-Feal, O., et al. (2022). Dualspysics: From fluid dynamics to multiphysics problems. *Computational Particle Mechanics, 9*(5), 867–895. https://doi.org/10.1007/s40571-021-00404-2

Feehley, C. E., Zheng, C., & Molz, F. J. (2000). A dual-domain mass transfer approach for modeling solute transport in heterogeneous aquifers: Application to the macrodispersion experiment (made) site. *Water Resources Research, 36*(9), 2501–2515. https://doi.org/10.1029/2000WR900148

Freeze, R. A., & Cherry, J. A. (1979). Groundwater prentice-hall (Vol. 176, pp. 161–177).

Guo, Z., Ma, R., Zhang, Y., & Zheng, C. (2021). Contaminant transport in heterogeneous aquifers: A critical review of mechanisms and numerical methods of non-Fickian dispersion. *Science China Earth Sciences, 64*(8), 1224–1241. https://doi.org/10.1007/s11430-020-9755-y

Holzner, M., Morales, V., Willmann, M., & Dentz, M. (2015). Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow. *Physical Review, 92*(1), 013015. https://doi.org/10.1103/PhysRevE.92.013015

Kang, P., Anna, P., Pereira Nunes, J. P., Bijeljic, B., Blunt, M., & Juanes, R. (2014). Pore-scale intermittent velocity structure underpinning anomalous transport through 3D porous media. *Geophysical Research Letters, 41*(17), 6184–6190. https://doi.org/10.1002/2014GL061475

Knight, C., O'Sullivan, C., van Wachem, B., & Dini, D. (2020). Computing drag and interactions between fluid and polydisperse particles in saturated granular materials. *Computers and Geotechnics, 117*, 103210. https://doi.org/10.1016/j.compgeo.2019.103210

Koltermann, C. E., & Gorelick, S. M. (1996). Heterogeneity in sedimentary deposits: A review of structure-imitating, process-imitating, and descriptive approaches. *Water Resources Research, 32*(9), 2617–2658. https://doi.org/10.1029/96WR00025

Le Borgne, T., Dentz, M., & Villermaux, E. (2015). The lamellar description of mixing in porous media. *Journal of Fluid Mechanics, 770*, 458–498. https://doi.org/10.1017/jfm.2015.117

Levy, M., & Berkowitz, B. (2003). Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. *Journal of Contaminant Hydrology, 64*(3), 203–226. https://doi.org/10.1016/S0169-7722(02)00204-8

Li, S., Peng, C., Wu, W., Wang, S., Chen, X., Chen, J., et al. (2020). Role of baffle shape on debris flow impact in step-pool channel: An SPH study. *Landslides, 17*(9), 2099–2111. https://doi.org/10.1007/s10346-020-01410-w
Lind, S., Stansby, P., & Rogers, B. (2015). Incompressible–compressible flows with a transient discontinuous interface using smoothed particle hydrodynamics (SPH). *Journal of Computational Physics*, **309**, 129–147. https://doi.org/10.1016/j.jcp.2015.12.005

Lista, L. (2023). Probability density functions. *Lecture Notes in Physics*, 37–73. https://doi.org/10.1007/978-3-031-19934-9_3

Lyu, H.-G., Deng, R., Sun, P.-N., & Miao, J.-M. (2021). Study on the wedge penetrating fluid interfaces characterized by different density-ratios: Numerical investigations with a multi-phase SPH model. *Ocean Engineering*, **237**, 109538. https://doi.org/10.1016/j.oceaneng.2021.109538

Manenti, S., Wang, D., Domínguez, J. M., Li, S., Amicarelli, A., & Albano, R. (2019). SPH modeling of water-related natural hazards. *Water*, **11**(9), 1875. https://doi.org/10.3390/w11091875

Meunier, P., & Villermaux, E. (2010). The diffusive strip method for scalar mixing in two dimensions. *Journal of Fluid Mechanics*, **662**, 134–172. https://doi.org/10.1017/s0022112010003162

Mokos, A., Rogers, B., & Stansby, P. (2016). A multi-phase particle shifting algorithm for SPH simulations of violent hydrodynamics with a large number of particles. *Journal of Hydraulic Research*, **55**(2), 1–20. https://doi.org/10.1080/00221686.2016.1212944

Monaghan, J. J. (1994). Simulating free surface flows with SPH. *Journal of Computational Physics*, **110**(2), 399–406. https://doi.org/10.1006/jcph.1994.1034

Morales, V. L., Dentz, M., Willmann, M., & Holzner, M. (2017). Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory. *Geophysical Research Letters*, **44**(18), 9361–9371. https://doi.org/10.1002/2017GL074326

Morris, J. P., Fox, P. J., & Zhu, Y. (1997). Modeling low Reynolds number incompressible flows using sph. *Journal of Computational Physics*, **136**(1), 214–226. https://doi.org/10.1006/jcph.1997.5776

Moukalled, F., Mangani, L., & Darwish, M. (2015). *The finite volume method in computational fluid dynamics: An advanced introduction with openfoam® and matlab* (1st ed.). Springer. https://doi.org/10.1007/978-3-319-16874-6

Nguyen, T. T., & Indraratna, B. (2020). The role of particle shape on hydraulic conductivity of granular soils captured through Kozeny–carman approach. *Geotechnique Letters*, **10**(3), 398–403. https://doi.org/10.1680/jgele.20.00032

Ottino, J., Ranz, W. E., & Macosko, C. W. (1979). A lamellar model for analysis of liquid-liquid mixing. *Chemical Engineering Science*, **34**(6), 877–890. https://doi.org/10.1016/0009-2509(79)85145-3

Ozgumus, T., Mobedi, M., & Ozkol, U. (2014). Determination of Kozeny constant based on porosity and pore to throat size ratio in porous medium with rectangular rods. *Engineering Applications of Computational Fluid Mechanics*, **8**(2), 308–318. https://doi.org/10.1080/19942060.2014.11015516

Peyrounette, M., Davit, Y., Quintard, M., & Lorthois, S. (2018). Multiscale modelling of blood flow in cerebral microcirculation: Details at capillary scale control accuracy at the level of the cortex. *PLoS One*, **13**(1), e0189474. https://doi.org/10.1371/journal.pone.0189474

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

Ranz, W. E. (1979). Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows. *AIChE Journal*, **25**(1), 41–47. https://doi.org/10.1002/aic.690250105

Rosenbom, A. E., Binning, P. J., Aamand, J., Dechesne, A., Smets, B. F., & Johnsen, A. R. (2014). Does microbial centimeter-scale heterogeneity impact MCPA degradation in and leaching from a loamy agricultural soil? *Science of the Total Environment*, **472**, 90–98. https://doi.org/10.1016/j.scitotenv.2013.11.009

Schulz, R., Ray, N., Zech, S., Rupp, A., & Knabner, P. (2019). Beyond Kozeny–carman: Predicting the permeability in porous media. *Transport in Porous Media*, **130**(2), 487–512. https://doi.org/10.1007/s11242-019-01321-y

Shao, X., Zhou, Z., & Wu, W. (2012). Particle-based simulation of bubbles in water–solid interaction. *Computer Animations and Virtual Worlds*, **23**(5), 477–487. https://doi.org/10.1002/cav.438

Skibitzke, H., Robinson, G., & Commission, U. A. E. (1963). *Dispersion in ground water flowing through heterogeneous materials*. U.S. Government Printing Office.

Souzy, M., Lhuissier, H., Méheust, Y., Borgne, T., & 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

Souzy, M., Zaier, I., Lhuissier, H., Borgne, T., & Metzger, B. (2017). Mixing lamellae in a shear flow. *Journal of Fluid Mechanics Rapids*, **838**, R3. https://doi.org/10.1017/jfm.2017.916

Stephen, S., & Simpson, E. (1987). Laboratory evidence of the scale effect in dispersion of solutes in porous media. *Water Resources Research*, **23**(8), 1667–1673. https://doi.org/10.1029/WR023i008p01667

Su, X., Jim Yeh, T.-C., Shu, L., Li, K., Brusseau, M. L., Wang, W., et al. (2020). Scale issues and the effects of heterogeneity on the dune-induced hyporheic mixing. *Journal of Hydrology*, **590**, 125429. https://doi.org/10.1016/j.jhydrol.2020.125429

Tanguy, L. B., Dentz, M., Bolster, D., Carrera, J., de Dreuzy, J. R., & Davy, P. (2010). Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media. *Advances in Water Resources*, **33**(12), 1468–1475. https://doi.org/10.1016/j.adwatres.2010.08.006

Villermaux, E. (2012). On dissipation in stirred mixtures. *Advances in Applied Mechanics*, **45**, 91–107. https://doi.org/10.1016/b978-0-12-380876-9.00003-3

Villermaux, E. (2019). Mixing versus stirring. *Annual Review of Fluid Mechanics*, **51**(1), 245–273. https://doi.org/10.1146/annurev-fluid-010518-040306

Wang, C., Han, L., Zhang, Y., Jiang, A., Wang, J., & Niu, X. (2023). Effects of physical properties and environmental conditions on the natural dispersion of oil. *Journal of Marine Science and Engineering*, **12**(1), 47. https://doi.org/10.3390/jmse12010047

Xiong, H., Hao, M., Zhao, D., Qiu, Y., & Chen, X. (2023). Study of the dynamics of water-enriched debris flow and its impact on slit-type barriers by a modified SPH–DEM coupling approach. *Acta Geotechnica*, **19**(2), 1019–1045. https://doi.org/10.1007/s11440-023-02106-w

Xu, W.-j., Yao, Z.-G., Luo, Y.-T., & Dong, X. (2020). Study on landslide-induced wave disasters using a 3D coupled SPH-DEM method. *Bulletin of Engineering Geology and the Environment*, **79**(1), 467–483. https://doi.org/10.1007/s10064-019-01558-3

Yeh, T.-C. J., & Liu, S. (2000). Hydraulic tomography: Development of a new aquifer test method. *Water Resources Research*, **36**(8), 2095–2105. https://doi.org/10.1029/2000WR090014

Yueting, L., Guo, N., Yang, Z., & Zhao, J. (2023). A fully resolved smoothed particle hydrodynamics-discrete element method study of the rheology of suspensions: The role of inertia and grain shape. *Physics of Fluids*, **35**(8), 083325. https://doi.org/10.1063/5.0161344

Zhan, L., Peng, C., Zhang, B., & Wu, W. (2019). A stabilized TL–WC SPH approach with GPU acceleration for three-dimensional fluid–structure interaction. *Journal of Fluids and Structures*, **86**, 329–353. https://doi.org/10.1016/j.jfluidstructs.2019.02.002

Zhang, C., Wang, J., Rezavand, M., Wu, D., & Hu, X. (2021). An integrative smoothed particle hydrodynamics method for modeling cardiac function. *Computer Methods in Applied Mechanics and Engineering*, **381**, 113847. https://doi.org/10.1016/j.cma.2021.113847