Modeling Longitudinal Dispersion in Variable Porosity Porous Media: Control of Velocity Distribution and Microstructures

Hydrodynamic dispersion process in relation with the geometrical properties of the porous media are studied in two sets of 6 porous media samples of porosity θ ranging from 0.1 to 0.25. These two sets of samples display distinctly different evolutions of the microstructures with porosity but share the same permeability trend with porosity. The methodology combines three approaches. First, numerical experiments are performed to measure pre-asymptotic to asymptotic dispersion from diffusion-controlled to advection-controlled regime using Time-Domain Random Walk solute transport simulations. Second, a porosity-equivalent network of bonds is extracted in order to measure the geometrical properties of the samples. Third, the results of the direct numerical simulations are interpreted as a Continuous Time Random Walk (CTRW) process controlled by the flow speed distribution and correlation. These complementary modeling approaches allow evaluating the relation between the parameters of the conceptual transport process embedded in the CTRW model, the flow field properties and the pore-scale geometrical properties. The results of the direct numerical simulations for all the 12 samples show the same scaling properties of the mean flow distribution, the first passage time distribution and the asymptotic dispersion vs. the Péclet number than those predicted by the CTRW model. It allows predicting the asymptotic dispersion coefficient D* from Pe = 1 to the largest values of Pe expected for laminar flow in natural environments (Pe≈ 4,000). D*∝Pe2−α for Pe≥Pecrit, where α can be inferred from the Eulerian flow distribution and Pecrit depends on porosity. The Eulerian flow distribution is controlled by the distribution of fractions of fluid flowing at each of the pore network nodes and thus is determined mainly by the distribution of the throat radius and the coordination number. The later scales with the number of throats per unit volume independently on the porosity. The asymptotic dispersion coefficient D* decreases when porosity increases for all Péclet values larger than 1 due to the increase with porosity of both α and the flow speed decorrelation length.

Hydrodynamic dispersion process in relation with the geometrical properties of the porous media are studied in two sets of 6 porous media samples of porosity θ ranging from 0.1 to 0.25. These two sets of samples display distinctly different evolutions of the microstructures with porosity but share the same permeability trend with porosity. The methodology combines three approaches. First, numerical experiments are performed to measure pre-asymptotic to asymptotic dispersion from diffusion-controlled to advectioncontrolled regime using Time-Domain Random Walk solute transport simulations. Second, a porosity-equivalent network of bonds is extracted in order to measure the geometrical properties of the samples. Third, the results of the direct numerical simulations are interpreted as a Continuous Time Random Walk (CTRW) process controlled by the flow speed distribution and correlation. These complementary modeling approaches allow evaluating the relation between the parameters of the conceptual transport process embedded in the CTRW model, the flow field properties and the porescale geometrical properties. The results of the direct numerical simulations for all the 12 samples show the same scaling properties of the mean flow distribution, the first passage time distribution and the asymptotic dispersion vs. the Péclet number than those predicted by the CTRW model. It allows predicting the asymptotic dispersion coefficient D * from Pe = 1 to the largest values of Pe expected for laminar flow in natural environments (Pe ≈ 4,000). D * ∝ Pe 2−α for Pe ≥ Pe crit , where α can be inferred from the Eulerian flow distribution and Pe crit depends on porosity. The Eulerian flow distribution is controlled by the distribution of fractions of fluid flowing at each of the pore network nodes and thus is determined mainly by the distribution of the throat radius and the coordination number. The later scales with the number of throats per unit volume independently on the porosity. The asymptotic dispersion coefficient D * decreases when porosity increases for all Péclet values larger than 1 due to the increase with porosity of both α and the flow speed decorrelation length.

