Temporal Evolution of Solute Dispersion in Three‐Dimensional Porous Rocks

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


Introduction
The transport of solutes in porous media is driven by the phenomenon of dispersion, which results from the interplay between advective spreading and diffusion.The former is triggered by the spatial variability of the fluid speed which is controlled by the geometry of the connected pore network (Alim et al., 2017;Datta et al., 2013;Puyguiraud et al., 2021;Valocchi et al., 2018) while the latter is ubiquitously controlled by the concentration gradients.The heterogeneity of the porous medium, which determines the flow distribution is therefore a primary parameter that controls dispersion from the pre-asymptotic to the Fickian regimes (Dentz et al., 2004;Sherman et al., 2021).Transport in porous media is considered in many fields of academic and industrial applications from material science, engineering and medicine to groundwater hydrology, environmental technologies and petroleum engineering, and at many scales from microfluidic applications to groundwater management.Beside being necessary for understanding and predicting the spreading of chemicals such as pollutants or bionutrients, modeling dispersion is required also to understand and predict solute-solute and solute-minerals reactions that can produce new solute species and trigger mineral dissolution and precipitation features, for instance.Dispersion in porous media has been extensively studied from the pore to the regional scale for decades (Dagan, 1990;Dentz et al., 2023;Gelhar & Axness, 1983;Saffman, 1959;Whitaker, 1967).Here we focus on hydrodynamic dispersion due to velocity fluctuations caused by the heterogeneity of the pore space.A main challenge concerns how continuum-scale solute transport can be modeled by macroscopic parameters, such as the dispersion coefficient, that can be inferred experimentally, by using direct pore scale simulations or upscaling methods such as volume averaging or stochastic modeling (Ahmadi et al., 1998;Bijeljic & Blunt, 2006;Brenner, 1980;Koch & Brady, 1985;Le Borgne et al., 2011;Puyguiraud et al., 2021;Scheven, 2013;Souzy et al., 2020).Dispersion and its pre-asymptotic behavior have been analyzed in terms of breakthrough curves, the time evolution of the spatial variance of concentration or particle distributions, or directly from particle velocities, using experiments and direct numerical pore scale simulations (Bijeljic et al., 2004;Gouze, Puyguiraud, Porcher, & Dentz, 2021;Gouze, Puyguiraud, Roubinet, & Dentz, 2021;Hulin & Plona, 1989;Khrapitchev & Callaghan, 2003;Puyguiraud et al., 2021;Sole-Mari et al., 2022).These studies show that the pore structure shapes the evolution of dispersion during the pre-asymptotic regime and then determines the asymptotic value.Hulin and Plona (1989) and Khrapitchev and Callaghan (2003) study the reversibility of pore-scale dispersion upon flow reversal, which addresses the issue of under which conditions hydrodynamic dispersion describes solute mixing or advective solute spreading.As mentioned above, the fundamental mechanisms of hydrodynamic dispersion are pore-scale velocity fluctuations and diffusion.The former mechanism is reversible in the Stokes regime, which holds for typical applications in groundwater resources.Irreversibility, or actual solute mixing is induced by the interaction of spatial velocity fluctuations and molecular diffusion (Dentz et al., 2023).Consider for example, a solute that evolves from an extended areal source.At early times, the solute front deforms due to velocity variability within the source distribution, which leads to a complex concentration distribution, which nevertheless is partially reversible.Hydrodynamic dispersion coefficients that are defined in terms of the spatial variance of the global solute distribution, measure at pre-asymptotic this advective spreading rather than actual solute mixing.This issue was recognized by Kitanidis (1988) in the context of solute dispersion in heterogeneous porous formations, and Bouchaud and Georges (1990) in the context of random walks in quenched disordered media.These authors propose to define dispersion coefficients from the second-centered moments of the solute or particle distributions that evolve from a point-like initial condition.In the absence of local scale dispersion or molecular diffusion, these dispersion coefficients are exactly zero.In the following, we refer to this concept as effective dispersion.The dispersion concept based on the spatial variance of the solute concentration evolving from an extended areal source, is called ensemble dispersion because it describes the dispersion of the ensemble of partial plumes (Green functions) of which the solute distribution is composed.As outlined above, at preasymptotic times ensemble dispersion measures advective solute spreading.In fact, it measures the center of mass fluctuations of the partial plumes that evolve from the point injections that constitute the spatially extended initial distribution (Bouchaud & Georges, 1990).Thus, the two different dispersion coefficients and their temporal evolution highlight different aspects of pore-scale transport, and allow to study the impact of medium heterogeneity and small scale mass transfer mechanisms on large scale dispersion.
Here, we use the longitudinal effective and ensemble dispersion coefficients to investigate pore-scale mixing and spreading, and to shed new light on the role of the medium structure and pore-scale diffusion for continuum scale solute transport.To this end, we perform three-dimensional direct numerical simulations of pore-scale flow and solute transport in a sand-pack medium and in a Berea sandstone of distinctly different heterogeneity levels, which can be measured, for instance, by the variance of the logarithm of the flow velocity distribution.While the time behaviors of both the effective and ensemble dispersion coefficients are governed by the interplay between pore-scale diffusion and velocity fluctuations, they evolve on different characteristic time scales, which reflect the relevant pore-scale mass transfer and heterogeneity mechanisms.These findings provide new insights into the pore-scale processes that govern continuum scale solute dispersion and mixing.
The paper is organized as follows.The methodology used to calculate flow and transport and measure dispersion are presented in Section 2. Section 3 provides a detailed study of the temporal evolution of the center of mass positions and the longitudinal effective and ensemble dispersion coefficients for the sand pack and Berea samples.It discusses the different time regimes and the implications for our understanding of the impact of pore-scale heterogeneity and mass transfer mechanisms on continuum scale transport.