INTRODUCTION
Modeling transport of solute in porous media is a prerequisite for many environmental and engineering applications, ranging from aquifers contaminant risk assessment to industrial reactors, filters and batteries design. The solutes can be pollutants, reactants and products involved in solute-solute or solute-mineral reactions, but also (bio-)nanoparticles or nutriments involved in the growth of bio-mass. The mechanism under consideration is the spatial dispersion which leads to the spreading and the mixing of dissolved chemicals, thus controlling the potential reactions in the flowing fluid and between the fluid and the porous media (Bear, 1972;Brenner and Edwards, 1993;Dentz et al., 2011). The dispersion process has been, and still is, a largely studied topic in the field of geosciences because rocks at depth are, as a general rule, porous media saturated with fluid(s) that move due to natural or artificial pressure gradients, and display a large spectrum of heterogeneities. In all these domains, reliable predictive models that can be parameterized by direct measurements are necessary, for example, to monitor and assess risks linked to the use of underground water resources, or in the course of industrial operations, such as hydrocarbon exploitation and CO 2 or underground nuclear waste storage.
Hydrodynamic dispersion is the macroscopic result of the mass transfers by diffusion and advection that occurs at the pore scale (Whitaker, 1967;Sahimi, 2011;De Anna et al., 2013). Together, diffusion and advection of solute produce a large spectrum of dispersion features because (natural) porous media display complex structures inducing a large diversity of velocity fields, and thus distinctly different speed distributions and spatial correlations. Probably the most obvious behavior that illustrates the complexity of dispersion mechanisms in porous media is the variably-lasting pre-asymptotic dispersion regime that cannot be modeled by a single Fickian dispersion coefficient. Preasymptotic, or non-Fickian, dispersion is commonly observed in laboratory experiments (Moroni and Cushman, 2001;Levy and Berkowitz, 2003;Seymour et al., 2004;Morales et al., 2017;Carrel et al., 2018;Souzy et al., 2020), and numerical simulations (Bijeljic et al., , 2013De Anna et al., 2013;Icardi et al., 2014;Kang et al., 2014;Li et al., 2018;Puyguiraud et al., 2019c). It is characterized by heavy-tailed arrival time distributions f t (t) and super-diffusive growth of the longitudinal displacement variance σ 2 (t). For a given porous medium, the duration of the non-Fickian regime is controlled by solute particles that move the slowest, which emphasizes the determinant role of both the regions where the velocity is low and the tortuosity of the flow paths. Asymptotically, dispersion converges toward Fickian behavior, characterized by the constant longitudinal dispersion coefficient D * (Bear, 1972;Brenner and Edwards, 1993).
Evaluating the longitudinal asymptotic dispersion coefficient D * is a fundamental issue, because most operational modeling tools have been constructed around the Fickian advection-dispersion equation that reads for transport in the direction of the mean flow, here the z-direction (Bear, 1972): where c is the solute concentration, θ is the connected porosity, u z = θ v z denotes Darcy's velocity, with v z being the mean pore velocity. Many experimental studies and mathematical developments on dispersion using mainly simple porous media have been performed since the pioneering works of Danckwerts (1953). The reader will find an exhaustive review of the different results and models of both longitudinal and transverse dispersion in Delgado (2006). A main well-observed feature of longitudinal dispersion D * is its non-linear increase with the mean flow velocity. It is recognized since the pioneering works of Saffman (1959) and then Bear (1972). It is generally expressed in terms of D * /d m vs. the Péclet number Pe = v e ℓ/d m , where ℓ is a characteristic length, d m is the molecular diffusion coefficient and v e is the mean Eulerian flow speed (v e = v 2 x + v 2 y + v 2 z , with v i denoting the flow velocity component i, see section 2.2). Simulations in networks of constant velocity tubes (Sahimi and Imdakm, 1988) of radius r following distributions such as P(r) ∝ re −r 2 (Chatzis and Dullien, 1985) indicated a relation of the form with β = 1.2 ∓ 0.1 (Sahimi, 2011), while for instance β = 2 in a single tube (Taylor, 1953). For infinite Pe, experimental particle tracking results (e.g, Souzy et al., 2020) give the relation D * /d m ≈ Pe, where the characteristic length ℓ is of the order of the pore length. However, it is worth noticing that in Souzy et al. (2020)'s experiments the lowest velocities cannot be measured because they use finite-size particles that cannot access to the vicinity of the solid. Interestingly, the behavior (Equation 2) with β ≃ 1.2 was cited in numerous studies concerning bead-packs and homogeneous sand-packs for intermediate Péclet numbers (Pfannkuch, 1963;Han et al., 1985;Sahimi et al., 1986;Seymour and Callaghan, 1997;Bijeljic et al., 2004). For instance, particle tracking simulations in pore-networks reported in Bijeljic and Blunt (2006) gave β = 1.2, for Pe < 400 and β = 1, for Pe > 400. Conversely, similar numerical simulations (using random walk particle tracking) performed by Puyguiraud et al. (2021) using digitized images of consolidated sandstone, gave a value of β = 1.65 for 10 ≤ Pe ≤ 10 5 . The few experimental data on rocks (obviously more heterogeneous than bead-packs) displayed a broader range of behaviors; for example Kinzel and Hill (1989) reported 1.30 ≤ β ≤ 1.33. However, it is worth noticing that evaluating dispersion in rocks, for a large range of Pe values, either at laboratory or field scale from tracer tests is challenging. For instance, controlling the boundary conditions and verifying that the tracer is conservative are some of the known issues that may introduce errors in the estimation. Yet, the main issue is probably linked to the fact that, by definition, the experimental results are interpreted using the Fickian model, whereas it is difficult to prove that dispersion is asymptotic without being able to measure the tracer breakthrough curves over several orders of magnitude in order to capture the low speed fraction of the solute transport (Gouze et al., 2008). We will show in section 3.4 that measuring asymptotic dispersion for large values of Pe in natural porous media is in fact virtually impossible using cm-or even meter-scale experiments. While measuring dispersion experimentally is burdensome, modeling approaches are now mature to perform numerical experiments. Direct numerical simulations (DNS) are unique tools for investigating both the pre-asymptotic and the asymptotic behavior in a common frame. They can be used to accurately measure D * , but also to study the mechanisms that produce dispersion in relation with the measurable (average) properties of the material, and to test upscaling theories. Recent works Blunt, 2006, 2007;De Anna et al., 2013;Puyguiraud et al., 2020Puyguiraud et al., , 2021 showed that hydrodynamic transport in porous media can be adequately conceptualized and modeled by a continuous time random walk (CTRW) that models streamwise transport through particle transitions over fixed spatial distance with a transition time given by the local flow speed and diffusion. The spatial distance at which particles speed changes corresponds to the decorrelation distance ℓ c of the mean flow speed. The CTRW integrates in a statistical framework parameters that are similar to the classical representation of porous media as a network of throats and pores. As such one can be tempted to investigate how ℓ c , which is a major ingredient of the CTRW model, is related to the topological and geometrical properties of the real 3-dimensional pore network. Moreover, the CTRW model predicts that asymptotic dispersion is controlled by the dispersion evolution during the pre-asymptotic regime which itself is controlled by the flow speed distribution. How the later is related to the properties of the pore network is a further issue that requires investigation.
The main objective of the present study, is to investigate the relation between the longitudinal dispersion D * (and its evolution with the mean flow rate) and the porous media microstructural properties in the frame of the theory proposed by Puyguiraud et al. (2021) which gives a generalized explanation of longitudinal dispersion (from pre-asymptotic to asymptotic regimes) and a formal relation between dispersion and the properties of the flow field (velocity distribution, velocity spatial decorrelation and flow path tortuosity).
The core of this study is a set of about 150 numerical experiments designed to measure pre-asymptotic to asymptotic dispersion from diffusion-controlled to advection-controlled regime in 12 sandstone-like samples of porosity ranging from 10 to 25%. For that, one first computes the steady-state Stokes flow field from which, the flow speed distribution and the decorrelation distance as well as advective tortuosity are derived. Then, the direct numerical simulation (DNS) of solute transport at pore scale, involving diffusion and advection, are performed using Time-Domain Random Walk (TDRW). The dispersion mechanisms are characterized from the time-resolved particles displacement variance and the first passage time distribution (FPT) given as outputs of the TDRW simulations. In parallel, the geometrical properties of the porous samples are evaluated from the computation of the bonds network model (BNM) for each of the samples, that is obtained from the medial axis transform, or squeletonization, of the connected porosity. This gives us the unique opportunity to characterize the topology of the connected porosity including the number of throats (bonds) and pores (network nodes) and the coordination (number of throats per pore), as well as the throat radius and length. Then, the results of the direct numerical simulations are analyzed in the light of the CTRW theory proposed by Puyguiraud et al. (2021) which provides quantitative links between the tail behaviors of the FPT distribution f t (t), the distribution of flow speeds v e , the particles displacement variance σ 2 (t) and the asymptotic dispersion coefficient D * scaling with the Pe value.
The methodology, including the conceptual and numerical tools used in this study are detailed in section 2. The geometrical and topological characteristics of the samples and the flow field properties are presented in section 3. The results of the direct numerical simulation of solute transport and the calculation of the dispersion coefficient for a large range of values of Pe are discussed in section 3.4. The conclusions of this study are exposed in section 4.

Porous Media Samples
The porous media are binary images made of 480 3 regular voxels (cubes) that are either void or solid. The first set of 6 samples, noted FSxx, were xx is replaced by the porosity value expressed in percent (ex: FS13 for the sample with θ = 0.13) was downloaded from the Digital Rocks Portal (Berg, 2016a). They were generated with the commercial software e-Core following a methodology described in Oren (2002), in order to mimic Fontainebleau sandstone at different porosity (Berg, 2016b). The op. cit. author indicated that they use identically parameterized silica grain sedimentation and compaction processes typical for Fontainebleau sandstones, the different porosity values (0.10, 0.13, 0.15, 0.21, and 0.25) being obtained by varying the amount of silica cement. As such, this process mimics the progressive diagenetic cementation by silica precipitation (from FS25 to FS10) of an initially poorly cemented sandstone. Conversely, we made the second set of samples by step-by-step homogeneous erosion of the solid phase starting from FS10. By removing 1 to 6 layers of solid at the solid-void interface we obtain 6 samples, denoted FSDxx of porosity 0.12, 0.15, 0.17, 0.20, 0.23, and 0.25. This process mimics homogeneous dissolution of the silica material. The top panel in Figure 1 displays the three-dimensional structure of the lowest porosity sample FS10, and the highest porosity samples FS25 and FSD25. It can be qualitatively appraised that the cement precipitation model used to construct FS25 increases the number of pores compared FS10, while the pore size is kept roughly similar. In contrast, the dissolution process producing FS25 from FS10 acts as increasing strongly the pore size, while the number of pores remains roughly unchanged. This set of sample is viewed as ideal for investigating dispersion of end-members of natural sandstones.
The top row displays the void space. The lines in the bottom panel show particle paths, the color scheme indicating the particle speed from white (u/ u ≤ 7 × 10 −4 ) to dark blue (u/ u = 10).

Flow
Flow simulations are performed on the three-dimensional binary images. The mesh used for solving the flow is obtained by dividing each of the image voxels by 2 in each of the directions so that 1 voxel of the raw image is represented by 8 cubic cells of size x = y = z = 2.85 × 10 −6 m. This procedure is applied for improving the resolution of the flow field in the smallest throats (Gjetvaj et al., 2015). The resulting discretization for the regular grid consists of 960 3 cubic cells. We are considering steady-state flow of an non-compressible Newtonian fluid at low Reynolds number so that the pore-scale flow velocity v(x) is given by the Stokes equation where p(x) is the fluid pressure. Stokes flow is solved using the finite volume SIMPLE (SemiImplicit Method for Pressure Linked Equations algorithm) scheme implemented in the SIMPLEFOAM solver of the OpenFOAM platform (Weller et al., 1998). Twenty layers are added at the inlet and outlet in order to minimize boundary effects (Guibert et al., 2016). The main flow direction is considered in the z-direction all over this study. We prescribe (1) a macroscopic pressure gradient ∇ * p between the inlet (z = 0) and the outlet (z = Lz) boundary conditions such that the Reynolds number Re is smaller than 10 −6 , i.e., laminar flow and (2) no-slip conditions at the voidsolid interfaces and at the remaining boundaries of the sample. After convergence, that is, once the normalized residual of the pressure and velocity components is below 10 −5 between two consecutive steps, we extract the components of the velocity at the voxel interfaces (v x , v y , v z ). The results of the flow simulations allow us to extract the three properties that control dispersion according to Puyguiraud et al. (2019b): (1) the Eulerian speed distribution p e (v) (2) the decorrelation distance ℓ c and (3) the advective tortuosity χ a . These fundamental flow properties are respectively displayed in Figures 5-7, and discussed in section 3.