Pore-Scale Flow and Transport
Flow in three-dimensional porous media, described as dual solid-void structures, is described by the Stokes equation together with the continuity equation (Leal, 2007), Water Resources Research where μ is the dynamic viscosity, u(x) is the Eulerian velocity and p(x) is the fluid pressure at position x = (x 1 , x 2 , x 3 ).Here, flow is driven by the macroscopic pressure gradient, which is aligned with the x-axis of the coordinated system.Zero-flux boundary conditions are set at the solid-void interface and at the lateral domain boundaries.
Transport of solutes is described by the advection-diffusion equation (ADE) for the solute concentration c(x, t) ∂c(x, t) ∂t where c(x, t) is the solute concentration at position x and time t, and D is the molecular diffusion coefficient.The advection-diffusion Equation 2 is equivalent to the Langevin equation (Risken, 1996) where ξ(t) is a Gaussian white noise with mean 〈ξ i 〉 = 0 and covariance 〈ξ j (t)ξ k (t)〉 = δ jk δ(t t′); δ jk is the Kronecker delta.
The characteristic pore length ℓ 0 , the mean velocity magnitude 〈v〉 = 〈|v(x)|〉 and the diffusion coefficient D set the advection time τ v = ℓ 0 /〈v〉 and the characteristic diffusion time τ D = ℓ 2 0 /D.The two time scales are compared by the Péclet number Pe = τ D /τ v = 〈v〉ℓ 0 /D.

Dispersion Coefficients
We analyze pore-scale solute dispersion and the underlying heterogeneity and transport mechanisms in terms of the temporal evolution of suitably defined longitudinal dispersion coefficients.To this end, we consider the concentration distribution c(x, t) that evolves from the normalized plane source where Ω f denotes the fluid domain.We choose a uniform injection across a plane, which mimics the mixing interface between two initially segregated constituents.The indicator function I( ⋅ ) is one if its argument is true and zero else, w and h denote the width and height of the medium and ϕ is porosity.The injection plane is large enough such that where Ω denotes the bulk domain, that is, the union of fluid domain and solid domain.The solute distribution can be decomposed into partial plumes (Green functions) g(x, t|x′) that satisfy Equation 2 for the initial conditions The concentration distribution c(x, t) can be written in terms of g(x, t|x′) as