Solute Transport
Pore-scale hydrodynamic transport is classically modeled by the advection-diffusion equation where c(x, t) is the solute concentration at position x and time t, d m is the molecular diffusion coefficient which is set equal to d m = 10 −9 m 2 /s, and v(x) is the flow velocity at position x which is obtained by solving the Stokes problem (see section 2.2). Here we use the time domain random walk (TDRW) method that is based on a finite volume discretization of Equation (4) (Delay et al., 2005). A detailed description of the TDRW method, its derivation and implementation using voxelized binary images can be found in Dentz (2012) and Russian et al. (2016); the main features of the method are given below. A study of the performance and accuracy of the TDRW method for a large range of values of the Péclet number can be found in Gouze et al. (2021). The domain discretization used for transport is that used for computing the flow, i.e., 960 3 cubic voxels .
The TDRW method is a grid-based method that models the displacement of particles in space and time according to the master equation that results from a finite volume discretization of the advection-diffusion equation. The ensemble average of the particle displacement gives the solution of the transport equation. A particle transition corresponds to a single transition of a constant length ξ = x from the center of a voxel j to the center of one of the 6 face-neighboring voxels i. The direction and the transition duration are random variables ruled by the local values of the fluid velocity at the voxel interface embedded into the local coefficients b ij (Russian et al., 2016).
where v ij is the velocity component of v j in the direction of voxel as a convention. The velocity at the solid-void interface is zero and d m = 0 if voxel i is a solid voxel. The recursive relations that describe the random walk from position x j to position x i of a given particle transition n are The probability p ij for a transition of length ξ from voxel j to voxel i is where [jk] denotes the summation over the nearest neighbors of voxel j. The transition time τ j is independent on the transition direction and is exponentially distributed ψ τ j (t) = τ j exp(−t/τ j ) with τ j the mean transition time from voxel j; Frontiers in Water | www.frontiersin.org The algorithm consists in computing once the probability p ij (Equation 7) and the mean transition time τ j (Equation 8) for each of the voxels belonging to the pore space and then solving the random walk (Equation 6) in which the direction for each particle transition is drawn from the p ij vector and the transition time is drawn from the exponential distribution of mean τ j .

Simulation Setup
For each sample, we performed simulations for different values of the Péclet number. The Péclet number is defined as Pe = v e λ/d m where λ is the mean throat length that is displayed for the 12 samples in Figure 2 and ranges from 6.5 × 10 −5 m to 8.8 × 10 −5 m. The different flow fields used for the TDRW simulations at different Péclet numbers are obtained by multiplying the raw flow field resulting from the Stokes simulation by a constant.
A pulse of constant concentration at the sample inlet (z = 0) is applied at t = 0 by locating particles in a flux weighted injection mode. Note that the pulse is formally an exponential distribution function of characteristic time τ j | z=0 whose mean value is negligible compared to the mean time required for the particles to move through the sample (Russian et al., 2016). Flux weighted injection means that the number of particles injected at a location is proportional to the local velocity. This corresponds to a constant concentration Dirichlet boundary condition. Particles that reach the sample outlet with a speed v out are reinjected randomly at the inlet plane at a position x satisfying the condition |v The distribution (PDF) of first passage times at a given distance Z from the injection location, that denotes the solute breakthrough curve (BTC) usually measured in laboratory or field tracer tests, is noted f t (t) (Figure 8). The apparent longitudinal dispersion coefficient D(t) is evaluated from the displacement variance σ 2 z (t) of the particles (Fischer, 1966): with σ 2 z (t) = (z(t) − z ) 2 − z(t) − z 2 . The asymptotic longitudinal dispersion coefficient D * = σ 2 z (t)/2t is obtained for t > t * , where t * is the time required for all the particles to sample the entire heterogeneity, i.e., when σ 2 z (t) ∼ t (see for example Figure 10).

The Equivalent Bond Network Model
We compute the bond network model (BNM) for each of the FS and FSD samples in order to extract the geometrical and topological characteristics of the connected porosity. The methodology to obtain the network representation of the connected porosity of the sample includes two main steps. The first one is the extraction of the void space skeleton which is the one-dimensional continuous object centrally located (and spatially referenced) inside the pore space. The skeleton can be computed using different approaches; here we used a thinning algorithm inspirited from the works of Lee et al. (1994) that provides the local medial-axis. The coordinate of the skeleton is known with a spatial resolution equivalent to that of the original 3D-image and associated with the local hydraulic radius r l normal to the local medial axis that is evaluated using a pondered 45 degree multi-ray method. Thus, the skeleton keeps the relevant geometrical and topological features of the pore space (Siddiqi and Pizer, 2008). The second step consists in transforming the skeleton into a network of bonds and nodes that connect three or more bonds. This yields an irregular lattice. The length λ of a given bond is the sum of the length of the skeleton components used to built this bond, so that the local tortuosity of the skeleton is embedded into λ. For each bond, the radius r h is obtained from the harmonic means (noted H ) of the local conductance, so that r h = ( r l H ) 1/4 . The algorithm is nonparametric; there is no assumption on any of the characteristics of the obtained lattice. Puyguiraud et al. (2021) proposed a continuous time random walk (CTRW) model that describes transport through particle transitions over the length ℓ c with a transition time that is given by the local flow speed and diffusion. The central assumption of this model is that transition times at subsequent CTRW steps are independent identically distributed random variables. Furthermore, it is assumed that particles move at the mean pore velocity, that is, it is assumed that during a transition particles are able to diffusively sample the velocities across pore conducts. The scale ℓ c is set equal to the decorrelation distance of particle speeds so that subsequent particle speeds can be considered statistically independent. The distribution of the Eulerian mean flow speeds p m (v) is obtained from the Eulerian speed PDF as p m (v) = −2v dp e (2v) dv .

Upscaled CTRW Model
As particles move at equidistant spatial steps, they sample flow speeds in a flux-weighted manner. This is due to the fact that particles are distributed at pore intersections according to the relative downstream fluxes. Thus, the distribution p v (v) of subsequent particle speeds are related to the distribution of Eulerian mean flow speeds through flux-weighting as .
At each turning point of the CTRW, particles are assigned a random speed from p v (v). The particle transition time distribution ψ(t) reflects both advection and diffusion. It is cutoff at times larger than τ D = ℓ 2 c /d m , the diffusion time over the decorrelation distance. For times small compared to the cut-off time, ψ(t) can be approximated by At times larger than τ D it is cut-off exponentially fast. The flow speed distribution is at the center of the transport process. In porous media, such as rocks, the mean flow speed can often be approximated by a Gamma-type distribution Puyguiraud et al., 2019b;Souzy et al., 2020) and displays a power-law scaling p e (v) ∼ v α−1 for v < v m . For sphere packs and simple structures such as sand-pack the linear flow profile close to the grains (due to the no-slip boundary condition) implies that p e (v) is flat at low velocities, so that α ≃ 1 . In more heterogeneous porous media, other values of α are expected. For example, Puyguiraud et al. (2021) found α ≈ 0.35 for a Berea sandstone sample. For such Gamma-type distributions, p e (v) ∼ v α−1 at small flow speeds, ψ(t) behaves for high Péclet numbers as ψ(t) ∼ t −2−α before the exponential cut-off at times larger than τ D . The tortuosity of particle trajectories in this framework is given by the ratio of the mean asymptotic particle speed ℓ c / τ ≡ v e (where τ denotes the particle mean travel time) and the mean streamwise flow velocity v z . Furthermore, for this type of flow speed distributions, the CTRW approach predicts some further interesting scaling laws that can be verified from direct numerical simulations. The behavior of particle breakthrough curves f (t, Z) at a control plane located at the streamwise location Z is analogous to the behavior of ψ(t). They show a powerlaw dependence as f (t, Z) ∼ t −2−α if Z/v z ≪ τ D (i.e., the peak time is much smaller than the cut-off time), and exponential decay for times larger than the cut-off time τ D . The predicted dependence of the asymptotic longitudinal dispersion coefficients on the Péclet number is for Pe ≫ 1 for 0 < α < 1 and for α = 1, see also Saffman (1959) and Koch and Brady (1985).
To sum-up, this upscaled model, constructed on the representation of the hydrodynamic transport as a CTRW process in a 1-dimensional network of bonds, is fully constrained, for any values of Pe > 1 by the knowledge of the distribution of Eulerian flow speeds p e (v) and the decorrelation distance ℓ c of particle speeds.
From here one can recognize on one hand the complementarity of the BNM and the DNS to explore the relation between the dispersion and the pore network characteristics, and on the other hand the conceptual framework that links the CTRW model and the BNM representation of the porous medium. This emphasizes the possibility of (1) relating the distribution of the Eulerian flow speed to the large scale transport behavior and (2) characterizing dispersion for different porous media based on the knowledge of the flow speed distribution. Indeed, the BNM gives us the information on the real topology of the pore network as well as the distribution and the average of bond properties (radius and length), while the DNS provides the information on the flow field (speed distribution and decorrelation distance as well as the advective tortuosity).

PORE NETWORK PROPERTIES, FLOW FIELDS, AND DISPERSION
The top row in Figure 1 illustrates the 3-dimensional structure of sample FS10 (θ = 0.1) and of both FS25 and FSD25 sharing the same porosity θ = 0.25. The bottom row in Figure 1 displays flow lines (and the local velocity) within the connected porosity for these three samples and gives a qualitative appraisal of the dissimilarities between the lowest porosity and the highest porosity samples on one hand, and on the other hand those occurring between the highest porosity sample of each of the two sets in relation with the pore network structures. In the following we will quantify these differences and their implications on dispersion.

Connected Porosity Geometrical Properties Retrieved From the BNM
As explained in section 2.1, we computed the Bond Network Model (BNM) for each of the 12 samples, in order to evaluated the topology and the geometry of the connected porosity and specifically how these characteristics change with the sample porosity for the FS and the FSD sets of samples. The main properties vs. porosity are summarized in Figures 2, 3. The topology of the connected porosity is characterized by the number of throats (network bonds) and pores (network nodes) per volume of rock (here the reference is the sample volume) as well as the coordination number κ that denotes the mean number of throats connected to a given pore. The bonds are characterized by the mean of the radius r h and length λ and by the radius r h distribution displayed in Figure 4.
For the FS set, decreasing porosity from the highest to the lowest porosity values is obtained by allocating increasing amounts of cement into localized clusters that acts as increasingly closing connections and thus decreasing the number of pores and throats and the coordination number. The fixed distribution of the cement clusters determines the length of the bonds  independently of the porosity (λ ≈ 65µm), but volume conservation imposes that the hydraulic radius r h increases with porosity. The distribution of r h / r h is wide, decreases almost monotonically from small to high r h and does not depends on porosity.
For the FSD set, increasing porosity from the lowest to the highest is obtained by homogeneous erosion of the solid phase, i.e., both the grains and the cement. The number of pores and throats as well as κ first decreases for θ ≤ 0.15 caused by merging of adjacent throats following a process which is roughly the opposite of that described for the FS set of samples. Then, the number of pores and throats stays almost constant for θ > 0.15. As a result, the increase of porosity is mainly due to the increase of the throat length λ and radius r h . The distribution of r h / r h is almost Gaussian around the mean value, and independent of the porosity for θ > 0.15. The transition from the original sample FS10 to the FSD12 and then FSD15 is well visible the r h distribution. Note that, as soon as the throats are widely distributed like for the FS set of samples, κ is an indicator of the potential local flow rate disorder at the network nodes because the probability of having upstream and downstream bonds of distinctly different flow rates is high.
Altogether, these results show that the two sets of samples are very different in terms of (1) the topology of the network; for the FSD set, the topology is almost similar for all the porosity range, while it is increasingly complex (with increasing tortuosity, see discussion below) as porosity decreases for the FS set of samples, and (2) the characteristic size of the throats which is almost independent of the porosity for the FD set whereas it increases with porosity for the FSD set.