Effective Dispersion
We define the longitudinal effective dispersion coefficient in terms of the second centered moments of the Green function g(x, t|x′) Water Resources Research which are measures for the spatial extension of the Green function.The ith longitudinal moment of g(x, t|x′) is defined by The first moments m 1 (t; x′) determine the center of mass position of g(x, t|x′).The average of κ L (t; x′) over all Green functions defines the effective second centered moment It is a measure for the average streamwise width of the Green function.Half the temporal rate of growth of κ eff L (t) defines the effective dispersion coefficients Thus, effective dispersion characterizes the typical spreading behavior for a solute evolving from a point-like injection.

Ensemble Dispersion
The longitudinal ensemble dispersion coefficients are defined in terms of the moments of the global concentration c(x, t), As per the second equality sign, the moments can be seen as averages over the ensemble of Green functions and as such are named the ensemble moments in the following.The second centered ensemble moment is defined by It is a measure for the streamwise spatial extension of the concentration distribution, or equivalently for the ensemble of Green functions.The temporal rate of growth of the second centered ensemble moments is measured by the ensemble dispersion coefficients

Center of Mass Fluctuations
The difference between the ensemble and effective variances, quantifies the variance of the center of mass fluctuations of the Green functions that constitute the solute plume.
Along the same lines, the difference between the ensemble and effective dispersion coefficients measures the dispersion of the center of mass positions of the Green functions Water Resources Research 16)

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

Porous Media
We study two three-dimensional porous media of different complexity, (a) a Berea sandstone sample and (b) a sand pack sample, which are illustrated in Figure 1.The Berea sample displays a complex pore structure with a porosity of ϕ = 0.18 and a permeability of 6.7 × 10 12 m 2 (Gouze, Puyguiraud, Roubinet, & Dentz, 2021).Sandstones are more or less diagenized (altered and cemented) buried sand structures and the Berea sandstone is often considered as a typical example of high permeability sedimentary reservoir (Churcher et al., 1991).The sand pack sample has a high porosity of ϕ = 0.37 with a more regular structure of the pore space and a higher permeability value of 2.9 × 10 11 m 2 .For both samples the characteristic pore diameter was evaluated from image analysis to be about 4.5 × 10 5 m and the characteristic length is ℓ 0 = 1.5 × 10 4 m (Gouze, Puyguiraud, Roubinet, & Dentz, 2021).The sand-pack image (Sand Pack LV60C) is obtained from the Imperial College image repository (Imperial College Consortium on Pore-scale Imaging and Modelling, 2014).It is a compact packing of irregular quartz grains of variable size that is a proxy of sub-surface aquifers (Di Palma et al., 2019).The different degrees of heterogeneity of the two samples can be quantified in terms of the distribution of flow speeds as discussed in the next section.