Permeability and Flow Field Properties
Permeability values k for the 12 samples computed using Darcy's law (k = v z µ/∇ * p) are plotted in the left panel of Figure 5. Permeability increases from 1.4×10 −13 m 2 for sample SF10 to 6.04×10 −12 m 2 (6.08×10 −12 m 2 ) for sample SF25 (SFD25) and are all aligned with the relation k ∼ θ 4 independently on geometrical characteristics of the pore space. The permeability computed on the BNM (solving a Kirchhoff problem) is also reported Figure 5 in order to evaluate the accuracy of the BNM.
The right panel of Figure 5 displays the advective tortuosity χ a , i.e., the mean tortuosity of the flow lines. The advective tortuosity is obtained from the ratio of the mean Eulerian speed v e to the mean velocity in the direction of the flow v z (Koponen et al., 1996;Ghanbarian et al., 2014;Puyguiraud et al., 2019c): For both the sets of samples, χ a decreases when porosity increases, but it is more pronounced for the FS set of samples. These trends seem to be mainly controlled by the increase of the throat radius as porosity increases, while the topological characteristic of the network plays a minor role which is probably resulting from a complex coupling of the geometrical and topologically parameters discussed above. This makes the advective tortuosity, which is one of the three parameters of the CTRW model proposed by Puyguiraud et al. (2021), an intrinsic characteristic of the hydrodynamic system that is essentially porosity-dependant.
The distributions of the Eulerian mean speed for the 12 samples are plotted in Figure 6. The dissimilarity of the p m (v) curves between the FS and the FSD sets is clearly visible. The FSD samples are displaying almost the same mean speed distributions with power-law trend p m (v) ∼ v α−1 for v < v m with α = 0.245 ± 0.05. For the FS set, the evolution of p m (v) with porosity includes two features. First, p m (v) gradually diverges from a Gamma distribution as porosity increases, with the occurrence of increasingly marked transition between the values of speed larger than the mean (v > v m ) and the power-law slope for the slower speed values. Second, the power-law slope for v ≪ v m increases when porosity decreases, ranging from β = α − 1 = 1.63 for θ = 0.25 to β = 1.75 for θ = 0.10.  These values are in agreement with the value of 1.65 found by Puyguiraud et al. (2021) for the Beara sandstone. As far as we know, they have been very few studies of the correlation between the flow speed distribution and the properties of the pore space microstructures (Siena et al., 2014;Matyka et al., 2016;Alim et al., 2017). For instance, Alim et al. (2017) investigated this issue using numerical simulations in 2-dimensional simple artificial porous media made of circular or elliptical discs placed on a square or triangular lattices with increasing disorder. By extracting and analyzing the corresponding network of tubes, following a procedure quite similar to that implemented for extracting the BNM (section 2.4), they concluded that the flow  distribution is mainly determined by the distribution of fractions of fluid flowing at each of the network node and not by the overall tube size distribution. Our results lead us to a similar conclusion for the complex 3-dimensional porous media studied here. The evolution of the mean flow speed with porosity for the FS set in comparison with the weak evolution of the mean flow speed with porosity for the FSD set appears to be correlated to the noticeable increase with porosity of the number of throats as well as the mean number of throats per pore κ (Figure 3) measured for the FS set, whereas both the number of throats and κ are almost constant for the FSD set of samples.

Speed Decorrelation Distance Length
The decorrelation distance ℓ c is evaluated from the Lagrangian flux weighted speed autocorrelation function ϒ vv (l) where l denotes the lag. The decorrelation distance ℓ c is given by the value of the lag corresponding to ϒ vv (l) = 1/e. The two panels at left of Figure 7 display the Lagrangian flux weighted speed autocorrelation function ϒ vv (l) for the two set of samples. The corresponding values of the decorrelation distance ℓ c vs. porosity are given in the third panel of Figure 7, and the ratio of the decorrelation distance to the mean throat length η = ℓ c /λ vs. porosity is given in the right panel.
For both the sample sets, the decorrelation distance ℓ c increases with porosity from about 150 µm at θ = 0.1 to about 240 µm for FS and 290 µm for FSD. The slight increase of ℓ c for the FSD set for θ > 0.15 compared to the FS set is caused by the increase of the throat radius and the decrease of tortuosity with porosity that are more important for FSD than for FS. The ratio η also displays an increase with porosity following a similar trend for both the FS and the FSD set of samples, the values for FSD being smaller of ∼ 0.5 unit than for FS. Thus, in average, the number of bond lengths traveled before losing the memory of the initial speed ranges from about 2 to 4. These values are in good agreement with the value of 4 obtained by Puyguiraud et al. (2021) by fitting DNS and CTRW for Berea sandstone of porosity 0.18.