Fluid Flow
As mentioned above, both pore geometries are based on X-Ray microtomography images.The geometries are meshed using regular hexahedron cells (voxels).This type of mesh has two major advantages.First, it perfectly fits the voxels of the X-Ray tomography images, and second, it allows for a simple and computationally efficient velocity interpolation scheme, which is required for the transport simulation based on random walk particle tracking (Mostaghimi et al., 2012).Each of the images is decomposed in 900 3 voxels of length l m = 1.060 × 10 6 m for the Berea sandstone and l m = 5.001 × 10 6 m for the sand pack.This implies that the throat (or pore constriction) diameters are resolved in average with 9 voxels for the sand pack and 45 voxels for the Berea sample, which ensures an accurate computation of the flow field (Gouze, Puyguiraud, Porcher, & Dentz, 2021).Note that the voxel size for the sand pack is comparable to the one chosen in Sole-Mari et al. (2022).Fluid flow in the pore space is solved numerically using the SIMPLE algorithm implemented in OpenFOAM (Weller et al., 1998).Pressure boundary conditions are set at the inlet (x = 0) and outlet (x = 900l m ) of the domains.No-slip boundary conditions are prescribed at the void-solid interface and at the lateral boundaries of the domain.Note that imposing impermeable lateral boundaries impacts transverse dispersion but hardly affects longitudinal dispersion and mixing which are controlled by the velocity variation along the trajectory and transverse diffusion within a pore (or throat).Once the solver has converged, the flow velocities are extracted at the centers of the interfaces of the mesh (i.e., at the six faces of each of the regular hexahedra that form the mesh) in the normal direction to the face.
The ratio between the mean flow speed 〈v〉 and the mean flow velocity 〈u〉 in streamwise direction defines the advective tortuosity χ = 〈v〉/〈u〉.For the Berea sample, we find χ = 1.64, and for the sand pack χ = 1.32.Since for Stokes flow, the flow velocities scale with the pressure gradient, the flow field is determined for a unit gradient and then scaled for the Péclet scenario under consideration.For example, for Pe = 200, the mean flow speeds are 〈v〉 = 2.67 × 10 3 m/s.The mean streamwise velocities can be obtained from the respective tortuosity values.Figure 1 shows the distribution of flow speeds for the two media, which can be used to quantify the flow heterogeneity of the two porous medium samples and highlight the differences between them (Alhashmi et al., 2016).

Random Walk Particle Tracking
Solute transport is modeled using random walk particle tracking (Noetinger et al., 2016).The numerical simulation is based on the discretized version of the Langevin Equation 3, The central limit theorem ensures that the sum of these uniform random variables is Gaussian distributed with zero mean and unit variance.The particle velocities u[x(t)] are interpolated from the velocities at the voxel faces using the algorithm of Mostaghimi et al. (2012), which implements a quadratic interpolation in the void voxels that are in contact with the solid and thus guarantees an accurate representation of the flow field in the vicinity of the solidvoid interface.The time step is variable and chosen such that the particle displacement at a given step is shorter than or equal to the side length of a voxel.The time step varies from Δt = 10 8 s at early times to get an accurate resolution of the moments to Δt = 10 3 s at late times to ensure faster simulations.The diffusion coefficient is set to D = 10 9 m 2 /s.
To investigate the effective and ensemble dispersion coefficients, particles are seeded at all voxels in the void space within a plane perpendicular to the mean flow direction, see Figure 2.There are N 0 = 1.5 × 10 5 voxels in the injection plane for the Berea sample, and N 0 = 3 × 10 5 for the sand pack.At each position, N p = 10 3 particles are injected, which makes a total of N = 1.5 × 10 8 particles for the Berea sandstone and N = 3 × 10 8 for the sand pack.We consider this scenario for Pe = 200 and Pe = 2,000.Each partial injection of particles defines a partial plume, or Green function.The first and second moments of a partial plume originating from x′ j are then obtained from the particle positions where i = 1, 2 and the subscript k denotes the kth particle.The effective and ensemble variances are given by

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

Center of Mass
Figure 3 shows the evolution of the streamwise center of mass position m 1 (t) of the global solute distribution c(x, t) in the top panels.The bottom panels show the rate of change δD m L (t) of the variance of the center of mass positions m 1 (t|x′) of partial plumes g(x, t|x′) defined by Equation 16.The center of mass of the global plume moves with the mean flow velocity 〈u〉, while the center of mass velocities of the partial plumes evolve from the velocities at the respective injection points toward the mean flow velocity.At short times t ≪ τ v , that is, travel distances shorter than the average pore size, the center of mass velocities are approximately constant, which implies m 1 (t; x′) = u 1 (x′)t and therefore where σ 2 0 denotes the initial velocity variability.The initial particle velocities persist until the plume starts sampling the flow field by transverse diffusion across streamlines, and by advection along the streamlines.This ballistic early time regime is observed for both the sand pack and Berea samples.