Dispersion
In this section, we are presenting the results of the transport DNS, discussing them in the frame of the scaling properties derived from the CTRW model proposed by Puyguiraud et al. (2021) and of the properties retrieved from the BNM (section 3.1).
The first passage time distributions f t (t) (or breakthrough curves) at a distance of 20 times the sample size are given in Figure 8 for Pe = 100 and also for purely advective transport (d m = 0; Pe = ∞). For the latter, all the curves display the power-law tailing that characterize pre-asymptotic (non-Fickian) regime over 3 to 4 orders of magnitude. The scaling f t (t) ∼ t −2−α predicted by Puyguiraud et al. (2021) with the values of α corresponding to those measured on the mean speed distribution is confirmed for all the samples. The comparison of the value of α (0.24 ≤ α ≤ 0.37) for the FS set of samples is given in Figure 9. For Pe = 100, even if it can be considered a quite large value for natural porous media, diffusion acts as increasing the rate at which f t (t) decreases with time and the α-dependent power-law trend is not present. Note that the beginning of the exponential decrease is visible for FSD25 at t ≈ 5τ D , where τ D = ℓ 2 c /d m ≈ 80s. We now focus on determining the asymptotic dispersion coefficient D * from the asymptotic regime of the displacement variance. Figure 10 displays, as an example, the displacement variance normalized to the throat length (σ 2 /λ 2 ) for the 12 samples in the case Pe = 100, but the following comments apply for all values of Pe larger than 1. All curves converge to the asymptotic regime (σ 2 /λ 2 ∼ t) for time t ≥ t a , where t a is independent of the value of Pe but depends on porosity; t a ≈ 10 3 s for θ = 0.1 and t a ≈ 10 4 s for θ = 0.25, i.e., about 40 and 120 times τ D , respectively. This point is important regarding the possibilities of measuring the asymptotic dispersion from laboratory experiments, deriving D * from the breakthrough curves, for instance. For Pe = 100, that corresponds to a mean flow speed of 1.5 × 10 −3 m/s for FS10, a sample of about 1.5 m long displaying the same properties of the mm-scale sample would be necessary to measure D * ; a distance of 60 m would be necessary for Pe = 4,000. This indicates that experimental measurement of D * can be performed only for low values of Pe, typically of the order Pe ≤ 10. However, for such low values of Pe it is not possible to measure α and thus determine the trend D * (Pe).
Conversely, the DNS allows us to perform numerical experiments over large range of Pe values; Figure 11 displays the value of D * vs. Pe for the 12 samples from diffusion-dominant regime (Pe = 10 −3 ) to advection-dominant (Pe = 2 × 10 4 ). These curves can be commented in terms of their slope and of their scaling with porosity, for Pe ≫ 1. Note that for Pe → 0 the ratio D * /d m is equal to the inverse of the diffusive tortuosity (D * /d m = χ −1 d ). For both the FS and FSD sets of samples, the relation D * /d m ∝ Pe 2−α predicted by the CTRW model for Pe ≫ 1 is observed. The values of α compared to those measured using the speed distribution and the tailing of f t (t) are given in Figure 9. The minimum value Pe c at which D * /d m ∝ Pe 2−α is effectively observed, is correlated with the shape of the mean speed distribution (Figure 6). For FS10, the trend p m (v) ∼ v α−1 with α = 0.24 extends up to 5 × 10 −3 v/ v m , while for FS25 the trend α = 0.37 extends up to 3 × 10 −4 v/ v m only. This gives values of Pe c ranging from 1,000 for FS10 to 50 for FS25. The same trend is observed for the FSD set of samples. These results demonstrate the clear control of the particle mean speed distribution on the evolution of D * with the Péclet number. However, both the two sets of samples display a scaling of D * with porosity, independently of the slope determined for Pe ≥ Pe c . The expected decrease of D * for all values of Pe > 1 when porosity increases, corresponding to a decrease of the slope of p m (v) for v ≪ v m is clearly visible for the FS set of samples. But, the results for the FSD set, that share the same mean speed distribution (Figure 6), show also a clear decrease of D * as porosity increases, which indicates that the dispersion scaling with porosity is not solely controlled by p m (v) for v ≪ v m . Indeed, the increase of D * with porosity is also related to the increase of the speed decorrelation distance ℓ c with porosity. In the frame of the CTRW model l c denotes the length at which a new velocity is drawn from the mean speed distribution, and as such ℓ c determines the rate at which the speed changes.
Furthermore, we observe in Figure 11 that D * shows different power-law behaviors for Pe < Pe c that can be related to the scaling behavior of the distribution of mean flow speeds and the transition time distribution. In the limit of infinite Pe, the transition time distribution is given by (Equation 12). For finite Pe, it is cut-off at the diffusion time τ D . The log-slope of ψ(t) at the cut-off time depends on the average flow speed v m . This is shown in Figure 12, which displays the distribution of purely advective transition times rescaled by τ v = ℓ c / v m for FS10 and FS25. The behavior of D * for Pe < Pe c corresponds to the powerlaw scaling of ψ(t) at dimensionless times equal to Pe. The slope FIGURE 12 | Distribution of advective transition times rescaled by τ v for FS10 and FS25. The dimensionless cut-off time is Pe = τ D /τ v . The vertical lines denote Pe = 10 (dashed lines), Pe = Pe c (solid lines) and Pe = 4,000 (dot line). The sloped lines denote the power-law behaviors t −2−α ′ for Pe < Pe c with α ′ = 0.38 and 0.79 for FS10 and FS25, respectively, and t −2−α for Pe ≥ Pe c with α = 0.23 and 0.37 for FS10 and FS25, respectively. of the ψ(t) curves display the power-law behaviors t −2−α ′ for Pe < Pe c with α ′ = 0.38 and 0.79 for FS10 and FS25, respectively. For Pe ≥ Pe c the values of α are similar to those reported in Figure 9 for f t (t), p m (v) and D * (Pe), i.e., α = 0.23 and 0.37 for FS10 and FS25, respectively.