Sand Pack Sample
The evolution of δD m L (t) for the sand pack sample is characterized by two regimes.The early time ballistic regime, and a sharp decay after a maximum that is assumed on the advective time scale τ v .This is at first counter-intuitive because transverse diffusion is the only mechanism that makes the partial plume sample the flow heterogeneity such that the differences between the center of mass positions of different partial plumes decrease.Thus, one would expect that the relevant time scale is set by the characteristic pore length and diffusion, that is, by the diffusion time τ D .Sampling occurs indeed by diffusion in transverse direction.However, the distance ℓ c to sample a new velocity depends on the flow rate because streamtubes in low velocity regions are wider than in high velocity regions.Since the flow rate is constant in a streamtube, Q c = ℓ 2 c 〈v〉, with Q c a characteristic flow rate, the decorrelation length becomes ℓ c = ̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ Q c /〈v〉 √ .Thus, the time scale at which particles decorrelate is From Figure 3, we observe that τ c ≈ τ v , which means that the characteristic flow rate is Q c ≈ Dℓ 0 .

Berea Sandstone Sample
For the Berea sample, we observe three different regimes for δD m L (t).The early time regime is ballistic as discussed above.The start of the second regime is marked by the advective time scale τ v as observed for the sand pack.Here, however, δD m L (t) does not assume a maximum on the advective time scale and then decays, but keeps increasing until the diffusion time τ D , where it reaches maximum and then shows a rapid decay.A possible explanation of this behavior in the second time regime is the contrast between particles that are initialized at moderate to high flow velocities, which decorrelate due to transverse velocity sampling on the one hand, and the persistence of particles in low velocity conducts on the other hand, which can give rise to the observed sub-linear increase of δD m L (t).These low velocities are eliminated on the time scale τ D , which sets the maximum transition time along a conduct.In other words, transition times of particles that move at low velocities along a conduct are cut-off at the diffusion time scale (Puyguiraud et al., 2021).
In summary, the evolution of the center of mass fluctuations is marked by the advection time scale for the sand pack sample, and by the advection and diffusion time scales for the Berea sample.The fact that the intermediate regime is not present for the sand pack sample can be explained by the spatial medium structures of the two samples shown in Figure 1.The structure of the Berea sample can be seen as a connected network of conducts, while the sand pack is more a connected network of pore bodies.These differences are also reflected in the evolutions of the effective and ensemble dispersion coefficients discussed in the next section.

Ensemble and Effective Dispersion
Figures 4 and 5 show the evolution of the effective and ensemble dispersion coefficients for the sand pack and Berea samples.One observes a marked difference between the ensemble and effective dispersion coefficients at short and intermediate times.At early times t < τ 0 = D/〈v〉 2 = Pe 1 τ v , diffusion dominates over advection, and both the ensemble and effective dispersion coefficients are equal to the molecular diffusion coefficient D. Note that the time scale τ 0 here is much smaller than the diffusion time across a pore throat.Thus, at t < τ 0 the diffusion process is not hindered by the confining porous medium and therefore both dispersion coefficients are equal to the molecular diffusion coefficient.
For τ 0 < t < τ v , advection starts dominating over diffusion.As outlined in the previous section, particles are transported at their initial velocities that persist over the characteristic length scale ℓ 0 .Thus, the ensemble dispersion coefficients evolve ballistically in this regime Water Resources Research where σ 2 0 is the initial velocity variance.It behaves in the same way as ΔD m L (t), see Equation 21.
This effect of the center of mass fluctuations between partial plumes is removed by the definition of the effective dispersion coefficients as the average dispersion coefficient of the partial plumes.For τ 0 < t < τ v , a partial plume is translated by its initial velocity.As its size increases by diffusion, the plume gets sheared by the transverse velocity contrast.Therefore, the effective dispersion coefficients D eff L (t) first remain at the value of the molecular diffusion coefficient and then increase steeply due to shear dispersion.Figures 4b and 5b show that the increase of the effective dispersion coefficients occurs for high Pe at earlier non-dimensional times than for low Pe.This indicates that the shear rate does not scale linearly with 〈u〉.In fact, a typical shear rate can be written as where ℓ γ is the scale of transverse velocity contrast.The latter is proportional to the typical streamtube size.That is, as ℓ 2 γ 〈v〉 = constant, we have ℓ γ ∼ 〈v〉 1/2 .The characteristic shear length scale decreases with increasing flow rate, and thus the shear rate scales as γ ∼ 〈u〉 3/2 .The characteristic shear time scale is then τ γ = γ 1 ∝ v /〈v〉 1/2 .This dependence explains the differences in the time behaviors of the effective dispersion coefficients for different Pe.The early time ballistic and shear dispersion behaviors for t < τ v are observed for both the sand pack and Berea samples.For t > τ v the dispersion behaviors are different as discussed in the following.

Sand Pack Sample
Figures 4a-4d show the evolution of the ensemble and effective dispersion coefficients for the sand pack sample.For times t > τ v , that is for mean travel distances larger than the average pore size, particles start sampling different flow velocities along their trajectories, and the ballistic behavior for the ensemble dispersion coefficients breaks down, see Figure 4a.
For purely advective transport, the ensemble dispersion coefficients continue growing non-linearly with time, which can be traced back to the broad distribution of transition time across pores (Puyguiraud et al., 2019).At finite Pe, the ensemble dispersion coefficients first follow the purely advective behavior and eventually cross over toward their asymptotic value on the time scale.The effective dispersion coefficients shown in Figure 4 cross over toward their asymptotic values, also on the time scale τ v .As shown in Figures 4c and 4d, they converge with D ens L (t).
As mentioned in Section 3.1, these behaviors are at first sight counter-intuitive because we expect the deviation from the purely advective behavior observed for D ens L (t) and the convergence of D eff L (t) toward D ens L (t) to be governed by diffusion.For ensemble dispersion, diffusion is the mechanism that decorrelates subsequent (low) velocities in time and thus leads to the separation of D ens L (t) from the (anomalous) purely advective behavior.Similarly, the mechanism by which the effective dispersion coefficients converge toward the ensemble dispersion coefficients is due to decorrelation of the velocities of the particles that start from the same point, which is due to diffusion in transverse direction.Thus one would expect that the dispersion coefficients evolve on the diffusion time scale τ D .
As discussed in Section 3.1.1,the decorrelation mechanism is indeed transverse diffusion across a length scale that is related to a typical streamtube width.Thus, the decorrelation time τ c is given by Equation 22, which is proportional to τ v .This observation explains the temporal evolution of the ensemble and effective dispersion coefficients on the time scale τ v .
Sole-Mari et al. ( 2022) study the longitudinal ensemble dispersion coefficients for transport in granular porous media.These authors find that the evolution time-scale depends on Pe and varies between the diffusion time scale for Pe smaller than a critical Pe c and the advection time scale for Pe ≫ Pe c .The estimate of these authors is based on the heuristic expression which we slightly modified such that D ens L (t = 0) = D is equal to the molecular diffusion coefficient.The asymptotic mechanical dispersion coefficient is denoted by δD ∞ L and β is a fitting parameter, which Sole-Mari et al. ( 2022) relate to Pe. Fits are performed using non-linear least squares.We find β ≈ 16 for Pe = 2,000 and β ≈ 13 for Pe = 200.As shown in Figure 4a, this expression can describe the crossover behavior from ballistic to asymptotic, but is not able to reproduce the ballistic behavior observed at times t ≪ τ v because it behaves as Here, we use the following heuristic expression to fit the data for the longitudinal ensemble dispersion coefficients with α a fitting parameter.At short times t ≪ τ v , Equation 26shows the correct ballistic behavior as D ens L (t) = D + σ 2 0 t. Figure 4a shows the best fits with α = 3.1 for Pe = 200 and α = 3.9 for Pe = 2,000.Note however, that unlike Equation 25, Equation 26 describes too fast a crossover from the ballistic to the asymptotic behavior.The effective dispersion coefficients show a more complex behavior and cannot be described by these simple regressions.