SUMMARY AND CONCLUSIONS
We performed numerical experiments of passive solute transport for two sets of porous media mimicking a large range of porosity and microstructures expected in sandstones. The aim was to test the validity of the CTRW model, to explore how the flow field characteristics are linked to the porous media geometrical properties and to determine the scaling of asymptotic dispersion coefficient D * with the Péclet number. The two sets of six samples share similar porosity, ranging from 0.1 to 0.25, and the same permeability-porosity trend k(θ ) but displays distinctly different microstructures and thus dispersion evolution.
The conceptual CTRW model of solute transport in porous media, as the one proposed by Puyguiraud et al. (2021), infers that solute spreading along particle paths is controlled by the transition time of the solute particles which is determined by the distribution of solute particle mean speeds p m (v), the velocity decorrelation distance ℓ c and diffusion. The effective tortuosity factor that depends on Pe and on the advective tortuosity χ a (that can be also easily evaluated form the flow field) allows mapping dispersion in the streamwise direction which is aligned with the mean pressure gradient. With decreasing Pe, the effective tortuosity of the solute particles increases and the control of p m (v) on dispersion decreases but remains important up to high values of Pe because of the wide distribution of the particles speeds toward low speed values. This means that for heterogeneous media, such as sandstones, the pre-asymptotic (non-Fickian) dispersion regime is likely to persist over long time scales.
We found that the scaling properties, measured by the coefficient α, predicted by Puyguiraud et al. (2021)'s model are effectively measurable for all the 12 studied samples. For instance, results shows that at high Pe, the tail of the breakthrough curves, that is controlled by the low flow speeds, scales as f t (t) ∼ t −2−α where α is given by the slope of the mean speed distribution p m (v) ∼ v α−1 , for v < v m . As Pe decreases, diffusion eventually dominates over low flow speeds, thus cuts off the power-law tail of the breakthrough curves and leads to Fickian behavior from which the asymptotic dispersion coefficient D * can be theoretically evaluated (Van Genuchten and Wierenga, 1986). However, the analysis of the displacement variance σ 2 (t) indicates that D * cannot be measured experimentally at laboratory scale, for high values of Pe, because the distance required for reaching the asymptotic regime is orders of magnitude larger than what is workable at laboratory scale. Thus, measuring experimentally the value of α, for determining how D * scales with Pe seems difficult.
The asymptotic dispersion coefficient D * was computed up to the largest values of Pe expected for laminar flow in natural environments. Results show that D * /d m ∝ Pe 2−α from Pe c up to the highest value of Pe (Pe = 4,000). Note that for the values of α expected for such heterogeneous rock samples, neither the trend D * ∼ Pe ln(Pe) (Saffman, 1959;Koch and Brady, 1985) assuming that the distribution of flow speeds is flat (α = 1), nor the trend D * ∼ Pe expected for α > 1 at high Pe are expected. For 1 < Pe < Pe c , D * /d m ∝ Pe 2−α ′ where α ′ > α depends on the mean speed distribution and the speed decorrelation distance ℓ c that are the parameters that determine the advective particle transition distribution and subsequently the value of Pe c . The mean particle speed remains correlated for longer distances in porous media with straighter and larger bonds (throats). As such ℓ c is a good indicator of the complexity of flow field, because it encompasses the effect of tortuosity that ubiquitously decreases with increasing porosity and the effect of the mean throat radius that ubiquitously increases with porosity, while the other structural parameters are distinctly different for the two sets of samples. Yet, when reported in term of number of bonds length traveled before speed decorrelates, it is observed that FS and FSD sets behave quite similarly; the equivalent number of pores (intersection nodes) crossed before losing the memory of the initial speed equals η − 1 and ranges from about 1 for θ = 0.1 to about 3 for θ = 0.25. We conjecture that the increase of the number nodes crossed before speed decorrelates is linked to the speed changes caused by the splitting of the flow at the network node and thus to both the mean radius of the bonds and the coordination number κ. Similar conjecture can be done for the distribution of the solute mean speed p m (v) which should be controlled by the speed changes caused by splitting of the flow where throats are connected, as it was anticipated by Alim et al. (2017) in numerical simulations in 2-dimensional simple artificial networks. The structural and hydrodynamic mechanisms that determine the flow distribution in 3-dimensional porous media, focusing on the impact pore size distribution, coordination number and local correlations on the speed distributions will be discussed in a forthcoming paper. Yet, from the results presented in this paper, one can conclude that the flow distribution, and thus the mean speed, are controlled by the distribution of fractions of fluid flowing at each of the network nodes which in turn is determined by the distribution of the throat radius (and not the mean) and the coordination number. At given porosity and mean bond radius the latter is controlled by the number of throats per unit volume that increases with porosity for the FS set and decrease with porosity for the FSD set of samples. We believe that these results give a first insight into both the mechanisms and the microstructural parameters that control dispersion in porous media.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.