Berea Sandstone Sample
Figures 5a-5d show the evolution of the ensemble and effective dispersion coefficients for the Berea sandstone sample.As seen in Figure 5a, the initial ballistic behavior for the ensemble dispersion coefficients breaks down on the time scale τ v when particles start sampling different flow velocities along their trajectories.For purely advective transport, we observe anomalous dispersion characterized by a super-linear growth of the ensemble dispersion coefficients, which can be traced back to broad distributions of advective particle transition times (Puyguiraud et al., 2019).Unlike for the sand pack, here the crossover toward the constant asymptotic long time values occurs on the diffusion time scale τ D .As discussed in Section 3.1.2,the temporal decorrelation of low velocities is due to diffusion along pore channels with the characteristic time scale τ D (Puyguiraud et al., 2021).Also the crossover of the effective dispersion coefficient toward the asymptotic value shown in Figure 5b occurs on the time scale τ D .
The convergence of the effective and ensemble dispersion coefficients shown in Figures 5c and 5d on the other hand, occurs on the decorrelation time scale τ c , see Equation 22.This time scale is set by transverse diffusion across streamtubes, which is the mechanism by which particles that originate from the same initial position start decorrelating and sampling different flow velocities.Thus, for t > τ c particle trajectories become independent, particles "forget" that they originate from the same initial condition, and dispersion is enhanced due to the increasing difference in travel distance between different particles.This is the same mechanism that causes longitudinal ensemble dispersion, and therefore the two dispersion coefficients converge.These complex behaviors cannot be captured by simple heuristic expressions as Equations 25 and 26, but require a more in-depth theoretical analysis, which is beyond the scope of this paper.

Conclusions
We investigate solute dispersion in three-dimensional porous rocks using detailed numerical simulations of porescale flow and transport.We consider a sand-like medium, and a Berea sandstone sample.The two media have quite distinct pore structures, which manifests in distinct pore-scale flow variability.The latter is quantified by the distribution of Eulerian flow speeds.The degree of flow heterogeneity is measured by the variance of the logarithm of the flow speed, which is significantly higher for the Berea sample than for the sand pack sample.
We use longitudinal effective and ensemble dispersion coefficients to probe pore-scale mixing and spreading behaviors and shed light on the impact of medium structure and diffusive mass transfer on solute and particle transport.Effective dispersion coefficients are defined in terms of the spatial average of the second-centered moments of the partial plumes (Green functions) that constitute the global solute distribution.Thus, they can be seen as a measure for the typical width of a mixing front.Ensemble dispersion coefficients on the other hand are defined in terms of the second centered moments of the global solute plume and as such are a measure for mixing front deformation due to the flow variability within the initial plume.The fundamental mechanisms that cause hydrodynamic dispersion are pore-scale flow variability and molecular diffusion, which govern the evolution of both the effective and ensemble dispersion coefficients.Their detailed time behaviors and the time scale on which they evolve provide information on the salient small scale mass transfer and heterogeneity mechanisms.
The early time behavior of the ensemble coefficient is ballistic as a result of the spatial persistence of flow velocities in the initial plume.The effective coefficients on the other hand are significantly smaller than their ensemble counterparts.Their early time evolution is dominated by shear dispersion, which results from the velocity gradients within the partial plumes, whose lateral extent initially increases by diffusion.The two dispersion coefficients start converging when the lateral extent of the partial plumes is large enough for the efficient sampling of the flow heterogeneity, and it is here, where dispersion in the sand pack and Berea sandstone behaves differently.For the sand pack, the evolution of effective dispersion is marked by the characteristic diffusion time across a streamtube, which sets the time for both convergence to ensemble dispersion and its asymptotic behavior.For the Berea sandstone, this time scale marks the time for convergence of effective and ensemble dispersion, which, however, still evolve non-linearly with time until they retrieve their asymptotic long time value on the time scale for diffusion over a typical pore length.The evolution of solute dispersion reflects the medium structure, which determines the microscopic mass transfer mechanisms.In fact, we hypothesize that the observed behaviors originate in the network-like medium structure in case of the Berea sample, and the strong connectivity of pores in the sand pack.While the behavior of ensemble dispersion can be captured by travel-time based approaches like the continuous time random walk in terms of flow variability and medium structure, it is still elusive how to quantify effective dispersion in these terms.Nevertheless, these findings give insight into how particles sample pore-scale velocity variability and how it depends on the medium structure.These are key elements for the upscaling of transport, mixing and reaction from the pore to the continuum scale.

Figure 1 .
Figure1.Eulerian velocity pdfs for the sand pack (blue circles) and the Berea sandstone (red squares).Inlay: The three-dimensional pore geometry of (left) the sand pack sample (5 mm 3 ) and of (right) the Berea sandstone (1 mm 3 ).The gray and blue colors represent the pore space and the solid phase, respectively.

Figure 2 .
Figure 2. Snapshots of the conservative simulation for the Berea sandstone for Pe = 2,000 at three different times t = 0.15τ v , t = 0.8τ v , and t = 5τ v .The density of particles represents the concentration.

Figure 3 .
Figure 3. Temporal evolution of the center of mass position for (a) sand pack and (b) Berea samples for Pe = 200 and Pe = 2,000.Temporal evolution of the variance of the center of mass position for (c) the sand pack and (d) Berea samples for Pe = 200 and Pe = 2,000.The dashed vertical lines denote (black) the advection time scale τ v , (yellow and orange) the respective diffusion time scales τ D .

Figure 4 .
Figure 4. Dispersion coefficients for the sand pack sample.Top panels: (a) Ensemble dispersion coefficients for (red solid line) Pe = 2,000 and (orange solid line) Pe = 200, and (b) corresponding effective dispersion coefficients.The vertical dashed lines denote the corresponding diffusion time scale τ D = τ v Pe.The dashed lines in panel (a) show expression (26) for (purple) α = 3.1 and (blue) α = 3.9.The dotted lines show expression (25) for (purple) β = 13.1 and (blue) β = 16.3.Bottom panels: (Black solid lines) Ensemble and (blue solid lines) effective dispersion coefficients for (c) Pe = 200 and (d) Pe = 2,000.The vertical black dashed lines denote the decorrelation time scale τ c = τ v , the blue dashed lines denotes the respective diffusion time scales.The horizontal dash-dotted lines denote the asymptotic short time and long time values.

Figure 5 .
Figure 5. Dispersion coefficients for the Berea sandstone sample.Top panels: (a) Ensemble dispersion coefficients for (red solid line) Pe = 2,000 and (orange solid line) Pe = 200, and (b) corresponding effective dispersion coefficients.The vertical dashed lines denote the corresponding diffusion time scale τ D = τ v Pe.Bottom panels: (Black solid lines) Ensemble and (blue solid lines) effective dispersion coefficients for (c) Pe = 200 and (d) Pe = 2,000.The vertical black dashed lines denote the decorrelation time scale τ c = τ v , the blue dashed lines denote the respective diffusion time scales.The horizontal dash-dotted lines denote the asymptotic short time and long time values.